###################################################
### chunk number 1: q1.1
###################################################
#line 59 "Lab1_2010_answers.Rnw"
options(useFancyQuotes=F)
library(foreign)
dat <- read.dta("http://www.quantoid.net/lab1_nes.dta")
dat$pidfac <- as.factor(dat$pid)
mod.pid <- lm(demtherm ~ age + racerec + pid, data=dat)
mod.pid2 <- lm(demtherm ~ age + racerec + pidfac, data=dat)
anova(mod.pid, mod.pid2, test="F")
 
 
###################################################
### chunk number 2: q1.2
###################################################
#line 79 "Lab1_2010_answers.Rnw"
mod.pid3 <- lm(demtherm ~ age + racerec + pid3, data=dat)
anova(mod.pid2, mod.pid3, test="F")
 
 
###################################################
### chunk number 3: q2.1
###################################################
#line 110 "Lab1_2010_answers.Rnw"
library(car)
mod <- lm(demtherm ~ age + racerec + income, data=dat)
Anova(mod)
dat$newinc<- as.numeric(dat$income)
mod2 <- lm(demtherm ~ age + racerec + newinc,data=dat)
anova(mod2, mod, test="F")
 
 
###################################################
### chunk number 4: q21fig1
###################################################
#line 123 "Lab1_2010_answers.Rnw"
library(effects)
print(
plot(effect("income", mod))
)
 
 
###################################################
### chunk number 5: q21fig2
###################################################
#line 130 "Lab1_2010_answers.Rnw"
print(
plot(effect("newinc", mod2))
)
 
 
###################################################
### chunk number 6: q22
###################################################
#line 140 "Lab1_2010_answers.Rnw"
library(splines)
AICs <- NULL
	for(i in 3:15){
		knot <- seq(2,24,length=i)[-c(1,i)]
		tmp <- lm(demtherm ~ age + racerec + 
			bs(newinc, knots=knot, degree=3), data=dat)
		AICs <- c(AICs, AIC(tmp))
	}
pdf("aic_plot.pdf", height=6, width=6)
plot(1:13, AICs, type="b", xlab="# Knots", ylab = "AIC")
invisible(dev.off())
	mod <- lm(demtherm ~ age + racerec + 
		bs(newinc, knots=13, degree=3), data=dat)
	summary(mod)
 
 
###################################################
### chunk number 7: q22fig1
###################################################
#line 159 "Lab1_2010_answers.Rnw"
print(plot(effect("bs(newinc, knots=13, degree=3)", mod, default.levels=100)))
 
 
###################################################
### chunk number 8: q22anova
###################################################
#line 166 "Lab1_2010_answers.Rnw"
anova(mod, mod2, test="F")
 
 
###################################################
### chunk number 9: op
###################################################
#line 175 "Lab1_2010_answers.Rnw"
tmp.dat <- na.omit(dat[,c("demtherm", "age", "racerec", "newinc")])
mod2 <- lm(demtherm ~ age + racerec + poly(newinc, 3),
data=tmp.dat)
summary(mod2)
 
 
###################################################
### chunk number 10: q23
###################################################
#line 197 "Lab1_2010_answers.Rnw"
library(mgcv)
gam.mod <- gam(demtherm ~ s(age) + racerec + newinc,data=tmp.dat)
gam.mod2 <- gam(demtherm ~ poly(age,2) + racerec + newinc,data=tmp.dat)
 
anova(gam.mod, gam.mod2, test='Chisq')
pred.dat <- data.frame(
	age = seq(18,90,by=1), 
	racerec = "white", 
	newinc = mean(dat$newinc, na.rm=T))
 
pred1 <- predict(gam.mod, newdata= pred.dat)
pred2 <- predict(gam.mod2, newdata= pred.dat)
 
pdf("gpredplot.pdf", height=6, width=6)
plot(pred.dat$age, pred1, xlab = "Age", ylab = "Predicted 
	Thermometer", type="l")
lines(pred.dat$age, pred2, lty=2)
legend("topleft", c("GAM", "Polynomial"), lty=c(1,2), inset=.01)
invisible(dev.off())
 
 
###################################################
### chunk number 11: q31
###################################################
#line 228 "Lab1_2010_answers.Rnw"
mod <- lm(demtherm ~ racerec + pid3 + libcon*age, data=dat)
 
 
###################################################
### chunk number 12: q31fig1
###################################################
#line 231 "Lab1_2010_answers.Rnw"
print(plot(effect("libcon*age", mod)))
 
 
###################################################
### chunk number 13: q31fig2
###################################################
#line 234 "Lab1_2010_answers.Rnw"
print(plot(effect("libcon*age", mod), x.var="age"))
 
 
###################################################
### chunk number 14: q32
###################################################
#line 242 "Lab1_2010_answers.Rnw"
library(DAMisc)
DAintfun2(mod, c("libcon","age"), name.stem="q32", plot.type="pdf")