[R] Random effect models

From: Jacques VESLOT <jacques.veslot_at_cirad.fr>
Date: Fri 28 Oct 2005 - 21:22:21 EST


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
Received on Fri Oct 28 21:29:39 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:40:52 EST