# Note: Some of these examples are marked as "don't run"
# to reduce the execution times of the examples.
mod.cowles <- glm(volunteer ~ sex + neuroticism*extraversion,
data=Cowles, family=binomial)
eff.cowles <- allEffects(mod.cowles, xlevels=list(extraversion=seq(0, 24, 6)),
given.values=c(sexmale=0.5))
eff.cowles
as.data.frame(eff.cowles[[2]])
# the following are equivalent:
eff.ne <- effect("neuroticism*extraversion", mod.cowles)
Eff.ne <- Effect(c("neuroticism", "extraversion"), mod.cowles)
all.equal(eff.ne$fit, Eff.ne$fit)
plot(eff.cowles, 'sex', ylab="Prob(Volunteer)")
plot(eff.cowles, 'neuroticism:extraversion', ylab="Prob(Volunteer)",
ticks=list(at=c(.1,.25,.5,.75,.9)))
plot(eff.cowles, 'neuroticism:extraversion', multiline=TRUE,
ylab="Prob(Volunteer)")
plot(effect('sex:neuroticism:extraversion', mod.cowles,
xlevels=list(extraversion=seq(0, 24, 6))), multiline=TRUE)
# a nested model:
mod <- lm(log(prestige) ~ income:type + education, data=Prestige)
# does not work: effect("income:type", mod, transformation=list(link=log, inverse=exp))
plot(Effect(c("income", "type"), mod, transformation=list(link=log, inverse=exp),
ylab="prestige")) # works
if (require(nnet)){
mod.beps <- multinom(vote ~ age + gender + economic.cond.national +
economic.cond.household + Blair + Hague + Kennedy +
Europe*political.knowledge, data=BEPS)
plot(effect("Europe*political.knowledge", mod.beps,
xlevels=list(political.knowledge=0:3)))
plot(Effect(c("Europe", "political.knowledge"), mod.beps,
xlevels=list(Europe=1:11, political.knowledge=0:3),
given.values=c(gendermale=0.5)),
style="stacked", colors=c("blue", "red", "orange"), rug=FALSE)
plot(effect("Europe*political.knowledge", mod.beps, # equivalent
xlevels=list(political.knowledge=0:3),
given.values=c(gendermale=0.5)),
style="stacked", colors=c("blue", "red", "orange"), rug=FALSE)
}
if (require(MASS)){
mod.wvs <- polr(poverty ~ gender + religion + degree + country*poly(age,3),
data=WVS)
plot(effect("country*poly(age, 3)", mod.wvs))
plot(Effect(c("country", "age"), mod.wvs), style="stacked")
plot(effect("country*poly(age, 3)", mod.wvs), style="stacked") # equivalent
plot(effect("country*poly(age, 3)", latent=TRUE, mod.wvs))
}
mod.pres <- lm(prestige ~ log(income, 10) + poly(education, 3) + poly(women, 2),
data=Prestige)
eff.pres <- allEffects(mod.pres, xlevels=50)
plot(eff.pres)
plot(eff.pres[1],
transform.x=list(income=list(trans=log10, inverse=function(x) 10^x)),
ticks.x=list(income=list(at=c(1000, 2000, 5000, 10000, 20000))))
# linear model with log-response and log-predictor
# to illustrate transforming axes and setting tick labels
mod.pres1 <- lm(log(prestige) ~ log(income) + poly(education, 3) + poly(women, 2),
data=Prestige)
# effect of the log-predictor
eff.log <- Effect(c("income"), mod.pres1)
# effect of the log-predictor transformed to the arithmetic scale
eff.trans <- Effect(c("income"), mod.pres1, transformation=list(link=log, inverse=exp))
#variations:
# y-axis: scale is log, tick labels are log
# x-axis: scale is arithmetic, tick labels are arithmetic
plot(eff.log)
# y-axis: scale is log, tick labels are log
# x-axis: scale is log, tick labels are arithmetic
plot(eff.log, transform.x=list(income=c(trans=log, inverse=exp)),
ticks.x=list(income=list(at=c(1000, 2000, 5000, 10000, 20000))),
xlab="income, log-scale")
# y-axis: scale is log, tick labels are airthmetic
# x-axis: scale is arithmetic, tick labels are arithmetic
plot(eff.trans, ylab="prestige")
# y-axis: scale is arithmetic, tick labels are airthmetic
# x-axis: scale is arithmetic, tick labels are arithmetic
plot(eff.trans, rescale.axis=FALSE, ylab="prestige")
# y-axis: scale is log, tick labels are arithmetic
# x-axis: scale is log, tick labels are arithmetic
plot(eff.trans, transform.x=list(income=c(trans=log, inverse=exp)),
ticks.x=list(income=list(at=c(1000, 2000, 5000, 10000, 20000))),
xlab="income, log-scale", ylab="prestige, log-scale",
main="Both effect and X in log-scale")
# y-axis: scale is arithmetic, tick labels are airthmetic
# x-axis: scale is log, tick labels are arithmetic
plot(eff.trans, transform.x=list(income=c(trans=log, inverse=exp)),
ticks.x=list(income=list(at=c(1000, 2000, 5000, 10000, 20000))),
rescale.axis=FALSE,
xlab="income, log-scale", ylab="prestige")
if (require(nlme)){ # for gls()
mod.hart <- gls(fconvict ~ mconvict + tfr + partic + degrees, data=Hartnagel,
correlation=corARMA(p=2, q=0), method="ML")
plot(allEffects(mod.hart))
detach(package:nlme)
}
if (require(lme4)){
data(cake, package="lme4")
fm1 <- lmer(angle ~ recipe * temperature + (1|recipe:replicate), cake,
REML = FALSE)
plot(Effect(c("recipe", "temperature"), fm1))
plot(effect("recipe:temperature", fm1), grid=TRUE) # equivalent
detach(package:lme4)
}
if (require(nlme) && length(find.package("lme4", quiet=TRUE)) > 0){
data(cake, package="lme4")
cake$rep <- with(cake, paste( as.character(recipe), as.character(replicate), sep=""))
fm2 <- lme(angle ~ recipe * temperature, data=cake,
random = ~ 1 | rep, method="ML")
plot(Effect(c("recipe", "temperature"), fm2))
plot(effect("recipe:temperature", fm2), grid=TRUE) # equivalent
}
detach(package:nlme)
if (require(poLCA)){
data(election)
f2a <- cbind(MORALG,CARESG,KNOWG,LEADG,DISHONG,INTELG,
MORALB,CARESB,KNOWB,LEADB,DISHONB,INTELB)~PARTY*AGE
nes2a <- poLCA(f2a,election,nclass=3,nrep=5)
plot(Effect(c("PARTY", "AGE"), nes2a), style="stacked")
}
# mlm example
if (require(heplots)) {
data(NLSY, package="heplots")
mod <- lm(cbind(read,math) ~ income+educ, data=NLSY)
eff.inc <- Effect("income", mod)
plot(eff.inc)
eff.edu <- Effect("educ", mod)
plot(eff.edu, rug=FALSE, grid=TRUE)
plot(Effect("educ", mod, response="read"))
detach(package:heplots)
}
# component + residual plot examples
Prestige$type <- factor(Prestige$type, levels=c("bc", "wc", "prof"))
mod.prestige.1 <- lm(prestige ~ income + education, data=Prestige)
plot(allEffects(mod.prestige.1, partial.residuals=TRUE)) # standard C+R plots
mod.prestige.2 <- lm(prestige ~ type*(income + education), data=Prestige)
plot(allEffects(mod.prestige.2, partial.residuals=TRUE))
mod.prestige.3 <- lm(prestige ~ type + income*education, data=Prestige)
plot(Effect(c("income", "education"), mod.prestige.3, partial.residuals=TRUE),
span=1)
# artificial data
set.seed(12345)
x1 <- runif(500, -75, 100)
x2 <- runif(500, -75, 100)
y <- 10 + 5*x1 + 5*x2 + x1^2 + x2^2 + x1*x2 + rnorm(500, 0, 1e3)
Data <- data.frame(y, x1, x2)
mod.1 <- lm(y ~ poly(x1, x2, degree=2, raw=TRUE), data=Data)
# raw=TRUE necessary for safe prediction
mod.2 <- lm(y ~ x1*x2, data=Data)
mod.3 <- lm(y ~ x1 + x2, data=Data)
plot(Effect(c("x1", "x2"), mod.1, partial.residuals=TRUE)) # correct model
plot(Effect(c("x1", "x2"), mod.2, partial.residuals=TRUE)) # wrong model
plot(Effect(c("x1", "x2"), mod.3, partial.residuals=TRUE)) # wrong modelRun the code above in your browser using DataLab