[R-sig-ME] Mixed-effects model with nested factors

[ONKELINX, Thierry]

Dear Manuel,

The fixed effect of sp IS estimated. You don't get a p-value because it has no df. You have no df because you used the factor in both the fixed effects as in the random effects. So remove sp from the random effect.

Furthermore, random effect with only a few levels is not a very good idea. You get unreliable estimates of the variance. Have a look at the plot below. It plots the confidence interval for the ratio of the sample variance and the population variance. Note that the interval is very wide with a small number of levels. 6 levels is a minimum minimorum. Ideally you should aim for 30 levels or more.

n <- 2:100
plot(n, qchisq(0.975, n - 1) / (n - 1), type = "l", ylim = c(0, 5.5))
lines(n, qchisq(0.025, n - 1) / (n - 1))
abline(h = 1, lty = 2, col = "red")
abline(v = c(2, 6, 10, 30), lty = 3, col = "blue")

HTH,

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 at 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-sig-mixed-models-bounces at r-project.org 
> [mailto:r-sig-mixed-models-bounces at r-project.org] Namens Manuel Ramon
> Verzonden: dinsdag 25 mei 2010 17:31
> Aan: r-sig-mixed-models at r-project.org
> Onderwerp: [R-sig-ME] Mixed-effects model with nested factors
> 
> Hello,
> I having problems on the definition of a mixed model with 
> nested factors. My data base is as follows:
>      y:      performance measure
>      sp:    species factor with two levels
>      ani:   individual with 6 levels but three of one belong 
> to the fist
> species
>              and the other three to the other. That is, sp 
> and ani are NESTED
>      v1:    treatment with two levels
> 
> The aim of my work is to study if there are significant 
> differences between both species (sp), and both treatments 
> (v1). I want to consider the ani effect as random effect.  My 
> attemp is the next:
> 
>      lme( y ~ sp + v1, data=mydata, random=~1|sp/ani )
> 
> My first question is if the fm1 model is well desing. When I 
> run this analysis, the factor sp is not estimated, but I do 
> not know why?
> 
> A simulated example could be:
> 
> ##    Simul data
>      set.seed(1234)
> 
>      N <- 24
>      error <- rnorm(N)   # error term
>      ani <- rep(1:(N/4),each=4)   # individual effect (random)
>      b1 <- c(0,2,4,1,2,5)
>      sp <- rep(1:2, each=N/2)   # species (fixed)
>      b2 <- c(1,5)
>      v1 <- rep(1:2,each=2,length=N)   # treatment (fixed)
>      b3 <- c(2,7)
> 
>      y <- b1[ani] + b2[sp] + b3[v1] + error
> 
>      dat <- data.frame(y,ani,sp,v1)
>           dat$ani <- as.factor(dat$ani)
>           dat$sp <- as.factor(dat$sp)
>           dat$v1 <- as.factor(dat$v1)
> 
>      library(nlme)
>      fm1 <- lme(y ~ sp + v1, data=dat, random=~1|sp/ani)
>      summary(fm1)
> 
> 
> ## Output
> Linear mixed-effects model fit by REML
> Data: dat
>       AIC      BIC   logLik
>   94.0726 100.3397 -41.0363
> 
> Random effects:
>  Formula: ~1 | sp
>         (Intercept)
> StdDev:   0.8539529
> 
>  Formula: ~1 | ani %in% sp
>         (Intercept) Residual
> StdDev:    2.226151 1.109317
> 
> Fixed effects: y ~ sp + v1
>                 Value         Std.Error     DF   t-value       p-value
> (Intercept)  5.124082   1.5921603   17    3.218320   0.005
> sp2           4.220951   2.2287665     0    1.893851     NaN
> v12            5.160433   0.4528769   17  11.394781   0.000
>  Correlation:
>     (Intr) sp2
> sp2 -0.700
> v12 -0.142  0.000
> 
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max
> -1.39622076 -0.59150658  0.06285846  0.41361144  1.75700646
> 
> Number of Observations: 24
> Number of Groups:
>          sp ani %in% sp
>           2           6
> Mensajes de aviso perdidos
> In pt(q, df, lower.tail, log.p) : Se han producido NaNs
> 
> Why the sp effect has not DF? It could be due to its presence 
> in the random term.
> Is my model rigth? I think that I'm missing something.
> 
> Thanks in advance
> 
> 
> --
> --
> Manuel Ramón
> m.ramon.fernandez at gmail.com
> 
> 	[[alternative HTML version deleted]]
> 
> 

Druk dit bericht a.u.b. niet onnodig af.
Please do not print this message unnecessarily.

Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message 
and any annex are purely those of the writer and may not be regarded as stating 
an official position of INBO, as long as the message is not confirmed by a duly 
signed document.