[R] Getting the correct factor level as Dunnett control in glht()
MacQueen, Don
macqueen1 at llnl.gov
Wed Feb 27 22:07:07 CET 2013
I think you can probably specify the base level using the contrMat()
function, which has a 'base=' argument.
Disclaimer:
This is a fragment from something I did recently, and is not intended to
be reproducible by anyone else. Rather, it is intended to provide a hint
as to how to use the contrMat() function for Andrew's key question. (so no
side comments about posting etiquette, please)
fita <- aov(vv ~ tfac, data=vsub)
ngrp <- table(vsub$tfac)
cm <- contrMat(ngrp)
dt.lt <- glht(fita, linfct=mcp(tfac=cm), alternative='less')
-Don
--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
On 2/26/13 3:53 PM, "Andrew Koeser" <arborkoeser at yahoo.com> wrote:
Hello all,
I would like to do a Dunnett test in glht(). However, the factor level I
want to use as the control is not the first.
dunn1<-glht(model3, linfct = mcp(Container = "Dunnett"), alternative =
"less")
The factor container has 8 levels, so it would be nice not to manually
enter in all of the contrasts. I originally discovered glht() when
working with a glm model and like the versatility it offers. I there a
place to enter a base level? I have not had any luck with this question
looking at online documentation. Is the book for multcomp a good reference?
Thanks to Dalgaard's book, I was able stumble through reordering the data.
>levels(petunia$Container)
[1] "bioplastic" "coir" "cow" "fertil" "net"
"peat" "plastic" "rice" "soilwrap"
[10] "straw"
>petunia$Container <-
>factor(petunia$Container,levels(petunia$Container)[c(7,1:6,8:10)])
>petunia2<-order(petunia$Container)
>petunia.sorted<-petunia[petunia2,]
>petunia.sorted
ID Container TotalWater LeafArea Distance Fresh Dry Bag Biomass
121 15 plastic 1961.00 332.0 109.10 42.700 11.00 7.90 3.10
122 38 plastic 2153.00 552.0 87.90 56.700 12.00 7.80 4.20
123 46 plastic 1880.00 394.0 81.10 48.100 11.60 7.90 3.70
124 66 plastic 2267.00 214.0 62.50 35.800 11.50 8.00 3.50
125 84 plastic 2150.50 362.0 46.50 40.400 10.60 7.80 2.80
126 86 plastic 2034.00 494.0 47.60 53.700 11.90 7.80 4.10
127 87 plastic 2342.00 383.5 47.05 44.700 12.80 7.90 4.90
128 96 plastic 1952.00 273.0 46.50 35.700 10.50 7.90 2.60
129 98 plastic 2411.00 532.0 37.10 58.400 12.20 7.90 4.30
130 106 plastic 2162.00 314.0 31.20 41.200 10.90 7.90 3.00
131 116 plastic 2531.00 389.0 19.40 51.000 12.10 7.90 4.20
132 119 plastic 2129.00 346.0 25.80 41.600 11.20 7.90 3.30
133 122 plastic 2429.00 464.0 21.30 51.600 11.90 7.90 4.00
134 123 plastic 2342.00 303.0 16.20 43.200 11.30 7.80 3.50
135 148 plastic 1709.00 331.0 135.30 40.800 10.80 7.80 3.00
136 152 plastic 2143.00 332.0 120.70 41.300 12.50 7.90 4.60
137 165 plastic 2213.50 440.5 122.80 51.300 12.80 7.90 4.90
138 178 plastic 2284.00 549.0 124.90 61.300 12.90 7.90 5.00
139 195 plastic 1994.00 314.0 111.10 39.300 10.80 7.90 2.90
140 199 plastic 2191.00 363.0 90.70 44.900 11.70 7.80 3.90
1 14 bioplastic 2120.00 433.0 108.20 48.800 11.60 7.90 3.70
2 20 bioplastic 1643.00 375.0 101.00 39.300 10.00 7.90 2.10
3 27 bioplastic 1735.00 432.0 94.70 43.500 10.70 7.90 2.80
After that, I reran the model and all worked out
model3b<-aov(TotalWater~Container, data=petunia.sorted)
dunn2<-glht(model3b, linfct = mcp(Container = "Dunnett"),alternative =
"greater")
summary(dunn2)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Dunnett Contrasts
Fit: aov(formula = TotalWater ~ Container, data = petunia.sorted)
Linear Hypotheses:
Estimate Std. Error t value Pr(>t)
bioplastic - plastic <= 0 -147.55 79.15 -1.864 1.000
coir - plastic <= 0 824.35 79.15 10.414 <1e-04 ***
cow - plastic <= 0 1380.28 79.15 17.438 <1e-04 ***
fertil - plastic <= 0 1572.60 79.15 19.868 <1e-04 ***
net - plastic <= 0 845.20 79.15 10.678 <1e-04 ***
peat - plastic <= 0 786.20 79.15 9.933 <1e-04 ***
rice - plastic <= 0 -36.62 79.15 -0.463 0.968
soilwrap - plastic <= 0 375.75 79.15 4.747 <1e-04 ***
straw - plastic <= 0 1358.20 79.15 17.159 <1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
Is there a better/easier way of doing this? I am showing my labmates R
as part of our weekly meetings (we are a SAS department). I think this
may scare them.
Thanks,
Andrew
--
Research and Teaching Assistant
Department of Crop Sciences
University of Illinois at Champaign-Urbana
[[alternative HTML version deleted]]
______________________________________________
R-help at 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