title: What is the most appropriate method for contrasts in lme (nlme)?
body: <p>I have a question regarding conducting post-hoc contrasts in a mixed glmm model. Essentially, my data is a series of response variables with fixed factors of 'Subject Type' (three types of subjects), 'Treatment' (four treatments), and 'Year' (two years), with 'Subject Number', and 'Location' as random effects.</p>
  
  <p>As I will be discussing estimates below, here are the means and standard deviations by Subject Type and Treatment:</p>
  
  <p>My model, which is working fine, is as follows:</p>
  
  <pre><code>library(nlme)
  
  ctrl &lt;- lmeControl(opt='optim') 
  
  myModel &lt;- lme(y ~ Treatment * SubjectType * Year, random=list(SubjectNumber=~1, Location=~1), control=ctrl, method="REML", data = dframe1, na.action="na.exclude")
  </code></pre>
  
  <p>I then get the following output, which all looks fine:</p>
  
  <pre><code>anova.lme(myModel)
                          numDF denDF  F-value p-value
  (Intercept)                 1   168 3199.405  &lt;.0001
  Treatment                   3   168    5.837  0.0008
  SubjectType                 2   168    7.502  0.0008
  Year                        1   168    9.410  0.0025
  Treatment:SubjectType       6   168    0.933  0.4725
  Treatment:Year              3   168    3.858  0.0106
  SubjectType:Year            2   168    6.322  0.0023
  Treatment:SubjectType:Year  6   168    1.383  0.2239
  </code></pre>
  
  <p>However, I now wish to run post-hoc contrasts to see where the differences lie within the main effects of 'Treatment' and 'SubjectType' (i.e. 'Is Subject Type 1 different from both Subject Types 2 and 3, or just 2/3?'.</p>
  
  <p>I have done a fair bit of reasearch and tried several methods, but I seem to be a bit lost and everything is giving me different results, so some clarification would be greatly appreciated. </p>
  
  <p>Initially, I used Tukey's with the following code, which I was happy with, but I have been advised that this doesn't take into account variation within the other main effects. Here is the code and output that I was using initially, and lsmeans are provided to reference when I speak of estimates in differences later in the question:</p>
  
  <pre><code>library(lsmeans)
  
  lsmeans(myModel, pairwise~SubjectType, adjust="tukey")
  
  NOTE: Results may be misleading due to involvement in interactions
  $lsmeans
   SubjectType     lsmean         SE  df  lower.CL  upper.CL
   SubjectType1    0.5531300 0.01853118 168 0.5165460 0.5897140
   SubjectType2    0.6078394 0.01853118 168 0.5712554 0.6444234
   SubjectType3    0.6545389 0.01853118 168 0.6179549 0.6911229
  
  Results are averaged over the levels of: Treatment 
  Confidence level used: 0.95 
  
  $contrasts
   contrast                        estimate         SE  df t.ratio p.value
   SubjectType1 - SubjectType2    -0.05470938 0.02620704 168  -2.088  0.0955
   SubjectType1 - SubjectType3    -0.10140890 0.02620704 168  -3.870  0.0005
   SubjectType2 - SubjectType3    -0.04669951 0.02620704 168  -1.782  0.1788
  
  Results are averaged over the levels of: Treatment 
  P value adjustment: tukey method for comparing a family of 3 estimates 
  
  
  
  lsmeans(myModel, pairwise~Treatment, adjust="tukey")
  
  NOTE: Results may be misleading due to involvement in interactions
  $lsmeans
   Treatment         lsmean         SE  df  lower.CL  upper.CL
   Treatment1        0.6379318 0.02139796 168 0.5956883 0.6801754
   Treatment2        0.5299945 0.02139796 168 0.4877510 0.5722380
   Treatment3        0.6123516 0.02139796 168 0.5701081 0.6545952
   Treatment4        0.6403998 0.02139796 168 0.5981563 0.6826434
  
  Results are averaged over the levels of: SubjectType 
  Confidence level used: 0.95 
  
  $contrasts
   contrast                           estimate         SE  df t.ratio p.value
   Treatment1 - Treatment2            0.107937332 0.03026128 168   3.567  0.0026
   Treatment1 - Treatment3            0.025580191 0.03026128 168   0.845  0.8327
   Treatment1 - Treatment4           -0.002468028 0.03026128 168  -0.082  0.9998
   Treatment2 - Treatment3           -0.082357141 0.03026128 168  -2.722  0.0358
   Treatment2 - Treatment4           -0.110405360 0.03026128 168  -3.648  0.0020
   Treatment3 - Treatment4           -0.028048219 0.03026128 168  -0.927  0.7905
  
  Results are averaged over the levels of: SubjectType 
  P value adjustment: tukey method for comparing a family of 4 estimates 
  </code></pre>
  
  <p>The estimates in this do look to be correct based on the data, and those flagged up as significant do appear to be different from each other when graphed. I have also done a similar comparison with the 'multcomp' package, which gives different results, and the estimated differences do not look to be correct:</p>
  
  <pre><code>library("multcomp")
  
  Contrasts2 &lt;- glht(myModel, linfct = mcp(SubjectType = "Tukey"))
  
  Warning message:
  In mcp2matrix(model, linfct = linfct) :
    covariate interactions found -- default contrast might be inappropriate
  
  summary(Contrasts2)
  
    Simultaneous Tests for General Linear Hypotheses
  
  Multiple Comparisons of Means: Tukey Contrasts
  
  
  Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1, 
      random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML", 
      na.action = "na.exclude", control = ctrl)
  
  Linear Hypotheses:
                                       Estimate Std. Error z value Pr(&gt;|z|)
  SubjectType2 - SubjectType1 == 0     0.26395    0.16575   1.593    0.249
  SubjectType3 - SubjectType1 == 0     0.31642    0.16575   1.909    0.136
  SubjectType3 - SubjectType2 == 0     0.05247    0.16575   0.317    0.946
  (Adjusted p values reported -- single-step method)
  
  ContrastsTreatment &lt;- glht(myModel, linfct = mcp(Treatment = "Tukey"))
  
  Warning message:
  In mcp2matrix(model, linfct = linfct) :
    covariate interactions found -- default contrast might be inappropriate
  
  summary(ContrastsTreatment)
  
       Simultaneous Tests for General Linear Hypotheses
  
  Multiple Comparisons of Means: Tukey Contrasts
  
  
  Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1, 
      random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML", 
      na.action = "na.exclude", control = ctrl)
  
  Linear Hypotheses:
                                      Estimate Std. Error z value Pr(&gt;|z|)   
  Treatment2 - Treatment1 == 0        -0.49321    0.16575  -2.976  0.01542 * 
  Treatment3 - Treatment1 == 0         0.08369    0.16575   0.505  0.95794   
  Treatment4 - Treatment1 == 0        -0.02263    0.16575  -0.137  0.99909   
  Treatment3 - Treatment2 == 0         0.57690    0.16575   3.481  0.00294 **
  Treatment4 - Treatment2 == 0         0.47059    0.16575   2.839  0.02373 * 
  Treatment4 - Treatment3 == 0        -0.10631    0.16575  -0.641  0.91856   
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  (Adjusted p values reported -- single-step method)
  </code></pre>
  
  <p>So then I tried doing a contrast matrix with glht, which again gives rather different results, and again the estimates do not look correct based on the data:</p>
  
  <pre><code>contrast.matrix.SubjectType &lt;- rbind(
   `SubjectType1 vs. Soanta` = c(1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
   `SubjectType1 vs. SubjectType3` = c(1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
   `SubjectType2 vs. SubjectType3` = c(0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
  
   compsSubjectType &lt;- glht(myModel, contrast.matrix.SubjectType)
  
   summary(compsSubjectType)
  
       Simultaneous Tests for General Linear Hypotheses
  
  Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1, 
      random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML", 
      na.action = "na.exclude", control = ctrl)
  
  Linear Hypotheses:
                                       Estimate Std. Error z value Pr(&gt;|z|)
  SubjectType1 vs. Soanta == 0         0.33426    0.26207   1.275    0.401
  SubjectType1 vs. SubjectType3 == 0   0.28180    0.26207   1.075    0.522
  SubjectType2 vs. SubjectType3 == 0  -0.05247    0.16575  -0.317    0.945
  (Adjusted p values reported -- single-step method)
  
   contrast.matrix.Treatment &lt;- rbind(
   `Treatment1 vs. Treatment 1` = c(1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
   `Treatment1 vs. Treatment 2` = c(1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
   `Treatment1 vs. Treatment4` = c(1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
   `Treatment 1 vs. Treatment 2` = c(0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
   `Treatment 1 vs. Treatment4` = c(0, 1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
   `Treatment 2 vs. Treatment4` = c(0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
  
   compsTreatment &lt;- glht(myModel, contrast.matrix.Treatment)
   summary(compsTreatment)
  
       Simultaneous Tests for General Linear Hypotheses
  
  Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1, 
      random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML", 
      na.action = "na.exclude", control = ctrl)
  
  Linear Hypotheses:
                                       Estimate Std. Error z value Pr(&gt;|z|)    
  Treatment1 vs. Treatment 1 == 0       1.0914     0.2621   4.165  &lt; 0.001 ***
  Treatment1 vs. Treatment 2 == 0       0.5145     0.2621   1.963  0.19465    
  Treatment1 vs. Treatment4 == 0        0.6208     0.2621   2.369  0.07935 .  
  Treatment 1 vs. Treatment 2 == 0     -0.5769     0.1657  -3.481  0.00268 ** 
  Treatment 1 vs. Treatment4 == 0      -0.4706     0.1657  -2.839  0.02195 *  
  Treatment 2 vs. Treatment4 == 0       0.1063     0.1657   0.641  0.91584    
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  (Adjusted p values reported -- single-step method)
  </code></pre>
  
  <p>Changing the correction obviously influences the results, but again, the estimates don't look right:</p>
  
  <pre><code>summary(glht(myModel, contrast.matrix.SubjectType), test = adjusted("none"))
  
       Simultaneous Tests for General Linear Hypotheses
  
  Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1, 
      random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML", 
      na.action = "na.exclude", control = ctrl)
  
  Linear Hypotheses:
                                        Estimate Std. Error z value Pr(&gt;|z|)
  SubjectType1 vs. Soanta == 0          0.33426    0.26207   1.275    0.202
  SubjectType1 vs. SubjectType3 == 0    0.28180    0.26207   1.075    0.282
  SubjectType2 vs. SubjectType3 == 0   -0.05247    0.16575  -0.317    0.752
  (Adjusted p values reported -- none method)
  
  summary(glht(myModel, contrast.matrix.Treatment), test = adjusted("none"))
  
       Simultaneous Tests for General Linear Hypotheses
  
  Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1, 
      random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML", 
      na.action = "na.exclude", control = ctrl)
  
  Linear Hypotheses:
                                       Estimate Std. Error z value Pr(&gt;|z|)    
  Treatment1 vs. Treatment 1 == 0       1.0914     0.2621   4.165 3.12e-05 ***
  Treatment1 vs. Treatment 2 == 0       0.5145     0.2621   1.963  0.04961 *  
  Treatment1 vs. Treatment4 == 0        0.6208     0.2621   2.369  0.01784 *  
  Treatment 1 vs. Treatment 2 == 0     -0.5769     0.1657  -3.481  0.00050 ***
  Treatment 1 vs. Treatment4 == 0      -0.4706     0.1657  -2.839  0.00452 ** 
  Treatment 2 vs. Treatment4 == 0       0.1063     0.1657   0.641  0.52125    
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  (Adjusted p values reported -- none method)
  </code></pre>
  
  <p>And finally, the model comparisons in the model summary do seem to fit what the data is showing, but again provide different p values and the estimates don't seem to match the data:</p>
  
  <pre><code>summary(myModel)
  
  Linear mixed-effects model fit by REML
   Data: dframe1 
          AIC      BIC  logLik
    0.9785978 104.0406 26.5107
  
  Random effects:
   Formula: ~1 | SubjectNumber
           (Intercept)
  StdDev: 0.0001623009
  
   Formula: ~1 | SubjectLocation %in% SubjectNumber
           (Intercept)  Residual
  StdDev: 0.0001623009 0.2029986
  
  Fixed effects: y ~ Treatment * SubjectType * Year 
                                                        Value  Std.Error  DF   t-value p-value
  (Intercept)                                       0.5982163 0.11720132 168  5.104177  0.0000
  TreatmentTreatment2                              -0.4932126 0.16574769 168 -2.975683  0.0034
  TreatmentTreatment3                               0.0836861 0.16574769 168  0.504900  0.6143
  TreatmentTreatment4                              -0.0226273 0.16574769 168 -0.136517  0.8916
  SubjectTypeSubjectType2                           0.2639539 0.16574769 168  1.592504  0.1132
  SubjectTypeSubjectType3                           0.3164198 0.16574769 168  1.909045  0.0580
  Year                                              0.0082773 0.07412461 168  0.111667  0.9112
  TreatmentTreatment2:SubjectTypeSubjectType2       0.0807378 0.23440263 168  0.344441  0.7309
  TreatmentTreatment3:SubjectTypeSubjectType2      -0.0983580 0.23440263 168 -0.419611  0.6753
  TreatmentTreatment4:SubjectTypeSubjectType2       0.1477391 0.23440263 168  0.630279  0.5294
  TreatmentTreatment2:SubjectTypeSubjectType3       0.3533741 0.23440263 168  1.507552  0.1335
  TreatmentTreatment3:SubjectTypeSubjectType3      -0.2917698 0.23440263 168 -1.244738  0.2150
  TreatmentTreatment4:SubjectTypeSubjectType3       0.0481905 0.23440263 168  0.205589  0.8374
  TreatmentTreatment2:Year                          0.2465353 0.10482803 168  2.351807  0.0198
  TreatmentTreatment3:Year                         -0.0791173 0.10482803 168 -0.754734  0.4515
  TreatmentTreatment4:Year                         -0.0326548 0.10482803 168 -0.311508  0.7558
  SubjectTypeSubjectType2:Year                     -0.1651260 0.10482803 168 -1.575208  0.1171
  SubjectTypeSubjectType3:Year                     -0.1671907 0.10482803 168 -1.594904  0.1126
  TreatmentTreatment2:SubjectTypeSubjectType2:Year -0.0516816 0.14824922 168 -0.348613  0.7278
  TreatmentTreatment3:SubjectTypeSubjectType2:Year  0.0718031 0.14824922 168  0.484341  0.6288
  TreatmentTreatment4:SubjectTypeSubjectType2:Year -0.0043490 0.14824922 168 -0.029336  0.9766
  TreatmentTreatment2:SubjectTypeSubjectType3:Year -0.2067818 0.14824922 168 -1.394826  0.1649
  TreatmentTreatment3:SubjectTypeSubjectType3:Year  0.2071016 0.14824922 168  1.396983  0.1643
  TreatmentTreatment4:SubjectTypeSubjectType3:Year  0.0218841 0.14824922 168  0.147617  0.8828
  </code></pre>
  
  <p>So my overall question is what is the most appropriate to use for determining comparisons between treatments / subject types? Why are the two Tukey's methods providing different results? Sorry for such a long-winded question, and thank you very much for your time and advice. If any answers could be provided in relatively simple terminology I would appreciate it! Thanks!</p>