Hi,

1/ I wonder why a anova.lme on a single lme object does not print the sum of squares 
(as expected from the help: "a data frame with the sums of squares, numerator degrees 
of freedom, denominator
degrees of freedom, F-values, and P-values").

Example:

>   fm2 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)
> anova(fm2)
            numDF denDF  F-value p-value
(Intercept)     1    80 4123.156  <.0001
age             1    80  114.838  <.0001
Sex             1    25    9.292  0.0054

How can one get the sum of squares easily?

2/  I have tried a home-made permutation test for comparison with results obtained 
above. I am a quite surprised that the intercept has a probability close to 1:

> lme.perm.test(fm2)
            p.value
(Intercept)   0.998
age           0.000
Sex           0.001

The home made function is grounded on 1000 permutations of the response variable 
(distance) and the comparison of the F-values obtained with the F-values observed with 
the original (non permuted)data.

Thanks in advance for any hint,

Patrick Giraudoux



PS: below the current script of the "draft" permutation test:


lme.perm.test<-function(lme.model,pn=1000){

    # Giraudoux 17.4.2004 permutation test for linear mixed effect models
    # lme.model = a lme object
    # pn = permutation number (default = 1000)
    # value: a data.frame of the p.values of each independent variable F
    # the permutation number is kept in the ouput attributes and thus
    # can be called from there
   
    an<-anova(lme.model)
    Fobs<-an[[3]]
    n<-rep(0,length(an[[3]]))
    
    for (i in 1:pn){
        lm1<-update(lme.model, sample(.)~.)
        an<-anova(lm1)
        Frdm<-an[[3]]
            for (j in 1:(length(an[[3]]))){
                if (Frdm[j]>=Fobs[j]) n[j]<-n[j]+1
            }
        }
    n<-n/pn
    names(n)<-row.names(an)
    n<-as.data.frame(n)
    names(n)<-"p.value"
    attributes(n)<-c(attributes(n),permutations=pn)
    return(n)
}



        [[alternative HTML version deleted]]

______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Reply via email to