library(lme4) library(lattice) ############################################################################################################ # Load a data set that includes some other predictors, both continuous and categorical. # This file will only show the use of one predictor at a time. # Note that I had to specify a "na.string" because the source file had some missing values. ############################################################################################################ myd<-read.csv("Example simple predictors VG.dat", na.string="") myd<-subset(myd, !is.na(Waited)) # screening missing values for DV for simplicity # I'm going to create new predictors for the transformed magnitude ratio and log delay. It makes things # easier when I start plotting. myd$logMagRatio <- log(myd$LLmagn/myd$SSmagn) # logDelay is a constant in this data, but still creating a column with the values myd$logDelayRatio <- log((myd$LLdelay+1)/(myd$SSdelay+1)) # "+ 1" wasn't necessary for this data, but still done for consistency # The first command excludes a predictor, the second one uses a between-subjects predictor, Sex. myr1<-glmer(Waited~logMagRatio+logDelayRatio-1+(logMagRatio+logDelayRatio-1|Subject), data=myd, family=binomial) contrasts(myd$Sex)=contr.sum(2) myr2<-glmer(Waited~logMagRatio+logDelayRatio + Sex:logMagRatio+Sex:logDelayRatio-1+(logMagRatio+logDelayRatio-1|Subject), data=myd, family=binomial) summary(myr1) summary(myr2) # Okay, some honesty here. I'm vexed by the output from the summary and the warning message. Something odd is # is going on with the coding of Sex for this dataset. So, I'm going to manually recode it using effect coding and try again. myd$c.Sex<-ifelse(myd$Sex=="f", 1, -1) myr2.1<-glmer(Waited~logMagRatio+logDelayRatio + c.Sex:logMagRatio+c.Sex:logDelayRatio-1+(logMagRatio+logDelayRatio-1|Subject), data=myd, family=binomial) # Okay, the manual version worked better, so let's compare # Look at the results summary(myr1) summary(myr2.1) AIC(myr1, myr2.1) # To compare the models # myr1 has the lower AIC and is thus the better model. Of course, the lack # statistical significance for the sex interactions in myr2.1 suggests the same thing. # If you're uncomfortable with interpreting the Fixed effects table in the summary, I'll unpack this one. # Average slope for magnitude ratio = 0.7692 Female slope = 0.7692+(1*-.1105) = 0.6587. Male slope = 0.7692 +(-1*-.1105) = 0.8797. # Ditto for delay slope: Average = -1.3872, Female = -1.3872+.4311 = -.9561, Male = -1.3872-.4311 = -1.8183 # Another way to do the above is using the lsmeans package library(lsmeans) lstrends(myr2.1, ~ c.Sex, var="logMagRatio", at = list(c.Sex = c(-1,1))) # to examine slope differences across conditions lstrends(myr2.1, ~ c.Sex, var="logDelayRatio", at = list(c.Sex = c(-1,1))) # to examine slope differences across conditions # Wait, the above didn't work! Oh... logDelayRatio is a constant so it's not so easy... lsmeans(myr2.1, ~ c.Sex, at = list(c.Sex = c(-1,1), logMagRatio=0)) # this grabs intercepts, but off my factor of log(8.5/3.5) -1.6133/log(8.5/3.5) = -1.818 -0.8483/log(8.5/3.5) = -0.956 ###################################################################################################### # Look at the results. Careful if you have missing values for a variable in xyplot - screws things up ###################################################################################################### xyplot(fitted(myr1)~logMagRatio|as.factor(Subject), data=myd, type="a", ylab="P(Waited)") xyplot(fitted(myr2.1)~logMagRatio|as.factor(Subject), group=Sex, data=myd, type="a", ylab="P(Waited)") xyplot(fitted(myr2.1)~logMagRatio|Sex, group=as.factor(Subject), data=myd, type="a", ylab="P(Waited)") ###################################################################################################### # Try the original fits with a continuous predictor ###################################################################################################### # Since I'm going to include it in an interaction, I'm going to center it first. myd$c.VGExp <- myd$VG.Experience - mean(myd$VG.Experience) myr3<-glmer(Waited~logMagRatio+logDelayRatio + c.VGExp:logMagRatio+c.VGExp:logDelayRatio-1+(logMagRatio+logDelayRatio-1|Subject), data=myd, family=binomial) # The standard regression weights are the top two. The modifications as a function of VG Experience are second. summary(myr3) # Nothing significant. # Let's try a more interesting example. This continuous variable is within subject so included as random effects. # More complex random effect structures take MUCH longer to fit. myd$c.level <- myd$Level - mean(myd$Level) myr4<-glmer(Waited~logMagRatio+logDelayRatio + c.level:logMagRatio+c.level:logDelayRatio-1+(logMagRatio+logDelayRatio + c.level:logMagRatio+c.level:logDelayRatio-1|Subject), data=myd, family=binomial) summary(myr4) # Control by magnitude is increasing quite a bit over the four levels (i.e., blocks of training). Control by # delay is also increasing, but not significantly so. lstrends(myr4, ~ c.level, var="logMagRatio", at = list(c.level = c(1,2,3,4)-2.408)) # to examine magnitude slope differences across levels coef(myr4) # Individual differences in all four effects. densityplot(~coef(myr4)$Subject[[3]]) # histogram of individual magnitude changes across levels; more positive than negative densityplot(~coef(myr4)$Subject[[4]]) # histogram of individual delay changes across levels; centered around zero