[R] LME and overall treatment effects
ONKELINX, Thierry
Thierry.ONKELINX at inbo.be
Fri Jul 15 11:06:49 CEST 2011
Dear Mark,
Interpreting one of the main effects when they are part of an interaction is, AFAIK, not possible.
Your statement about comparing treatments when Year is continuous is not correct. The parameters of treatment assume that Year == 0! Which might lead to very strange effect when year is not centered to a year close to the observed years. Have a look at the example below
set.seed(123456)
dataset <- expand.grid(cYear = 0:20, Treatment = factor(c("A", "B", "C")), Obs = 1:3)
dataset$Year <- dataset$cYear + 2000
Trend <- c(A = 1, B = 0.1, C = -0.5)
TreatmentEffect <- c(A = 2, B = -1, C = 0.5)
sdNoise <- 1
dataset$Value <- with(dataset, TreatmentEffect[Treatment] + Trend[Treatment] * cYear) + rnorm(nrow(dataset), sd = sdNoise)
lm(Value ~ Year * Treatment, data = dataset)
lm(Value ~ cYear * Treatment, data = dataset)
If you want to focus on the treatment effect alone but take the year effect into account, then add Year as a random effect.
library(lme4)
lmer(Value ~ 0 + Treatment + (0 + Treatment|Year), data = dataset)
In your case you want to cross the random effect of year with those of plot. Crossed random effects are hard to do with the nlme package but easy with the lme4 package.
Model <- lmer(Species~ 0 + Treatment + (0 + Treatment|Year) + (1|Plot/Quadrat) ,na.action =na.omit,data=UDD)
Best regards,
Thierry
----------------------------------------------------------------------------
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium
Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
Thierry.Onkelinx op inbo.be
www.inbo.be
To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
> -----Oorspronkelijk bericht-----
> Van: r-help-bounces op r-project.org [mailto:r-help-bounces op r-project.org]
> Namens Mark Bilton
> Verzonden: donderdag 14 juli 2011 22:05
> Aan: Bert Gunter
> CC: r-help op r-project.org
> Onderwerp: Re: [R] LME and overall treatment effects
>
> Ok...lets try again with some code...
>
> ---------------------------------------------------------------------------
> Hello fellow R users,
>
> I am having a problem finding the estimates for some overall treatment effects
> for my mixed models using 'lme' (package nlme). I hope someone can help.
>
> Firstly then, the model:
> The data: Plant biomass (log transformed) Fixed Factors: Treatment(x3 Dry, Wet,
> Control) Year(x8 2002-2009) Random Factors: 5 plots per treatment, 10 quadrats
> per plot (N=1200 (3*5*10*8)).
>
> I am modelling this in two ways, firstly with year as a continuous variable
> (interested in the difference in estimated slope over time in each treatment
> 'year*treatment'), and secondly with year as a categorical variable (interested in
> difference between 'treatments').
>
> --------------
> ie: (with Year as either numeric or factor)
>
> Model<-lme(Species~Year*Treatment,random=~1|Plot/Quadrat,na.action =
> na.omit,data=UDD)
> ---------------------
>
> When using Year as a continuous variable, the output of the lme means that I
> can compare the 3 treatments within my model...
> i.e. it takes one of the Treatment*year interactions as the baseline and
> compares (contrasts) the other two to that.
> -----------------
> ie
>
> Fixed effects: Species ~ Year * Treatment
> Value Std.Error DF t-value p-value
> (Intercept) 1514.3700 352.7552 1047 4.292978 0.0000
> Year -0.7519 0.1759 1047 -4.274786 0.0000
> Treatment0 -461.9500 498.8711 12 -0.925991 0.3727
> Treatment1 -1355.0450 498.8711 12 -2.716222 0.0187
> Year:Treatment0 0.2305 0.2488 1047 0.926537 0.3544
> Year:Treatment1 0.6776 0.2488 1047 2.724094 0.0066
>
> so Year:Treatment0 differs from baseline Year:Treatment-1 by 0.2305 and
> Year:Treatment1 is significantly different (p=0.0066) from
> Year:Treatment-1
> ----------------------
>
> I can then calculate the overall treatment*year effect using 'anova.lme(Model)'.
> -------------------------
> > anova.lme(Model1)
> numDF denDF F-value p-value
> (Intercept) 1 1047 143.15245 <.0001
> Year 1 1047 19.56663 <.0001
> Treatment 2 12 3.73890 0.0547
> Year:Treatment 2 1047 3.83679 0.0219
>
> so there is an overall difference in slope between treatments (Year:Treatment
> interaction) p=0.0219
> --------------------------
>
> However, the problem comes when I use Year as a categorical variable.
> Here, I am interested solely in the Treatment effect (not the interaction with
> year). However, the output for the two labelled 'Treatment's against the third
> comparison, are not the overall effect but are a comparison within a year
> (2002).
> --------------------------
> Fixed effects: Species ~ Year * Treatment
> Value Std.Error DF t-value p-value
> (Intercept) 6.42 1.528179 1029 4.201079 0.0000
> Year2003 4.10 1.551578 1029 2.642471 0.0084
> Year2004 5.00 1.551578 1029 3.222526 0.0013
> Year2005 -1.52 1.551578 1029 -0.979648 0.3275
> Year2006 -3.08 1.551578 1029 -1.985076 0.0474
> Year2007 -2.40 1.551578 1029 -1.546813 0.1222
> Year2008 2.24 1.551578 1029 1.443692 0.1491
> Year2009 -4.30 1.551578 1029 -2.771372 0.0057
> Treatment0 0.46 2.161171 12 0.212848 0.8350
> Treatment1 0.50 2.161171 12 0.231356 0.8209
> Year2003:Treatment0 -2.46 2.194262 1029 -1.121106 0.2625
> Year2004:Treatment0 -1.34 2.194262 1029 -0.610684 0.5415
> Year2005:Treatment0 0.34 2.194262 1029 0.154950 0.8769
> Year2006:Treatment0 1.60 2.194262 1029 0.729174 0.4661
> Year2007:Treatment0 1.76 2.194262 1029 0.802092 0.4227
> Year2008:Treatment0 -3.22 2.194262 1029 -1.467464 0.1426
> Year2009:Treatment0 1.80 2.194262 1029 0.820321 0.4122
> Year2003:Treatment1 0.22 2.194262 1029 0.100261 0.9202
> Year2004:Treatment1 3.48 2.194262 1029 1.585954 0.1131
> Year2005:Treatment1 5.00 2.194262 1029 2.278670 0.0229
> Year2006:Treatment1 4.70 2.194262 1029 2.141950 0.0324
> Year2007:Treatment1 6.08 2.194262 1029 2.770863 0.0057
> Year2008:Treatment1 2.32 2.194262 1029 1.057303 0.2906
> Year2009:Treatment1 5.56 2.194262 1029 2.533881 0.0114
>
> so Treatment0 (in year 2002) is different to baseline Treatment-1 (in year 2002)
> by 0.46 p=0.8350
> ----------------------------------------
>
> I can still get my overall effect (using anova.lme) but how do I calculate the
> estimates (with P-Values if possible) for each seperate overall treatment
> comparison (not within year). I can do this in SAS using 'estimate' or 'lsmeans',
> but for various reasons I want to do it in R as well.
> I tried 'glht' (package 'multicomp') but this only works if there are no interactions
> present, otherwise again it gives a comparison for one particular year.
>
> --------------------------------------
> ie
> require(multcomp)
> summary(glht(Model, linfct = mcp(Treatment = "Tukey")))
>
> (sorry, I can't get this to work at home but trust me that it's the same as for the
> contrasts in the fixed effects outputs:
> ie Treatment0 (in year 2002) is different to baseline Treatment-1 (in year 2002)
> by 0.46 p=0.8350)
> ---------------------------------------
>
>
> Very grateful for any assistance,
> Mark
>
>
> ----------------------------------
> Hope that helps a bit more m
>
> ______________________________________________
> R-help op 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.
More information about the R-help
mailing list