An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20100525/71fe2510/attachment.pl>
Mixed-effects model with nested factors
2 messages · Manuel Ramon, 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.