[R] LME and overall treatment effects

[Mark Bilton]

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