Testing zero variance component in a mixed model

Question

I find this post relevant. The mock data-set is constructed the same way as in the post but with a different seed to produce zero variance in the mixed model. Think this as a three way anova, with A as temperature factor, B as day factor and C as replicate. The output of fit1 gives a SD of 0.455 for subject, and anova with p-value 0.02483 suggests a significant non-zero SD. I did some research and find this chi-square test is a conservative test considering it is on the boundary, so p-value should have been even lower. My question is whenever we see zero SD for some factors (in fit1, SD for factors A and B are zero), does this imply that the variance estimates for other factors are not reliable? The residual SD is correctly estimated to be 1 though.

> set.seed(8)
> d <- data.frame(
+   Y = rnorm(96),
+   subject = factor(rep(1:12, 4)),
+   A = factor(rep(1:2, each=24)),
+   B = factor(rep(rep(1:2, each=12))),
+   C = factor(rep(rep(1:2, each=48))))
> 
> fit1 <- lmer(Y ~  (1|subject) + (1|A) + (1|B), data=d)
> 
> fit2 <- lmer(Y ~   (1|A) + (1|B), data=d)
> 
> fit1
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: Y ~ (1 | subject) + (1 | A) + (1 | B)
   Data: d
REML criterion at convergence: 286.2724
Random effects:
 Groups   Name        Std.Dev.
 subject  (Intercept) 0.455   
 A        (Intercept) 0.000   
 B        (Intercept) 0.000   
 Residual             1.008   
Number of obs: 96, groups:  subject, 12; A, 2; B, 2
Fixed Effects:
(Intercept)  
   -0.09281  
> 
> anova(fit1,fit2)
refitting model(s) with ML (instead of REML)
Data: d
Models:
fit2: Y ~ (1 | A) + (1 | B)
fit1: Y ~ (1 | subject) + (1 | A) + (1 | B)
     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
fit2  4 297.52 307.78 -144.76   289.52                           
fit1  5 294.48 307.31 -142.24   284.48 5.0357      1    0.02483 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Is there a reason why you think that 0 SD for one random effect variance make another random effect variance unreliable? — Erik Ruzek, Feb 21 '20 at 03:07
No. I guess I am wrong that this significant non-zero result is actually a small probability event based on the simulations I did after I post this. — Statisfun, Feb 21 '20 at 16:33

score 2 · Answer 1 · answered Feb 22 '20 at 19:45

Note that the reported p-value from the anova() method when testing for a variance component is not totally correct. In particular, classic hypothesis testing assumes that the null hypothesis is an interior point of the parameter space. However, this does not hold in the case you test a variance component. This is because the null hypothesis that the variance is zero is on the boundary of the parameter space. In some simple settings, not the one you are considering, it has been shown that the asymptotic distribution is a mixture of $\chi^2$ distributions.

With regard to your question, the fact that you are fitting a model for which the true values of the variance components is zero could affect the numerical stability and hence the behavior of the likelihood ratio test.

Testing zero variance component in a mixed model

1 Answers1