[R] F tests for random effect models

[Jacques VESLOT]

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
>>
>>
>>
>>
> 
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: mca2.txt
Url: https://stat.ethz.ch/pipermail/r-help/attachments/20051027/ce7d42e4/mca2.txt