[R] Wide confidence intervals or Error message in a mixed effects model (nlme)

[Bert Gunter]

Inline below. -- Bert

On Mon, Jul 25, 2011 at 7:42 AM, Dieter Menne
<dieter.menne at menne-biomed.de> wrote:
>
> Menelaos Stavrinides-2 wrote:
>>
>> I am analyzing a dataset on the effects of six pesticides on population
>> growth rate of a predatory mite. The response variable is the population
>> growth rate of the mite (ranges from negative to positive) and the
>> exploratory variable is a categorical variable (treatment). The
>> experiment was blocked in time (3 blocks / replicates per block) and it
>> is unbalanced - at least 1 replicate per block. I am analyzing the data
>> in nlme using model<-lme(growth.rate~treatment,random=~1|block).
>>
>>
>> In another study, I am investigating the interactions between pesticides
>> in a two-way design: (pesticideA x no pesticide A) crossed with
>> (pesticideB x no pesticide B). The blocking is as above, and the data
>> are unbalanced again. The model is defined as
>> model<-lme(growth.rate~pestA*pestB,random=~1|block). When I run
>> intervals (model), I usually get the following error message: "Error in
>> intervals.lme(model) : Cannot get confidence intervals on var-cov
>> components: Non-positive definite approximate variance-covariance".
>>
>>
>
> It depends on your data which you did not show. At least show str() or
> summary(), it could also be a forgotten factor(). With simulated data, it is
> quite easy to get reasonable-looking cases where var-cov is degenerate.

Yes, especially when you must estimate the the b(b+1)/2 = 6 covariance
parameters for b=3 blocks  of the unstructured full variance
covariance matrix. You might try using a structured matrix (e.g.
?pdDiag ) to reduce the number of random effects parameters you are
estimating. Remember, this is a nonlinear estimation problem, and you
need a "lot" of data to get good estimates in nonlinear estimation.

I would also address future questions of this sort to the
R-sig-mixed-models list, as this is not really an r-help kind of
question.

-- Bert

>
> Dieter
>
> library(nlme)
> set.seed(47)
> d = expand.grid(block=LETTERS[1:3],treatment=letters[1:6])
> d$growth.rate = as.integer(d$treatment)*0.2+rnorm(nrow(d))
> #d.lme = lme(growth.rate~treatment, random=~1|block,data=d)
>
> for (i in 1:100){
>  d1= d[sample(nrow(d) ,nrow(d)-3)  ,]
>  d.lme = lme(growth.rate~treatment, random=~1|block,data=d1)
>  iv = try(intervals(d.lme), TRUE)
>  if (inherits(iv, "try-error")){
>    print(iv)
>    print(table(d1[,c("block","treatment")]))
>  }
> }
>
>
>
>
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/Wide-confidence-intervals-or-Error-message-in-a-mixed-effects-model-nlme-tp3692637p3692986.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
"Men by nature long to get on to the ultimate truths, and will often
be impatient with elementary studies or fight shy of them. If it were
possible to reach the ultimate truths without the elementary studies
usually prefixed to them, these would not be preparatory studies but
superfluous diversions."

-- Maimonides (1135-1204)

Bert Gunter
Genentech Nonclinical Biostatistics