Hi,
I have a few questions about glht() and the interpretation of output from
Tukey's in multcomp package for lme() model.
The main issue is that I noticed that a plot that I produced with code
letters seem to contradict the graph itself. I provide data and code
below. I end with my questions.
A few things about data set. "LMA.vcp" is continuous response variable.
"Canopy.position.vcp" is independent categorical variable and
"Leaf.height.vcp" is independent covariate.
"plots.vcp" and "tree.vcp" are random effects.
> data.vcp ##data-set##
LMA.vcp Canopy.position.vcp Leaf.height.vcp plots.vcp tree.vcp
1 96.8 BOTTOM 14.5 2 1
2 104.4 BOTTOM 14.7 2 2
3 95.0 BOTTOM 14.7 2 3
4 105.5 BOTTOM 17.4 4 1
5 105.1 BOTTOM 17.6 4 2
6 94.2 BOTTOM 12.6 11 1
7 85.1 BOTTOM 13.0 11 2
8 101.6 BOTTOM 13.5 11 3
9 95.8 BOTTOM 14.0 12 4
10 91.8 BOTTOM 14.7 12 5
11 101.8 BOTTOM 16.7 17 5
12 108.7 BOTTOM 16.9 17 6
13 83.5 BOTTOM 14.5 19 6
14 86.7 BOTTOM 14.7 19 7
15 88.5 BOTTOM 13.1 21 4
16 103.9 BOTTOM 13.0 21 5
17 98.2 BOTTOM 13.7 21 6
18 121.5 BOTTOM 17.1 23 3
19 112.6 BOTTOM 16.7 23 4
20 93.0 BOTTOM 13.1 28 7
21 84.6 BOTTOM 12.1 28 8
22 86.6 BOTTOM 14.3 29 8
23 108.3 BOTTOM 17.0 31 7
24 78.1 BOTTOM 17.4 31 8
25 100.0 MIDDLE 16.7 2 1
26 99.8 MIDDLE 16.6 2 2
27 82.8 MIDDLE 16.7 2 3
28 103.2 MIDDLE 19.4 4 1
29 104.7 MIDDLE 19.9 4 2
30 80.7 MIDDLE 14.9 11 1
31 87.8 MIDDLE 15.0 11 2
32 88.1 MIDDLE 15.6 11 3
33 85.0 MIDDLE 16.4 12 4
34 98.9 MIDDLE 17.0 12 5
35 69.4 MIDDLE 19.1 17 5
36 106.6 MIDDLE 18.9 17 6
37 96.8 MIDDLE 16.6 19 6
38 96.3 MIDDLE 16.7 19 7
39 83.3 MIDDLE 15.2 21 4
40 96.1 MIDDLE 15.3 21 5
41 91.2 MIDDLE 15.6 21 6
42 118.9 MIDDLE 19.2 23 3
43 111.6 MIDDLE 18.6 23 4
44 88.6 MIDDLE 15.2 28 7
45 96.0 MIDDLE 14.4 28 8
46 92.2 MIDDLE 16.3 29 8
47 114.9 MIDDLE 19.1 31 7
48 106.1 MIDDLE 19.3 31 8
49 95.4 TOP 18.9 2 1
50 92.9 TOP 18.4 2 2
51 93.4 TOP 18.7 2 3
52 90.3 TOP 21.4 4 1
53 89.8 TOP 22.2 4 2
54 82.8 TOP 17.2 11 1
55 75.2 TOP 17.0 11 2
56 85.4 TOP 17.7 11 3
57 99.9 TOP 18.8 12 4
58 90.4 TOP 19.3 12 5
59 93.5 TOP 21.4 17 5
60 98.7 TOP 20.9 17 6
61 85.5 TOP 18.7 19 6
62 97.7 TOP 18.7 19 7
63 72.2 TOP 17.3 21 4
64 93.8 TOP 17.5 21 5
65 79.2 TOP 17.4 21 6
66 118.3 TOP 21.3 23 3
67 101.6 TOP 20.5 23 4
68 83.9 TOP 17.3 28 7
69 85.7 TOP 16.7 28 8
70 94.4 TOP 18.3 29 8
71 117.0 TOP 21.2 31 7
72 101.1 TOP 21.2 31 8
##Model##
>modelVCP1<-lme(LMA.vcp~Canopy.position.vcp+Leaf.height.vcp,random=~1 |
plots.vcp/tree.vcp,data=data.vcp,na.action=na.omit) ##lme model##
>modelVCP1_glht<-glht(modelVCP1,linfct=mcp(Canopy.position.vcp="Tukey"))
>summary(modelVCP1_glht) ##summary output of Tukey's test##
*Simultaneous Tests for General Linear Hypotheses*
*
*
*Multiple Comparisons of Means: Tukey Contrasts*
*
*
*
*
*Fit: lme.formula(fixed = LMA.vcp ~ Canopy.position.vcp + Leaf.height.vcp, *
* data = data.vcp, random = ~1 | plots.vcp/tree.vcp, na.action = na.omit)
*
*
*
*Linear Hypotheses:*
* Estimate Std. Error z value Pr(>|z|) *
*MIDDLE - BOTTOM == 0 -8.926 2.777 -3.214 0.00383 ** *
*TOP - BOTTOM == 0 -19.817 4.080 -4.857 < 0.001 ****
*TOP - MIDDLE == 0 -10.891 2.769 -3.934 < 0.001 ****
*---*
*Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 *
*(Adjusted p values reported -- single-step method)*
*
*
*
*
Based on the output it appears that the mean (LMA) is greatest at the
bottom and lowest at the bottom. However, when I receive the code letters
(a,b,c, etc) from
line and from a graph, the output suggests otherwise.
>modelVCP1_glht_cld<-cld(modelVCP1_glht,level=0.05,decreasing=FALSE)
>modelVCP1_glht_cld
>old.par<-par(mai=c(1,1,1.25,1))
>plot(modelVCP1_glht_cld,col=c("yellow","red","blue"))
>par(old.par)
First Question:
As you can see, the code letters in the boxplot graph seem to contradict
the plot itself as well as the Tukey output that I displayed above.
Can anyone explain why this may be happening?
Second Question:
Are the estimate values in the output the difference between the mean of
one level and another level?
Third Question:
If results from the Tukey's test are reported, does one typically report
the difference between levels or the means with code letters? If so, how
are the means extracted from
the model "modelVCP1_glht"?
Thanks!
--
Adam P. Coble
Ph.D. Student
Michigan Technological University
School of Forest Resources and Environmental Science
Houghton, MI 49931
[[alternative HTML version deleted]]
______________________________________________
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.
