R Notebook

In the datafile: j = stim and i = participant, c = condition, y = dependent variable The packages below are already on my R installation, but for you to install them, you just need to type: install.packages(“lme4”), then the library command to load them

library(readr) #library to read in all sorts of datafiles
mixedexample<- read.csv(url("http://psych.colorado.edu/~cjudd/mixedexample.csv")) #read the csv file in directly from the website
mixedexample <- within(mixedexample, {
  i <- factor(i)
  j <- factor(j)
  }) # making sure the items and participant level variables are factors 

Note on one command below: the control option here is optional, the default optimization function for lme4 is fine, but bobyqa is recommended by ben bolker (one of the package creators and a god), so I go with that - things also tend to converge nicer when you have a ton of data with bobyqa

library(lme4) #best mixed effect library ever
library(lmerTest) #add on to lme4 to get p-values because people will die if you don't have them
m.a<-lmer(y ~ c + (1+c|i) , data = mixedexample, REML = FALSE,control=lmerControl(optimizer="bobyqa")) # model fit without stimulus effect
summary(m.a) # get a print out of the model results

Linear mixed model fit by maximum likelihood  
t-tests use  Satterthwaite approximations to degrees of freedom ['lmerMod']
Formula: y ~ c + (1 + c | i)
   Data: mixedexample
Control: lmerControl(optimizer = "bobyqa")

     AIC      BIC   logLik deviance df.resid 
  5310.7   5339.5  -2649.3   5298.7      894 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3712 -0.6962  0.0061  0.6767  2.7851 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 i        (Intercept)  4.011   2.003        
          c            3.484   1.867    0.30
 Residual             19.228   4.385        
Number of obs: 900, groups:  i, 30

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)  -0.1804     0.3938 30.0000  -0.458     0.65    
c             2.5211     0.4490 30.0000   5.615 4.11e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
  (Intr)
c 0.209 

m.a.5<-lmer(y ~ c + (1|i) , data = mixedexample, REML = FALSE,control=lmerControl(optimizer="bobyqa")) # do we need that random slope for condition?
summary(m.a.5) 

Linear mixed model fit by maximum likelihood  
t-tests use  Satterthwaite approximations to degrees of freedom ['lmerMod']
Formula: y ~ c + (1 | i)
   Data: mixedexample
Control: lmerControl(optimizer = "bobyqa")

     AIC      BIC   logLik deviance df.resid 
  5322.1   5341.3  -2657.1   5314.1      896 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4306 -0.6594  0.0316  0.6727  2.7564 

Random effects:
 Groups   Name        Variance Std.Dev.
 i        (Intercept)  3.981   1.995   
 Residual             20.129   4.487   
Number of obs: 900, groups:  i, 30

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  -0.1804     0.3938  30.0000  -0.458     0.65    
c             2.5211     0.2991 870.0000   8.429   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
  (Intr)
c 0.000 

m.a1<-lmer(y ~ c + (1+c|i) + (1|j), data = mixedexample, REML = FALSE,control=lmerControl(optimizer="bobyqa"))# adding in the stimuli effect
summary(m.a1)

Linear mixed model fit by maximum likelihood  
t-tests use  Satterthwaite approximations to degrees of freedom ['lmerMod']
Formula: y ~ c + (1 + c | i) + (1 | j)
   Data: mixedexample
Control: lmerControl(optimizer = "bobyqa")

     AIC      BIC   logLik deviance df.resid 
  5195.1   5228.7  -2590.5   5181.1      893 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.11699 -0.66211 -0.01678  0.62057  2.80807 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 i        (Intercept)  4.200   2.049        
          c            4.114   2.028    0.27
 j        (Intercept)  3.485   1.867        
 Residual             15.561   3.945        
Number of obs: 900, groups:  i, 30; j, 30

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)   
(Intercept)  -0.1804     0.5229 51.8200  -0.345  0.73147   
c             2.5211     0.8191 40.5400   3.078  0.00373 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
  (Intr)
c 0.087 

anova(m.a,m.a.5) #use anova to run a chi-sqaured test on model deviance, m.a - model with random condition slope, fits better

Data: mixedexample
Models:
..1: y ~ c + (1 | i)
object: y ~ c + (1 + c | i)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
..1     4 5322.1 5341.3 -2657.1   5314.1                             
object  6 5310.7 5339.5 -2649.3   5298.7 15.436      2  0.0004447 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(m.a,m.a1) #use anova to run a chi-sqaured test on model deviance, m.a1 is much much better

Data: mixedexample
Models:
object: y ~ c + (1 + c | i)
..1: y ~ c + (1 + c | i) + (1 | j)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
object  6 5310.7 5339.5 -2649.3   5298.7                             
..1     7 5195.1 5228.7 -2590.5   5181.1 117.62      1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

randomeffects<-ranef(m.a1) #get the random effects
randomeffects$j #RE for j

randomeffects$i # RE for i

fixedeffects<-fixef(m.a1) #grab the fixed effects
fixedeffects

(Intercept)           c 
 -0.1804222   2.5210667 

library(sjPlot) ###couple of cool packages that make nice figures directly from model
library(sjstats)
library(sjmisc)
icc(m.a1) #get icc for both random effects

Linear mixed model
 Family: gaussian (identity)
Formula: y ~ c + (1 + c | i) + (1 | j)

  ICC (i): 0.180668
  ICC (j): 0.149932

sjp.lmer(m.a1,type ="re",sort.est = "sort.all",facet.grid=FALSE) #nice way to plot random effects

sjp.lmer(m.a1, type = "fe") # can do the same for FE