# Binomial GLM for attracting large mammal from presence/absence data

This is my first post here so I apologise if it is formatted incorrectly or missing key info.

I need a little bit of help/direction with interpreting R results for a project. Basic premise of project is testing the success of different attractants (independent variable) in attracting bears (presence/absence, dependent variable). Other independent variable in the model is trap presence (Y/N).

I have created a binomial GLM and have found no significant p-values (which is fine, of course - you can probably tell the data wasn't fab from the massive confidence intervals and it's also relatively small n=36) but would appreciate a sanity check. Ideally the results would have indicated a statistically sig. difference in bear presence, but I am happy to accept the null hypothesis that bear presence is not affected by trap presence or attractant type in this experiment.

#Extract deviance residuals and calculated theta value
devresid <- resid(m1, type = "deviance")
hist(devresid)
theta <- 31.858/24

#create new model removing interaction effect of trap/attractant, checked by AIC
m2<- update(m1,~.-Attractant:Trap)
AIC(m1,m2)

#AIC indicates Model 2 is better so find overall effect of interaction, insignificant p-value
anova(m2,m1)

#find the individual effects
m3a<- update(m2,~.-Attractant)
m3b<- update(m2,~.-Trap)

#Find the best model for explaining, m3a has lowest AIC but only 2d.f and m3b did not reduce enough so continued on with model 2
AIC(m2,m3a,m3b)

#find significance of each effect, neither statistically significant
anova(m3a, m2, test = "Chi")
anova(m3b, m2, test = "Chi" )

#Visualise model 2 as it retains the two main effects
eff <- Effect("Attractant", m2)
plot(eff, ci.style = "bars", xlab = "Attractant", ylab = "Predicted Bear Presence")


The big questions I have are essentially -

1. Am I on the right track with choice of statistical test?
2. Is there a better/more concise way of doing this?
3. In terms of interpreting results - I have chosen to show the effect of attractant on bear presence as this is what is of most interest. Ignoring the massive confidence intervals for now, could you interpret results as 'predicted bear presence of 0.5 per week with a bait choice of wine?

I know this is a very low-level question but I would really appreciate any help, even if it is to tell me I'm entirely wrong and have to start from scratch!

Many thanks

You shouldn't do model selection (ie. use AIC or ANOVA F-test) to decide whether to keep or drop the interaction term.

From a science point of view: You started the investigation with the hypothesis that there may be an interaction between attractant type and trap presence. The result you obtained is that the interaction is insignificant; you should report it in the analysis summary explicitly, not "by omission".

From a statistics point of view: You know that the interaction term is insignificant because you performed a hypothesis test. If you remove the term and then continue with the smaller model as if you never evaluated the bigger model, the statistical properties — p-values, confidence intervals — are no longer correct. You don't account properly for all the tests you've run.

Use the full model, m1; it represents your entire study, with correct p-values. Unfortunately, you didn't learn much about bears but you acknowledge that: the sample is small and the signal to noise ratio is low.

And you can visualize the (lack of) interaction by plotting the probability of bear presence for all combinations of attractant type and trap presence. This plot is more revealing than omitting the interaction. Note: The figure is generated with fake bear data.

Here is R code to make an interaction effect plot with emmeans.

library("emmeans")

attractant_type <- c(
"Blank", "Coffee", "Creosote", "Honey", "Vanilla", "Wine"
)

n <- 120

set.seed(1234)

data <- data.frame(
attractant = rep(attractant_type, n / length(attractant_type)),
trap = sample(c("yes", "no"), n, replace = TRUE),
y = rbinom(n, size = 1, prob = 0.5)
)

model <- glm(y ~ attractant * trap, data = data, family = binomial)

emmip(model, trap ~ attractant, type = "response", CIs = TRUE)