As a supplement to Bens reply: The pbkrtest package allows calculation of degrees of freedom by the Kenward-Roger method; something like
library(lme4) model1 <- lmer(value~group + (1|animal), data=bip) summary(model1) anova(model1) model0 <- update(model1, .~.-group) anova(model1, model0) library(pbkrtest) KRmodcomp(model1, model0) PBmodcomp(model1, model0) The doBy package provides facilities for calculating LSmeans (se the vignette) for groupwise comparisons (and so does the lsmeans) package. See also facilities in the multcomp packages. Cheers Søren -----Original Message----- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Ben Bolker Sent: 23. januar 2014 23:35 To: r-h...@stat.math.ethz.ch Subject: Re: [R] degrees of freedom 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-mod...@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 ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
