Andrew Koeser
2013-Feb-26 23:53 UTC
[R] Getting the correct factor level as Dunnett control in glht()
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.sortedID 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]]
MacQueen, Don
2013-Feb-27 21:07 UTC
[R] Getting the correct factor level as Dunnett control in glht()
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.sortedID 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.