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