Skip to content

multcomp and glm

2 messages · Jesse.Whittington@pc.gc.ca, Torsten Hothorn

Jesse.Whittington@pc.gc.ca
# 
I have run the following logistic regression model:

options(contrasts=c("contr.treatment", "contr.poly"))
m <- glm(wolf.cross ~ null.cross + feature, family = "binomial")

where:
wolf.cross = likelihood of wolves crossing a linear feature
null.cross = proportion of random paths that crossed a linear feature
feature = CATEGORY of linear feature with 5 levels: high-use road, low-use
road, high-use trail, low-use trail, and railway line

I would like to determine whether wolves are more likely to cross some
features than others and am using csimtest in the package multcomp to do
so.

x <- coef(model)
var.cov <- vcov(model)
df <- model$df.residual
contrast.matrix <- contrMat(rep(2,length(x), type = "Tukey")
post.hoc <- csimtest(estpar = x, df=df, covm = var.cov, cmatrix =
contrast.matrix)
summary(post.hoc)

My questions are:
(1) Have I entered my parameters for the posthoc analysis correctly and
does it matter which categorical variable I use as my reference category?

(2) How do I omit the continuous variable "null.cross" from the list of
multiple comparisons?

I would be most grateful for any assistance you can give me,
Jesse Whittington


Jasper National Park
Box 10,
Jasper, Alberta
Canada
(780) 852-6187
Torsten Hothorn
# 
0) you should use asympt = TRUE for a logistic regression model, that is
using the asymptotic join multivariate normal distribution of the
estimated parameters instead of the t-distribution. This will NOT work
with `csimtest' because of a bug that will be fixed in the next release
(but `csimint' will do).
I think there are some problems here.
will tell you about ALL estimated parameters, that is: intercept,
null.cross and the 4 treatment contrasts for factor feature.
will compute the contrasts for ALL parameters and I guess you would like
to have it for the levels of `feature' only.

As I pointed out in the last email thread on "Multiple comparison and
lme (again, sorry)" (r-help, May 14th), I can't see any `high-level' way
to compute Tukey contrasts and the corresponding convariance matrix
except doing the (model dependent) calculations by `hand' because you
can't estimate all parameters due to a design matrix with reduced rank.

The basic problem is that the S `contrasts' function (and all
model.fit functions) implicitly assume that there are at most k-1
contrasts for a factor at k levels, something that does not
cover Tukey contrasts.

Anyway, Tukey contrasts would implement all-pair comparisons, that is
comparing every combination of levels of `feature'. Is this really your
question?
just omit it from x and var.cov (maybe you should omit the intercept,
too):

csimint(estpar = x[-(1:2)], covm = var.cov[-(1:2),-(1:2)], cmatrix =
        diag(4), asympt = TRUE)

will give you simultaneous CI's for the treatment contrasts of `feature'.

There will be high-level functions for those problems in the next
release with some examples (and I can send you the current state,
if you like).

Best,

Torsten