A zero estimate of a variance possibly indicates the model is wrong.

Consider e.g. the following dataset :
> dat
   tube position          y
1     1   bottom  0.6743106
2     1   bottom  1.6816419
3     1   bottom -0.2978006
4     1   bottom  3.2026162
5     1      top 12.7948831
6     1      top  4.9543361
7     1      top 10.2241050
8     1      top  5.0544860
9     2   bottom  6.8313433
10    2   bottom  6.4114389
11    2   bottom  8.2796931
12    2   bottom  7.2528927
13    2      top  1.0694828
14    2      top  3.8067290
15    2      top  0.8501896
16    2      top -0.1316989
17    3   bottom  4.7749865
18    3   bottom  1.6121819
19    3   bottom  2.2984276
20    3   bottom -0.7272648
21    3      top  4.9126798
22    3      top  3.6382901
23    3      top  5.0407385
24    3      top  4.4799782


> str(dat)
'data.frame':   24 obs. of  3 variables:
 $ tube    : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 2 2 ...
 $ position: Factor w/ 2 levels "bottom","top": 1 1 1 1 2 2 2 2 1 1 ...
 $ y       : num  0.674 1.682 -0.298 3.203 12.795 ...




The covariance between the means of the "top" level and the "bottom" level of the tube factor is negative:

> means.bottom <- with(subset(dat, subset= position=="bottom"), 
+ sapply(split(y, tube), mean)
+ )
> means.top <- with(subset(dat, subset= position=="top"), 
+ sapply(split(y, tube), mean)
+ ) 
> cov(means.bottom, means.top)
[1] -9.84539


Therefore we get a negative estimate of the correlation with the following lme() model:

> lme(y ~ position, data=dat, random= list(tube = pdCompSymm(~ 0+position )))
Linear mixed-effects model fit by REML
  Data: dat 
  Log-restricted-likelihood: -53.07383
  Fixed: y ~ position 
(Intercept) positiontop 
   3.499539    1.224978 

Random effects:
 Formula: ~0 + position | tube
 Structure: Compound Symmetry
               StdDev   Corr  
positionbottom 3.157148       
positiontop    3.157148 -0.988
Residual       2.100170       

Number of Observations: 24
Number of Groups: 3 

> # covariance estimate (the slight difference between the result and -9.84539is due to the rounding):
> 3.157148^2*-0.988
[1] -9.847972


Now,you can check that the above lme() model is equivalent to the following lmer() model, but with different interpretations of the variance components:

> lmer(y ~ position + (1|tube) + (1|tube:position), data=dat)
Linear mixed model fit by REML 
Formula: y ~ position + (1 | tube) + (1 | tube:position) 
   Data: dat 
   AIC   BIC logLik deviance REMLdev
 119.3 125.2 -54.64    115.5   109.3
Random effects:
 Groups        Name        Variance   Std.Dev.  
 tube:position (Intercept) 9.9678e+00 3.1572e+00
 tube          (Intercept) 3.0907e-20 1.7580e-10
 Residual                  4.4107e+00 2.1002e+00
Number of obs: 24, groups: tube:position, 6; tube, 3




I claim that these two models are equivalent because they yield the same marginal models, but with one caveat: the covariance in the first model corresponds to the ratio of the tube variance over the sum of the tube variance and the interaction variance in the second model. Therefore the covariance is constrained to be positive in the second model, and this is why we get a (almost) zero estimate of the tube variance. Using PROC MIXED in SAS we would get the same result with the NOBOUND option in the second model.


We can check that the two models yield exactly the same results when the covariance is positive:

> dat$y[dat$position=="top"] <- -dat$y[dat$position=="top"]
> lmer(y ~ position + (1|tube) + (1|tube:position), data=dat)
Linear mixed model fit by REML 
Formula: y ~ position + (1 | tube) + (1 | tube:position) 
   Data: dat 
   AIC BIC logLik deviance REMLdev
 116.1 122 -53.07    110.8   106.1
Random effects:
 Groups        Name        Variance Std.Dev.
 tube:position (Intercept) 0.12224  0.34963 
 tube          (Intercept) 9.84539  3.13774 
 Residual                  4.41071  2.10017 
Number of obs: 24, groups: tube:position, 6; tube, 3

> # correlation estimate:
> 9.84539/(9.84539+0.12224)
[1] 0.9877363



> lme(y ~ position, data=dat, random= list(tube = pdCompSymm(~ 0+position )))
Linear mixed-effects model fit by REML
  Data: dat 
  Log-restricted-likelihood: -53.07383
  Fixed: y ~ position 
(Intercept) positiontop 
   3.499539   -8.224056 

Random effects:
 Formula: ~0 + position | tube
 Structure: Compound Symmetry
               StdDev   Corr 
positionbottom 3.157148      
positiontop    3.157148 0.988
Residual       2.100170      

Number of Observations: 24
Number of Groups: 3 



If you want to reproduce my example, my dataset was simulated:


set.seed(666) 
I <- 2 # positions (fixed factor)
J <- 3 # tubes (random factor)
K <- 4 # repeats 
n <- I*J*K
tube <- rep(1:J, each=I)
position <- rep(c("bottom","top"), times=J)
Mu_i <- 3*(1:I)
    tube <- rep(tube, each=K)
    position <- rep(position, each=K)
    sigmaw <- 2 

Mu_ij <- c(t(rmvnorm(J, mean=Mu_i, sigma=3*diag(I)-1)) )  
    Mu_ij <- rep(Mu_ij, each=K)
    dat <- data.frame(tube, position)
    dat$y <- rnorm(n, Mu_ij, sigmaw)
    dat$tube <- factor(dat$tube)



Kind regards,
Stéphane



________________________________
 De : Douglas Bates <bates@stat.wisc.edu>

Cc : "r-sig-mixed-models@r-project.org" <r-sig-mixed-models@r-project.org> 
Envoyé le : Vendredi 11 mai 2012 21h14
Objet : Re: [R-sig-ME] Variance estimate at 0 with lmer - not with lme - issue or rounding ?



> Hello all,
>
> I'm getting puzzled by the following:
>
> I'm fitting a simple hierarchical model with individual identity as random factor. Fitted with lmer (lme4 package) the result look like this:
>
> Formula: response ~ some fixed effects + (1 | indiv)
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  indiv     (Intercept) 0.0000   0.0000
>  Residual             1.4802   1.2166
> Number of obs: 116, groups: name, 69
>
> The variance estimate for individuals is a perfect 0 (which is always suspicious...and maybe an indication of problem in data distribution across effects)

No, an estimate of zero is not suspicious.  It is simply an indication
that the variability between individuals is not significantly larger
than what one would expect from the random variability in the
response.  (By the way, is that really the output from the model fit
because the name "name" in the last line should be "indiv".)

> I then fitted the model with lme (in nlme package) and got:
> Formula: lme(fixed=response ~ some fixed effects,random= ~ 1 | indiv)
> Random effects:
>  Formula: ~1 | indiv
>          (Intercept) Residual
> StdDev: 9.796623e-05 1.2166
> which I agree is a (very) small variance estimate but is less suspicious.

One of the internal differences between lme and lmer is that the
variance component estimates in lmer can go to zero gracefully whereas
in lme they can't because lme works on the scale of the logarithm of
the relative standard deviation.  It took a long time to figure out
how to solve the "going to zero gracefully" problem.

> I'm hitting this "issue" regularly with this dataset, as under some combination of fixed effect variables the variance estimate for individuals suddenly drops to a perfect 0.
>
> So to help me figure out if I should look for a problem in the distribution of data across variables, my question is:
>
> Is the difference in variance estimated by lmer and lme function such as the one observed here likely just caused by a rounding issue (and thus I should assume that with lmer the estimate is correct and simply 0), or is it sign of a real "problem" in my dataset?

Just a rounding issue.

> As ever, thanks for your help (here and in other posts that provide answers so we don't have to ask)
>
> Jehol
>        [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

_______________________________________________
R-sig-mixed-models@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
	[[alternative HTML version deleted]]

