# NOT RUN {
# }
# NOT RUN {
###########################################################################
# we will stay as close to the older optimizer of lme4 as possible -
# this requires the optimx package and using the control option of lmer()
###########################################################################
require(optimx)
###########################################################################
# fitting a cosine with a spline (simulated data)
###########################################################################
require("rms", quietly=TRUE, character=TRUE)
require("lme4", quietly=TRUE, character=TRUE)
dfr = makeSplineData.fnc()
table(dfr$Subject)
xylowess.fnc(Y ~ X | Subject, data = dfr)
# the smoother doesn't recognize the cosine function implemented in makeSplineData.fnc()
dev.off()
dfr.lmer = lmer(Y ~ rcs(X, 5) + (1|Subject), data = dfr,
control = lmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
plotLMER.fnc(dfr.lmer)
# comparison with ols from Design package
dfr.lm = lm(Y~Subject+rcs(X), data=dfr, x=T, y=T)
dfr$fittedOLS = fitted(dfr.lm)
dfr$fittedLMER = as.vector(dfr.lmer@pp$X %*% fixef(dfr.lmer))
# we plot the lmer() fit in blue, the ols() fit in red (both adjusted for
# subject S1), and plot the underlying model in green
plot(dfr[dfr$Subject=="S1",]$X,
dfr[dfr$Subject=="S1",]$fittedLMER + ranef(dfr.lmer)[[1]]["S1",],
col="blue", ylim = c(24,30), xlab="X", ylab="Y", type="n")
lines(dfr[dfr$Subject=="S1",]$X, dfr[dfr$Subject=="S1",]$fittedOLS, col="red")
lines(dfr[dfr$Subject=="S1",]$X, dfr[dfr$Subject=="S1",]$fittedLMER, col="blue")
lines(dfr[dfr$Subject=="S1",]$X, dfr[dfr$Subject=="S1",]$y+
ranef(dfr.lmer)[[1]]["S1",], col="green")
legend(2,30,c("30+cos(x)", "lmer (S1)", "ols (S1)"), lty=rep(1,3),
col=c("green", "blue", "red"))
#############################################################
# a model with a raw polynomial
#############################################################
bg.lmer = lmer(LogRT ~ PC1+PC2+PC3 + ReadingScore +
poly(OrthLength, 2, raw=TRUE) + LogFrequency + LogFamilySize +
(1|Word) + (1|Subject)+(0+OrthLength|Subject) +
(0+LogFrequency|Subject), data = beginningReaders,
control=lmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
pars = par()
par(mfrow=c(3,3), mar=c(5,5,1,1))
plotLMER.fnc(bg.lmer, fun=exp, ylabel = "RT (ms)")
#############################################################
# a model with an interaction involving numeric predictors
#############################################################
danish.lmer = lmer(LogRT ~ PC1 + PC2 + PrevError + Rank +
ResidSemRating + ResidFamSize + LogWordFreq*LogAffixFreq*Sex +
poly(LogCUP, 2, raw=TRUE) + LogUP + LogCUPtoEnd +
(1|Subject) + (1|Word) + (1|Affix), data = danish,
control=lmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
danish.lmerA = lmer(LogRT ~ PC1 + PC2 + PrevError + Rank +
ResidSemRating + ResidFamSize + LogWordFreq*LogAffixFreq*Sex +
poly(LogCUP, 2, raw=TRUE) + LogUP + LogCUPtoEnd +
(1|Subject) + (1|Word) + (1|Affix), data = danish,
subset=abs(scale(resid(danish.lmer)))<2.5,
control=lmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
# plot for reference level of Sex
plotLMER.fnc(danish.lmerA, pred = "LogAffixFreq",
intr=list("LogWordFreq", round(quantile(danish$LogWordFreq),3), "beg",
list(c("red", "green", "blue", "yellow", "purple"), rep(1,5))),
ylimit=c(6.5,7.0))
# this model has a significant three-way interaction
# for visualization, we can either relevel Sex and refit,
# or make use of the control option. First releveling:
danish$Sex=relevel(danish$Sex, "F")
danish.lmerF = lmer(LogRT ~ PC1 + PC2 + PrevError + Rank +
ResidSemRating + ResidFamSize + LogWordFreq*LogAffixFreq*Sex +
poly(LogCUP, 2, raw=TRUE) + LogUP + LogCUPtoEnd +
(1|Subject) + (1|Word) + (1|Affix), data = danish,
control=lmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
danish$Sex=relevel(danish$Sex, "M")
danish.lmerM = lmer(LogRT ~ PC1 + PC2 + PrevError + Rank +
ResidSemRating + ResidFamSize + LogWordFreq*LogAffixFreq*Sex +
poly(LogCUP, 2, raw=TRUE) + LogUP + LogCUPtoEnd +
(1|Subject) + (1|Word) + (1|Affix), data = danish,
control=lmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
# Next preparing for using the control option:
#
# names(fixef(danish.lmer))[10] # SexM
# unique(danish.lmer@pp$X[,10]) # 1 0
par(mfrow=c(2,2))
plotLMER.fnc(danish.lmer, pred="LogWordFreq", ylimit=c(6.5,7.0),
intr=list("LogAffixFreq", round(quantile(danish$LogAffixFreq),2), "end"),
control=list("SexM", 0))
mtext("females", line=1.5, cex=0.9)
plotLMER.fnc(danish.lmer, pred="LogWordFreq", ylimit=c(6.5,7.0),
intr=list("LogAffixFreq", round(quantile(danish$LogAffixFreq),2), "end"),
control=list("SexM", 1))
mtext("males", line=1.5, cex=0.9)
plotLMER.fnc(danish.lmerF, pred="LogWordFreq", ylimit=c(6.5,7.0),
intr=list("LogAffixFreq", round(quantile(danish$LogAffixFreq),2), "end"))
mtext("females", line=1.5, cex=0.9)
plotLMER.fnc(danish.lmerM, pred="LogWordFreq", ylimit=c(6.5, 7.0),
intr=list("LogAffixFreq", round(quantile(danish$LogAffixFreq),2), "end"))
mtext("males", line=1.5, cex=0.9)
par(mfrow=c(1,1))
#############################################################
# calculating effect sizes, defined as max - min
#############################################################
# effect size for a covariate
dfr = plotLMER.fnc(danish.lmerA, pred = "LogCUP", withList=TRUE)
max(dfr$LogCUP$Y)-min(dfr$LogCUP$Y)
# effect size for a factor
dfr = plotLMER.fnc(danish.lmerA, pred = "PrevError", withList=TRUE)
max(dfr$PrevError$Y)-min(dfr$PrevError$Y)
# effect sizes for the quantiles in an interaction plot
dfr = plotLMER.fnc(danish.lmerA, pred = "LogAffixFreq",
withList=TRUE,
intr=list("LogWordFreq", round(quantile(danish$LogWordFreq),3), "beg"))
unlist(lapply(dfr$LogAffixFreq, FUN=function(X)return(max(X$Y)-min(X$Y))))
#############################################################
# plotting an interaction between two factors
#############################################################
danish$WordFreqFac = danish$LogWordFreq > median(danish$LogWordFreq)
danish.lmer2 = lmer(LogRT ~ WordFreqFac*Sex +
(1|Subject) + (1|Word) + (1|Affix), data = danish,
control=lmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
plotLMER.fnc(danish.lmer2, pred = "Sex",
intr=list("WordFreqFac", c("TRUE", "FALSE"), "end",
list(c("red", "blue"), rep(1,2))),
ylimit=c(6.7,6.9), cexsize=1.0, addlines=TRUE)
#############################################################
# a generalized linear mixed-effects model
#############################################################
dative.lmer = glmer(RealizationOfRecipient ~
AccessOfTheme + AccessOfRec + LengthOfRecipient + AnimacyOfRec +
AnimacyOfTheme + PronomOfTheme + DefinOfTheme + LengthOfTheme +
SemanticClass + Modality + (1|Verb),
data = dative, family = "binomial",
control=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
par(mfrow=c(3,4),mar=c(5,5,1,1))
plotLMER.fnc(dative.lmer, fun=plogis, addlines=TRUE)
# with user-specified labels for the x-axis
par(mfrow=c(3,4),mar=c(5,5,1,1))
plotLMER.fnc(dative.lmer, fun=plogis, addlines=TRUE,
xlabs=unlist(strsplit("abcdefghij","")))
par(pars)
# }
Run the code above in your browser using DataLab