[R-sig-ME] Random effects confidence intervals too large: infinite variance or blocking factor estimate is really 0?

Fernandez de la pradilla, J. I. n.fernandez at abdn.ac.uk
Fri Sep 12 13:22:53 CEST 2008


Dear mixed-up R world,

I am fitting a lmm model in "nlme" package in R version 2.6.0 to my population growth data with the following structure:

Response:       population growth between two sequential time steps (a continuous variable)
Fixed effects (full model version):      Treat + Season + Year + Clust + Treat:Season + Treat:Year +Treat:Clust + Season:Year + Treat:Season:Year (all categorical variables)
Random effects:         1|Block or 1|Block/Plot (all categorical)

The experimental design consists of population growth rates at the plot level taken at different 3 different seasons nested within 2 different years. The spatial experimental design consists of 3 areas (Clusters) each one of them with 2 replicate Blocks each one with 4 plots (1 plot for every Treatment). So we have: 3 areas * 2 Blocks * 4 Treatments (Plots) * 6Time steps (2 years*3 seasons).

I have explored different structures of the random effects to take into account the spatial and temporal design and the best fits seems to be the ones stated above, which are a quite a good description of the experiment. I have considered Block as the basic level of replication given that the experimental set up is randomised for Blocks by design. I plan to incorporate some autoregressive component (AR1) later on and see how it performs.

The following analysis is just a subsample of the original one, with just 2 areas(Clusters) * 2 Treatments(Plots). The subsampled database looks like this:

Time TimeN Year Season Clust Block Treat Plot      Dens      Growth
1    t1     1 2006     Sp    AB     A     I   AI 0.2400000 -0.12000000
2    t1     1 2006     Sp    AB     A    IV  AIV 0.2916667 -0.01166667
3    t1     1 2006     Sp    AB     B     I   BI 0.0000000  0.12000000
4    t1     1 2006     Sp    AB     B    IV  BIV 0.2800000  0.00000000
5    t1     1 2006     Sp    CD     C     I   CI 0.0000000  0.04000000
6    t1     1 2006     Sp    CD     C    IV  CIV 0.0000000  0.04000000
7    t1     1 2006     Sp    CD     D     I   DI 0.0400000 -0.04000000
8    t1     1 2006     Sp    CD     D    IV  DIV 0.1250000  0.03500000
9    t2     2 2006     Su    AB     A     I   AI 0.1200000  0.28000000
10   t2     2 2006     Su    AB     A    IV  AIV 0.2800000  0.42833333
11   t2     2 2006     Su    AB     B     I   BI 0.1200000  0.60000000
12   t2     2 2006     Su    AB     B    IV  BIV 0.2800000  0.36000000
13   t2     2 2006     Su    CD     C     I   CI 0.0400000 -0.04000000
14   t2     2 2006     Su    CD     C    IV  CIV 0.0400000  0.08000000
15   t2     2 2006     Su    CD     D     I   DI 0.0000000  0.00000000
16   t2     2 2006     Su    CD     D    IV  DIV 0.1600000  0.12000000
17   t3     3 2006     Fa    AB     A     I   AI 0.4000000 -0.08000000
18   t3     3 2006     Fa    AB     A    IV  AIV 0.7083333  0.05166667
19   t3     3 2006     Fa    AB     B     I   BI 0.7200000 -0.36000000
20   t3     3 2006     Fa    AB     B    IV  BIV 0.6400000 -0.04000000
21   t3     3 2006     Fa    CD     C     I   CI 0.0000000  0.00000000
22   t3     3 2006     Fa    CD     C    IV  CIV 0.1200000  0.12000000
23   t3     3 2006     Fa    CD     D     I   DI 0.0000000  0.08000000
24   t3     3 2006     Fa    CD     D    IV  DIV 0.2800000  0.28000000
25   t4     4 2007     Sp    AB     A     I   AI 0.3200000  0.04000000
26   t4     4 2007     Sp    AB     A    IV  AIV 0.7600000 -0.23619048
27   t4     4 2007     Sp    AB     B     I   BI 0.3600000 -0.20000000
28   t4     4 2007     Sp    AB     B    IV  BIV 0.6000000  0.04000000
29   t4     4 2007     Sp    CD     C     I   CI 0.0000000  0.00000000
30   t4     4 2007     Sp    CD     C    IV  CIV 0.2400000  0.08000000
31   t4     4 2007     Sp    CD     D     I   DI 0.0800000  0.08000000
32   t4     4 2007     Sp    CD     D    IV  DIV 0.5600000  0.20000000
33   t5     5 2007     Su    AB     A     I   AI 0.3600000  0.28000000
34   t5     5 2007     Su    AB     A    IV  AIV 0.5238095  0.39285714
35   t5     5 2007     Su    AB     B     I   BI 0.1600000  0.60000000
36   t5     5 2007     Su    AB     B    IV  BIV 0.6400000  0.20000000
37   t5     5 2007     Su    CD     C     I   CI 0.0000000  0.28000000
38   t5     5 2007     Su    CD     C    IV  CIV 0.3200000  0.52000000
39   t5     5 2007     Su    CD     D     I   DI 0.1600000  0.28000000
40   t5     5 2007     Su    CD     D    IV  DIV 0.7600000  0.16000000
41   t6     6 2007     Fa    AB     A     I   AI 0.6400000 -0.24000000
42   t6     6 2007     Fa    AB     A    IV  AIV 0.9166667 -0.31666667
43   t6     6 2007     Fa    AB     B     I   BI 0.7600000 -0.40000000
44   t6     6 2007     Fa    AB     B    IV  BIV 0.8400000 -0.36000000
45   t6     6 2007     Fa    CD     C     I   CI 0.2800000 -0.24000000
46   t6     6 2007     Fa    CD     C    IV  CIV 0.8400000 -0.56000000
47   t6     6 2007     Fa    CD     D     I   DI 0.4400000 -0.24000000
48   t6     6 2007     Fa    CD     D    IV  DIV 0.9200000 -0.24000000


And the structure:

'data.frame':   48 obs. of  10 variables:
 $ Time  : Factor w/ 7 levels "t1","t2","t3",..: 1 1 1 1 1 1 1 1 2 2 ...
 $ TimeN : int  1 1 1 1 1 1 1 1 2 2 ...
 $ Year  : Factor w/ 3 levels "2006","2007",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Season: Factor w/ 3 levels "Fa","Sp","Su": 2 2 2 2 2 2 2 2 3 3 ...
 $ Clust : Factor w/ 2 levels "AB","CD": 1 1 1 1 2 2 2 2 1 1 ...
 $ Block : Factor w/ 4 levels "A","B","C","D": 1 1 2 2 3 3 4 4 1 1 ...
 $ Treat : Factor w/ 2 levels "I","IV": 1 2 1 2 1 2 1 2 1 2 ...
 $ Plot  : Factor w/ 8 levels "AI","AIV","BI",..: 1 2 3 4 5 6 7 8 1 2 ...
 $ Dens  : num  0.240 0.292 0.000 0.280 0.000 ...
 $ Growth: num  -0.1200 -0.0117  0.1200  0.0000  0.0400 ...


So first we go for the most parameterized version with random effects specified as (1|Block/Plot), therefore allowing for 6 temporal replicates of every "Plot" level.

The full model output fitting a looks fine:
> Gm4c<-lme(fixed=Growth~Season+Year+Season:Year, random=~1|Block/Plot,method="ML")
> summary(Gm4c)

Linear mixed-effects model fit by maximum likelihood
 Data: NULL
      AIC       BIC logLik
  -29.808 -12.96719 23.904

Random effects:
 Formula: ~1 | Block
         (Intercept)
StdDev: 4.096235e-07

 Formula: ~1 | Plot %in% Block
        (Intercept)  Residual
StdDev: 8.33304e-08 0.1470565

 (Fixed effects omitted)

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-2.4919563 -0.5316126  0.1190019  0.5751758  2.5259568

Number of Observations: 48
Number of Groups:
          Block Plot %in% Block
              4               8

But I don't get any confidence intervals:

> intervals(Gm4c)
Error in intervals.lme(Gm4c) :
  Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance


So I reduce the random effects specification to allow just estimation of effects given by the grouping "Block" factor, therefore allowing temporal replication and increasing df. However, I still obtain ridiculously large confidence intervals for this grouping factor:

> Gm3c<-lme(fixed=Growth~Season+Year+Season:Year, random=~1|Block,method="ML")
> summary(Gm3c)

Linear mixed-effects model fit by maximum likelihood
 Data: NULL
      AIC       BIC logLik
  -31.808 -16.83839 23.904

Random effects:
 Formula: ~1 | Block
         (Intercept)  Residual
StdDev: 3.352947e-07 0.1470565

(Fixed effects omitted)

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-2.4919563 -0.5316126  0.1190019  0.5751758  2.5259568

Number of Observations: 48
Number of Groups: 4

Intervals:

> intervals(Gm3c)

Approximate 95% confidence intervals

 Fixed effects:
                        lower         est.      upper
(Intercept)       -0.09870606  0.006458333  0.1116227
SeasonSp          -0.14726657  0.001458333  0.1501832
SeasonSu           0.07335843  0.222083333  0.3708082
Year2007          -0.47976657 -0.331041667 -0.1823168
SeasonSp:Year2007  0.11327241  0.323601191  0.5339300
SeasonSu:Year2007  0.23127836  0.441607143  0.6519359
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: Block
                lower         est. upper
sd((Intercept))     0 3.352947e-07   Inf

 Within-group standard error:
    lower      est.     upper
0.1203951 0.1470565 0.1796221

So my question is two-fold:

1) Are the large confidence intervals obtained in the model Gm3c realistic (i.e. the variance due to Block effects is simply = 0 and Block is not explaining any variance at all) or just a computing problem (i.e. R is unable to calculate the confidence intervals due to limited replication)?

1)  When I allow temporal replication (6 time steps) for estimating variance due to "Plot" (model Gm4c), I would expect to get at least some confidence intervals (even as large as in Gm3c ), but I don't get anything. Does anybody has any idea why is that?

Thanks in advance. I look forward to your suggestions and welcome your views,

Nacho


Nacho Villar
PhD Student
ACES - Aberdeen Centre for Environmental Sustainability
University of Aberdeen-Zoology Building
Tillydrone Avenue - AB24 2TZ Aberdeen
n.fernandez at abdn.ac.uk
http://www.aces.ac.uk/people/




The University of Aberdeen is a charity registered in Scotland, No SC013683.




More information about the R-sig-mixed-models mailing list