Re: [R] Random effect models

From: Christoph Buser <buser_at_stat.math.ethz.ch>
Date: Sat 29 Oct 2005 - 00:01:05 EST

Dear Jacques

You may be interested in the answer of Doug Bates to another post. The question is different, but in the answer there is a part about testing if a variance component may be zero:

http://finzi.psych.upenn.edu/R/Rhelp02a/archive/51080.html

Hope this helps.

Best regards,

Christoph Buser



Christoph Buser <buser@stat.math.ethz.ch> Seminar fuer Statistik, LEO C13
ETH (Federal Inst. Technology)	8092 Zurich	 SWITZERLAND
phone: x-41-44-632-4673		fax: 632-1228

http://stat.ethz.ch/~buser/

Jacques VESLOT writes:
> Dear R-users,
>
> Sorry for reposting. I put it in another way :
>
> I want to test random effects in this random effect model :
> Rendement ~ Pollinisateur (random) + Lignee (random) + Pollinisateur:Lignee (random)
>
> Of course :
> summary(aov(Rendement ~ Pollinisateur * Lignee, data = mca2))
> gives wrong tests for random effects.
>
> But :
> summary(aov1 <- aov(Rendement ~ Error(Pollinisateur * Lignee), data = mca2))
> gives no test at all, and I have to do it with mean squares lying in summary(aov1).
>
> With "lme4" package (I did'nt succeed in writing a working formula with lme from "nlme" package),
> I can "see" standard deviations of random effects (I don't know how to find them) but I can't find F
> tests for random effects.
>
> I only want to know if there is an easy way (a function ?) to do F tests for random effects in
> random effect models.
>
> Thanks in advance.
>
> Best regards,
>
> Jacques VESLOT
>
>
> Data and output are as follows :
>
> > head(mca2)
> Lignee Pollinisateur Rendement
> 1 L1 P1 13.4
> 2 L1 P1 13.3
> 3 L2 P1 12.4
> 4 L2 P1 12.6
> 5 L3 P1 12.7
> 6 L3 P1 13.0
>
>
> > names(mca2)
> [1] "Lignee" "Pollinisateur" "Rendement"
>
> > dim(mca2)
> [1] 100 3
>
> > replications(Rendement ~ Lignee * Pollinisateur, data = mca2)
> Lignee Pollinisateur Lignee:Pollinisateur
> 20 10 2
>
> > summary(aov1 <- aov(Rendement ~ Error(Pollinisateur * Lignee), data = mca2)
>
> Error: Pollinisateur
> Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 9 11.9729 1.3303
>
> Error: Lignee
> Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 4 18.0294 4.5074
>
> Error: Pollinisateur:Lignee
> Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 36 5.1726 0.1437
>
> Error: Within
> Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 50 3.7950 0.0759
>
>
> # F tests :
>
> > Femp <- c(tab1[1:3, 3]/tab1[c(3,3,4), 3])
> > names(Femp) <- c("Pollinisateur", "Lignee", "Interaction")
> > Femp
> Pollinisateur Lignee Interaction
> 9.258709 31.370027 1.893061
>
> > 1 - pf(Femp, tab1[1:3,1], tab1[c(3,3,4),1])
> Pollinisateur Lignee Interaction
> 4.230265e-07 2.773448e-11 1.841028e-02
>
> # Standard deviation :
>
> > variances <- c(c(tab1[1:3, 3] - tab1[c(3,3,4), 3]) / c(2*5, 2*10, 2), tab1[4,3])
> > names(variances) <- c(names(Femp), "Residuelle")
> > variances
> Pollinisateur Lignee Interaction Residuelle
> 0.11866389 0.21818333 0.03389167 0.07590000
>
> # Using lmer :
>
> > library(lme4)
> > lme1 <- lmer(Rendement ~ (1|Pollinisateur) + (1|Lignee) + (1|Pollinisateur:Lignee), data = mca2))
> > summary(lme1)
> Linear mixed-effects model fit by REML
> Formula: Rendement ~ (1 | Pollinisateur) + (1 | Lignee) + (1 | Pollinisateur:Lignee)
> Data: mca2
> AIC BIC logLik MLdeviance REMLdeviance
> 105.3845 118.4104 -47.69227 94.35162 95.38453
> Random effects:
> Groups Name Variance Std.Dev.
> Pollinisateur:Lignee (Intercept) 0.033892 0.18410
> Pollinisateur (Intercept) 0.118664 0.34448
> Lignee (Intercept) 0.218183 0.46710
> Residual 0.075900 0.27550
> # of obs: 100, groups: Pollinisateur:Lignee, 50; Pollinisateur, 10; Lignee, 5
>
> Fixed effects:
> Estimate Std. Error DF t value Pr(>|t|)
> (Intercept) 12.60100 0.23862 99 52.808 < 2.2e-16 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> > anova(lme1)
> Analysis of Variance Table
> Erreur dans ok[, -nc] : nombre de dimensions incorrect
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sat Oct 29 00:13:37 2005

This archive was generated by hypermail 2.1.8 : Sat 29 Oct 2005 - 03:44:34 EST