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