Test for effects of categorical variables on a binary response variable considering their interactions?

Question

i am looking for a test similar to a 2-way ANOVA that would work on a binary response variable. My response variable is presence/absence of plant species. My explanatory variables are Treatment (3 treatment groups) and habitat types(2 sites). First, I would like to know whether Treatment, Site and their interaction have a significant effect on P/A. Second, if either Treatment or Site is significant, I would like to test all pairs of treatment groups or sites to know which pairs of levels are significantly different, as I would normally do with an ANOVA.

1. Can I use the logistic regression to construct a model and then the anova function to extract p-values as mentioned by KTwildginger (Test for effects of two categorical variables on a binary response variable?)?

mymodel <- glm(Survival ~ Treatment*Site, data=survivaldata, family="binomial") anova(mymodel, test="Chisq")

2. Can I use lsmeans to compute the pairwise comparisons from the contrasts?

Answer 1

You seem to be on the right track with your thinking. If you call your response variable Presence and define it such that:

Presence = 1 if species is present and 0 if it is not 

then your binary logistic regression model can be specified as:

mymodel <- glm(Presence ~ Treatment*Site, 
               family=binomial(link="logit"), 
               data = mydata) 

This model actually models the log odds of species presence as a function of treatment, site and their interaction. As a consequence, the model will report coefficients on the log odds scale:

summary(mymodel)

coefficients(mymodel) 

To report coefficients on a more interpretable scale, you can simply exponentiate the reported coefficients:

exp(coefficients(mymodel))

In terms of testing the significance of the main effects and interaction, you have several options, as explained at https://rcompanion.org/rcompanion/d_04.html:

### --------------------------------------------------------------
### Example of Type I, II, III tests
### --------------------------------------------------------------

### needed for type III tests
options(contrasts = c("contr.sum", "contr.poly"))

anova(mymodel, test="Chisq")              # Type I tests

install.packages("car")
library(car)

Anova(mymodel, type="II")                 # Type II tests

Anova(mymodel, type="III")                # Type III tests

See this thread for insights on when to use each type of test: Choice between Type-I, Type-II, or Type-III ANOVA.

The lrm() function in the rms package can also be used to fit the model and it reports several quantities which will help you understand the performance of your model (see the item stats in the help file for this function: https://www.rdocumentation.org/packages/rms/versions/5.1-2/topics/lrm):

install.packages("rms")

library(rms)

mymodel <- lrm(Presence ~ Treatment*Site, 
               x = TRUE, y = TRUE,
               data = mydata)

mymodel

Answer 2

Thanks for your quick answer. I think I should use anova (mymodel, type="III") instead of anova (mymodel, test="Chisq"), because we have an unbalanced model (sample sizes are different). I run these codes to test for effects of categorical variables on a binary response variable,

mymodel <- glm(Stellaria.media ~ Treatment *site, family=binomial(link="logit"), data = my data)

library(car)

Anova(mymodel, type="III")

to see whether both factors and interactions are supported:

library(MuMIn)

options(na.action = "na.fail")

BESRGN1 <- dredge(mymodel)

to compute the pairwise comparisons from the contrasts and to see the significance of the effects of the supported parameters:

library(lsmeans)

rg = regrid(ref.grid(mymodel), transform = TRUE)

GN1LSPAIR<-lsmeans(rg, pairwise~Treatment *site, type="response")

I will be really grateful if you let me know your opinion. You think is it true?