[R-sig-ME] Gmodels, estimable() question

eirintan at stud.ntnu.no eirintan at stud.ntnu.no
Fri May 13 11:06:57 CEST 2011


Hi,

I'm a master student in statistics at NTNU in Norway, writing my thesis on
Linear mixed effect models. I was recommended to try this mailing list  
by Douglas Bates, because I have some questions regarding estimable  
(gmodels) in R and  denominator degrees of freedom for conditional  
tests for fixed-effect terms, which I hope you can be able to help me  
with.

I'm testing different contrasts, using estimable (gmodels) in R. This method
uses the denominator degrees of freedom from Pinheiro & Bates (2000) page
91.

My data are repeated measurements (4 for each individual) on 32  
individuals, across factors "time" and "diet" which both contains two  
levels (respectively, Day "0" and "6" and diet "A" and "B"). I also  
have a factor sex, which contains the gender for the 32 individuals.  
My model is coded as follows:

model <- lme(resp ~ factor(sex)*factor(time)*factor(diet), random = ~1 | id,
data = this.data)

Hence,
N=32*4=128
Q=1 (only grouped by id)
m_0=1, m_1=32, m_2=128
p_0=1, m_1=1, p_2= 6
Which results in:
denDF_1 = 32-(1+1)=30
denDF_2 = 128-(32+6)=90

I want to find the mean contrasts for male and female, and so I use a weight
which is: (number of male participants)/(number of participants)=0.59375.
And then my contrasts of interest is:
a0 <- c(1,0.59375,0,0,0,0,0,0)
a6 <- c(1,0.59375,1,0,0.59375,0,0,0)
b0 <- c(1,0.59375,0,1,0,0.59375,0,0)
b6 <- c(1,0.59375,1,1,0.59375,0.59375,1,0.59375)
A <- a6-a0
B <- b6-b0
AB <- B-A
ContrastTable <- rbind(estimable(model, a0), estimable(model, a6),
estimable(model, b0), estimable(model, b6), estimable(model, A),
estimable(model, B), estimable(model, AB))

The results are:
                 Estimate Std. Error    t value      DF     Pr(>|t|)
A0              607.50000   55.37051 10.9715437 17.8125 2.360399e-09
A6              810.43750   55.37051 14.6366262 17.8125 2.274803e-11
B0              566.65625   55.37051 10.2338992 17.8125 6.939256e-09
B6              706.43750   55.37051 12.7583702 17.8125 2.144311e-10
A6-A0           202.93750   47.15961  4.3032052 53.4375 7.210324e-05
B6-B0           139.78125   47.15961  2.9640032 53.4375 4.528098e-03
(B6-B0)-(A6-A0) -63.15625   66.69377 -0.9469588 53.4375 3.479224e-01

With the warning: Degrees of freedom vary among parameters used to construct
linear contrast(s): 1. Using the smallest df among the set of parameters.

Ok, so we use denDF_1 in A0, A6, B0 and B6, since sex(estimated at level 1)
is present in the linear functions, and denDF_2 in (A6-A0), (B6-B0) and
(B6-B0)-(A6-A0), since only factors estimated at level 2 is present in these
linear functions.

So here are my questions:

1) Why do the interaction terms between sex (between-subject factor/level 1)
and time/diet/time*diet (within-subject factors/level 2) have degrees of
freedom equal to the within-subject factors? Shouldn't the smallest df in
the interaction be used?

2) Why is the contrasts degrees of freedom weighted in the estimable-output
from R? (30*0.59375=17.8125 and 90*0.59375=53.4375) Is this correct? Should
we be "punished" for wanting to use the mean?

3) If I'm estimating the same contrasts for male and female separately
(male=1, female=0), is it correct to do like estimable() and let A0, A6, B0
and B6 have df=90 for female and df=30 for male, just because of the
parametrization? If I had chosen to let male=0 and female=1 the results
would be different! Or even if I use "sum to zero" contrast coding (male=-1,
female=1), instead of "treatment" coding, A0, A6, B0 and B6 both for female
and male have df=30. But the linear functions estimated for the mean (0)
will have df=90.

I hope one of you you have the time to answer my questions, or refer  
to literature in these areas. Thank you so much for you're time!

Kind regards,
Eirin Ostgard.




More information about the R-sig-mixed-models mailing list