[R] Random effect models

[Jacques VESLOT]

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