Re: [R] Random effect models

From: Doran, Harold <HDoran_at_air.org>
Date: Fri 28 Oct 2005 - 22:01:04 EST


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@stat.math.ethz.ch
[mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Jacques VESLOT Sent: Friday, October 28, 2005 7:22 AM
To: R-help@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@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 Fri Oct 28 22:06:08 2005

This archive was generated by hypermail 2.1.8 : Sat 29 Oct 2005 - 00:13:02 EST