Post-hoc test after a GLM with binary data

Question

I have run various logistic regressions (GLM) from the binomial family and they have produced some very interesting results. I would now like to run post-hoc tests to find out which levels of the explanatory variables are significant, but I am finding it very difficult to find a post-hoc test that is compatible with my data.

I am looking at a population of skeletons. I have recorded the presence/absence of infectious lesions at the level of the individual. I have a response variable (Periostitis, which is a type of lesion and which has two levels: present and absent) and I have four explanatory variables (Age.Group, Sex, Cemetery, and Time.Period). These represent the age of the individual (A, B, C, and D), the sex of the individual (Male or Female), the cemetery the individual was recovered from (Cem1, Cem2, and Cem3), and the time period the individual is dated to (TP1 and TP2), respectivly.

I have run logistic regressions because my data is binomial. An example of which can be seen here:

Call:
glm(formula = Periosna$Periostitis ~ Periosna$Cemetery, family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9103   0.5931   0.6841   0.6841   0.9537  

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)             1.3332     0.2087   6.387 1.69e-10 ***
Periosna$CemeteryCem2   0.3155     0.5311   0.594   0.5525    
    Periosna$CemeteryCem3  -0.7811     0.3557  -2.196   0.0281 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 244.03  on 221  degrees of freedom
Residual deviance: 238.04  on 219  degrees of freedom
AIC: 244.04

Number of Fisher Scoring iterations: 4

I now want to run a post-hoc test to see which levels (Cem1, Cem2, or Cem3) are significant. I have read everything on this site and on others and have tried to run the glht function from the multcomp package. I see that there are different types of tests (matrix, character, expression, mcp, and mlf), but I do not know which to use. There is a lot of literature online about mcp, but I do not think that I can use this as my data is binomial. When I try to use it, I get an error code:

> library(multcomp)
> 
> results<-glht(modelA,linfct=mcp(Cemetery="Tukey"))
Error in mcp2matrix(model, linfct = linfct) : 
  Variable(s) ‘Cemetery’ have been specified in ‘linfct’ but cannot be     found in ‘model’! 

What would you suggest I use to analyze the levels of my explanatory variables?

Thank you very much for your help. I am stumped!

Answer 1

The problem you haveg is that are abusing (misusing) the formula notation. You never use a formula like this:

Periosna$Periostitis ~ Periosna$Cemetery

because, as you've found out it, it totally breaks the logic that you want to then use for subsequent operations.

You asked mcp() to set up Tukey contrasts for the covariate Cemetery, yet you fitted a model with covariate name Periosna$Cemetery, hence it quite rightly told you that a variable of the stated name was not involved in the model fit.

Instead, the call should have been

glm(Periostitis ~ Cemetery, data = Periosna, family = binomial)

noting that we use the variable names directly, as they exist in Periosna, and we tell R where to locate the variables in the formula via the data argument.

Now when you call glht() it should be able to find the Cemetery variable.