## ----setopts,echo=FALSE--------------------------------------------- library("knitr") opts_chunk$set(fig.width=5,fig.height=5,#tidy=TRUE, out.width="0.8\\textwidth",echo=TRUE) #options(prompt=" ") options(continue="+ ", prompt="R> ", width=70) options(show.signif.stars=FALSE, scipen=3) ## ----setup, include=FALSE, cache=FALSE, results='hide', echo=FALSE------------ library(car) library(effects) render_sweave() options(width=80, digits=5, str=list(strict.width="cut")) strOptions(strict.width="cut") ## ----------------------------------------------------------------------------- library("car") # also loads the carData package Prestige$type <- factor(Prestige$type, levels=c("bc", "wc", "prof")) lm1 <- lm(prestige ~ education + poly(women, 2) + log(income)*type, data=Prestige) ## ----------------------------------------------------------------------------- S(lm1) ## ----fig11,include=TRUE,fig.width=5,fig.height=4,fig.show='hide'-------------- library("effects") e1.lm1 <- predictorEffect("education", lm1) plot(e1.lm1) ## ----------------------------------------------------------------------------- brief(e1.lm1$model.matrix) ## ----------------------------------------------------------------------------- e1a.lm1 <- predictorEffect("education", lm1, focal.levels=5) e1a.lm1 summary(e1a.lm1) as.data.frame(e1a.lm1) ## ----------------------------------------------------------------------------- e2.lm1 <- predictorEffect("income", lm1, focal.levels=5) as.data.frame(e2.lm1) ## ----fig12,include=TRUE,fig.width=5,fig.height=5,fig.show='hide'-------------- plot(predictorEffect("income", lm1), lines=list(multiline=TRUE)) ## ----fig13,include=TRUE,fig.width=5,fig.height=5,fig.show='hide'-------------- plot(predictorEffect("type", lm1), lines=list(multiline=TRUE)) ## ----fig14,include=TRUE,fig.width=7,fig.height=8,fig.show='hide'-------------- eall.lm1 <- predictorEffects(lm1) plot(eall.lm1) ## ----eval=FALSE--------------------------------------------------------------- # plot(eall.lm1) # plot(predictorEffects(lm1)) # plot(predictorEffects(lm1, ~ income + education + women + type)) ## ----eval=FALSE--------------------------------------------------------------- # plot(predictorEffects(lm1, ~ type + education)) ## ----eval=FALSE--------------------------------------------------------------- # plot(predictorEffects(lm1, ~ women)) # plot(predictorEffects(lm1)[[2]]) # plot(predictorEffect("women", lm1)) ## ----fig21a,include=TRUE,fig.width=5,fig.height=4.5,fig.show='hide'----------- e3.lm1 <- predictorEffect("type", lm1) plot(e3.lm1, lines=list(multiline=TRUE)) ## ----fig21b,include=TRUE,fig.width=6,fig.height=5,fig.show='hide'------------- plot(e3.lm1, lines=list(multiline=FALSE)) # the default ## ----fig22a,include=TRUE,fig.width=5,fig.height=4.5,fig.show='hide'----------- e3.lm1 <- predictorEffect("type", lm1, xlevels=list(income=c(5000, 15000, 25000))) plot(e3.lm1, lines=list(multiline=TRUE), confint=list(style="bars")) ## ----fig22b,include=TRUE,fig.width=5.5,fig.height=5,fig.show='hide'----------- plot(e3.lm1, lines=list(multiline=FALSE), # the default lattice=list(layout=c(3, 1))) ## ----fig23,include=TRUE,fig.width=5,fig.height=4,fig.show='hide'-------------- e4.lm1 <- predictorEffect("education", lm1, se=list(type="scheffe", level=.99), vcov.=hccm) plot(e4.lm1) ## ----------------------------------------------------------------------------- lm2 <- lm(log(prestige) ~ log(income) + education + type, Prestige) ## ----fig30,include=TRUE,fig.width=5,fig.height=4,fig.show='hide'-------------- plot(predictorEffects(lm2, ~ income)) ## ----fig31,include=TRUE,fig.width=5,fig.height=4,fig.show='hide'-------------- plot(predictorEffects(lm2, ~ income), axes=list( x=list(income=list(transform=list(trans=log, inverse=exp))) )) ## ----fig32,include=TRUE,fig.width=5,fig.height=5,fig.show='hide'-------------- plot(predictorEffects(lm2, ~ income), main="Transformed Plot", axes=list( grid=TRUE, x=list(rotate=30, rug=FALSE, income=list(transform=list(trans=log, inverse=exp), lab="income, log-scale", ticks=list(at=c(2000, 5000, 10000, 20000)), lim=c(1900, 21000)) ))) ## ----fig33,include=TRUE,fig.width=4,fig.height=4,fig.show='hide'-------------- # default: plot(predictorEffects(lm2, ~ education), main="Default log(prestige)") # Change only tick-mark labels to arithmetic scale: plot(predictorEffects(lm2, ~ education), main="log(prestige), Arithmetic Ticks", axes=list(y=list(transform=list(trans=log, inverse=exp), lab="prestige", type="rescale"))) # Replace log(presige) by prestige: plot(predictorEffects(lm2, ~ education), main="Prestige in Arithmethic Scale", axes=list(y=list(transform=exp, lab="prestige"))) ## ----------------------------------------------------------------------------- library("lme4") # for lmer() Blackmore$tran.exercise <- bcnPower(Blackmore$exercise, lambda=0.25, gamma=0.1) mm1 <- lmer(tran.exercise ~ I(age - 8)*group + (I(age - 8) | subject), data=Blackmore) ## ----fig33a,include=TRUE,fig.width=5,fig.height=5,fig.show='hide'------------- e1.mm1 <- predictorEffect("age", mm1) plot(e1.mm1, lines=list(multiline=TRUE), confint=list(style="auto")) ## ----fig33b,include=TRUE,fig.width=5,fig.height=5,fig.show='hide'------------- f.trans <- function(x) bcnPower(x, lambda=0.25, gamma=0.1) f.inverse <- function(x) bcnPowerInverse(x, lambda=0.25, gamma=0.1) plot(e1.mm1, lines=list(multiline=TRUE), confint=list(style="auto"), axes=list(x=list(age=list(lab="Age (years)")), y=list(transform=list(trans=f.trans, inverse=f.inverse), type="response", lab="Exercise (hours/week)")), lattice=list(key.args=list(x=.20, y=.75, corner=c(0, 0), padding.text=1.25)), main="" ) ## ----------------------------------------------------------------------------- data("Blowdown", package="alr4") gm1 <- glm(y ~ log(d) + s + spp, family=binomial, data=Blowdown) ## ----fig34,include=TRUE,fig.width=6.5,fig.height=6.5,fig.show='hide'---------- plot(predictorEffects(gm1), axes=list(grid=TRUE, x=list(rug=FALSE, rotate=35))) ## ----fig35,include=TRUE,fig.width=3.5,fig.height=3.5,fig.show='hide'---------- e1.gm1 <- predictorEffect("spp", gm1) plot(e1.gm1, main="type='rescale'", axes=list(y=list(type="rescale", lab="logit scale, probability labels"), x=list(rotate=30), grid=TRUE)) plot(e1.gm1, main="type='link'", axes=list(y=list(type="link", lab="logit scale, logit labels"), x=list(rotate=30), grid=TRUE)) plot(e1.gm1, main="type='response'", axes=list(y=list(type="response", grid=TRUE, lab="probabilty scale, probability labels"), x=list(rotate=30), grid=TRUE)) ## ----fig36,include=TRUE,fig.width=5.5,fig.height=4.5,fig.show='hide'---------- or <- order(as.data.frame(e1.gm1)$fit) # order smallest to largest Blowdown$spp1 <- factor(Blowdown$spp, # reorder levels of spp levels=levels(Blowdown$spp)[or]) gm2 <- update(gm1, ~ . - spp + spp1) # refit model plot(predictorEffects(gm2, ~ spp1), main="type='response', ordered", axes=list(y=list(type="response", lab="probabilty scale, probability labels"), x=list(rotate=30, spp=list(lab="Species")), grid=TRUE)) ## ----fig37,include=TRUE,fig.width=9,fig.height=12,fig.show='hide'------------- gm3 <- update(gm2, ~ . + s:log(d)) # add an interaction plot(predictorEffects(gm3, ~ s + d), axes=list(x=list(rug=FALSE, rotate=90), y=list(type="response", lab="Blowdown Probability")), lattice=list(layout=c(1, 5))) ## ----fig38,include=TRUE,fig.width=9,fig.height=5,fig.show='hide'-------------- plot(predictorEffects(gm3, ~ s + d, xlevels=list(d=c(5, 40, 80), s=c(0.1, 0.5, 0.9))), axes=list(grid=TRUE, x=list(rug=FALSE), y=list(type="response", lab="Blowdown probability")), lines=list(multiline=TRUE)) ## ----fig39,include=TRUE,fig.width=7,fig.height=7,fig.show='hide'-------------- gm4 <- update(gm3, ~ . + spp:log(d)) plot(predictorEffects(gm4, ~ d, xlevels=list(s=c(0.1, 0.5, 0.9))), axes=list(grid=TRUE, y=list(type="response"), x=list(rug=FALSE)), lines=list(multiline=TRUE)) ## ----fig310,include=TRUE,fig.width=7,fig.height=5,fig.show='hide'------------- plot(predictorEffects(gm4, ~ d, xlevels=list(s=c(0.1, 0.5, 0.9))), axes=list(grid=TRUE, y=list(type="response"), x=list(rug=FALSE)), lines=list(multiline=TRUE, z.var="spp", lty=1:9), lattice=list(layout=c(3, 1))) ## ----fig311,include=TRUE,fig.width=5.5,fig.height=5.5,fig.show='hide'--------- plot(predictorEffects(gm3, ~ d, xlevels=list(s=c(0.1, 0.5, 0.9))), axes=list(grid=TRUE, x=list(rug=FALSE), y=list(type="response")), lines=list(multiline=TRUE), confint=list(style="auto")) ## ----fig312,include=TRUE,fig.width=7,fig.height=6,fig.show='hide'------------- gm5 <- update(gm2, ~ . + spp:s) plot(predictorEffects(gm5, ~ spp, xlevels=list(s=c(0.1, 0.5, 0.9))), axes=list(grid=TRUE, y=list(type="response"), x=list(rug=FALSE, rotate=30)), lines=list(multiline=TRUE), confint=list(style="auto")) ## ----fig314,include=TRUE,fig.width=8,fig.height=6,fig.show='hide'------------- plot(predictorEffects(gm5, ~ spp, xlevels=list(s=c(0.1, 0.5, 0.9))), rug=FALSE, axes=list(grid=TRUE, y=list(type="response"), x=list(rotate=30)), lines=list(multiline=TRUE), confint=list(style="auto"), lattice=list(key.args=list(space="right", columns=1, border=TRUE, fontfamily="serif", cex=1.25, cex.title=1.5))) ## ----fig313,include=TRUE,fig.width=13,fig.height=5.5,fig.show='hide'---------- plot(predictorEffects(gm3, ~ s + d, xlevels=list(s=6, d=6)), axes=list(x=list(rug=FALSE, rotate=90), y=list(ticks=list(at=c(.999, .99, .95, .8, .5, .2, .05)))), lattice=list(layout=c(3, 2))) ## ----fig313b,include=TRUE,fig.width=6,fig.height=10,fig.show='hide'----------- plot(predictorEffect("s", gm3, xlevels=list(d=6)), axes=list(x=list(rug=FALSE, rotate=90), y=list(ticks=list(at=c(.999, .99, .95, .8, .5, .2, .05)))), lattice=list(layout=c(3, 2), array=list(row=1, col=1, nrow=2, ncol=1, more=TRUE))) plot(predictorEffect("d", gm3, xlevels=list(s=6)), axes=list(x=list(rug=FALSE, rotate=90), y=list(ticks=list(at=c(.999, .99, .95, .8, .5, .2, .05)))), lattice=list(layout=c(3, 2), array=list(row=2, col=1, nrow=2, ncol=1, more=FALSE))) ## ----fig316,include=TRUE,fig.width=7,fig.height=5,fig.show='hide'------------- plot(predictorEffects(gm4, ~ d, xlevels=list(s=c(0.1, 0.5, 0.9))), axes=list(grid=TRUE, x=list(rug=FALSE), y=list(type="response")), lines=list(multiline=TRUE, z.var="spp", lty=1:9), lattice=list(layout=c(3, 1), strip=list(factor.names=TRUE, values=TRUE, cex=1.5))) ## ----fig315,include=TRUE,fig.width=7,fig.height=6,fig.show='hide'------------- gm5 <- update(gm2, ~ . + spp:s) plot(predictorEffects(gm5, ~ spp, xlevels=list(s=c(0.1, 0.5, 0.9))), symbols=list(pch=15:17, cex=1.5), axes=list(grid=TRUE, y=list(type="response"), x=list(rotate=30)), lines=list(multiline=TRUE), confint=list(style="auto"), lattice=list(key.args=list(cex=1.5, cex.title=1.5))) ## ----fig51,include=TRUE,fig.width=10,fig.height=9,fig.show='hide'------------- lm5 <- lm(prestige ~ log(income) + education + women + type, Prestige) plot(predictorEffects(lm5, residuals=TRUE), axes=list(grid=TRUE, x=list(rotate=30)), partial.residuals=list(smooth=TRUE, span=0.75, lty="dashed")) ## ----fig52,include=TRUE,fig.width=10,fig.height=5,fig.show='hide'------------- options(scipen=10) # suppress scientific notation lm6 <- lm(infantMortality ~ group*ppgdp, data=UN) plot(predictorEffects(lm6, ~ ppgdp, partial.residuals=TRUE), axes=list(x=list(rotate=25), y=list(lim=c(0, 150))), id=list(n=1), lattice=list(layout=c(3, 1))) ## ----fig53,include=TRUE,fig.width=10,fig.height=5,fig.show='hide'------------- lm7 <- lm(log(infantMortality) ~ group*log(ppgdp), data=UN) plot(predictorEffects(lm7, ~ ppgdp, partial.residuals=TRUE), axes=list(x=list(rotate=25)), id=list(n=1), lattice=list(layout=c(3, 1))) ## ----fig54,include=TRUE,fig.width=10,fig.height=5,fig.show='hide'------------- plot(predictorEffects(lm7, ~ ppgdp, partial.residuals=TRUE), axes=list(x=list(rotate=25), y=list(transform=list(trans=log, inverse=exp), type="response", lab="Infant Mortality")), id=list(n=1), lattice=list(layout=c(3, 1))) ## ----------------------------------------------------------------------------- S(lm2) ## ----fig55,include=TRUE,fig.width=8,fig.height=4,fig.show='hide'-------------- plot(Effect(c("income", "type"), lm2, residuals=TRUE), axes=list(x=list(rotate=30)), partial.residuals=list(span=0.9), layout=c(3, 1)) ## ----------------------------------------------------------------------------- library("MASS") # for polr() Womenlf$partic <- factor(Womenlf$partic, levels=c("not.work", "parttime", "fulltime")) # order response levels or1 <- polr(partic ~ log(hincome) + children, data=Womenlf) S(or1) ## ----fig41,include=TRUE,fig.width=6.5,fig.height=6.5,fig.show='hide'---------- plot(predictorEffects(or1), axes=list(grid=TRUE), lattice=list(key.args=list(columns=1))) ## ----fig62,include=TRUE,fig.width=6,fig.height=4,fig.show='hide'-------------- plot(predictorEffects(or1), axes=list(grid=TRUE, y=list(style="stacked")), lattice=list(key.args=list(columns=1))) ## ----------------------------------------------------------------------------- library("nnet") # for multinom() mr1 <- multinom(vote ~ age + gender + economic.cond.national + economic.cond.household + Blair + Hague + Kennedy + Europe*political.knowledge, data=BEPS) ## ----fig42,include=TRUE,fig.width=6.5,fig.height=6.5,fig.show='hide'---------- plot(predictorEffects(mr1, ~ age + Blair + Hague + Kennedy), axes=list(grid=TRUE, x=list(rug=FALSE)), lattice=list(key.args=list(columns=1)), lines=list(multiline=TRUE, col=c("blue", "red", "orange"))) ## ----fig43,include=TRUE,fig.width=10,fig.height=5,fig.show='hide'------------- plot(predictorEffects(mr1, ~ Europe + political.knowledge, xlevels=list(political.knowledge=0:3, Europe=c(1, 6, 11))), axes=list(grid=TRUE, x=list(rug=FALSE, Europe=list(ticks=list(at=c(1, 6, 11))), political.knowledge=list(ticks=list(at=0:3))), y=list(style="stacked")), lines=list(col=c("blue", "red", "orange")), lattice=list(key.args=list(columns=1), strip=list(factor.names=FALSE))) ## ----eval=FALSE--------------------------------------------------------------- # effectsTheme()