[R] multcomp and glm

[Torsten Hothorn]

> 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)
>

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).

> 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?
>

I think there are some problems here.

> x <- coef(model)
> var.cov <- vcov(model)

will tell you about ALL estimated parameters, that is: intercept,
null.cross and the 4 treatment contrasts for factor feature.

> contrast.matrix <- contrMat(rep(2,length(x), type = "Tukey")

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?

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

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

> 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
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>
>