[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