Loading libraries

library(lme4)    # loads the library containing the lmer and glmer commands for multilevel modeling
## Loading required package: Matrix
library(lattice) # use for basic graphing

Load 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.424
summary(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.338
AIC(myr1, myr2) # To compare the models
##      df      AIC
## myr1  5 343.8445
## myr2  5 395.5878

Note 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)")