[R] degrees of freedom

[Ben Bolker]

Iain Gallagher <iaingallagher <at> btopenworld.com> writes:

> 
> Hello List
> 
> I have been asked to analyse some data for a colleague. 
> The design consists of a two sets of animals 
> 
> First set of three - one leg is treated and the
> other is not under two different conditions (control &
> overload are the same animals - control leg is control
> (!) for treated leg; 
> 
> Second set of three - one leg is treated and the other
> is not under two different conditions (high_fat and
> high_fat_overload are the same animals with high_fat 
> being control leg for high_fat_overload).
> 
> Ideally I'd like to find differences between the treatments.
 
bip <- structure(list(group = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 
3L, 3L, 4L, 4L, 4L), .Label = c("control", "overload", "high_fat", 
"high_fat_overload"), class = "factor"), variable = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "BiP", class = "factor"), 
animal = structure(c(1L, 3L, 5L, 1L, 3L, 5L, 2L, 4L, 6L, 
2L, 4L, 6L), .Label = c("rat1_c", "rat1_hf", "rat2_c", "rat2_hf", 
"rat3_c", "rat3_hf"), class = "factor"), value = c(404979.65625, 
783511.8125, 677277.625, 1576900.375, 1460101.875, 1591022, 
581313.75, 992724.1875, 1106941.5, 996600.375, 1101696.5, 
1171004.375)), .Names = c("group", "variable", "animal", 
"value"), row.names = c(NA, 12L), class = "data.frame")

> 
> I chose to analyse this as a mixed effects model with treatment
> as a fixed effect and animal as random.
>

library(lme4)
model1 <- lmer(value~group + (1|animal), data=bip)
summary(model1)
 
> And then compare this to no treatment with:
 
anova(model1)
 
> From this I wanted to work out whether 'treatment' 
> was significantly affecting BiP levels by calculating
> the critical value of F for this design. I have 2
>  groups of animals and 3 animals per group. My calculation
> for the degrees of freedom for treatment is 4-1=3.
> 
> I'm not sure about the degrees of freedom for the denominator 
> though. Since I'm comparing a model with
> treatment to one without (i.e. the grand mean) would the 
> df for my denominator be 6-1=5?
> 
> So I'd then have:
> 
> qf(0.95,3,5)
> 
> for my critical F value?
> 
> Best
> 
> iain


  I started to answer this, but then realized I'd really recommend that
you re-post this to r-sig-mixed-models at r-project.org.  I have a couple
of points for you to think about that might help:

  * since you only have two treatments, I think you can analyze this
as a _paired_ model, that is, reduce the data to (treatment-control,
i.e. overload - non_overload) for each animal.  Then you'll have 6 data
points, 3 in each group, and you can just do a regular 1-way ANOVA on
them.
  * I *think* you've only got 1 df for treatment
  * You might also be able to handle this problem via aov() with
an Error stratum

  Ben Bolker