[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