R help - Jun 2003 - degrees of freedom in a LME model

Dear All,

I am analysing some data for a colleague (not my data, gotta be published
so I cannot divulge).

My response variable is the number of matings observed per day for some
fruitlies.

My factors are:
Day: the observations were taken on 9 days
Regime: 3 selection regimes
Line: 3 replicates per selection regime.

I have 81 observations in total

The lines are coded A to I, so I do not need to do any extra grouping.

my model is:

anova(lme(Matings ~ Day * Regime, random = ~1| Line/Day, mydata))

I would expect to have:
1 df per Day
2 df per Regime
2 df per Day * Regime
6 df per Line %in% Regime
6 df per Day * Line %in% Regime,


so my anova would have:

	numDF	denDF
int	1	63
Day	1	6
Regime	2	6
D*R	2	6

what I get is:

	numDF	denDF
int	1	69
Day	1	69
Regime	2	6
D*R	2	69

why is lme not calculating correctly the Line/Day interation ?

I am using R 1.7.0 under W2K, although I updated the packages and I get the
warning "nlme lib built under R1.7.1..."

Regards,

Federico 


========================
Federico C.F. Calboli

Department of Biology
University College London
Room 327
Darwin Building
Gower Street
London
WClE 6BT

Tel: (+44) 020 7679 4395 
Fax (+44) 020 7679 7096
f.calboli at ucl.ac.uk

Federico Calboli <f.calboli at ucl.ac.uk> writes:
> Dear All,
> 
> I am analysing some data for a colleague (not my data, gotta be published
> so I cannot divulge).
> 
> My response variable is the number of matings observed per day for some
> fruitlies.
> 
> My factors are:
> Day: the observations were taken on 9 days
> Regime: 3 selection regimes
> Line: 3 replicates per selection regime.
> 
> I have 81 observations in total
> 
> The lines are coded A to I, so I do not need to do any extra grouping.
> 
> my model is:
> 
> anova(lme(Matings ~ Day * Regime, random = ~1| Line/Day, mydata))
> 
> I would expect to have:
> 1 df per Day
> 2 df per Regime
> 2 df per Day * Regime
> 6 df per Line %in% Regime
> 6 df per Day * Line %in% Regime,
> 
> 
> so my anova would have:
> 
> 	numDF	denDF
> int	1	63
> Day	1	6
> Regime	2	6
> D*R	2	6
> 
> what I get is:
> 
> 	numDF	denDF
> int	1	69
> Day	1	69
> Regime	2	6
> D*R	2	69
> 
> why is lme not calculating correctly the Line/Day interation ?
I think your calculation is based on using only within-strata
information whereas lme uses both within-strata and between-strata
information for estimates of effects.

The way that we calculate denominator degrees of freedom is described
on pp. 90-91 of Pinheiro and Bates (2000).  For each term in the
fixed-effects we determine the innermost level of the random effects
at which is it changing.  Because Regime is constant for each Line it
has the fewest degrees of freedom but Day is changing within Line so
terms in Day have more degrees of freedom.  Is this what you intended?

I must admit I am having difficulty understanding the structure of the
experiment but it is still Monday morning for me so perhaps that is
not surprising.
> I am using R 1.7.0 under W2K, although I updated the packages and I get the
> warning "nlme lib built under R1.7.1..."
> 
> Regards,
> 
> Federico 

-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/

Dear Prof. Bates,

Thank you for your reply. I did actually check the book at pages 89-92, but
I have to say I found it a bit "cryptic", if not downright confusing,
for a
genetist like me. Any day of the week. 

To use the example in the book, I cannot see why MACHINE is inner to
WORKER. Or where the Pi = sum of df for the term estimatet at level i... 

To me it's just a two way anova with one random effect + one interaction.
My df correspond to those in the book, but I calculated them in the "bog
standard" way of 2 df for 3 machines, 5 df for 6 workers and 2 * 5 = 10 df
for the interaction...oh well...

I attach the datasest (I did a sample with replace = TRUE, should have
thought of this earlier...). I hope thing would be clearer. It is obvious
that if I consider day a "factor", as I have just one datapoint per
day per
line within regime, I end up using all my degrees of freedom:

2 df for 3 regimes
6 df for line within regime
8 df for day (9 days, day as factor)
16 df for regime * day
48 df for line * day
total 80 df, with 0 df remaining for the error. A bit of a problem I
daresay, but I did not collect the data!

BUT I did calculate my anova considering day as a continuous variable,

2 df for 3 regimes
6 df for line within regime
1 df for day (day as a number)
2 df for regime * day
6 df for line * day
total 17 df, with 63 df remaining for the error. 

I still do not get why the interaction term Regime*Day is not tested on the
interaction Line-in-regime*Day... In my analysis LME is clumping Error and
Line-in-regime*Day together, judging by the df I was talking about in my
previuos email. To me it should be just another "bog standard"
situation...But if I were smart enough to be a statistician I would not be
here doing genetics ;)

Regards,
Federico Calboli



Day	Line	Regime	Matings
2	501	es	0.4
4	501	es	0.32989691
9	501	es	0.48484848
11	501	es	0.72727273
16	501	es	0.34042553
18	501	es	0.56470588
25	501	es	0.37509377
30	501	es	0.22222222
32	501	es	0.77777778
2	502	es	0.57142857
4	502	es	1.06666667
9	502	es	0.16
11	502	es	0.4
16	502	es	0.4
18	502	es	0.33333333
25	502	es	0.80808081
30	502	es	0.48
32	502	es	0.25531915
2	503	es	0
4	503	es	0.72727273
9	503	es	0.3
11	503	es	0.77777778
16	503	es	0.34042553
18	503	es	1.06666667
25	503	es	0.32989691
30	503	es	1.24444444
32	503	es	0.19153725
2	fb1	fb	0.72727273
4	fb1	fb	0.42105263
9	fb1	fb	0.32989691
11	fb1	fb	0.32323232
16	fb1	fb	0.37509377
18	fb1	fb	0.51612903
25	fb1	fb	0.4
30	fb1	fb	0.24742268
32	fb1	fb	0.24742268
2	fb2	fb	1.33333333
4	fb2	fb	0.25263158
9	fb2	fb	0.66666667
11	fb2	fb	1.06666667
16	fb2	fb	0.97959184
18	fb2	fb	0.42105263
25	fb2	fb	0.57142857
30	fb2	fb	1.15555556
32	fb2	fb	0.80808081
2	fb3	fb	0.18952106
4	fb3	fb	0.66666667
9	fb3	fb	1.22033898
11	fb3	fb	0.35955056
16	fb3	fb	1.68421053
18	fb3	fb	0.57461174
25	fb3	fb	0.93506493
30	fb3	fb	0.80808081
32	fb3	fb	0.57461174
2	mb1	mb	0.22222222
4	mb1	mb	0.3902439
9	mb1	mb	0.42105263
11	mb1	mb	0.48484848
16	mb1	mb	0.55555556
18	mb1	mb	1.68421053
25	mb1	mb	0.4
30	mb1	mb	0.16
32	mb1	mb	0.09523809
2	mb2	mb	0.77777778
4	mb2	mb	0.68817204
9	mb2	mb	0.25263158
11	mb2	mb	1.33333333
16	mb2	mb	0.5
18	mb2	mb	0.37509377
25	mb2	mb	0.16494845
30	mb2	mb	2.46153846
32	mb2	mb	0.47058824
2	mb3	mb	1.06666667
4	mb3	mb	0.5
9	mb3	mb	0.5
11	mb3	mb	0.38095238
16	mb3	mb	0.26373626
18	mb3	mb	0.87912088
25	mb3	mb	0.55555556
30	mb3	mb	0.74418605
32	mb3	mb	1.33333333
>
>I think your calculation is based on using only within-strata
>information whereas lme uses both within-strata and between-strata
>information for estimates of effects.
>
>The way that we calculate denominator degrees of freedom is described
>on pp. 90-91 of Pinheiro and Bates (2000).  For each term in the
>fixed-effects we determine the innermost level of the random effects
>at which is it changing.  Because Regime is constant for each Line it
>has the fewest degrees of freedom but Day is changing within Line so
>terms in Day have more degrees of freedom.  Is this what you intended?
>
>I must admit I am having difficulty understanding the structure of the
>experiment but it is still Monday morning for me so perhaps that is
>not surprising.

========================
Federico C.F. Calboli

Department of Biology
University College London
Room 327
Darwin Building
Gower Street
London
WClE 6BT

Tel: (+44) 020 7679 4395 
Fax (+44) 020 7679 7096
f.calboli at ucl.ac.uk

Dear All,

I know I am quite obsessive and downright annoying (I apologize about that,
but it's the way I am), but I would like to get my understanding of the way
nlme calculates degrees of freedom straight.

For instance, on page 91 in Pinheiro and Bates (2000), on the examle of
anova(fm2Machie), how is the sum of DF *Pi* corresponding to the terms
estimated at level *i* calculated? In the example presenting
anova(fm2Machine), *P1 = 0*. I just fail to see why. Same thing for *P2 2*
(although this seems intuitive, but intuitive could be miles off the
real reason) and *P3 = 0*.

Incidentally, should I change the grouping, putting *machine* outside and
*worker* inside, would anything change?

A second thing: is there any substantial difference between the classical
decomposition of DF for an ANOVA and the method used by lme for the
interaction between a fixed and a random effect, in case the random
variable is nominal and the fixed one continuous?

Best regards,
Federico Calboli


========================
Federico C.F. Calboli

Department of Biology
University College London
Room 327
Darwin Building
Gower Street
London
WClE 6BT

Tel: (+44) 020 7679 4395 
Fax (+44) 020 7679 7096
f.calboli at ucl.ac.uk