Solved – Testing the variance component in a mixed effects model

asymptoticshypothesis testingmixed modelp-value

Say $y=X\beta+ Zu +\epsilon$ is our mixed effects model where $u=(u_1,..,u_r)$ and $u_{j} \stackrel{i.i.d.}{\sim} N(0, \sigma^2_{a})$ for $j=1,…,r$ and $\epsilon=(\epsilon_1,…,\epsilon_n)$ are i.i.d. $N(0, \sigma^2_{b})$, furthermore $\epsilon_j$ and $u_i$ are also assumed to be independent for all $j$'s and all $i$'s.

I am interested in testing the hypothesis $H_{0}:\sigma_{a}^2=0$ vs $H_{1}: H_{0}$ is not true. The ${\bf lmer}$ package in R does provide a bootstrap p-value for this problem, but my question is: Is there is any non-bootstrap way of obtaining this p-value (asymptotic or otherwise)?

Best Answer

This is usually done with maximum likelihood ratio between original model and a model omitting the variance coefficient to be estimate (random intercept/random slope/random co-variance between slope and intercept).

A good example is in these tutorials:

When model has more than one random coefficient: http://www.bodowinter.com/tutorial/bw_LME_tutorial.pdf (p.12)
When model has one random coefficient: http://www.stat.wisc.edu/~ane/st572/notes/lec21.pdf (p.13)

Sample R code:

> model1 = lmer(resp ˜ fixed1 + (1 | random1))
> model2 = lm(resp ˜ fixed1)
> chi2 = -2*logLik(model2, REML=T) +2*logLik(model1, REML=T)
> chi2
[1] 5.011
> pchisq(chi2, df=1, lower.tail=F)
[1] 0.02518675

Related Solutions

Solved – In practice how is the random effects covariance matrix calculated in a mixed effects model

The Goldstein .pdf @probabilityislogic linked is a great document. Here's a list of some references that discuss your particular question:

Harville, 1976: Extension of the Gauss-Markov Theorem to include the estimation of random effects.

Harville, 1977: Maximum likelihood approaches to variance component estimation and to related problems.

Laird and Ware, 1982: Random-effects models for longitudinal data.

McCulloch, 1997: Maximum likelihood algorithms for generalized linear mixed models.

The SAS User Guide excerpt for the MIXED procedure has some great information about covariance estimation and many more sources (starting on page 3968).

There are numerous quality textbooks on longitudinal/repeated measures data analysis, but here's one that goes into some detail about implementation in R (from the authors of lme4 and nlme):

Pinheiro and Bates, 2000: Mixed-Effects Models in S and S-PLUS.

EDIT: One more relevant paper: Lindstrom and Bates, 1988: Newton-Raphson and EM Algorithms for linear mixed-effects models for repeated-measures data.

EDIT 2: And another: Jennrich and Schluchter, 1986: Unbalanced Repeated-Measures Models with Structured Covariance Matrices.

Solved – Mixed model specification with nlme in R

If I understand your question correctly, you can specify your model with nested random effects like this:

fit.1 <- lme(Change ~ Dose*Time, random=~1|ID/Dose, data=mydata)

To specify the covariance structure, e.g. a simple compound symmetry form, try this:

fit.2 <- lme(Change ~ Dose*Time, random=~1|ID/Dose, data=mydata, cor=corCompSymm())

To look at the estimated parameters try:

summary(fit.1)

To get all estimated coefficients try:

coef(fit.1)

To get the p-values then use:

anova(fit.1)

Notice that if you need to specify the covariance structure of the residuals, you'll have to use nlme as although lme4 (i.e. the lmer function) is a more advanced package, currently it does not support that feature.

Best Answer

Related Solutions

Solved – In practice how is the random effects covariance matrix calculated in a mixed effects model

Solved – Mixed model specification with nlme in R

Related Question