# What post-hoc test should be used for a glmer model with a binary response, and a continuous and categorical predictor?

I'm a bit of a newbie with stats and R, so need a bit of direction to find a suitable post-hoc test for my glmer model. I'm trying to find if presence is affected by environmental factors for each species. 24 surveys were completed at each site.

The model has a binary dependent variable (absent/present) and the predictor variables are interactive terms between a multiple continuous variables(eg temp and pH) and a categorical variable (species, n=3). The random effect is the survey site ID data was collected from. Only interactive terms, rather than the continuous factor in isolation, produce significant results when an anova is run on the model. Species by itself has a large effect because one species is much rarer than the others.

I'm trying to tease apart how the presence of these species varies across pH and between species. I've tried lsmeans test with Tukey, and Firth's Bias-Reduced Logistic Regression, emmeans based on some other posts I read where people had similar questions. I ran the effects function on the interactive terms, so had a rough expectation of what a post hoc could show, but the results logistf (firth's) have produced I was not expecting. Emmeans and tukey both gave the same results and ignored the continuous variable I assume because it's not a factor.

When I run firth's regression it produces chi-squared and p values that are either infinity for chi values, some with infinite degrees of freedom, or the p values astronomically small, even though what I saw through effects suggested no significant difference. I can't tell with the interactive term if there truly is an effect of the environmental variable or if the significant effect is because of the difference in species.

If I wasn't clear enough about something please let me know and if anyone has any suggestions or advice, they would be massively appreciated. Thanks!

Glmer code, anova output and the logistf and effects output for pH are below. Based on effects, I was not expecting significant difference but there was one in logistif.

 ###glmer model
> Large<-glmer(Abs.Pres~ Species:Q.Depth+Species:Conductivity+Species:Temp+Species:pH+Species:DO.P+(1|QID),
+              nAGQ=0,
+              family=binomial,
+              data=Stacked_Pref)
####Anova output
> anova(Large)
Analysis of Variance Table
npar  Sum Sq Mean Sq F value
Species:Q.Depth         3 234.904  78.301 78.3014
Species:Conductivity    3  32.991  10.997 10.9970
Species:Temp            3  39.001  13.000 13.0004
Species:pH              3  25.369   8.456  8.4562
Species:DO.P            3  34.930  11.643 11.6434

####logistf run on pH
> Lp<-logistf(Abs.Pres~Species:pH, data=Stacked_Pref, contrasts.arg=list(pH="contr.treatment", Species="contr.sum"))
> Lp
logistf(formula = Abs.Pres ~ Species:pH, data = Stacked_Pref,
contrasts.arg = list(pH = "contr.treatment", Species = "contr.sum"))
Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood

coef   se(coef) lower 0.95 upper 0.95    Chisq            p
(Intercept)         1.9711411 0.57309880  0.8552342  3.1015114 12.09107 5.066380e-04
SpeciesGoby:pH     -0.3393185 0.07146049 -0.4804047 -0.2003108 23.31954 1.371993e-06
SpeciesMosquito:pH -0.3001385 0.07127771 -0.4408186 -0.1614419 18.24981 1.937453e-05
SpeciesRFBE:pH     -0.4771393 0.07232469 -0.6200179 -0.3365343 45.73750 1.352096e-11

Likelihood ratio test=267.0212 on 3 df, p=0, n=3945

###effect function output on pH
> SpE<-effect("Species:pH", Large)
> summary(SpE)

Species*pH effect
pH
Species             7        7.7        8.5       9.3         10
Goby     0.22239538 0.23898961 0.25896972 0.2800056 0.29924424
Mosquito 0.36689425 0.34004541 0.31057990 0.2825744 0.25936811
RFBE     0.09393222 0.09413637 0.09437017 0.0946045 0.09480996

Lower 95 Percent Confidence Limits
pH
Species             7        7.7       8.5        9.3         10
Goby     0.13722030 0.16103685 0.1753282 0.17341519 0.16408392
Mosquito 0.24476920 0.23994376 0.2148559 0.17474573 0.13820850
RFBE     0.05387189 0.05905686 0.0588516 0.05251263 0.04504883

Upper 95 Percent Confidence Limits
pH
Species            7      7.7       8.5       9.3        10
Goby     0.3396283 0.339411 0.3648593 0.4189090 0.4815962
Mosquito 0.5088941 0.456809 0.4258216 0.4228475 0.4333341
RFBE     0.1587829 0.146802 0.1479552 0.1645751 0.1886773
$$$$

• Can you add the R formula for your glmer model so that we can tell what fixed and random effects it includes? Also, can you explain where pH comes in - you initially mentioned temperature (temp) as your continuous predictor, but then phrased your research questions in terms of pH. Some details about your study design would also be helpful (e.g., what was measured and how often?). Finally, how does logistf come into play here - as far as I know, it can't accommodate random effects? – Isabella Ghement May 16 '20 at 15:24
• Maybe I don't understand effect() correctly, but what I think I see in your outputs are estimates of the predicted P(Abs.Pres), whereas the tests in the model summary have to do with comparisons of those. Note that for Goby, those estimates seem to increase with pH, whereas they decrease for Mosquito and kind of stay the same for RFBE – Russ Lenth May 19 '20 at 16:40

My suggestion is to do something like this:

library(emmeans)

emt <- emtrends(Large, "Species", var = "pH")
emt          # estimated slopes for each species
pairs(emt)   # pairwise comparisons of slopes


... and similar with var = "Q.Depth", var = "Conductivity", etc. Note the slopes are on the logit scale, i.e., (change in logit(p)) / (change in pH).

You can visualize these trends using

emmip(Large, Species ~ pH, cov.reduce = range)


or

emmip(Large, Species ~ pH, type = "response",
at = list(pH = c(... several pH values ...))


The latter would show the trends on the probability scale (they will be curves), while the former shows the trends on the logit scale (which are straight lines).

## Caution

However, I seriously question the appropriateness of the model Large. It includes only the terms Species:Q.Depth, Species:Conductivity, etc., without main effects. This model forces all the fitted lines to go through the origin, which in this context means that your estimated probability is 0.5 when Q.Dept, Conductivity, etc. are all equal to zero. This seems highly unrealistic. It is almost always a mistake to leave the intercept and main effects out of a model. I suggest you fit the model with this formula:

glmer(Abs.Pres ~ Species * (Q.Depth + Conductivity + Temp + pH +
DO.P) + (1 | QID), ...)
`

(I am also a big fan of whitespace in one's code...). I have every reason to believe that your estimates and tests will be radically different.