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 ?
