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@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.
