[R] Random effect models

[Doran, Harold]

Dear Jacques

I think your question is a little confusing, but let me take a stab at
what I think you're getting at. It seems you are trying to find the
statistical significance of a variance component in your lme model, and
not the significance of a fixed effect. If this is what you are looking
for you will not find this as standard output in lme or lmer. Not in
summary() or anova()

I don't know what you mean by "I can see the sd of ranodm effects but
you do not know how to find them"

Other software programs for mixed linear models do indeed provide this
as ouput, but not the mixed model functions in R. This topic has been
discussed often on this list and the package developer, Doug Bates, has
noted numerous times, with good empirical examples, why such a test may
potentially be misleading. Use the archives to study this rationale and
see the numerous threads on the topic.

Note that the standard deviations of the random effects are *not* the
standard errors of those variance components.

I also noticed you are loading the lme4 library and using lmer. You
should load Matrix, which has the most recent lmer function in that
package.

If I am mistaken, then please clarify and maybe someone else can help.

-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Jacques VESLOT
Sent: Friday, October 28, 2005 7:22 AM
To: R-help at stat.math.ethz.ch
Subject: [R] Random effect models

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