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.