Quick question about the usage of glht. I'm working with a data set
from an experiment where the response is bounded at 0 whose variance
increases with the mean, and is continuous. A Gamma error
distribution with a log link seemed like the logical choice, and so
I've modeled it as such.
However, when I use glht to look for differences between groups, I get
significant differences where there are none. Now, I'm all for
eyeballing means +/- 95% CIs. However, I've had reviewers and
committee members all tell me that I needed them. Oy. Here's the
code and some of the sample data that, when visualized, is clearly not
different in the comparisons I'm making, and, yet, glht (at least, how
I'm using it, which might be improper) says that the differences are
there.
Hrm.
I'm guessing I'm just using glht improperly, but, any help would be
appreciated!
trt<-c("d", "b", "c", "a",
"a", "d", "b", "c", "c",
"d", "b", "a")
trt<-as.factor(trt)
resp<-c(0.432368576, 0.265148862, 0.140761439, 0.218506998,
0.105017007, 0.140137615, 0.205552589, 0.081970097, 0.24352179,
0.158875904, 0.150195422, 0.187526698)
#take a gander at the lack of differences
boxplot(resp ~ trt)
#model it
a.glm<-glm(resp ~ trt, family=Gamma(link="log"))
summary(a.glm)
#set up the contrast matrix
contra<-rbind("A v. B" = c(-1,1,0,0),
"A v. C" = c(-1,0,1,0),
"A v. D" = c(-1,0,0,1))
library(multcomp)
summary(glht(a.glm, linfct=contra))
---
Yields:
Linear Hypotheses:
Estimate Std. Error z value p value
A v. B == 0 1.9646 0.6201 3.168 0.00314 **
A v. C == 0 1.6782 0.6201 2.706 0.01545 *
A v. D == 0 2.1284 0.6201 3.433 0.00137 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported)
-Jarrett
----------------------------------------
Jarrett Byrnes
Population Biology Graduate Group, UC Davis
Bodega Marine Lab
707-875-1969
http://www-eve.ucdavis.edu/stachowicz/byrnes.shtml
[[alternative HTML version deleted]]
On Mon, 14 Apr 2008, Jarrett Byrnes wrote:> Quick question about the usage of glht. I'm working with a data set > from an experiment where the response is bounded at 0 whose variance > increases with the mean, and is continuous. A Gamma error > distribution with a log link seemed like the logical choice, and so > I've modeled it as such. > > However, when I use glht to look for differences between groups, I get > significant differences where there are none. Now, I'm all for > eyeballing means +/- 95% CIs. However, I've had reviewers and > committee members all tell me that I needed them. Oy. Here's the > code and some of the sample data that, when visualized, is clearly not > different in the comparisons I'm making, and, yet, glht (at least, how > I'm using it, which might be improper) says that the differences are > there. > > Hrm. > > I'm guessing I'm just using glht improperly, but, any help would be > appreciated!Good guess. Compare this to below:> coef(a.glm)-coef(a.glm)[1](Intercept) trtb trtc trtd 0.000000 1.964595 1.678159 2.128366>You are testing the hypothesis that B - 2 * A == 0, etc. HTH, Chuck> > trt<-c("d", "b", "c", "a", "a", "d", "b", "c", "c", "d", "b", "a") > trt<-as.factor(trt) > > resp<-c(0.432368576, 0.265148862, 0.140761439, 0.218506998, > 0.105017007, 0.140137615, 0.205552589, 0.081970097, 0.24352179, > 0.158875904, 0.150195422, 0.187526698) > > #take a gander at the lack of differences > boxplot(resp ~ trt) > > #model it > a.glm<-glm(resp ~ trt, family=Gamma(link="log")) > > summary(a.glm) > > #set up the contrast matrix > contra<-rbind("A v. B" = c(-1,1,0,0), > "A v. C" = c(-1,0,1,0), > "A v. D" = c(-1,0,0,1)) > library(multcomp) > summary(glht(a.glm, linfct=contra)) > --- > Yields: > > Linear Hypotheses: > Estimate Std. Error z value p value > A v. B == 0 1.9646 0.6201 3.168 0.00314 ** > A v. C == 0 1.6782 0.6201 2.706 0.01545 * > A v. D == 0 2.1284 0.6201 3.433 0.00137 ** > --- > Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 > (Adjusted p values reported) > > > -Jarrett > > > > > ---------------------------------------- > Jarrett Byrnes > Population Biology Graduate Group, UC Davis > Bodega Marine Lab > 707-875-1969 > http://www-eve.ucdavis.edu/stachowicz/byrnes.shtml > > > [[alternative HTML version deleted]] > >Charles C. Berry (858) 534-2098 Dept of Family/Preventive Medicine E mailto:cberry at tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
To reply, in case anyone reads this, the problem was of course in the setup
of the matrix, and there are two possible solutions.
The first, for a model with only a single set of groupings, is to use mcp.
So, even with this contrast matrix
contra<-rbind("A v. B" = c(-1,1,0,0),
"A v. C" = c(-1,0,1,0),
"A v. D" = c(-1,0,0,1))
The following will produce the desired analysis:
summary(glht(a.glm, linfct=mcp(trt=contra)))
If one wants to use the coefficients, however, one would need something like
as follows
contrb<-rbind("A v. B" = c(0,1,0,0),
"A v. C" = c(0,0,1,0),
"A v. D" = c(0,0,0,1))
Note, that coefficients are important in the case of factorial designs (e.g.
if there was a block and block:trt effect included in the model). In this
case, one needs to look at the coefficients and setup the appropriate
contrast as on contr.
At least, I think so. I have not yet found a way to use mcp for factorial
designs. Does one exist?
-Jarrett
jebyrnes wrote:>
> Quick question about the usage of glht. I'm working with a data set
> from an experiment where the response is bounded at 0 whose variance
> increases with the mean, and is continuous. A Gamma error
> distribution with a log link seemed like the logical choice, and so
> I've modeled it as such.
>
> I'm guessing I'm just using glht improperly, but, any help would be
> appreciated!
>
> trt<-c("d", "b", "c", "a",
"a", "d", "b", "c", "c",
"d", "b", "a")
> trt<-as.factor(trt)
>
> resp<-c(0.432368576, 0.265148862, 0.140761439, 0.218506998,
> 0.105017007, 0.140137615, 0.205552589, 0.081970097, 0.24352179,
> 0.158875904, 0.150195422, 0.187526698)
>
> #take a gander at the lack of differences
> boxplot(resp ~ trt)
>
> #model it
> a.glm<-glm(resp ~ trt, family=Gamma(link="log"))
>
> summary(a.glm)
>
> #set up the contrast matrix
> contra<-rbind("A v. B" = c(-1,1,0,0),
> "A v. C" = c(-1,0,1,0),
> "A v. D" = c(-1,0,0,1))
> library(multcomp)
> summary(glht(a.glm, linfct=contra))
> ---
> Yields:
>
> Linear Hypotheses:
> Estimate Std. Error z value p value
> A v. B == 0 1.9646 0.6201 3.168 0.00314 **
> A v. C == 0 1.6782 0.6201 2.706 0.01545 *
> A v. D == 0 2.1284 0.6201 3.433 0.00137 **
> ---
> Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> (Adjusted p values reported)
>
>
>
--
View this message in context:
http://www.nabble.com/glht-with-a-glm-using-a-Gamma-distribution-tp16694729p16709118.html
Sent from the R help mailing list archive at Nabble.com.