# Comparing different treatments using logistic regression

I have some data on testing different pesticides on insects. I have applied different concentrations and different pesticides (total 11 treatments) for 24 hours, and have 6 different data points. I'm interested in comparing the mortality and function motor capacity between each treatment. My output is binary (dead or alive, functional or non functional).

One of the 11 treatments is just exposure to H2O (control). My problem arises when I try to compare treatments where all the insects are alive or all are dead, therefore having a variable with one level. For instance, when I try to compare H2O (all alive) with any treatment.

As an example, my H2O variable will be, for a particular time point:

alive_DI_hour4_H2O = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1


I've tried to use Chi-squared and Fischer tests and it did not work since each variable needs to have at least 2 levels. I'm now trying to use logistic regression, since I read it's a better approach and I could do a post hoc analysis with pairwise comparisons between treatments. Again, I face the same problem when I try to compare H2O with any treatment and it only gives me p-values close to 1, when in some cases it's clearly not well estimated.

For instance if I try to look at differences in function motor capacity at hour 4

treatment_hour4=treatment[time==4] ## Specific for hour 4
a=glm(FMS[time==4] ~treatment_hour4, family="binomial")


H2O is the first "treatment", so when I do summary (a) it compares H2O with all the other tests using Wald test:

Call:
glm(formula = FMS[time == 4] ~ treatment_hour4, family = "binomial")

Deviance Residuals:
Min        1Q    Median        3Q       Max
-1.36144  -1.22740   0.00013   1.05892   1.29325

Coefficients:
Estimate Std.   Error z   value   Pr(>|z|)
(Intercept)               1.857e+01  4.210e+02   0.044    0.965
treatment2               -7.022e-10  5.954e+02   0.000    1.000
treatment3               -3.783e-09  5.954e+02   0.000    1.000
treatment4               -1.825e+01  4.210e+02  -0.043    0.965
treatment5               -1.828e+01  4.210e+02  -0.043    0.965
treatment6               -1.814e+01  4.210e+02  -0.043    0.966
treatment7               -1.840e+01  4.210e+02  -0.044    0.965
treatment8               -1.843e+01  4.210e+02  -0.044    0.965
treatment9               -1.883e+01  4.210e+02  -0.045    0.964
treatment10              -1.845e+01  4.210e+02  -0.044    0.965
treatment11              -1.842e+01  4.210e+02  -0.044    0.965

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 3362.2  on 2639  degrees of freedom
Residual deviance: 2631.1  on 2629  degrees of freedom
AIC: 2653.1

Number of Fisher Scoring iterations: 17


Most p-values are close to 1, when in some of the treatments more than half of the insects died. Also when I try to look at multiple comparisons by doing:

library(multcomp)
summary(glht(a, mcp(treatment_hour4="Tukey")))


I get this warning:

Warning messages:
1: In RET$pfunction("adjusted", ...) : lower == upper 2: In RET$pfunction("adjusted", ...) : Completion with error > abseps


And again, p-values that do not make sense.

I have done other type of analysis, by calculating the percentages of alive or funcional insects for each time point and doing ANOVA and TukeyHSD after arcsine square root transformation. Here I can see differences between treatments when having into account all time points but I would like to analyse separately for each time point as well.

I've also looked at each comparison individually, by defining a label and using anova this way:

label=c(rep(0,40),rep(1,40)) # 40 insects in each group
alive_DI_dose1_hour2_H2O_treatment1=c(alive_DI_dose1_hour2_H2O,alive_DI_dose1_hour2__treatment1)
anova(glm(alive_DI_dose1_hour2_H2O_treatment1~label,family="binomial"),test="LRT")


But I'm not sure if this is a good way to do it, but the p-values make more sense. Could this be right if at the end I correct the p-values for multiple comparisons?

Basically, what I would like to know is why Wald test is not appropriate and fails to estimate the p-values for this kind of data, and if I can use anova to test the model significance and how to interpret that outcome as treatment mortality.

The last two methods assume that nothing truly results in an event exactly 0 or 100% of the time, while the first allows confidence intervals for the odds ratio to include 0 or $\infty$.