3
$\begingroup$

I am quite new to R and also to LME statistics and would greatly appreciate any help. I am trying to figure out whether I have changed the variables in my study accordingly, and how to read the results I got from an lme using package lme4.

My data look like:

Subject Age sex Session Match Learned Average.Rating 
 1      20   2      1     1  Target            4.4        
 1      20   2      2     1  Target            3.3        
 1      20   2      3     1  Target            3.3        
 1      20   2      4     1  Target            3.1        
 1      20   2      5     1  Target            2.0        
 1      20   2      1     1 Control            3.0        

I wanted to use an lme as over the Session variable there was a high attrition rate and I didn't want to have to discount entire individuals (as you would have to for a standard repeated measures ANOVA) but rather discount specific data points.

Firstly I changed my Subject variable to a factor, my Learned condition to a 2 level factor ('control' and 'target'), and Session to an ordered factor as this is time over the course of weeks (weeks 1, 2, 3, 4, 5).

#changing subject to factor
MatchedData$Subject<- as.factor(MatchedData$Subject)

#changing the Learned condition to a factor with 2 levels - target or 
control
MatchedData$Learned<- factor(MatchedData$Learned, 
                         levels = c('Target', 'Control'))


#changing block to be an ordered factor - 1 goes before 2, etc.
MatchedData$Session<- ordered(MatchedData$Session, 
                       levels = c(1, 2, 3, 4, 5))

I would like to know whether changing Session from an integer to an ordered factor makes sense - as it is a time course over weeks where a subject can't possibly do Session 1 before Session 2, etc.

A figure of my data:

Violin Plot of my data

Below is my model, which assumes a random intercept for every Subject, and also a random slope - the effect of time (Session). What I am interested in is the effect of Session and Learned on Average.Rating.

m.allsess <- lmer(Average.Rating~Session*Learned + (1+Session|Subject), REML = FALSE, MatchedData)

My summary(m.allsess) looks like this

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Average.Rating ~ Session * Learned + (1 + Session | Subject)
Data: MatchedData

 AIC      BIC   logLik deviance df.resid 
750.4    852.9   -349.2    698.4      356 

Scaled residuals: 
 Min      1Q  Median      3Q     Max 
 -2.4520 -0.6333 -0.1035  0.6021  4.3905 

  Random effects:
  Groups   Name        Variance Std.Dev. Corr                   
 Subject  (Intercept) 0.231059 0.48069                         
          Session.L   0.115409 0.33972  -0.43                  
          Session.Q   0.008404 0.09167  -0.95  0.41            
          Session.C   0.001357 0.03684   0.92 -0.23 -0.76      
          Session^4   0.004692 0.06850   0.40  0.60 -0.28  0.64
          Residual             0.269142 0.51879                         
      Number of obs: 382, groups:  Subject, 44

Fixed effects:
                          Estimate Std. Error t value
(Intercept)               2.04127    0.08217  24.842
Session.L                -0.97504    0.09967  -9.783
Session.Q                 0.41282    0.08661   4.766
Session.C                -0.21320    0.08491  -2.511
Session^4                 0.03327    0.08687   0.383
LearnedControl            0.20778    0.05343   3.889
Session.L:LearnedControl  0.16551    0.11858   1.396
Session.Q:LearnedControl -0.22804    0.11965  -1.906
Session.C:LearnedControl  0.15812    0.11894   1.329
Session^4:LearnedControl  0.00510    0.12071   0.042

Correlation of Fixed Effects:
            (Intr) Sssn.L Sssn.Q Sssn.C Sssn^4 LrndCn S.L:LC S.Q:LC S.C:LC
Session.L   -0.140                                                        
Session.Q   -0.142  0.097                                                 
Session.C    0.037 -0.012  0.082                                          
Session^4    0.056  0.030 -0.032  0.086                                   
LearndCntrl -0.325 -0.062  0.012  0.021 -0.011                            
Sssn.L:LrnC -0.034 -0.595 -0.046  0.003  0.011  0.105                     
Sssn.Q:LrnC  0.006 -0.039 -0.691 -0.057  0.015 -0.018  0.066              
Sssn.C:LrnC  0.010  0.002 -0.057 -0.700 -0.055 -0.030 -0.004  0.082       
Sssn^4:LrnC -0.005  0.010  0.015 -0.056 -0.695  0.016 -0.016 -0.022  0.079

I am not sure how to interpret this output. I believe the L, Q, C, and ^4 are giving me fits for a linear, quadratic, cubic, etc. models, but are these necessary? What would I report when writing this up for a journal article? I am used to ANOVA output which gives you main effects of each of your variables of interest (Learned and Session) but the interaction term.

When I use anova(m.allsess) I get something like what I am looking for but I am not sure that this is an appropriate way to look at the lme results.

Analysis of Variance Table
                Df Sum Sq Mean Sq F value
Session          4 46.856 11.7139 43.5232
Learned          1  3.807  3.8067 14.1437
Session:Learned  4  2.216  0.5540  2.0585

I have also tried to look at the results with the Anova(m.allsess) function, using type 3 sums of squares as I am interested in whether there is an interaction term. See below:

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: Average.Rating
                   Chisq Df Pr(>Chisq)    
(Intercept)     617.1297  1  < 2.2e-16 ***
Session         139.7322  4  < 2.2e-16 ***
Learned          15.1234  1  0.0001007 ***
Session:Learned   8.2342  4  0.0833664 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I have also ran the above model with Session as an integer and I get results which I am much more accustomed to, when using the summary() function.

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Average.Rating ~ Session * Learned + (1 + Session | Subject)
   Data: MatchedData

     AIC      BIC   logLik deviance df.resid 
   753.0    784.6   -368.5    737.0      374 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9964 -0.6378 -0.0560  0.6187  4.3088 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Subject  (Intercept) 0.448218 0.6695        
          Session     0.009999 0.1000   -0.78
 Residual             0.306934 0.5540        
Number of obs: 382, groups:  Subject, 44

Fixed effects:
                       Estimate Std. Error t value
(Intercept)             2.99655    0.13544  22.125
Session                -0.31912    0.03244  -9.836
LearnedControl          0.03582    0.12713   0.282
Session:LearnedControl  0.05752    0.03995   1.440

Correlation of Fixed Effects:
            (Intr) Sessin LrndCn
Session     -0.797              
LearndCntrl -0.469  0.551       
Sssn:LrndCn  0.420 -0.616 -0.895

The results are in the format which I am used to from a repeated-measures ANOVA with two main effects and one interaction term.

Any help with my above questions would be extremely appreciated.

$\endgroup$
  • $\begingroup$ Google "polynomial contrasts". The decision between ordered factor and integer is the decision between polynomial and strictly linear relationship. $\endgroup$ – Roland Feb 14 '18 at 12:57
  • $\begingroup$ Hi Roland, Thanks for your comment! I took a look and I think I understand what you mean. When I changed the contrasts with contrasts(MatchedData$Session) <- contr.treatment(base = 1, n = 5) it then gave the exact same results in the anova() as converting the Session data into a regular factor. But I am still unsure whether I should be using a factor or an integer for time. Any help would be appreciated. $\endgroup$ – V Mileva Feb 21 '18 at 11:44
  • $\begingroup$ For mixed models, you can treat the time-variable as numeric, and you can indicate different spaces between time-points with related values, e.g. time = 1, 3, and 12 for baseline, after 3 and after 12 months or so. Usually, you may get warnings when treating time as factor, because you then have equal or more groups of random effects as observations in your model. This doesn't happen if time is numeric. $\endgroup$ – Daniel Aug 12 '18 at 10:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.