###################################################
### chunk number 1: makenonlin
###################################################
#line 210 "Lecture5_2010.Rnw"
library(apsrtable)
set.seed(123)
x <- rep(1:5, 100)
x <- x[order(x)]
 
y <- 2 + x + rnorm(500,0,2)
x <- as.factor(x)
 
mod <- lm(y ~ x)
apsrtable(mod, model.names="")
 
 
###################################################
### chunk number 2: glht1
###################################################
#line 235 "Lecture5_2010.Rnw"
L <- matrix(c(
0, 2, -1, 0, 0, 
0, 3, 0, -1, 0, 
0, 4, 0, 0, -1), ncol=5, byrow=T)
c <- rep(0, 3)
 
e <- matrix(residuals(mod), ncol=1)
s2e <- (t(e)%*% e)/with(mod, df.residual)
b <- coef(mod)
X <- model.matrix(mod)
 
F0 <- (t(L%*%b-c) %*% solve(L%*%solve(t(X) %*% X)%*% 
    t(L))%*%(L%*%b - c))/(3*s2e)
 
cat("F = ",F0, "\n", sep="")
cat("Pr(>F) = ", 1-pf(F0, 3, with(mod, df.residual)), "\n", sep="")
 
 
###################################################
### chunk number 3: nlsimfig
###################################################
#line 260 "Lecture5_2010.Rnw"
plot(seq(1,5), c(0, b[2:5]), pch=19, type="o", xlab="x", ylab="predicted")
abline(a=-1, b=1, lty=2)
 
 
###################################################
### chunk number 4: anovatest
###################################################
#line 278 "Lecture5_2010.Rnw"
library(xtable)
restricted.mod <- lm(y ~ as.numeric(x))
unrestricted.mod <- lm(y ~ x)
xtable(anova(restricted.mod, unrestricted.mod, test="F"))
 
 
###################################################
### chunk number 5: glmlin
###################################################
#line 292 "Lecture5_2010.Rnw"
library(foreign)
anes <- read.dta("anes1992.dta")
anes[["pidfac"]] <- as.factor(anes[["pid"]])
unrestricted.mod <- glm(votedem ~ retnat + pidfac + age + male + educ + 
	black + south, data=anes, family=binomial)
restricted.mod <- glm(votedem ~ retnat + pid + age + male + educ + 
	black + south, data=anes, family=binomial)
 
 
###################################################
### chunk number 6: glmout
###################################################
#line 307 "Lecture5_2010.Rnw"
apsrtable(unrestricted.mod, restricted.mod, model.names="", Sweave=T)
 
 
###################################################
### chunk number 7: plotglm1
###################################################
#line 318 "Lecture5_2010.Rnw"
library(effects)
unres.eff <- effect("pidfac", unrestricted.mod)
res.eff <- effect("pid", restricted.mod, default.levels=7)
plot(c(1,7), c(0,1), type="n", xlab="Party ID", ylab="Pr(Vote Dem)")
polygon(x=c(1:7,7:1,1), y=plogis(c(unres.eff[["lower"]], 		
	rev(unres.eff[["upper"]]), unres.eff[["lower"]][1])), 
	col=rgb(0,1,0,.25,1), border=NA)
polygon(x=c(1:7,7:1,1), y=plogis(c(res.eff[["lower"]], 		
	rev(res.eff[["upper"]]), res.eff[["lower"]][1])), 
	col=rgb(0,0,1,.25,1), border=NA)
lines(c(1:7), plogis(res.eff[["fit"]]), lty=1, col="blue", lwd=2)
lines(c(1:7), plogis(unres.eff[["fit"]]), lty=1, col="green", lwd=2)
 
legend("topright", c("Linear","Non-linear"), fill=c(rgb(0,0,1,.25,1),
	rgb(0,1,0,.25,1)), inset=.01)
 
 
###################################################
### chunk number 8: loess
###################################################
#line 519 "Lecture5_2010.Rnw"
dat <- read.table("weakliem.txt", header=T)
out.loess <- loess(secpay ~ gini, data=dat, span=.75, degree=1,
    family="symmetric")
plot(secpay ~ gini, data=dat)
lines(out.loess$fitted[order(dat$gini)] ~
    dat$gini[order(dat$gini)])
out.loess <- loess(secpay ~ gini, data=dat, span=.75, degree=2,
    family="symmetric")
lines(out.loess$fitted[order(dat$gini)] ~
    dat$gini[order(dat$gini)], lty=2)
legend("topright", c("Degree = 1", "Degree = 2"),
    lty=c(1,2), inset=.01)
 
 
dat <- read.table("weakliem.txt", header=T)
out.loess <- loess(secpay ~ gini, data=dat, span=.6, degree=2,
    family="symmetric")
plot(secpay ~ gini, data=dat)
lines(out.loess$fitted[order(dat$gini)] ~
    dat$gini[order(dat$gini)])
out.loess <- loess(secpay ~ gini, data=dat, span=.75, degree=2,
    family="symmetric")
lines(out.loess$fitted[order(dat$gini)] ~
    dat$gini[order(dat$gini)], lty=2)
out.loess <- loess(secpay ~ gini, data=dat, span=.5, degree=2,
    family="symmetric")
lines(out.loess$fitted[order(dat$gini)] ~
    dat$gini[order(dat$gini)], lty=3)
out.loess <- loess(secpay ~ gini, data=dat, span=1, degree=2,
    family="symmetric")
lines(out.loess$fitted[order(dat$gini)] ~
    dat$gini[order(dat$gini)], lty=4)
legend("topright", c("Degree = 1", "Degree = 2"),
    lty=c(1,2), inset=.01)
 
 
## Added in class - Loess Confidence Bounds
 
dat <- read.table("weakliem.txt", header=T)
out.loess <- loess(secpay ~ gini, data=dat, span=.75, degree=2,
    family="symmetric")
 
s <- unique(dat$gini)
s <- s[order(s)]
 
preds <- predict(out.loess, newdata=data.frame(gini=s), se=T)
pred.lower <- preds$fit - 2*preds$se
pred.upper <- preds$fit + 2*preds$se
 
plot(s, preds$fit, type="l", ylim=range(c(pred.lower, pred.upper)))
lines(s, pred.lower, lty=2)
lines(s, pred.upper, lty=2)
 
###################################################
### chunk number 9: cpr1
###################################################
#line 652 "Lecture5_2010.Rnw"
data(Prestige)
Prestige.model<-lm(prestige ~ income + education + 
	women, data=Prestige)
library(car)
crPlot(Prestige.model, "income")
 
 
###################################################
### chunk number 10: cpr2
###################################################
#line 659 "Lecture5_2010.Rnw"
cr.plot(Prestige.model, "education")
 
 
###################################################
### chunk number 11: cpr3
###################################################
#line 662 "Lecture5_2010.Rnw"
cr.plot(Prestige.model, "women")
 
 
###################################################
### chunk number 12: orthpoly
###################################################
#line 814 "Lecture5_2010.Rnw"
mod <- lm(prestige ~ poly(income, 3) +
    poly(education, 3), data=Prestige)
apsrtable(mod, model.names="", Sweave=T)
Prestige$wom.dum <- as.numeric(Prestige$women > mean(Prestige$women))
mod.a <- mod <- lm(prestige ~ poly(income, 3) + wom.dum + 
 
 
## Added in class - Setting values of Other variables in Effects
poly(education, 3), data=Prestige)
pol <- poly(Duncan$education, 3)
pol2 <- poly(1:100, 3, coefs=attr(pol, "coef"))
colnames(pol2) <- NULL
 
eff2 <- effect("poly(income, 3)", mod, 
	given.values=c("wom.dum"=0, 
	"poly(education, 3)1" = pol2[16,1],
	"poly(education, 3)2" = pol2[16,2],
	"poly(education, 3)3" = pol2[16,3]))
 
pol <- poly(Duncan$education, 3)
 
mod.b<- lm(prestige ~ income + I(income^2) + I(income^3) + wom.dum + 
    poly(education, 3), data=Prestige)
 
 
 
 
###################################################
### chunk number 13: polyeff1
###################################################
#line 833 "Lecture5_2010.Rnw"
library(effects)
print(
plot(effect("poly(income, 3)", mod,
    default.levels=100, se=T))
)
 
 
###################################################
### chunk number 14: polyeff2
###################################################
#line 840 "Lecture5_2010.Rnw"
print(
plot(effect("poly(education, 2)", mod,
    default.levels=100, se=T))
)
 
 
###################################################
### chunk number 15: ornbox
###################################################
#line 917 "Lecture5_2010.Rnw"
data(Ornstein)
mod <- lm(interlocks ~ box.cox.var(interlocks+1) + assets + 
	sector + nation, data=Ornstein)
summary(mod)
 
 
## Added in class - Interpreting Nonlinear transformations of Y
 
Ornstein$new_inter <- (Ornstein$interlocks+1)^.31
mod <- lm(new_inter ~ box.cox.var(new_inter) + assets + 
	sector + nation, data=Ornstein)
summary(mod)
 
Ornstein$new_inter <- (Ornstein$interlocks+1)^.3
mod <- lm(new_inter ~ assets + 
	sector + nation, data=Ornstein)
 
eff <- effect("assets", mod)
newfit <- eff$fit^(1/.3)
newlower <- eff$lower^(1/.3)
newupper <- eff$upper^(1/.3)
plot(eff$x$assets, newfit, type="l", 
	ylim=range(c(newlower, newupper)))
points(Ornstein$assets, Ornstein$interlocks)
lines(eff$x$assets, newlower, lty=2)
lines(eff$x$assets, newupper, lty=2)
 
 
 
 
###################################################
### chunk number 16: ornboxfig
###################################################
#line 936 "Lecture5_2010.Rnw"
av.plots(mod,'box.cox.var(interlocks + 1)',
    col='black', identify.points=F)
 
 
###################################################
### chunk number 17: boxtid
###################################################
#line 981 "Lecture5_2010.Rnw"
Prestige$loginc <- log(Prestige$income)
 
box.tidwell(prestige ~ loginc + education,
    ~poly(women, 2), data=Prestige)
 
 
###################################################
### chunk number 18: avinc
###################################################
#line 996 "Lecture5_2010.Rnw"
newinc <- with(Prestige, log(income))
newed <- with(Prestige, education^2)
mod <- lm(prestige ~ newinc + newed + poly(women, 2), data=Prestige)
av.plots(mod,"newinc", identify.points=F, col="black")
 
 
###################################################
### chunk number 19: aved
###################################################
#line 1002 "Lecture5_2010.Rnw"
av.plots(mod,"newed", identify.points=F, col="black")
 
 
###################################################
### chunk number 20: incdiag1
###################################################
#line 1042 "Lecture5_2010.Rnw"
inc.mod1 <- lm(prestige ~ log(income) + education + 
    women, data=Prestige)
apsrtable(inc.mod1, model.names="", Sweave=T)
 
 
###################################################
### chunk number 21: incdiag2
###################################################
#line 1059 "Lecture5_2010.Rnw"
box.tidwell(prestige ~ income, ~ education + women, data=Prestige)
 
 
###################################################
### chunk number 22: incfix
###################################################
#line 1074 "Lecture5_2010.Rnw"
newinc <- with(Prestige, income^0.08073)
inc.mod2 <- lm(prestige ~ newinc + education +
	women, data=Prestige)
apsrtable(inc.mod2, model.names="", Sweave=T)
 
 
###################################################
### chunk number 23: eddiag1
###################################################
#line 1093 "Lecture5_2010.Rnw"
ed.mod1 <- lm(prestige ~ newinc + 
	poly(education, 3) + women, data=Prestige)
apsrtable(ed.mod1, model.names="", Sweave=T)
 
 
###################################################
### chunk number 24: edfig
###################################################
#line 1107 "Lecture5_2010.Rnw"
print(
plot(effect("poly(education, 3)", ed.mod1))
)
 
 
###################################################
### chunk number 25: wompostinc
###################################################
#line 1121 "Lecture5_2010.Rnw"
cr.plot(inc.mod2,"women", col="black")
 
 
###################################################
### chunk number 26: wom2
###################################################
#line 1137 "Lecture5_2010.Rnw"
wom.mod1 <- lm(prestige ~ newinc + education + 
	poly(women, 2), data=Prestige)
apsrtable(wom.mod1, model.names="", Sweave=T)
 
 
###################################################
### chunk number 27: womeff
###################################################
#line 1154 "Lecture5_2010.Rnw"
print(
plot(effect("poly(women, 2)", wom.mod1))
)
 
 
###################################################
### chunk number 28: boxinc2
###################################################
#line 1173 "Lecture5_2010.Rnw"
box.tidwell(prestige ~ income, ~education + 
    poly(women, 2), data=Prestige)
 
 
###################################################
### chunk number 29: finalmod
###################################################
#line 1190 "Lecture5_2010.Rnw"
final.mod <- lm(prestige ~ log(income) + education + 
	poly(women, 2), data=Prestige)
apsrtable(final.mod, model.names="", Sweave=T)
 
 
###################################################
### chunk number 30: effincfinal
###################################################
#line 1204 "Lecture5_2010.Rnw"
print(
plot(effect("log(income)", final.mod))
)
 
 
###################################################
### chunk number 31: effwomfinal
###################################################
#line 1209 "Lecture5_2010.Rnw"
print(
plot(effect("poly(women, 2)", final.mod))
)