### R code from vignette source '/home/jimrc/classes/stat505/notes/ARMChapter5/chapter5-C-Notes.Rnw' ################################################### ### code chunk number 1: fit8 ################################################### require(arm) options(width=55,digits=3) attach(wells) wells.fit8 <- glm (switch ~ (c.dist100 + c.arsenic + c.educ4)^2, family=binomial(link="logit")) pred.8 <- fitted(wells.fit8) display (wells.fit8) binned.resids <- function (x, y, nclass=sqrt(length(x))){ shinglex <- co.intervals(x, number=nclass, overlap=0) break.x <- cut(x, c(shinglex[,1], shinglex[nrow(shinglex),2])) n <- tapply(x, break.x, length) twoSE <- 2*tapply(y, break.x, sd)/sqrt(n) output <- cbind( tapply(x, break.x, mean), tapply(y, break.x, mean), n, shinglex[,1], shinglex[,2], twoSE ) colnames(output) <- c ("xbar", "ybar", "n", "x.lo", "x.hi", "twoSE") output } ################################################### ### code chunk number 2: fit8plot ################################################### ## Residual Plot (Figure 5.13 (a)) par(mfrow=c(1,2)) plot( pred.8, switch - pred.8, ylab="Observed - estimated", xlab="P(switch)",main="Raw Residual Plot") br.8 <- binned.resids(pred.8, switch-pred.8, nclass=40) plot(range(br.8[,1]), range(br.8[,2], br.8[,6],-br.8[,6]), xlab="Estimated Pr (switching)", ylab="Average residual", type="n", main="Binned residual plot", mgp=c(2,.5,0)) abline (0,0, col="gray", lwd=.5) lines (br.8[,1], br.8[,6], col="gray", lwd=.5) lines (br.8[,1], -br.8[,6], col="gray", lwd=.5) points (br.8[,1], br.8[,2], pch=19, cex=.5) ################################################### ### code chunk number 3: dist-As-Fit ################################################### par(mfrow=c(1,2)) br.dist <- binned.resids (dist, switch-pred.8, nclass=40) plot(range(br.dist[,1]), range(br.dist[,2],br.dist[,6],-br.dist[,6]), xlab="Distance to nearest safe well", ylab="Average residual", type="n", main="Binned residual plot", mgp=c(2,.5,0)) abline (0,0, col="gray", lwd=.5) lines (br.dist[,1], br.dist[,6], col="gray", lwd=.5) lines (br.dist[,1], -br.dist[,6], col="gray", lwd=.5) points (br.dist[,1], br.dist[,2], pch=19, cex=.5) # arsenic (Figure 5.13 (b)) br.arsenic <- binned.resids (arsenic, switch-pred.8, nclass=40) plot(range(0,br.arsenic[,1]), range(br.arsenic[,2],br.arsenic[,6],-br.arsenic[,6]), xlab="Arsenic level", ylab="Average residual", type="n", main="Binned residual plot", mgp=c(2,.5,0)) abline (0,0, col="gray", lwd=.5) lines (br.arsenic[,1], br.arsenic[,6], col="gray", lwd=.5) lines (br.arsenic[,1], -br.arsenic[,6], col="gray", lwd=.5) points (br.arsenic[,1], br.arsenic[,2], pch=19, cex=.5) ################################################### ### code chunk number 4: fit9 ################################################### c.log.arsenic <- log(arsenic) - mean (log(arsenic)) wells.fit9 <- update( wells.fit8, . ~ (c.dist100 + c.log.arsenic + c.educ4)^2) display (wells.fit9) ################################################### ### code chunk number 5: fit9plots ################################################### binary.jitter <- function(a, jitt=.05){ jitter <- runif (length(a), 0, jitt) a + (a==0)*jitter - (a==1)* jitter } par(mfrow=c(1,2)) plot(arsenic,binary.jitter(switch), xlim=c(0,10), xlab="Arsenic concentration in well water", ylab="Pr (switching)", xaxs="i", yaxs="i", mgp=c(2,.5,0), pch=20, cex=.1) ## dist = 48 so c.dist=0 curve (invlogit(cbind(1, 0, log(x)-.314, 0, 0*(log(x)-.314),0,0) %*% coef(wells.fit9)), from=.50, lwd=.5, add=TRUE) curve (invlogit(cbind(1, -.4833, log(x)-.314, 0, -.4833*(log(x)-.314),0,0) %*% coef(wells.fit9)), from=.50, lwd=.5, add=TRUE) text (1.2, .8, "if dist = 0", adj=0, cex=.8) text (1.8, .6, "if dist = 48.3", adj=0, cex=.8) ## Graph of binned residuals for log model fit.9 (Figure 5.15 (b)) pred.9 <- fitted(wells.fit9) br.fit.9 <- binned.resids (arsenic, switch-pred.9, nclass=40) plot(range(0,br.fit.9[,1]), range(br.fit.9[,2],br.fit.9[,6],-br.fit.9[,6]), xlab="Arsenic level", ylab="Average residual", type="n", main="Binned residual plot\n for model with log (arsenic)", mgp=c(2,.5,0)) abline (0,0, col="gray", lwd=.5) lines (br.fit.9[,1], br.fit.9[,6], col="gray", lwd=.5) lines (br.fit.9[,1], -br.fit.9[,6], col="gray", lwd=.5) points (br.fit.9[,1], br.fit.9[,2], pch=19, cex=.5) ################################################### ### code chunk number 6: errorRate ################################################### errorRate <- function(y,pred,c){ sum( y != (pred >= c))/length(y) } c( errorRate(switch, pred.9,.5), errorRate(switch, 1, .5)) ################################################### ### code chunk number 7: anova ################################################### anova(wells.fit5,update(wells.fit6,.~.-assoc), wells.fit6) ################################################### ### code chunk number 8: delta1 ################################################### wells.fit10 <- update(wells.fit3, .~ . + educ4) display(wells.fit10) delta <- invlogit( cbind(1, 1, arsenic, educ4) %*% coef(wells.fit10)) - invlogit( cbind(1, 0, arsenic, educ4) %*% coef(wells.fit10)) c(mean(delta), sd(delta)) ################################################### ### code chunk number 9: delta2 ################################################### delta <- invlogit( cbind(1, dist/100, 1, educ4) %*% coef(wells.fit10)) - invlogit( cbind(1, dist/100, .5, educ4) %*% coef(wells.fit10)) c(mean(delta), sd(delta)) delta <- invlogit( cbind(1, dist/100, arsenic,3 ) %*% coef(wells.fit10)) - invlogit( cbind(1, dist/100, arsenic, 0) %*% coef(wells.fit10)) c(mean(delta), sd(delta)) ################################################### ### code chunk number 10: delta3 ################################################### wells.fit11 <- update(wells.fit10, .~.+I(dist/100):arsenic) delta <- invlogit( cbind(1, 1, arsenic, educ4, 1*arsenic) %*% coef(wells.fit11)) - invlogit( cbind(1,0 ,arsenic , educ4,0*arsenic) %*% coef(wells.fit11)) c(mean(delta), sd(delta)) ################################################### ### code chunk number 11: badplot ################################################### par(mfrow=c(1,1)) x <- rnorm(40) y <- as.numeric(x < .5) plot(x,y) abline(v=.5,col=2) curve(invlogit(-500*(x-.502)),-3,2,add=TRUE,col=3)