##------------------------------------------------------------------------## ## Script for Chapter 7, An R and S-PLUS Companion to Applied Regression ## ## John Fox ## ## Sage Publications, 2002 ## ##------------------------------------------------------------------------## # preliminaries options(width=65) options(digits=5) # Basics args(plot.default) plot(c(0,1), c(0,1), type='n', xlab='', ylab='') par('col') par() plot(1:25, xlab='Symbol Number', ylab='', type='n') for (pch in 1:25) points(pch, pch, pch=pch) lines(1:25, type='h', lty=2) plot(1:26, xlab='letters', ylab='', type='n', axes=F) box() for (letter in 1:26) points(letter, 27 - letter, pch=letters[letter]) ## R alternative plot(1:25, pch=1:25, xlab='Symbol Number', ylab='') lines(1:25, type='h', lty='dashed') plot(26:1, xlab='letters', ylab='', pch=letters, axes=F, frame=T) plot(c(1,7), c(0,1), type='n', axes=F, xlab='Line Type (lty)', ylab='') box() axis(1, at=1:6) for (lty in 1:6) lines(c(lty, lty, lty + 1), c(0, 0.5, 1), lty=lty) ## R alternative: plot(c(1,7), c(0,1), type='n', axes=F, frame=T, xlab='Line Type (lty)', ylab='') axis(1, at=1:6) for (lty in 1:6) lines(c(lty, lty, lty + 1), c(0, 0.5, 1), lty=lty) # text par(mfrow=c(1,2)) plot(c(0,1), c(0,1), axes=F, type='n', xlab='', ylab='') box() text(x=c(.2, .5), y=c(.2, .7), c('example text', 'another string')) title('(a)') plot(c(0,1), c(0,1), axes=F, type='n', xlab='', ylab='') box() text(locator(3), c('one','two','three')) title('(b)') locator() # arrows plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='') arrows(x0=1:5, y0=rep(0.1, 5), x1=1:5, y1=seq(0.3, 0.9, len=5), code=3) title('(a) arrows') plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='') segments(x0=1:5, y0=rep(0.1, 5), x1=1:5, y1=seq(0.3, 0.9, len=5)) title('(b) segments') ## S-PLUS alternative: plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='') arrows(x1=1:5, y1=rep(0.1, 5), x2=1:5, y2=seq(0.3, 0.9, len=5)) title('(a) arrows') plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='') segments(x1=1:5, y1=rep(0.1, 5), x2=1:5, y2=seq(0.3, 0.9, len=5)) title('(b) segments') # polygons plot(c(0,1), c(0,1), type='n', xlab='', ylab='') polygon(c(.2,.8,.8), c(.2,.2,.8), col=1) polygon(c(.2,.2,.8), c(.2,.8,.8)) # legend plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='', frame.plot=T) legend(locator(1), legend=c('group A', 'group B', 'group C'), lty=1:3, pch=1:3) # S-Plus alternative: plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='') legend(locator(1), legend=c('group A', 'group B', 'group C'), lty=1:3, marks=1:3) # colors piechart(rep(1, length(palette())), col=palette()) palette() rainbow(10) gray(0:8/8) piechart(rep(1,100), col=rainbow(100), labels=rep('',100)) piechart(rep(1,100), col=gray(0:100/100), labels=rep('',100)) # S-PLUS alternative: piechart(rep(1,16), col=1:16) # Extended example: effect plots library(car) data(Cowles) dim(Cowles) Cowles[sort(sample(1421, 10)),] mod.cowles <- glm(volunteer ~ neuroticism * extraversion + sex, data=Cowles, family=binomial) summary(mod.cowles) # show full response surface neuroticism <- 0:24 extraversion <- 0:24 sex <- c('male', 'female') graph.data <- expand.grid(neuroticism=neuroticism, extraversion=extraversion, sex=sex) graph.data\$fit <- predict(mod.cowles, newdata=graph.data, type='response') graph.data par(mfrow=c(2,1)) prob <- matrix(graph.data\$fit[graph.data\$sex=='male'], 25, 25) persp(neuroticism, extraversion, prob, phi=30, theta=45, expand=0.65, d=2, shade=0.75, ticktype='detailed', zlab='Probability(Volunteer)', main='Males') prob <- matrix(graph.data\$fit[graph.data\$sex=='female'], 25, 25) persp(neuroticism, extraversion, prob, phi=30, theta=45, expand=0.65, d=2, shade=0.75, ticktype='detailed', zlab='Probability(Volunteer)', main='Females') # make effect plot extraversion <- 0:24 neuroticism <- seq(0, 24, by=6) graph.data <- expand.grid(neuroticism=neuroticism, extraversion=extraversion) X <- cbind(constant=1, as.matrix(graph.data), sex=0.5, neuro.extra=graph.data\$neuroticism * graph.data\$extraversion) X logit <- X %*% coefficients(mod.cowles) se <- sqrt(diag(X %*% Var(mod.cowles) %*% t(X))) prob <- matrix(1/(1+exp(-logit)), 5, 25,) low <- matrix(1/(1+exp(-(logit - se))), 5, 25) high <- matrix(1/(1+exp(-(logit + se))), 5, 25) prob plot(range(extraversion), range(c(low, high + 0.05)), type='n', xlab='Introversion-Extraversion', ylab='Probability of Volunteering', xaxt='n', xlim=c(0,30)) axis(1, at=seq(0, 24, by=6)) for (neuro in 1:5){ lines(extraversion, prob[neuro,]) text(25, prob[neuro,25], paste('N = ', neuroticism[neuro]), adj=0) } extra <- seq(1, 25, by=6) for (neuro in c(1, 3, 5)){ arrows(extraversion[extra], low[neuro, extra], extraversion[extra], high[neuro, extra], angle=90, code=3, lty=2, length=0.05) } ## S-PLUS alternative: for (neuro in c(1, 3, 5)){ segments(extraversion[extra], low[neuro, extra], extraversion[extra], high[neuro, extra]) } text(locator(1), 'Stability-Neuroticism', adj=1, cex=0.75)