Example one group glmer
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. The data is in CSV format, so feel free to open it to examine the structure. R will look for the data in the working directory. If your file isn’t there, you either need to move it there (use “getwd()” to find out what it is) or change the working directory (it’s in your menu).
I’m assuming that your data is in the long format, in comma-delimited form, and has no missing values.
If you have missing values, you may need to use na.strings in the read.csv call. This data set always has SS delays that are zero. The data were simulated to include people with different k values from a Kirby questionnaire - not ideal for the approach being advocated.
myd<-read.csv("Example one group data glmer.dat") The following is the basic multilevel model without an intercept with the recommended log transformations of the magnitude and delay predictors.
myr1<-glmer(Waited~log(LLmagn/SSmagn)+log(LLdelay+1)-1+(log(LLmagn/SSmagn)+log(LLdelay+1)-1|Subject), data=myd, family=binomial)This version assumes a linear relationship for both predictors. The “I” is necessary because we’re doing math on the predictors (dividing them or adding 1).
myr2<-glmer(Waited~I(LLmagn/SSmagn)+I(LLdelay+1)-1+(I(LLmagn/SSmagn)+I(LLdelay+1)-1|Subject), data=myd, family=binomial)Look at the results for both the logarithmic and linear models.
summary(myr1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Waited ~ log(LLmagn/SSmagn) + log(LLdelay + 1) - 1 + (log(LLmagn/SSmagn) +
## log(LLdelay + 1) - 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 343.8 368.2 -166.9 333.8 967
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.3087 -0.0065 0.0000 0.0024 2.6672
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject log(LLmagn/SSmagn) 148.4 12.182
## log(LLdelay + 1) 2.2 1.483 1.00
## Number of obs: 972, groups: Subject, 36
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## log(LLmagn/SSmagn) 35.3427 5.2162 6.776 1.24e-11 ***
## log(LLdelay + 1) -2.3050 0.4419 -5.216 1.83e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## l(LL/S
## lg(LLdly+1) -0.424summary(myr2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Waited ~ I(LLmagn/SSmagn) + I(LLdelay + 1) - 1 + (I(LLmagn/SSmagn) +
## I(LLdelay + 1) - 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 395.6 420.0 -192.8 385.6 967
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3047 -0.0018 0.0002 0.0334 3.4062
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject I(LLmagn/SSmagn) 3.256767 1.80465
## I(LLdelay + 1) 0.009065 0.09521 1.00
## Number of obs: 972, groups: Subject, 36
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## I(LLmagn/SSmagn) 5.3815 0.7037 7.647 2.06e-14 ***
## I(LLdelay + 1) -0.1452 0.0271 -5.357 8.48e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## I(LL/S
## I(LLdely+1) -0.338AIC(myr1, myr2) # To compare the models## df AIC
## myr1 5 343.8445
## myr2 5 395.5878Note that a likelihood ratio test could normally be performed using “anova(myr1, myr2),” but an LR test is only valid when the two models are “nested” meaning that one contains a subset of the terms used in the other model. This is not the case for these two models.
Look at the results. Note that plotting binomial data isn’t very interesting, but I’ll try it anyway in the second graph. The second graph will jitter the 0/1 values to make them easier to see. These graphs aren’t super pretty. Better control can be obtained using ggplot which will be demonstrated in a separate file.
xyplot(fitted(myr1)~log(LLmagn/SSmagn)|as.factor(Subject), data=myd, type="a", ylab="P(Waited)")
xyplot(fitted(myr1)+jitter(Waited)~log(LLmagn/SSmagn)|as.factor(Subject), data=myd, type=c("a", "p"), distribute.type=T, ylab="P(Waited)")
Individual differences
Let’s examine individual differences. I’ll focus on the myr1 model.
An aside - for this simulated data, the Kirby survey was used assuming a hyperbolic discounting function. This means that the delay and magnitude predictors are highly correlated (r = -.97). Thus, the individual subject beta weights are also highly correlated (r = 1.0). This confirms the problems with the Kirby for identifying the independent contributions of delay and magnitude.
coef(myr1) # shows the individual regression weights## $Subject
## log(LLmagn/SSmagn) log(LLdelay + 1)
## 1 50.49909 -0.4597151
## 2 51.76085 -0.3060965
## 3 50.64251 -0.4422536
## 4 49.37580 -0.5964746
## 5 51.76085 -0.3060965
## 6 50.64251 -0.4422536
## 7 49.41046 -0.5922547
## 8 49.33829 -0.6010413
## 9 45.70386 -1.0435301
## 10 46.98487 -0.8875674
## 11 45.67258 -1.0473379
## 12 42.35730 -1.4509699
## 13 44.24235 -1.2214668
## 14 44.24235 -1.2214668
## 15 38.88699 -1.8734768
## 16 38.86653 -1.8759679
## 17 36.76002 -2.1324337
## 18 36.76002 -2.1324337
## 19 36.76002 -2.1324337
## 20 33.03570 -2.5858671
## 21 32.92600 -2.5992222
## 22 31.04058 -2.8287705
## 23 31.04058 -2.8287705
## 24 31.04058 -2.8287705
## 25 29.03532 -3.0729094
## 26 27.44849 -3.2661041
## 27 24.66269 -3.6052724
## 28 27.40570 -3.2713142
## 29 24.66269 -3.6052724
## 30 21.53772 -3.9857349
## 31 20.12675 -4.1575194
## 32 18.00711 -4.4155833
## 33 18.00711 -4.4155833
## 34 18.00711 -4.4155833
## 35 14.40431 -4.8542208
## 36 11.53868 -5.2031091
##
## attr(,"class")
## [1] "coef.mer"Plots the combinations; each dot is an individual. Not very interesting for this simulated data.
plot(coef(myr1)$Subject[,1] ~ coef(myr1)$Subject[,2], ylab="log(LLmagn/SSmagn)", xlab="log(LLdelay+1)") 
A histogram of each of the regression weights.
hist(coef(myr1)$Subject[,1], xlab="log(LLmagn/SSmagn) slope") 
hist(coef(myr1)$Subject[,2], xlab="log(LLdelay+1) slope") 
Non-zero SS delays
Now, let’s use data with non-zero SS delays.
myd<-read.csv("Simulated NonZero data.dat") The following model uses the log ratios discussed in the manuscript.
myrNZ<-glmer(Waited~log(LLmagn/SSmagn)+log((LLdelay+1)/(SSdelay+1))-1+(log(LLmagn/SSmagn)+log((LLdelay+1)/(SSdelay+1))-1|Subject), data=myd, family=binomial)
summary(myrNZ)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ log(LLmagn/SSmagn) + log((LLdelay + 1)/(SSdelay + 1)) -
## 1 + (log(LLmagn/SSmagn) + log((LLdelay + 1)/(SSdelay + 1)) -
## 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 498.4 526.4 -244.2 488.4 1995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2739 -0.0020 0.0003 0.0216 6.6002
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject log(LLmagn/SSmagn) 1.532 1.238
## log((LLdelay + 1)/(SSdelay + 1)) 4.308 2.076 -0.55
## Number of obs: 2000, groups: Subject, 20
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## log(LLmagn/SSmagn) 8.3268 0.8471 9.83 < 2e-16 ***
## log((LLdelay + 1)/(SSdelay + 1)) -4.7046 0.6516 -7.22 5.21e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## l(LL/S
## l((LL+1)/+1 -0.753Picture of the results.
xyplot(fitted(myrNZ)~log(LLmagn/SSmagn)|as.factor(Subject), group=cut(log((LLdelay+1)/(SSdelay+1)), 3), data=myd, type="a", ylab="P(Waited)", auto.key=T)
The above plot mixes lots of combinations of continous values of delay and magnitude, so it’s particularly ugly. This analysis really would be better served by a ggplot of the model fit. However, I’m going to demonstrate another alternative using a bubble plot where the color of each dot indicates the likelihood of waiting for each combination of the two ratios. In the example below, green = wait, red = don’t.
myd$fit<-fitted(myrNZ) # grabs the P(Waiting) predicted by the model
xyplot(log(LLmagn/SSmagn)~log((LLdelay+1)/(SSdelay+1))|as.factor(Subject), data=myd,
col=rgb(1-sqrt(myd$fit), sqrt(myd$fit), 0), # This sets the color in RGB notation.
panel=function(x,y,..., col, subscripts) {
panel.xyplot(x,y,col=col[subscripts])
}
)
Here I’m plotting the same data but where the size of the dot corresponds to the P(Waiting). The problem with this command is that the small dots are too small for most purposes. Dot sizes are designated using the “cex” parameter. First I’m changing the range using a power function.
xyplot(log(LLmagn/SSmagn)~log((LLdelay+1)/(SSdelay+1))|as.factor(Subject), data=myd,
cex=myd$fit^.25,
panel=function(x,y,..., cex, subscripts) {
panel.xyplot(x,y,cex=cex[subscripts], col="black")
}
)
Still too small, so next I’m adding a constant to the point size (designed with “cex”) to increase the minimum size and changing the power to .10 to shrink the range of circle sizes.
xyplot(log(LLmagn/SSmagn)~log((LLdelay+1)/(SSdelay+1))|as.factor(Subject), data=myd,
cex=myd$fit^.1+.2,
panel=function(x,y,..., cex, subscripts) {
panel.xyplot(x,y,cex=cex[subscripts], col="black")
}
)