R help - Mar 2009 - Multiple Comparisons for (multicomp - glht) for glm negative binomial (glm.nb)

Hi
 
I have some experimental data where I have counts of the number of
insects collected to different trap types rotated through 5 different
location (variable -location), 4 different chemical attractants [A, B,
C, D]  were applied to the traps (variable - semio) and all were
trialled at two different CO2 release rates [1, 2] (variable CO2) I also
have a selection of continuous variables measuring meteorological
conditions to account for any bias cause by changing weather conditions,
the data is over dispersed so I have fitted a negative binomial glm
(glm.nb) and simplified using stepAIC from the MASS package etc.
 
There are significant differences in the number of insects attracted to
the different chemical (semio) and to the two different CO2 (release
rates) I have then used the glht function from the multcomp package to
do multiple comparisons to see what the specific differences between the
levels are for semio and CO2 using the code below which works great but
what I would like to do is to do comparisons combining the factors e.g a
comparison for semioA at CO at level 1 vs Semio A at CO2 level 2  etc to
see which is the best combination, is this possible or should I have
started of with my counts already split up into this e.g. a treatment
variable(semioA at CO2 level 1 = A1, semioA at CO2 level 2 = A2 etc), I
started with them this way as we have no prior knowledge that increasing
co2 will have any effect.
 
I have had a quick try with the data split into a treatment factor
(instead of semio and CO2 level) but I can not get convergence with
glm.nb I think this may be to do with to many zeros in the data set, do
you know if glht or another multiple comparison will work or
zeroinflated negative binomial regression(zeroinfl() from the pscl
library)?
 
Any help or ideas will be gratefully appreciated. Many thanks in
advance.
 
Lara
 
 
semiochemical <-read.csv("G:/semiochemical_data.csv", header=T)

semiochemical$location2<-factor(semiochemical$location)
levels(semiochemical$location2)<-c("1","2","3","4","5")
 
semiochemical$semio2<-factor(semiochemical$semio)
levels(semiochemical$semio2)<-c("A","B","C","D")
 
semiochemical$CO22<-factor(semiochemical$CO2)
levels(semiochemical$CO22)<-c("1","2")
 
model1<-glm.nb(total ~ semio + CO2 + location + temp + mean.wind.speed)
 
model.glht.Semio <- glht(model1, linfct=mcp(semio="Tukey"))
model.glht.CO2 <- glht(model1, linfct=mcp(CO2="Tukey"))
 

	[[alternative HTML version deleted]]

On Sun, 22 Mar 2009, lara harrup (IAH-P) wrote:
> Hi
>
> I have some experimental data where I have counts of the number of
> insects collected to different trap types rotated through 5 different
> location (variable -location), 4 different chemical attractants [A, B,
> C, D]  were applied to the traps (variable - semio) and all were
> trialled at two different CO2 release rates [1, 2] (variable CO2) I also
> have a selection of continuous variables measuring meteorological
> conditions to account for any bias cause by changing weather conditions,
> the data is over dispersed so I have fitted a negative binomial glm
> (glm.nb) and simplified using stepAIC from the MASS package etc.
>
> There are significant differences in the number of insects attracted to
> the different chemical (semio) and to the two different CO2 (release
> rates) I have then used the glht function from the multcomp package to
> do multiple comparisons to see what the specific differences between the
> levels are for semio and CO2 using the code below which works great but
> what I would like to do is to do comparisons combining the factors e.g a
> comparison for semioA at CO at level 1 vs Semio A at CO2 level 2  etc to
> see which is the best combination, is this possible or should I have
> started of with my counts already split up into this e.g. a treatment
> variable(semioA at CO2 level 1 = A1, semioA at CO2 level 2 = A2 etc), I
> started with them this way as we have no prior knowledge that increasing
> co2 will have any effect.
>
> I have had a quick try with the data split into a treatment factor
> (instead of semio and CO2 level) but I can not get convergence with
> glm.nb I think this may be to do with to many zeros in the data set, do
> you know if glht or another multiple comparison will work or
> zeroinflated negative binomial regression(zeroinfl() from the pscl
> library)?
Yes, it does, because zeroinfl() (and hurdle()) both return models with 
the usual extractor methods (coef, vcov, et al). The only nuisance is that 
you can't use the mcp() interface because it will be confused about the 
two parts of the model (count vs. zero-inflation part) which might both 
contain the same variables. What you have to do is set up your contrast 
matrix by hand. An example using the NMES1988 data from the "AER"
package
is included below.

hth,
Z

## packages
library("pscl")
library("multcomp")

## data
data("NMES1988", package = "AER")
nmes <- NMES1988[, c(1, 6:8)]

## models
fm_pois <- glm(visits ~ ., data = nmes, family = poisson)
fm_nbin <- glm.nb(visits ~ ., data = nmes)
fm_zinb <- zeroinfl(visits ~ . | ., data = nmes, dist = "negbin")

## generalized linear hypotheses
glht_pois <- glht(fm_pois, linfct = mcp(health = "Tukey"))
glht_nbin <- glht(fm_nbin, linfct = mcp(health = "Tukey"))

## compute contrasts "by hand" for ZINB
nr <- length(levels(nmes$health))
contr <- matrix(0, nrow = nr, ncol = length(coef(fm_zinb)))
colnames(contr) <- names(coef(fm_nbin))
rownames(contr) <- paste(levels(nmes$health)[c(2, 3, 3)], 
levels(nmes$health)[c(1, 1, 2)], sep = " - ")
contr[,3:4] <- contrMat(numeric(nrow(contr)), type = "Tukey")[,-2]
glht_zinb <- glht(fm_zinb, linfct = contr)

## multiple comparisons
summary(glht_pois)
summary(glht_nbin)
summary(glht_zinb)

> Any help or ideas will be gratefully appreciated. Many thanks in
> advance.
>
> Lara
>
>
> semiochemical <-read.csv("G:/semiochemical_data.csv",
header=T)
>
> semiochemical$location2<-factor(semiochemical$location)
>
levels(semiochemical$location2)<-c("1","2","3","4","5")
>
> semiochemical$semio2<-factor(semiochemical$semio)
>
levels(semiochemical$semio2)<-c("A","B","C","D")
>
> semiochemical$CO22<-factor(semiochemical$CO2)
> levels(semiochemical$CO22)<-c("1","2")
>
> model1<-glm.nb(total ~ semio + CO2 + location + temp + mean.wind.speed)
>
> model.glht.Semio <- glht(model1, linfct=mcp(semio="Tukey"))
> model.glht.CO2 <- glht(model1, linfct=mcp(CO2="Tukey"))
>
>
> 	[[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.
>
>