Sorry,

Actually I gave my data in an image file (.RData) - I've just checked my send 
emails.
Am I to give data in another format, such as text ? Here are they in text 
(.txt).

The output are :

> 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)
> summary(lmer(Rendement ~ (1 |Pollinisateur) + (1 | Lignee) + (1 | Pollinisateur:Lignee), data=mca2))

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


Thanks,

Jacques VESLOT


Prof Brian Ripley a écrit :
Nothing was enclosed, nor was the output from summary.aov, so we are left guessing.

On Thu, 27 Oct 2005, Jacques VESLOT wrote:

Dear R-users,

My question is how to get right F tests for random effects in random effect models (I hope this question has not been answered too many times yet - I didn't find an answer in rhelp archives).

My data are in mca2 (enc.) :

names(mca2)
[1] "Lignee"        "Pollinisateur" "Rendement"

dim(mca2)
[1] 100   3

replications(Rendement ~ Lignee * Pollinisateur, data = mca2)
             Lignee        Pollinisateur Lignee:Pollinisateur
                 20                   10                    2

Of course, summary(aov(Rendement ~ Pollinisateur * Lignee, data = mca2)) gives wrong tests of random effects. But, summary(aov1 <- aov(Rendement ~ Error(Pollinisateur * Lignee), data = mca2)) gives no test at all, and I have to do it like this :

tab1 <- matrix(unlist(summary(aov1)), nc=5, byrow=T)[,1:3]

Femp <- c(tab1[1:3, 3]/tab1[c(3,3,4), 3])

names(Femp) <- c("Pollinisateur", "Lignee", "Interaction")

1 - pf(Femp, tab1[1:3,1], tab1[c(3,3,4),1])

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 (but don't know how to find them) with :

library(lme4)
summary(lmer(Rendement ~ (1 |Pollinisateur) + (1 | Lignee) + (1 | Pollinisateur:Lignee), data=mca2))

but I can't get F tests.

Thanks in advance.

Best regards,

Jacques VESLOT





Lignee  Pollinisateur   Rendement
L1      P1      13.4
L1      P1      13.3
L2      P1      12.4
L2      P1      12.6
L3      P1      12.7
L3      P1      13
L4      P1      12.6
L4      P1      12.6
L5      P1      11.9
L5      P1      11.6
L1      P2      12.6
L1      P2      12.7
L2      P2      12.1
L2      P2      11.3
L3      P2      12.4
L3      P2      11.9
L4      P2      12.2
L4      P2      12.1
L5      P2      11.3
L5      P2      11.5
L1      P3      13.3
L1      P3      13.2
L2      P3      12.9
L2      P3      12.3
L3      P3      12.1
L3      P3      12.9
L4      P3      12.8
L4      P3      13.4
L5      P3      11.7
L5      P3      11.8
L1      P4      13.5
L1      P4      14
L2      P4      13.5
L2      P4      12.7
L3      P4      13.3
L3      P4      13.5
L4      P4      13.4
L4      P4      13.4
L5      P4      13
L5      P4      13.1
L1      P5      13.7
L1      P5      13.8
L2      P5      12.5
L2      P5      13.1
L3      P5      12.8
L3      P5      12.5
L4      P5      13.7
L4      P5      13.8
L5      P5      12.2
L5      P5      12.1
L1      P6      12.8
L1      P6      13.1
L2      P6      12
L2      P6      11.8
L3      P6      12.4
L3      P6      12.2
L4      P6      13.3
L4      P6      13.5
L5      P6      12.4
L5      P6      11.5
L1      P7      13.2
L1      P7      13.2
L2      P7      12.4
L2      P7      12.2
L3      P7      12.1
L3      P7      12.1
L4      P7      12.9
L4      P7      13.3
L5      P7      12.6
L5      P7      12.2
L1      P8      13.4
L1      P8      12.8
L2      P8      12.2
L2      P8      12.4
L3      P8      12
L3      P8      12.2
L4      P8      12.7
L4      P8      13.4
L5      P8      11.8
L5      P8      11.8
L1      P9      12.6
L1      P9      13
L2      P9      11.9
L2      P9      12.3
L3      P9      12.6
L3      P9      12.5
L4      P9      13.5
L4      P9      12.9
L5      P9      11.5
L5      P9      12.1
L1      P10     12.7
L1      P10     12.5
L2      P10     11.8
L2      P10     12.1
L3      P10     12.3
L3      P10     12.1
L4      P10     12.5
L4      P10     12.7
L5      P10     12.3
L5      P10     12.1
______________________________________________
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

Reply via email to