[R] NaNs in Nested Mixed Model

[ONKELINX, Thierry]

Dear Johan,

A few remarks.

- R-sig-mixed models is a better list for asking questions about mixed model.
- I presume that Nymphs is the number of insects? In that case you need a generalised linear (mixed) model with poisson family
- What are you interessed in? The variability among genotypes or the effect of each genotype. 
	You can achieve the first with a glmm like glmer(Nymphs ~ Species + (1|Genotype), family = "poisson"). Genotype will be implicitly nested in Species.  Note that since you have only 4 genotypes, you will not get very reliable estimates of the genotype variance.
	For the latter you cannot use a mixed model so you need a simple glm(Nymphs ~ Species/Genotype, family = "poisson"). Note that several coefficients will be NaN, because you cannot estimate them.

Best regards,

Thierry

----------------------------------------------------------------------------
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium

Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium

tel. + 32 54/436 185
Thierry.Onkelinx op inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
  

> -----Oorspronkelijk bericht-----
> Van: r-help-bounces op r-project.org 
> [mailto:r-help-bounces op r-project.org] Namens Johan Stenberg
> Verzonden: dinsdag 8 maart 2011 16:52
> Aan: r-help op r-project.org
> Onderwerp: [R] NaNs in Nested Mixed Model
> 
> Dear R users,
> 
> I have a problem with something called "NaNs" in a nested mixed model.
> 
> The background is that I have studied the number of insect 
> nymphs emerging from replicated Willow genotypes in the 
> field. I have 15 replicates each of 4 Willow genotypes 
> belonging two 2 Willow species.
> Now I want to elucidate the effect of Willow genotype on the 
> number of emerging nymphs. Previously I performed a simple 
> one-way anova with "genotype" as explanatory factor and 
> "number of nymphs emerging" as dependent variable, but the 
> editor of the journal I've submitted this piece to wants me 
> to nest Willow genotype within Willow species before he 
> accepts the paper for publication [Species*Genotype(Species)].
> 
> The fact that I didn't include "Willow species" as a factor 
> in my initial analysis reflects that I am not very interested 
> in the species factor per se - I am just interested in if 
> genetic variation in the host plant is important, but 
> "species" is of course a factor that structures genetic diversity.
> 
> I thought the below model would be appropriate:
> 
> > model<-lme(Nymphs~Species*Genotype,random=~1|Species/Genotype)
> 
> ...but I then get the error message "Error in MEEM(object, conLin,
> control$niterEM) : Singularity in backsolve at level 0, block 1"
> 
> I then tried to remove "Genotype" from the fixed factors, but 
> then I get the error message "NaNs produced".
> 
> > model<-lme(Nymphs~Species,random=~1|Species/Genotype)
> > summary(model)
> Linear mixed-effects model fit by REML
>  Data: NULL
>        AIC      BIC    logLik
>   259.5054 269.8077 -124.7527
> 
> Random effects:
>  Formula: ~1 | Species
>         (Intercept)
> StdDev:   0.9481812
> 
>  Formula: ~1 | Genotype %in% Species
>         (Intercept) Residual
> StdDev:   0.3486937 1.947526
> 
> Fixed effects: Nymphs ~ Species
>                  Value Std.Error DF   t-value p-value
> (Intercept)   2.666667  1.042243 56  2.558585  0.0132
> Speciesviminalis -2.033333  1.473954  0 -1.379510     NaN
>  Correlation:
>              (Intr)
> Speciesviminalis -0.707
> 
> Standardized Within-Group Residuals:
>        Min         Q1        Med         Q3        Max
> -1.4581821 -0.3892233 -0.2751795  0.3439871  3.1630658
> 
> Number of Observations: 60
> Number of Groups:
>           Species Genotype %in% Species
>             2             4
> Warning message:
> In pt(q, df, lower.tail, log.p) : NaNs produced
> ***********
> 
> Do you have any idea what these error messages mean in my 
> case and how I can get around them?
> 
> Thank you on beforehand! (data set attached).
> 
> Johan
>