Extract the different components of variance in a linear mixed model in R

Question

Consider a mixed model as follows.

library(lme4)
# Load data
data <- structure(list(blk = c(1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3L),
                       gent = c(1, 2, 3, 4, 7, 11, 12, 1, 2, 3, 4, 5, 9, 1, 2, 3, 4, 8, 6, 10L),
                       yld = c(83, 77, 78, 78, 70, 75, 74, 79, 81, 81, 91, 79, 78, 92, 79, 87, 81, 96, 89, 82L),
                       syld = c(250, 240, 268, 287, 226, 395, 450, 260, 220, 237, 227, 281, 311, 258, 224, 238, 278, 347, 300, 289L)),
                  .Names = c("blk", "gent", "yld", "syld"), class = "data.frame", row.names = c(NA, -20L))
data$blk <- as.factor(data$blk)
data$gent <- as.factor(data$gent)

The data is unbalanced.

# Mixed effect model
frmla <- "syld ~ 1 + gent + (1|blk)"
library(lme4)
model <- lmer(formula(frmla), data = data)

model
Linear mixed model fit by REML ['merModLmerTest']
Formula: syld ~ 1 + gent + (1 | blk)
   Data: data
REML criterion at convergence: 73.9572
Random effects:
 Groups   Name        Std.Dev.
 blk      (Intercept)  9.385  
 Residual             16.919  
Number of obs: 20, groups:  blk, 3
Fixed Effects:
(Intercept)        gent2        gent3        gent4        gent5        gent6        gent7        gent8        gent9  
    256.000      -28.000       -8.333        8.000       32.127       43.678      -36.805       90.678       62.127  
     gent10       gent11       gent12  
     32.678      132.195      187.195  

Primarily I want to compare the gent levels by LS means.

library("lmerTest")
lsmeans(model)
Least Squares Means table:
         gent Estimate Standard Error   DF t-value Lower CI Upper CI p-value    
gent  1   1.0    256.0           11.2  6.9    22.9      229      283  <2e-16 ***
gent  2   5.0    228.0           11.2  6.9    20.4      201      255  <2e-16 ***
gent  3   6.0    247.7           11.2  6.9    22.2      221      274  <2e-16 ***
gent  4   7.0    264.0           11.2  6.9    23.6      237      291  <2e-16 ***
gent  5   8.0    288.1           18.5  8.0    15.6      245      331  <2e-16 ***
gent  6   9.0    299.7           18.5  8.0    16.2      257      342  <2e-16 ***
gent  7  10.0    219.2           18.5  8.0    11.8      177      262  <2e-16 ***
gent  8  11.0    346.7           18.5  8.0    18.8      304      389  <2e-16 ***
gent  9  12.0    318.1           18.5  8.0    17.2      275      361  <2e-16 ***
gent  10  2.0    288.7           18.5  8.0    15.6      246      331  <2e-16 ***
gent  11  3.0    388.2           18.5  8.0    21.0      346      431  <2e-16 ***
gent  12  4.0    443.2           18.5  8.0    24.0      401      486  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

In addition I am interested in variance partitioning.

The variance component due to random effect and residual can be estimated as follows.

VCrandom <- VarCorr(model)
print(VCrandom, comp = "Variance")
 Groups   Name        Variance
 blk      (Intercept)  88.083 
 Residual             286.250

How to partition the total variance into components due to each of the factors gent and blk along with the residual ? Something similar to the output given by PROC MIXED of SAS, where MSE is computed even when estimation is by ML or REML instead of least squares.

Should I treat the fixed effect as random just for the purpouse of getting variance component ?

frmla2 <- "syld ~ 1 + (1|gent) + (1|blk)"
model2 <- lmer(formula(frmla2), data = data)
model2

VCrandom2 <- VarCorr(model2)
print(VCrandom2, comp = "Variance")
 Groups   Name        Variance
 gent     (Intercept) 4152.08 
 blk      (Intercept)  116.11 
 Residual              274.92 

If there is no random effect, variance components can be estimated using the least squares approach (ANOVA, Sum of squares, MSE).

The package mixlm has provision for variance partitioning using SS in case of mixed models.

library(mixlm)

mixlm <- lm(syld ~ 1 + r(gent) + r(blk), data)

Anova(mixlm, type="III")

Analysis of variance (unrestricted model)
Response: syld
          Mean Sq   Sum Sq Df F value Pr(>F)
gent      5360.49 58965.36 11   18.73 0.0009
blk        638.58  1277.17  2    2.23 0.1886
Residuals  286.25  1717.50  6       -      -

            Err.term(s) Err.df VC(SS)
1 gent              (3)      6 3044.5
2 blk               (3)      6   52.8
3 Residuals           -      -  286.3
(VC = variance component)

               Expected mean squares
gent      (3) + 1.66666666666667 (1)
blk       (3) + 6.66666666666667 (2)
Residuals (3)                       

WARNING: Unbalanced data may lead to poor estimates

The estimates are different

# Total variance
var(data$syld)

|source   |  model1|  model2|  mixlm|
|:--------|-------:|-------:|------:|
|gent     |      NA| 4152.08| 3044.5|
|blk      |  88.083|  116.11|   52.8|
|Residual | 286.250|  274.92|  286.3|

Can fixed effect variance be extracted using predict function as suggested here In R: How to extract the different components of variance in a linear mixed model! ?

var(predict(model))

Which is the most appropriate method compatible with (RE)ML estimates in lme4 ?

Answer 1

This isn't exactly an answer, but whenever I see a question about "explained variance" in mixed models, I always think of this email from Douglas Bates, the original author of lme4 and co-author of nlme, on the R-Sig-ME mailing list on Feb 26, 2010, in response to:

I have written some code to implement Gelman & Pardoe's Rsq for an lmer object. It gives some believable results, but it's difficult to be confident because of the translation from Bayesian into frequentist paradigms.

If anyone is interested then I'd be really happy to discuss this off-list and share/develop the code.

I think it is useful to have this email here on CV. This is Bates' reply:

Assuming that one wants to define an R^2 measure, I think an argument could be made for treating the penalized residual sum of squares from a linear mixed model in the same way that we consider the residual sum of squares from a linear model. Or one could use just the residual sum of squares without the penalty or the minimum residual sum of squares obtainable from a given set of terms, which corresponds to an infinite precision matrix. I don't know, really. It depends on what you are trying to characterize.

In other words, what's the purpose? What aspect of the R^2 for a linear model are you trying to generalize?

I'm sorry if I sound argumentative but discussions like this sometimes frustrate me. A linear mixed model does not behave exactly like a linear model without random effects so a measure that may be appropriate for the linear model does not necessarily generalize. I'm not saying that this is the case but if the request is "I don't care what the number means or if indeed it means anything at all, just give me a number I can report", that's not the style of statistics I practice.

I regard Bill Venables' wonderful unpublished paper "Exegeses on Linear Models" (just put the name in a search engine to find a copy - there is only one paper with "Exegeses" and "Linear Models" in the title) as required reading for statisticians. As Bill emphasizes in that paper, statistics is not just a collection of formulas (many of which are based on approximations). It's about models and comparing how well different models fit the observed data. If we start with a formula and only ask ourselves "How do we generalize this formula?" we're missing the point. We should start at the model.

In a linear model the R^2 statistic is a dimensionless comparison of the quality of the current model fit, as measured by the residual sum of squares, to the fit one would obtain from a trivial model. When the current model can be shown to contain a model with an intercept term only (and whose coefficient will be estimated by the mean response) then that model fit is the trivial model. Otherwise the trivial model is a prediction of zero for each response. We know that the trivial model will produce a greater residual sum of squares than the current model fit because the models are nested. The R^2 is the proportion of variability not accounted for by the trivial model but accounted for by the current model (my apologies to my grammar teachers for having juxtaposed prepositions).

The interesting point there is that when you think of the relationships between models you can determine how you handle the case of a model that does not have an intercept term. If you start from the formula instead you can end up calculating a negative R^2 because you compare models that are not nested. Such nonsensical results are often reported. (I think it was the Mathematica documentation that gave a careful explanation of why you get a negative R^2 instead of recognizing that the formula they were using did not apply in certain cases.)

It may be that there is a sensible measure of the quality of fit from a linear mixed model that generalizes the R^2 from a linear model. I don't see an obvious candidate but I will freely admit that I haven't thought much about the problem. I would ask others who are thinking about this to consider both the "what" and the "why". George Mallory's justification of "because it's there" for attempting to climb Everest is perhaps a good justification for such endeavors (Mallory may have questioned his rationale as he lay freezing to death on the mountain). I don't think it is a good justification for manipulating formulas.