###################################################
### chunk number 1: asfac
###################################################
#line 194 "Lecture4_2010.Rnw"
library(car)
data(Duncan)
type2 <- with(Duncan, as.factor(type))
contrasts(type2)
 
 
###################################################
### chunk number 2: relev
###################################################
#line 203 "Lecture4_2010.Rnw"
type2 <- relevel(type2, ref="bc")
contrasts(type2)
 
 
###################################################
### chunk number 3: lm1
###################################################
#line 245 "Lecture4_2010.Rnw"
mod1<-lm(prestige~income+education+type, data=Duncan)
summary(mod1)
 
 
###################################################
### chunk number 4: xtanova
###################################################
#line 260 "Lecture4_2010.Rnw"
library(xtable)
xtable(Anova(mod1))
 
 
###################################################
### chunk number 5: qv
###################################################
#line 395 "Lecture4_2010.Rnw"
library(qvcalc)
qvtype<-qvcalc(mod1,"type")
summary(qvtype)
 
 
###################################################
### chunk number 6: plotqv
###################################################
#line 409 "Lecture4_2010.Rnw"
plot(qvtype)
 
 
###################################################
### chunk number 7: plotfloat
###################################################
#line 431 "Lecture4_2010.Rnw"
library(Epi)
float1 <- float(mod1, "type")
cis <- float1$coef + 
	t(sapply(sqrt(float1$var), 
	function(x)x*1.96*c(-1,1)))
plot(1:3, float1$coef, pch=16, 
	ylim = range(c(cis)), axes=F,
	xlab = "", ylab = "Coefficients")
segments(1:3, cis[,1], 
	1:3, cis[,2])
axis(2)
axis(1, at=1:3, 
	labels=c("BC", "Prof","WC"))
box()
 
 
###################################################
### chunk number 8: ornmod
###################################################
#line 487 "Lecture4_2010.Rnw"
data(Ornstein)
omod <- lm(interlocks ~ 
	nation + sector + log2(assets), 
	data=Ornstein)
oqv <- qvcalc(omod, "sector")
pdf("ornqv.pdf", height=6, width=6)
plot(oqv)
invisible(dev.off())
library(apsrtable)
apsrtable(omod, Sweave=T)
 
 
###################################################
### chunk number 9: ornfp
###################################################
#line 512 "Lecture4_2010.Rnw"
library(factorplot)
ofp <- factorplot(
    omod, 
    "sector")
pdf("omod.pdf", 
    height=6, 
    width=6)
plot(ofp, 
	bonferroni=F)
invisible(dev.off())
 
 
###################################################
### chunk number 10: bonomod
###################################################
#line 544 "Lecture4_2010.Rnw"
pdf("omod2.pdf", 
    height=6, 
    width=6)
plot(ofp, 
	bonferroni=T)
invisible(dev.off())
 
 
###################################################
### chunk number 11: fpsum
###################################################
#line 560 "Lecture4_2010.Rnw"
summary(ofp)
 
 
###################################################
### chunk number 12: fpp
###################################################
#line 569 "Lecture4_2010.Rnw"
print(ofp, sig=T)
 
 
###################################################
### chunk number 13: fpp2
###################################################
#line 579 "Lecture4_2010.Rnw"
print(ofp, sig="bon")
 
 
###################################################
### chunk number 14: contrsum
###################################################
#line 665 "Lecture4_2010.Rnw"
contrasts(Duncan[["type"]])<-'contr.sum'
contrasts(Duncan[["type"]])
 
 
###################################################
### chunk number 15: duncanmod
###################################################
#line 836 "Lecture4_2010.Rnw"
Duncan.mod <- lm(prestige ~ income*type, data=Duncan)
 
 
###################################################
### chunk number 16: duncantab
###################################################
#line 847 "Lecture4_2010.Rnw"
library(apsrtable)
apsrtable(Duncan.mod, model.names="", Sweave=T)
 
 
###################################################
### chunk number 17: xtanduncan
###################################################
#line 870 "Lecture4_2010.Rnw"
xtable(Anova(Duncan.mod))
 
 
###################################################
### chunk number 18: inccat
###################################################
#line 892 "Lecture4_2010.Rnw"
library(car)
library(xtable)
 
Duncan$inc.cat <- cut(Duncan$income, 3)
 
mod <- lm(prestige~ inc.cat * type + education,
  data=Duncan)
 
 
###################################################
### chunk number 19: summod
###################################################
#line 907 "Lecture4_2010.Rnw"
summary(mod)
 
 
###################################################
### chunk number 20: calceff
###################################################
#line 920 "Lecture4_2010.Rnw"
eg <- expand.grid(levels(Duncan$inc.cat), 
  levels(Duncan$type))
newdat <- data.frame(
  inc.cat = eg[,1], 
  type = eg[,2], 
  education = mean(Duncan$education), 
  prestige=50)
 
preds <- predict(mod, newdata=newdat, 
  interval="confidence")
 
 
###################################################
### chunk number 21: cat2mat
###################################################
#line 941 "Lecture4_2010.Rnw"
eff.mat <- matrix(round(preds[,1],3), 
  ncol=nlevels(Duncan$inc.cat), 
  nrow = nlevels(Duncan$type))
rownames(eff.mat) <- levels(Duncan$inc.cat)
colnames(eff.mat) <- levels(Duncan$type)
 
cis <- apply(round(preds[,2:3], 3), 1, 
  function(x)paste("(", x[1], ", ", 
    x[2], ")", sep=""))
 
cis.mat <- matrix(cis, ncol=nlevels(Duncan$type))
out.mat <- rbind(eff.mat[1,], 
  cis.mat[1,], 
eff.mat[2,], cis.mat[2,], eff.mat[3,], 
  cis.mat[3,])
 
levs <- levels(Duncan$inc.cat)
rownames(out.mat) <- c(levs[1], " ", 
  levs[2], "  ", levs[3], "   ")
 
 
###################################################
### chunk number 22: ltable
###################################################
#line 972 "Lecture4_2010.Rnw"
xtable(out.mat, align="lccc")
 
 
###################################################
### chunk number 23: catcon
###################################################
#line 1005 "Lecture4_2010.Rnw"
mod <- lm(prestige ~ income*type + education, 
  data=Prestige)
 
 
###################################################
### chunk number 24: summod2
###################################################
#line 1015 "Lecture4_2010.Rnw"
summary(mod)
 
 
###################################################
### chunk number 25: anova
###################################################
#line 1026 "Lecture4_2010.Rnw"
Anova(mod)
 
 
###################################################
### chunk number 26: condinc
###################################################
#line 1037 "Lecture4_2010.Rnw"
A <- cbind(c(0,1,0,0,0,0,0), 
           c(0,1,0,0,0,1,0),
           c(0,1,0,0,0,0,1))
b <- coef(mod)
v <- vcov(mod)
 
cond.eff <- t(A) %*% b
cond.var <- t(A) %*% v %*% A 
cond.se <- sqrt(diag(cond.var))
cond.t <- cond.eff/cond.se
cond.p <- 2*pt(cond.t, mod$df.residual, 
  lower.tail=F)
mat <- cbind(cond.eff, cond.se, cond.t, cond.p)
colnames(mat) <- c("Effect", "SE", "t", "p-value")
rownames(mat) <- c("Blue Collar", "Professional", 
  "White Collar")
 
 
###################################################
### chunk number 27: xtab1
###################################################
#line 1063 "Lecture4_2010.Rnw"
xtable(mat, align = "lcccc", digits=3)
 
 
###################################################
### chunk number 28: condtype
###################################################
#line 1074 "Lecture4_2010.Rnw"
s <- seq(min(Prestige$income), 
  max(Prestige$income), length=25)
A1 <- rbind(0,0,1,0,0,s,0)
A2 <- rbind(0,0,0,1,0,0,s)
A3 <- rbind(0,0,1,-1,0,s,-s)
 
bcprof.eff <- t(A1) %*% b
bcwc.eff <- t(A2) %*% b
profwc.eff <- t(A3) %*% b
 
bcprof.se <- sqrt(diag(t(A1) %*% v %*% A1))
bcwc.se <- sqrt(diag(t(A2) %*% v %*% A2))
profwc.se <- sqrt(diag(t(A3) %*% v %*% A3))
 
 
###################################################
### chunk number 29: prepdata
###################################################
#line 1097 "Lecture4_2010.Rnw"
 
plot.dat <- data.frame(
  income = rep(s,3), 
  eff = c(bcprof.eff, bcwc.eff, profwc.eff), 
  se = c(bcprof.se, bcwc.se, profwc.se), 
  comp = rep(1:3, each=25))
 
plot.dat$lower <- plot.dat$eff - qt(.975, 
  mod$df.residual, lower.tail=F)*plot.dat$se
plot.dat$upper <- plot.dat$eff + qt(.975, 
  mod$df.residual, lower.tail=F)*plot.dat$se
 
plot.dat$comp <- factor(plot.dat$comp, 
  levels=1:3, labels=c("PROF-BC", "WC-BC", "PROF-WC"))
 
 
###################################################
### chunk number 30: makegraph
###################################################
#line 1120 "Lecture4_2010.Rnw"
library(lattice)
print(
xyplot(eff ~ income | comp, data=plot.dat, 
  xlab = "Income", ylab = "Conditional Effect", 
  as.table=T, 
  ylim=range(c(plot.dat$lower, plot.dat$upper)), 
  panel = function(x,y,subscripts){
    panel.lines(x,y, col="black")
    panel.lines(x, plot.dat$lower[subscripts], 
      lty=2, col="black")
    panel.lines(x, plot.dat$upper[subscripts], 
      lty=2, col="black")
    panel.abline(h=0, col="gray65")}, 
  key = list(space = "top", 
    lines = list(col=c("black", "black"), lty=c(1,2)), 
    text= list(c("Effect", "95% CI"))))
)
 
 
###################################################
### chunk number 31: quant2
###################################################
#line 1172 "Lecture4_2010.Rnw"
mod <- lm(prestige ~ income*education + type, data=Prestige)
summary(mod)
 
 
###################################################
### chunk number 32: DAMplot
###################################################
#line 1205 "Lecture4_2010.Rnw"
 
library(DAMisc)
DAintfun2(mod, c("income", "education"), plot.type="pdf")
 
 
###################################################
### chunk number 33: intfun2
###################################################
#line 1260 "Lecture4_2010.Rnw"
DAintfun(mod, c("income","education"), theta=-45, phi=20)
 
 
###################################################
### chunk number 34: cent_data
###################################################
#line 1324 "Lecture4_2010.Rnw"
set.seed(123)
X <- mvrnorm(1000, c(10,10), matrix(c(1,.4,.4,1), ncol=2), empirical=T)
X <- cbind(1, X, apply(X, 1, prod))
b <- matrix(c(2,3,-4,3), ncol=1)
Y <- X %*% b + rnorm(1000, 0, 4)
df <- data.frame(Y=Y, x1=X[,2], x2=X[,3])
mod1 <- lm(Y ~ x1*x2, data=df)
df2 <- as.data.frame(apply(df, 2, function(x)x-mean(x)))
mod2 <- lm(Y ~ x1*x2, data=df2)
 
 
###################################################
### chunk number 35: tabcent
###################################################
#line 1337 "Lecture4_2010.Rnw"
apsrtable(mod1, mod2, model.names=c("Not Cent","Cent"), Sweave=T)
 
 
###################################################
### chunk number 36: viftab
###################################################
#line 1342 "Lecture4_2010.Rnw"
tab <- cbind(vif(mod1), vif(mod2))
colnames(tab) <- c("No Cent","Cent")
xtable(tab, caption="VIF Statistics")
 
 
###################################################
### chunk number 37: pctiletab
###################################################
#line 1358 "Lecture4_2010.Rnw"
library(xtable)
tab <- cbind(
	with(df, quantile(x1, probs=c(.25,.5,.75))), 
	with(df2, quantile(x1, probs=c(.25,.5,.75))), 
	with(df, quantile(x2, probs=c(.25,.5,.75))), 
	with(df2, quantile(x2, probs=c(.25,.5,.75))))
rownames(tab) <- c("25th","50th","75th")
colnames(tab) <- c("x1","x1 (cent)","x2","x2 (cent)")
xtable(tab)
 
 
###################################################
### chunk number 38: condeff
###################################################
#line 1373 "Lecture4_2010.Rnw"
b1 <- coef(mod1)
b2 <- coef(mod2)
v1 <- vcov(mod1)
v2 <- vcov(mod2)
eff1.2 <- round(cbind(1, tab[,1]) %*% b1[c(3,4)], 3)
eff2.2 <- round(cbind(1, tab[,2]) %*% b2[c(3,4)], 3)
eff1.1 <- round(cbind(1, tab[,3]) %*% b1[c(2,4)], 3)
eff2.1 <- round(cbind(1, tab[,4]) %*% b2[c(2,4)], 3)
 
se1.2 <- round(sqrt(diag(cbind(1, tab[,1]) %*% v1[c(3,4), c(3,4)] %*% t(cbind(1, tab[,1])))), 3)
se2.2 <- round(sqrt(diag(cbind(1, tab[,2]) %*% v2[c(3,4), c(3,4)] %*% t(cbind(1, tab[,2])))), 3)
se1.1 <- round(sqrt(diag(cbind(1, tab[,3]) %*% v1[c(2,4), c(2,4)] %*% t(cbind(1, tab[,3])))), 3)
se2.1 <- round(sqrt(diag(cbind(1, tab[,4]) %*% v2[c(2,4), c(2,4)] %*% t(cbind(1, tab[,4])))), 3)
 
tab2 <- cbind(c(rbind(c(eff1.1), c(se1.1))), c(rbind(c(eff2.1), c(se2.1))),c(rbind(c(eff1.2), c(se1.2))), c(rbind(c(eff2.2), c(se2.2))))
 
rownames(tab2) <- c("25th","","50th"," ","75th","  ")
colnames(tab2) <- c("eff x1",  "eff x1 (cent)", "eff x2", "eff x2 (cent)")
xtable(tab2, caption = "Conditional Effects of x1 and x2")
 
 
###################################################
### chunk number 39: modrelimp
###################################################
#line 1528 "Lecture4_2010.Rnw"
mod1<-lm(prestige~income+education+type, data=Duncan)
apsrtable(mod1, model.names="", Sweave=T)
 
 
###################################################
### chunk number 40: relimpsum
###################################################
#line 1539 "Lecture4_2010.Rnw"
library(relimp)
relimp(mod1, set1=2, set2=4:5)