More REML exercise

August 25, 2013
By

(This article was first published on Wiekvoet, and kindly contributed to R-bloggers)

Last week I tried exercise 1 of the SAS(R) proc mixed with R libraries lme4 and MCMCglm. So this week I aimed for exercise 2 but ended up redoing exercise 1 with nlme.
Exercise 2 results gave me problems with library lme4 and latter parts of the exercise had special correlation structures which seem not feasible with lme4. So library nlme came in the picture. However, while reading up on it, I concluded that Pinheiro and Bates (2000) used this same data in their book. I am not about to reproduce the book in a blog, so I have to check in more detail what I can blog. So, for now I am redoing exercise 1 with library nlme.

Running the analysis

The first part of the analysis is actually clear sailing. The random effect specification is quite different from the other libraries, but it was pretty simple.
m1 <- lme(Y ~ A + B + A*B,random=~1|Block/A,data=sp)
summary(m1)
Linear mixed-effects model fit by REML
Data: sp
AIC      BIC    logLik
137.7618 145.7752 -59.88092

Random effects:
Formula: ~1 | Block
(Intercept)
StdDev:    7.899104

Formula: ~1 | A %in% Block
(Intercept) Residual
StdDev:    3.921982 3.059593

Fixed effects: Y ~ A + B + A * B
Value Std.Error DF   t-value p-value
(Intercept)  37.00  4.667411  9  7.927307  0.0000
A2            1.00  3.517318  6  0.284308  0.7857
A3          -11.00  3.517318  6 -3.127383  0.0204
B2           -8.25  2.163459  9 -3.813338  0.0041
A2:B2         0.50  3.059593  9  0.163420  0.8738
A3:B2         7.75  3.059593  9  2.533016  0.0321
Correlation:
(Intr) A2     A3     B2     A2:B2
A2    -0.377
A3    -0.377  0.500
B2    -0.232  0.308  0.308
A2:B2  0.164 -0.435 -0.217 -0.707
A3:B2  0.164 -0.217 -0.435 -0.707  0.500

Standardized Within-Group Residuals:
Min           Q1          Med           Q3          Max
-1.165347655 -0.544427403 -0.007998877  0.526246640  1.473815718

Number of Observations: 24
Number of Groups:
Block A %in% Block
4           12

Covariance

VarCorr(m1)

Variance     StdDev
Block =     pdLogChol(1)
(Intercept) 62.395850    7.899104
A =         pdLogChol(1)
(Intercept) 15.381944    3.921982
Residual     9.361111    3.059593

Likelihood

-2*logLik(m1)
‘log Lik.’ 119.7618 (df=9)

ANOVA

anova(m1)

numDF denDF  F-value p-value

(Intercept)     1     9 55.34418  <.0001
A               2     6  4.06957  0.0764
B               1     9 19.38873  0.0017
A:B             2     9  4.01929  0.0566

Predict A level 1

This is a problem. I went through the book, browsed the web, only to conclude that getting the mean is not difficult, but getting the associated standard error is not supported. Not because of a hate for standard errors, but because it is an intrinsically difficult problem. I seem to remember reading years ago that lme4 had the MCMC functionality precisely for this reason. Which does not help me if I try stubbornly to get the standard error anyway.

Not a solution

After reading http://glmm.wikidot.com/faq, https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q2/003585.html and  http://stackoverflow.com/questions/14358811/extract-prediction-band-from-lme-fit I drafted this, which does not really work for my problem as it does not provide level A1 without specifying a level for B.
newdat <- unique(sp[sp\$A==1,c(‘A’,’B’,’Block’)])
newdat\$pred <- predict(m1, newdat, level = 0)
Designmat <- model.matrix(eval(eval(m1\$call\$fixed)[-2]), newdat[-3])
predvar <- diag(Designmat %*% m1\$varFix %*% t(Designmat))
newdat\$SE <- sqrt(predvar)
newdat\$SE2 <- sqrt(predvar+m1\$sigma^2)
newdat
A B Block  pred       SE      SE2
1  1 1     1 37.00 4.667411 5.580846
4  1 1     2 37.00 4.667411 5.580846
7  1 1     3 37.00 4.667411 5.580846
10 1 1     4 37.00 4.667411 5.580846
13 1 2     1 28.75 4.667411 5.580846
16 1 2     2 28.75 4.667411 5.580846
19 1 2     3 28.75 4.667411 5.580846
22 1 2     4 28.75 4.667411 5.580846

Combining parameters

Since level A1 is just a weighted sum of known parameters with their own variances. But it is so unclear which level of variance is used. It does provide the 4.54 which SAS has labelled ‘A1 mean broad’.
fixef(m1)%*%c(1,0,0,.5,0,0)
[,1]
[1,] 32.875
sqrt(vcov(m1)[1,1]+ .25 *vcov(m1)[4,4]+vcov(m1)[1,4])
[1] 4.540329

Contrasts

With a different parameterization of the model and suitable contrasts. This also provides the ‘A1 mean broad’ of SAS.
contrasts(sp\$B) <- contr.sum(levels(sp\$B))
m2 <- lme(Y ~ A + B + A*B -1 ,random=~1|Block/A,data=sp)
summary(m2)

Linear mixed-effects model fit by REML

Data: sp
AIC      BIC    logLik
141.9207 149.9341 -61.96036
Random effects:
Formula: ~1 | Block
(Intercept)
StdDev:    7.899104
Formula: ~1 | A %in% Block
(Intercept) Residual
StdDev:    3.921982 3.059593
Fixed effects: Y ~ A + B + A * B – 1
Value Std.Error DF   t-value p-value
A1    32.875  4.540329  6  7.240665  0.0004
A2    34.125  4.540329  6  7.515975  0.0003
A3    25.750  4.540329  6  5.671395  0.0013
B1     4.125  1.081730 10  3.813338  0.0034
A2:B1 -0.250  1.529797 10 -0.163420  0.8734
A3:B1 -3.875  1.529797 10 -2.533016  0.0297
Correlation:
A1     A2     A3     B1     A2:B1
A2     0.757
A3     0.757  0.757
B1     0.000  0.000  0.000
A2:B1  0.000  0.000  0.000 -0.707
A3:B1  0.000  0.000  0.000 -0.707  0.500
Standardized Within-Group Residuals:
Min           Q1          Med           Q3          Max
-1.165347654 -0.544427403 -0.007998877  0.526246640  1.473815718
Number of Observations: 24
Number of Groups:
Block A %in% Block
4           12

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...