I am having a problem conducting post hoc analyses on a binomial glm using the data linked above.
Experiment: I am looking at the detectability (presence/absence) of an insect marker given 3 explanatory variables; application method (apptreat - factor, 2 levels), marker (factor, 2 levels), and exposure time (exp - numerical, 3 levels). Response variable (detect) is recorded as 1 for present and 0 for absent. After running the model several times to determine whether or not any interactions were significant, I came up with the following model.
Model:
id.glm2 <- glm(detec~apptreat+marker+exp+apptreat*exp, family=binomial, data=indiv_detec2)
anova(id.glm2, test="Chisq")
Output:
Analysis of Deviance Table
Model: binomial, link: logit
Response: detec
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 906 1182.37
apptreat 1 238.847 905 943.52 < 2.2e-16 ***
marker 1 156.844 904 786.68 < 2.2e-16 ***
exp 1 138.098 903 648.58 < 2.2e-16 ***
apptreat:exp 1 9.414 902 639.17 0.002153 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
I would like to do a post hoc multiple comparison analysis to compare the differences within each variable including the interaction using the glht function in multcomp CRAN. I am able to get comparisons for the main effects with a warning about interactions being present (as expected from other posts). In order to perform posthoc analysis on the interaction I added a column to the data for the interaction using the following code and included it as a main effect in the glm model
indiv_detec2$AE <- interaction(indiv_detec2$apptreat, indiv_detec2$exp)
id.glm2 <- glm(detec~apptreat+marker+exp+AE, family=binomial, data=indiv_detec2)
anova(id.glm2, test="Chisq")
When I use the glht() function to do the posthoc analysis on the interaction variable (AE) it results with the following error message
glht(m1, linfct = mcp(AE = "Tukey"))
Error in modelparm.default(model, ...) : dimensions of coefficients and covariance matrix don't match
The data is unbalanced but I don't see how that would prevent the posthoc on the interaction variable and not the other variables. I know this problem has been brought up on other posts but I haven't been able to get the solutions mentioned in those posts to work. I'm assuming it's something I am doing wrong/simple mistake and just can't figure it out.
UPDATE I have tested the model with the all data as factors and where "exp" (exposure time) and "detec" (detection) are continuous. Both ways result in the same error messages being generated, however, different interactions are detected as being significant depending on how the data is classified. Knowing this, I will continue to use the model above for simplicity since the problem remains the same
When attempting the glht() function using the model
id.glm3<-glm(detec~apptreat+marker+exp+apptreat*exp+marker*exp, family=binomial, data=indiv_detec2)
summary(glht(id.glm3, linfct = mcp(app = "Tukey")))
I am unable to run the test on the interaction (which is normal for glht() and mcp() functions) and the warning message below appears when running the test on the main effects "apptreat" and "exp".
Warning message: In mcp2matrix(model, linfct = linfct) : covariate interactions found -- default contrast might be inappropriate
Error in modelparm.default(model, ...) : dimensions of coefficients and covariance matrix don't match
The summary output looks like this
Call:
glm(formula = detec ~ marker + apptreat + exp + apptreat * exp,
family = binomial, data = indiv_detec2)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.0545 -0.5493 -0.1415 0.6164 3.1055
Coefficients: (1 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.8141 0.3672 -13.111 < 2e-16 ***
markerMilk 2.7846 0.2365 11.774 < 2e-16 ***
apptreatLitter 3.6469 0.4933 7.394 1.43e-13 ***
exp24 0.2145 0.4038 0.531 0.595
exp48 3.0381 0.3948 7.696 1.41e-14 ***
apptreatLitter:exp24 -0.2675 0.5697 -0.469 0.639
apptreatLitter:exp48 NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1182.37 on 906 degrees of freedom
Residual deviance: 623.15 on 901 degrees of freedom
AIC: 635.15
Number of Fisher Scoring iterations: 5
When I amend the data set to include a column representing the interaction as suggested in other posts and run the model with the interaction as a main effect
indiv_detec2$AE<-interaction(indiv_detec$apptreat, indiv_detec$exp)
id.glm3<-glm(detec~apptreat+marker+exp+AE, family=binomial, data=indiv_detec2)
I get the following error message when running the glht() function
summary(glht(id.glm3, mcp(AE="Tukey")))
Error in modelparm.default(model, ...) : dimensions of coefficients and covariance matrix don't match
The summary output looks like this
Call:
glm(formula = detec ~ marker + apptreat + exp + AE, family = binomial,
data = indiv_detec2)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.0545 -0.5493 -0.1415 0.6164 3.1055
Coefficients: (3 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.81408 0.36719 -13.111 < 2e-16 ***
markerMilk 2.78462 0.23650 11.774 < 2e-16 ***
apptreatLitter 3.64689 0.49325 7.394 1.43e-13 ***
exp24 -0.05294 0.57806 -0.092 0.927
exp48 3.03810 0.39479 7.696 1.41e-14 ***
AEDirect.24 0.26748 0.56973 0.469 0.639
AELitter.24 NA NA NA NA
AEDirect.48 NA NA NA NA
AELitter.48 NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1182.37 on 906 degrees of freedom
Residual deviance: 623.15 on 901 degrees of freedom
AIC: 635.15
Number of Fisher Scoring iterations: 5
I have read, and suspect, that whatever is causing the NA to appear in the summary output is the problem, but I have not seen an answer that clearly explains what's going on so I can fix the problem.
I hope this clarifies
