Skip to content

Why am I getting a Variance of 0 for my random effect

11 messages · Daniel Ezra Johnson, Gustavo Betini, Kevin E. Thorpe +3 more

Original
Kevin E. Thorpe
# 
Hello.

I'm getting a variance of 0 on a random effect and I don't know why.
I suspect I've not set the model up correctly.  My transcript is below
with my own comments sprinkled in for time to time.

A little bit about the data (which I will provide off-list if 
requested).  We have nurses managing an aspect of patient care
according to different algorithms.  Interest focuses on of the
algorithms result in different outcomes.  I have restricted this
to only nurses who did each algorithm twice (in case my problem
was being caused by some nurses doing only one algorithm, possibly
only one time).

I figured that since I have multiple observations per nurse, I
should treat nurse as a random effect, but maybe I confused myself
again.


R version 2.11.1 Patched (2010-07-21 r52598)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

   Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

 > library(lattice)
 > library(lme4)

 > str(data1)
'data.frame':   72 obs. of  3 variables:
  $ RN        : int  1 1 2 3 7 7 9 9 15 15 ...
  $ Assignment: Factor w/ 2 levels "E","N": 1 1 1 1 1 1 1 1 1 1 ...
  $ AUChr     : num  12.26 7.23 9.26 4.04 10.31 ...
 > tmp1 <- 
with(data1,aggregate(AUChr,list(RN=RN,Assigment=Assignment),mean))
 > names(tmp1)[3] <- "Mean"
 >
 > tmp2 <- 
with(data1,aggregate(AUChr,list(RN=RN,Assignment=Assignment),var))
 > names(tmp2)[3] <- "Variance"
 >
 > meanvar <- merge(tmp1,tmp2)

The point of this is to show that the means are not all the same,
nor are the variances.

 > meanvar
    RN Assignment   Mean  Variance
1   1          E  9.745  12.65045
2   1          N  7.185   1.36125
3  15          E 10.605  15.07005
4  15          N 10.385   4.41045
5  16          E  8.175   0.00845
6  16          N  8.420   1.03680
7   2          E  7.300   7.68320
8   2          N  6.950   1.00820
9  21          E  9.670   9.41780
10 21          N 10.535   2.44205
11 22          E  7.720   2.04020
12 22          N  7.930   1.21680
13 24          E  9.555  10.35125
14 24          N  9.330   0.38720
15 25          E  8.240   0.92480
16 25          N  9.485   0.00125
17 27          E  8.635   0.08405
18 27          N  7.745   3.72645
19 28          E  9.635   8.61125
20 28          N  8.315  10.35125
21  3          E  6.005   7.72245
22  3          N 11.435  55.44045
23 31          E  9.590   9.94580
24 31          N 10.570  16.70420
25 35          E  9.055   0.32805
26 35          N  9.925  14.41845
27 36          E  9.040   2.08080
28 36          N  7.395   1.14005
29  5          E  8.430   3.38000
30  5          N 17.385 139.94645
31  6          E  6.930   0.24500
32  6          N  8.330   1.72980
33  7          E 10.650   0.23120
34  7          N  7.375   0.09245
35  9          E  8.885   7.56605
36  9          N  8.405   0.73205

Model with "Assignment" (algorithm).

 > lmer(AUChr~Assignment+(1|RN),data=data1,REML=FALSE)
Linear mixed model fit by maximum likelihood
Formula: AUChr ~ Assignment + (1 | RN)
    Data: data1
    AIC   BIC logLik deviance REMLdev
  365.7 374.8 -178.8    357.7   356.9
Random effects:
  Groups   Name        Variance Std.Dev.
  RN       (Intercept) 0.0000   0.0000
  Residual             8.4152   2.9009
Number of obs: 72, groups: RN, 18

Fixed effects:
             Estimate Std. Error t value
(Intercept)   8.7703     0.4835   18.14
AssignmentN   0.5131     0.6837    0.75

Correlation of Fixed Effects:
             (Intr)
AssignmentN -0.707


Model without the algorithm variable.

 > lmer(AUChr~(1|RN),data=data1,REML=FALSE)
Linear mixed model fit by maximum likelihood
Formula: AUChr ~ (1 | RN)
    Data: data1
    AIC   BIC logLik deviance REMLdev
  364.3 371.1 -179.1    358.3   358.5
Random effects:
  Groups   Name        Variance Std.Dev.
  RN       (Intercept) 0.000    0.0000
  Residual             8.481    2.9122
Number of obs: 72, groups: RN, 18

Fixed effects:
             Estimate Std. Error t value
(Intercept)   9.0268     0.3432    26.3
 >
 > sessionInfo()
R version 2.11.1 Patched (2010-07-21 r52598)
Platform: i686-pc-linux-gnu (32-bit)

locale:
  [1] LC_CTYPE=en_US       LC_NUMERIC=C         LC_TIME=en_US
  [4] LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=en_US
  [7] LC_PAPER=en_US       LC_NAME=C            LC_ADDRESS=C
[10] LC_TELEPHONE=C       LC_MEASUREMENT=en_US LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] lme4_0.999375-34   Matrix_0.999375-42 lattice_0.18-8

loaded via a namespace (and not attached):
[1] grid_2.11.1   nlme_3.1-96   stats4_2.11.1
 >
 > proc.time()
    user  system elapsed
   3.488   0.056   3.536
Daniel Ezra Johnson
# 
Try data1$RN <- as.factor(data1$RN).

On Wed, Aug 11, 2010 at 2:13 PM, Kevin E. Thorpe
<kevin.thorpe at utoronto.ca> wrote:
Kevin E. Thorpe
#
On 08/11/2010 02:16 PM, Daniel Ezra Johnson wrote:
Thanks, but that has no effect.  That is I get the same results.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
On Wed, Aug 11, 2010 at 1:13 PM, Kevin E. Thorpe

        
<kevin.thorpe at utoronto.ca> wrote:

            

      
      
      
    

        
It's not a bug - it's a feature.  ML estimates or REML estimates of
variance components can be zero.  This simply indicates that the
variability in the response associated with the factor, RN in your
case, is not sufficient to warrant the additional complexity in the
model.

            

      
      
      
    

        
You are quite correct that it is appropriate to allow for the
possibility of RN having an effect on the response and that it should
be incorporated as a random effect if it were in the model but the
results indicate that RN does not have sufficient effect on the
response.

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
On Wed, Aug 11, 2010 at 1:16 PM, Daniel Ezra Johnson

        
<danielezrajohnson at gmail.com> wrote:

            

      
      
      
    

        
A reasonable suggestion but it wouldn't make a difference.  Note that
lmer is already picking up on there being 18 levels of the RN factor.

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
On Wed, Aug 11, 2010 at 1:40 PM, Gustavo Betini <betinig at uoguelph.ca> wrote:

            

      
      
      
    

        
Yes.  For example,

            

      
      
      
    

        
npt = 7 , n =  3
rhobeg =  0.2 , rhoend =  2e-07
   0.020:  11:      2368.50; 1.09296 -0.173139 0.0953204
  0.0020:  30:      2364.50; 1.48770 -0.374305 0.0138819
 0.00020:  42:      2364.50; 1.48462 -0.372458 0.00762182
 2.0e-05:  58:      2364.50; 1.48417 -0.372319 0.00114305
 2.0e-06:  74:      2364.50; 1.48420 -0.372480  0.00000
 2.0e-07:  80:      2364.50; 1.48420 -0.372481  0.00000
At return
 85:     2364.5016:  1.48420 -0.372481 2.77475e-07

            

      
      
      
    

        
Linear mixed model fit by REML ['merMod']
Formula: cog ~ tos + trt:tos + (tos | id)
   Data: Early
REML criterion at convergence: 2364.502

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 id       (Intercept) 166.40   12.900
          tos          10.48    3.237   -1.000
 Residual              75.54    8.691
Number of obs: 309, groups: id, 103

Fixed effects:
            Estimate Std. Error t value
(Intercept)  120.783      1.824   66.22
tos          -22.470      1.494  -15.04
tos:trtY       7.646      1.447    5.28


The resulting model no longer fulfills the technical definition of a
linear mixed-effects model.

  
  


  


      

    

    

      
      

        


  

        

      Kevin E. Thorpe
    

    

      
      #
    

  


  

      
On 08/11/2010 02:25 PM, Douglas Bates wrote:

            

      
      
      
    

        
Thanks Doug.  Your response is helpful, as always.  If RN does not
contribute a random effect, would it be appropriate to revert to a
standard regression model, or is it best to leave the unimportant
random effect?  In the case of ordinary regression, dropping
variables based on their p-values compromises inference.  Does the
same apply with dropping a random effect with no variance?

Kevin

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
On Wed, Aug 11, 2010 at 2:18 PM, Kevin E. Thorpe

        
<kevin.thorpe at utoronto.ca> wrote:

            

      
      
      
    

        
When the random effects variance is zero the model reverts to the
linear regression model in the sense that the two models give the same
predicted values, the same log-likelihood, coefficient estimates for
the fixed effects and standard errors.  That is, there is no need to
continue to represent the model as a linear mixed model if the only
variance component parameter's estimated value is zero.

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Andrew Robinson
    

    

      
      #
    

  


  

      
On Wed, Aug 11, 2010 at 03:05:06PM -0500, Douglas Bates wrote:

            

      
      
      
    

        
It's worth noting that this behaviour will not always be what you want
to happen.  There are times when you may want the inference and
estimation that arise from a model to reflect the inclusion of random
effects, even if the mode of the density of those effects is zero.

This is true, for example, in design-based inference.  For example,
you may feel strongly that the experimental design (or sample design)
has elements of clustering, and that to fit a model that ignores the
clustering will result in negatively biased estimates of the standard
errors of the fixed effects.

If your inference is model based, then this behaviour should be
perfectly fine.

Best wishes,

Andrew

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      John Maindonald
    

    

      
      #
    

  


  

      
Surely you do want to treat nurses as a random effect, by analysing 
summary data at the nurse level, if not in a multi-level model.

The zero variance may (if it really would prefer to be negative) be telling 
you that there is a systematic difference between the two times, for which 
your model needs to account.  Maybe there is a learning effect -- 2nd time 
is systematically different from the first.  Did your model account for such 
an effect?

Or (requires more thought to model), those who do badly the first time
may learn rather more from their experience than those who did
moderately well, doing better than average next time?  It appears that
the data have the information needed to get insight on these questions.

The most insightful approach might well be separate regressions
for 2-1 differences and 2+1 averages.  I'd do those analyses whatever 
else you do.  

John Maindonald             email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.
http://www.maths.anu.edu.au/~johnm

        
On 12/08/2010, at 4:24 AM, Kevin E. Thorpe wrote: