[R] Random effect models

[Christoph Buser]

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 at 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 at 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