Solved – How to estimate variance components with lmer for models with random effects and compare them with lme results

anovalme4-nlmervariance

I performed an experiment where I raised different families coming from two different source populations. Each family was assigned one of two treatments. After the experiment I measured several traits on each individual.
To test for an effect of either treatment or source as well as their interaction, I used a linear mixed effect model with family as random factor, i.e.

lme(fixed=Trait~Treatment*Source,random=~1|Family,method="ML")

so far so good,
Now I have to calculate the relative variance components, i.e. the percentage of variation that is explained by either treatment or source as well as the interaction.

Without a random effect, I could easily use the sums of squares (SS) to calculate the variance explained by each factor. But for a mixed model (with ML estimation), there are no SS, hence I thought I could use Treatment and Source as random effects too to estimate the variance, i.e.

lme(fixed=Trait~1,random=~(Treatment*Source)|Family, method="REML")

However, in some cases, lme does not converge, hence I used lmer from the lme4 package:

lmer(Trait~1+(Treatment*Source|Family),data=DATA)

Where I extract the variances from the model using the summary function:

model<-lmer(Trait~1+(Treatment*Source|Family),data=regrexpdat)
results<-VarCorr(model)
variances<-results[,3]

I get the same values as with the VarCorr function. I use then these values to calculate the actual percentage of variation taking the sum as the total variation.

Where I am struggling is with the interpretation of the results from the initial lme model (with treatment and source as fixed effects) and the random model to estimate the variance components (with treatment and source as random effect). I find in most cases that the percentage of variance explained by each factor does not correspond to the significance of the fixed effect.

For example for the trait HD,
The initial lme suggests a tendency for the interaction as well as a significance for Treatment. Using a backward procedure, I find that Treatment has a close to significant tendency. However, estimating variance components, I find that Source has the highest variance, making up to 26.7% of the total variance.

The lme:

anova(lme(fixed=HD~as.factor(Treatment)*as.factor(Source),random=~1|as.factor(Family),method="ML",data=test),type="m")
                                      numDF denDF  F-value p-value
(Intercept)                                1   426 0.044523  0.8330
as.factor(Treatment)                       1   426 5.935189  0.0153
as.factor(Source)                          1    11 0.042662  0.8401
as.factor(Treatment):as.factor(Source)     1   426 3.754112  0.0533

And the lmer:

summary(lmer(HD~1+(as.factor(Treatment)*as.factor(Source)|Family),data=regrexpdat))
Linear mixed model fit by REML 
Formula: HD ~ 1 + (as.factor(Treatment) * as.factor(Source) | Family) 
   Data: regrexpdat 
    AIC    BIC logLik deviance REMLdev
 -103.5 -54.43  63.75   -132.5  -127.5
Random effects:
 Groups   Name                                      Variance  Std.Dev. Corr                 
 Family   (Intercept)                               0.0113276 0.106431                      
          as.factor(Treatment)                      0.0063710 0.079819  0.405               
          as.factor(Source)                         0.0235294 0.153393 -0.134 -0.157        
          as.factor(Treatment)L:as.factor(Source)   0.0076353 0.087380 -0.578 -0.589 -0.585 
 Residual                                           0.0394610 0.198648                      
Number of obs: 441, groups: Family, 13

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.02740    0.03237  -0.846

Hence my question is, is it correct what I am doing? Or should I use another way to estimate the amount of variance explained by each factor (i.e. Treatment, Source and their interaction). For example, would the effect sizes be a more appropriate way to go?

Best Answer

One common way to determine the relative contribution of each factor to a model is to remove the factor and compare the relative likelihood with something like a chi-squared test:

pchisq(logLik(model1) - logLik(model2), 1)

As the way that likelihoods are calculated between functions may be slightly different, I typically will only compare them between the same method.

Related Solutions

Maximum Likelihood – Using REML or ML to Compare Mixed Effects Models With Different Fixed Effects

Zuur et al., and Faraway (from @janhove's comment above) are right; using likelihood-based methods (including AIC) to compare two models with different fixed effects that are fitted by REML will generally lead to nonsense.

Faraway (2006) Extending the linear model with R (p. 156):

The reason is that REML estimates the random effects by considering linear combinations of the data that remove the fixed effects. If these fixed effects are changed, the likelihoods of the two models will not be directly comparable

These two questions discuss the issue further: Allowed comparisons of mixed effects models (random effects primarily) ; REML vs ML stepAIC

Solved – Different results obtained with lmer() and aov() for three-way repeated-measures experiment

Using Jake Westfall's EMSfunction documented HERE, I estimated the variance components from the mixed model. The random variances above that were estimated to be 0 are indeed negative as derived by EMS.

cbind(rev(ans[c(grep("e",names(ans)),grep("s",names(ans)))])/
    c(1,2,2,2,4,4,4,5,1))


s        69.838600
a:s       7.711213
b:s      -1.423696
c:s      -1.098723
a:b:s    -1.498620
a:c:s    -1.212590
b:c:s    -1.244919
a:b:c:s  -2.650683
e       922.447967

When running lmer.4 on the aggregated data-set, lmergives identical results (except for the 3-way interaction, but it is still very similar).

> lmer.4 <- lmer(rating ~ timepoint_n*att_cond_n*gng_cond_n +(1|subjectID) + (0 + timepoint_n | subjectID) 
                                  (0 + att_cond_n | subjectID) + (0+gng_cond_n |subjectID) +  (0+timepoint_n:att_cond_n|subjectID) + 
                                  (0+timepoint_n:gng_cond_n|subjectID) + (0 + att_cond_n:gng_cond_n|subjectID), 
                                 data = d_aggr, control=lmerControl(optCtrl=list(maxfun=1e9)))
> summary(lmer.4)
Linear mixed model fit by REML 
t-tests use  Satterthwaite approximations to degrees of freedom ['lmerMod']
Formula: rating ~ timepoint_n * att_cond_n * gng_cond_n + (1 | subjectID) +  
  (0 + timepoint_n | subjectID) + (0 + att_cond_n | subjectID) +  
  (0 + gng_cond_n | subjectID) + (0 + timepoint_n:att_cond_n |  
                                    subjectID) + (0 + timepoint_n:gng_cond_n | subjectID) + (0 +      att_cond_n:gng_cond_n | subjectID)
Data: d_aggr
Control: lmerControl(optCtrl = list(maxfun = 1e+09))

REML criterion at convergence: 2438.7

Scaled residuals: 
  Min      1Q  Median      3Q     Max 
-4.0726 -0.1788 -0.0138  0.2457  4.3407 

Random effects:
  Groups      Name                   Variance  Std.Dev. 
subjectID   (Intercept)            2.850e+02 1.688e+01
subjectID.1 timepoint_n            3.651e+01 6.042e+00
subjectID.2 att_cond_n             2.257e-11 4.750e-06
subjectID.3 gng_cond_n             1.269e+00 1.126e+00
subjectID.4 timepoint_n:att_cond_n 0.000e+00 0.000e+00
subjectID.5 timepoint_n:gng_cond_n 8.132e-01 9.018e-01
subjectID.6 att_cond_n:gng_cond_n  6.839e-01 8.270e-01
Residual                           4.694e+01 6.851e+00
Number of obs: 328, groups:  subjectID, 41

Fixed effects:
                                   Estimate   Std. Error      df  t value Pr(>|t|)    
(Intercept)                       140.88293    2.66360  40.00000  52.892  < 2e-16 ***
timepoint_n                        -6.92561    1.01664  40.00000  -6.812 3.43e-08 ***
att_cond_n                         -0.13354    0.37828 120.00000  -0.353  0.72470    
gng_cond_n                          1.35854    0.41718  40.00000   3.256  0.00230 ** 
timepoint_n:att_cond_n             -0.02744    0.37828 120.00000  -0.073  0.94230    
timepoint_n:gng_cond_n              1.36220    0.40365  40.00000   3.375  0.00165 ** 
att_cond_n:gng_cond_n               1.00549    0.39972  40.00000   2.515  0.01601 *  
timepoint_n:att_cond_n:gng_cond_n   0.98963    0.37828 120.00000   2.616  0.01004 *  
  ---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) tmpnt_ att_c_ gng_c_ tmpnt_n:t__ tmpnt_n:g__ a__:__
timepoint_n 0.000                                                     
att_cond_n  0.000  0.000                                              
gng_cond_n  0.000  0.000  0.000                                       
tmpnt_n:t__ 0.000  0.000  0.000  0.000                                
tmpnt_n:g__ 0.000  0.000  0.000  0.000  0.000                         
att_cnd_:__ 0.000  0.000  0.000  0.000  0.000       0.000             
tmpn_:__:__ 0.000  0.000  0.000  0.000  0.000       0.000       0.000

Therefore, it seems that as lmer cannot achieve the fit of aov by not being able/allowed to estimate negative variance components, the results differ, while the difference is gone when running the lmer representation of the model on the aggregated data, eliminating the negatively estimated variance components.