Example simple predictors
Loading libraries
library(lme4) # loads the library containing the lmer and glmer commands for multilevel modeling## Loading required package: Matrixlibrary(lattice) # use for basic graphingLoad 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 simplicityI’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 consistencyThe 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)## Warning: Some predictor variables are on very different scales: consider
## rescalingcontrasts(myd$Sex)=contr.sum(2)
myr2<-glmer(Waited~logMagRatio+logDelayRatio
+ Sex:logMagRatio+Sex:logDelayRatio-1+(logMagRatio+logDelayRatio-1|Subject), data=myd, family=binomial)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(myr1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Waited ~ logMagRatio + logDelayRatio - 1 + (logMagRatio + logDelayRatio -
## 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 24116.7 24157.2 -12053.4 24106.7 24219
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8037 -0.6393 -0.2516 0.6747 7.3727
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject logMagRatio 1.003 1.002
## logDelayRatio 3.480 1.865 -0.75
## Number of obs: 24224, groups: Subject, 70
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## logMagRatio 0.7476 0.1259 5.938 2.89e-09 ***
## logDelayRatio -1.3020 0.2309 -5.640 1.70e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## lgMgRt
## logDelayRat -0.763
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingsummary(myr2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Waited ~ logMagRatio + logDelayRatio + Sex:logMagRatio + Sex:logDelayRatio -
## 1 + (logMagRatio + logDelayRatio - 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 24116.5 24173.2 -12051.3 24102.5 24217
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.7542 -0.6402 -0.2519 0.6725 7.3231
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject logMagRatio 0.9905 0.9952
## logDelayRatio 3.3034 1.8175 -0.75
## Number of obs: 24224, groups: Subject, 70
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## logMagRatio 0.8797 0.1951 4.509 6.50e-06 ***
## logDelayRatio -1.3872 0.2278 -6.089 1.14e-09 ***
## logMagRatio:Sexf -0.2210 0.2532 -0.873 0.3827
## logDelayRatio:Sex1 0.4311 0.2277 1.893 0.0583 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## lgMgRt lgDlyR lgMR:S
## logDelayRat -0.586
## logMgRt:Sxf -0.770 0.142
## lgDlyRt:Sx1 0.585 -0.191 -0.760
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Some predictor variables are on very different scales: consider rescalingOkay, some honesty here. I’m vexed by the output from the summary and the warning message about dropping one column. 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)## Warning: Some predictor variables are on very different scales: consider
## rescaling*Okay, the manual version worked better, so let’s compare by looking at the results
summary(myr1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Waited ~ logMagRatio + logDelayRatio - 1 + (logMagRatio + logDelayRatio -
## 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 24116.7 24157.2 -12053.4 24106.7 24219
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8037 -0.6393 -0.2516 0.6747 7.3727
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject logMagRatio 1.003 1.002
## logDelayRatio 3.480 1.865 -0.75
## Number of obs: 24224, groups: Subject, 70
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## logMagRatio 0.7476 0.1259 5.938 2.89e-09 ***
## logDelayRatio -1.3020 0.2309 -5.640 1.70e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## lgMgRt
## logDelayRat -0.763
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingsummary(myr2.1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Waited ~ logMagRatio + logDelayRatio + c.Sex:logMagRatio + c.Sex:logDelayRatio -
## 1 + (logMagRatio + logDelayRatio - 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 24116.5 24173.2 -12051.3 24102.5 24217
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.7542 -0.6402 -0.2519 0.6725 7.3231
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject logMagRatio 0.9905 0.9952
## logDelayRatio 3.3034 1.8175 -0.75
## Number of obs: 24224, groups: Subject, 70
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## logMagRatio 0.7692 0.1273 6.045 1.5e-09 ***
## logDelayRatio -1.3872 0.2291 -6.056 1.4e-09 ***
## logMagRatio:c.Sex -0.1105 0.1271 -0.869 0.3847
## logDelayRatio:c.Sex 0.4311 0.2290 1.883 0.0597 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## lgMgRt lgDlyR lMR:.S
## logDelayRat -0.762
## lgMgRt:c.Sx -0.190 0.146
## lgDlyRt:c.S 0.146 -0.196 -0.762
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingAIC(myr1, myr2.1) # To compare the models## df AIC
## myr1 5 24116.71
## myr2.1 7 24116.52There is no preference for either, so I’ll go with the simpler model. Of course, the lack of 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 emmeans package.
library(emmeans)
emtrends(myr2.1, ~ c.Sex, var="logMagRatio", at = list(c.Sex = c(-1,1))) # to examine slope differences across conditions## c.Sex logMagRatio.trend SE df asymp.LCL asymp.UCL
## -1 0.8797453 0.1962302 Inf 0.495141 1.2643495
## 1 0.6587276 0.1619023 Inf 0.341405 0.9760502
##
## Trends are based on the logit (transformed) scale
## Confidence level used: 0.95emtrends(myr2.1, ~ c.Sex, var="logDelayRatio", at = list(c.Sex = c(-1,1))) # to examine slope differences across conditions## Error in emtrends(myr2.1, ~c.Sex, var = "logDelayRatio", at = list(c.Sex = c(-1, : Provide a nonzero value of 'delta.var'Wait, the second one didn’t work! Oh… logDelayRatio is a constant so it’s not so easy…
emmeans(myr2.1, ~ c.Sex, at = list(c.Sex = c(-1,1), logMagRatio=0)) # this grabs intercepts, but result is off by factor of log(8.5/3.5) due to the constant issue## NOTE: Results may be misleading due to involvement in interactions## c.Sex emmean SE df asymp.LCL asymp.UCL
## -1 -1.6133467 0.3142335 Inf -2.229233 -0.9974603
## 1 -0.8483411 0.2577410 Inf -1.353504 -0.3431781
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95-1.6133/log(8.5/3.5) = -1.818 and -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)## Warning: Some predictor variables are on very different scales: consider
## rescalingThe standard regression weights are the top two. The modifications as a function of VG Experience are second.
summary(myr3) ## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ logMagRatio + logDelayRatio + c.VGExp:logMagRatio +
## c.VGExp:logDelayRatio - 1 + (logMagRatio + logDelayRatio -
## 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 24116.8 24173.5 -12051.4 24102.8 24217
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8423 -0.6383 -0.2515 0.6751 7.3702
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject logMagRatio 0.9844 0.9921
## logDelayRatio 3.4758 1.8643 -0.76
## Number of obs: 24224, groups: Subject, 70
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## logMagRatio 0.729987 0.125399 5.821 5.84e-09 ***
## logDelayRatio -1.308894 0.231759 -5.648 1.63e-08 ***
## logMagRatio:c.VGExp -0.013974 0.013429 -1.041 0.298
## logDelayRatio:c.VGExp -0.006951 0.024789 -0.280 0.779
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## lgMgRt lgDlyR lMR:.V
## logDelayRat -0.773
## lgMgRt:.VGE 0.123 -0.098
## lgDlyR:.VGE -0.098 0.128 -0.773
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingNothing 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)## Warning: Some predictor variables are on very different scales: consider
## rescaling## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00264452
## (tol = 0.001, component 1)summary(myr4)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ logMagRatio + logDelayRatio + c.level:logMagRatio +
## c.level:logDelayRatio - 1 + (logMagRatio + logDelayRatio +
## c.level:logMagRatio + c.level:logDelayRatio - 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 23012.5 23125.8 -11492.3 22984.5 24210
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2263 -0.5651 -0.2568 0.5798 8.4366
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject logMagRatio 0.6982 0.8356
## logDelayRatio 3.2433 1.8009 -0.66
## logMagRatio:c.level 0.3634 0.6028 -0.06 0.16
## logDelayRatio:c.level 1.0238 1.0118 -0.06 0.14 -0.80
## Number of obs: 24224, groups: Subject, 70
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## logMagRatio 0.86999 0.10953 7.943 1.97e-15 ***
## logDelayRatio -1.30818 0.22474 -5.821 5.85e-09 ***
## logMagRatio:c.level 0.32787 0.08282 3.959 7.53e-05 ***
## logDelayRatio:c.level -0.21289 0.13447 -1.583 0.113
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## lgMgRt lgDlyR lgMR:.
## logDelayRat -0.686
## lgMgRt:c.lv 0.003 0.107
## lgDlyRt:c.l -0.083 0.139 -0.824
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## convergence code: 0
## Model failed to converge with max|grad| = 0.00264452 (tol = 0.001, component 1)Control by magnitude is increasing quite a bit over the four levels (i.e., blocks of training). Control by delay is also increasing (becoming more negative), but not significantly so.
emtrends(myr4, ~ c.level, var="logMagRatio", at = list(c.level = c(1,2,3,4)-2.408)) # to examine magnitude slope differences across the four game levels, 1 through 4.## c.level logMagRatio.trend SE df asymp.LCL asymp.UCL
## -1.408 0.4083495 0.1597115 Inf 0.09532065 0.7213783
## -0.408 0.7362215 0.1145115 Inf 0.51178314 0.9606598
## 0.592 1.0640935 0.1201529 Inf 0.82859818 1.2995888
## 1.592 1.3919655 0.1716953 Inf 1.05544884 1.7284821
##
## Trends are based on the logit (transformed) scale
## Confidence level used: 0.95coef(myr4) # Individual differences in all four effects.## $Subject
## logMagRatio logDelayRatio logMagRatio:c.level logDelayRatio:c.level
## 401a 0.8959205 -2.71909249 -0.276995792 0.415906683
## 401b 1.0655225 -2.77842158 0.533139387 -0.741979676
## 401c -0.0456459 1.86076406 0.140839153 0.620301026
## 401d 0.6995972 -0.90666668 1.105947046 -0.675759835
## 402a 0.2234378 -0.19436439 -0.123541070 -0.961888170
## 402b 0.8819489 -2.34828385 0.621950763 -1.062186704
## 402c 1.5292786 -2.08688181 0.538142823 -1.938275758
## 402d 1.9855286 -5.79097373 -0.047243512 -0.060137610
## 403a 1.3741556 -2.96211633 -0.154224859 0.137447250
## 403b 0.2562086 0.65567358 -0.103382397 0.950852884
## 403c 0.3628135 -3.29423747 0.233347186 -0.637589317
## 403d 1.4259015 -2.35044620 0.992891877 -0.921586626
## 404a 1.3513720 -2.76010368 0.490185170 -0.240362609
## 404b 0.9852560 -2.53002082 -0.415293568 1.195865732
## 404c 0.3237980 0.62970877 1.145007547 -1.238778278
## 404d -1.0414860 4.98424351 1.766673657 -0.785420291
## 405a 0.1885768 -0.83126591 -0.315880173 1.079496950
## 405c -0.2700104 0.76696586 1.015444302 -1.416576451
## 405d 0.7026481 -0.73127926 0.834851830 -1.024112294
## 406a 0.4781498 -0.28863812 1.050142755 -1.959380800
## 406b -0.3374578 0.07961507 0.791917697 -1.157756708
## 406c -0.3029449 1.64327484 -0.515422895 1.428743529
## 406d 0.1601692 0.38505833 0.280168153 0.005841234
## 407a 1.5782594 -2.81544742 -0.173004264 -0.005298535
## 407b 0.5781053 1.56663259 1.030284299 -0.971982781
## 407c 2.6221696 -3.28158856 0.724232687 -1.804592647
## 407d 0.1888997 -0.82815546 0.783503086 -0.215051424
## 408a 1.9861068 -0.95048971 -0.279148499 1.128624770
## 408b 0.9673073 -1.95723209 -0.290886475 0.226514571
## 408c 0.9255218 0.93832677 0.518257145 0.836218157
## 408d 0.2302907 -1.92906572 0.529136358 -1.426580688
## 409a 2.0644727 -2.41385935 -0.226875866 0.312059891
## 409b -0.4412231 1.14883857 0.729339004 -0.337180901
## 409c 2.2411058 -3.25433091 1.014397483 -1.070970905
## 409d 0.2876285 -1.68273695 -1.054779437 1.966594507
## 410a 0.4901790 -2.75194742 -0.024229512 -0.856345092
## 410b 1.0142285 -0.89429895 0.764093385 -1.316714805
## 410d 1.5699623 0.01567417 -0.291172532 1.697137619
## 411a 0.5368710 -1.62866319 -0.592807123 0.359129866
## 411b 0.1762327 -0.26913168 -0.090389660 0.080867491
## 411c 1.6566554 -1.86316840 -0.093902630 0.320107701
## 411d 1.6954835 -2.29749814 -0.140341062 0.807456128
## 412a 1.1654914 -3.80665755 0.141760056 0.049017545
## 412b 1.4073541 -2.08298582 -0.083804679 0.041146429
## 412c 0.3459633 -1.02562901 -0.191035654 0.871317418
## 412d 0.8397674 -1.27259285 0.893208077 -1.053159479
## 413a 0.8538666 1.30942729 -0.098576579 1.476272004
## 413b 0.6404707 -0.13025751 -0.007158202 -0.086260931
## 413c 0.2409353 -0.58883193 0.186577893 -0.087484986
## 413d -0.4371907 0.78226170 0.517503129 -0.210903499
## 414a 0.3866139 0.27889884 0.268227457 0.246587284
## 414b 2.0258127 -1.65758762 0.504085913 -0.693553345
## 414c 0.6159960 0.34055605 -0.317995089 1.229640389
## 414d 0.1499799 -3.48085839 0.228941066 -0.678952806
## 415a 0.7029899 -2.27456283 1.152320824 -1.647093470
## 415b 0.8945547 -1.73205931 0.761562058 -0.325768607
## 415c 1.4278927 -2.03477588 0.546140046 -0.244675979
## 415d 0.1861072 -0.53019443 0.527037191 0.297579366
## 416a 1.0021836 -0.50807857 0.920562824 -1.699246770
## 416b 0.7518099 -1.33686276 -0.148076502 0.004169621
## 416c 0.9847632 -1.29178212 0.180526314 -0.027872267
## 416d 1.5307151 -2.35909663 0.272400281 0.156284605
## 417a 1.9970638 -1.45601604 0.675373762 -0.508488383
## 417b -0.3935161 -2.71978720 0.559699576 0.129047670
## 417c 1.5332204 -1.93163901 0.381242448 -0.164964155
## 417d 1.2182209 -2.57656430 0.504303906 -1.023035366
## 418a 1.3059376 -2.24598313 0.720134573 -0.905620309
## 418b 2.2096059 -3.55891523 0.464523263 -0.680713907
## 418d 0.9951792 -2.18065744 -0.212110810 0.646130532
## 419c 1.6122793 -3.44294476 0.699961590 -0.522882930
##
## attr(,"class")
## [1] "coef.mer"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