# Unexpected post-hoc result using binomial GLM

I have a dataset that has a categorical factor and numerical response variables as proportions. A simplified deput() included at the end of the body. But, here is what the data looks like:

> head(df)
ID treatment day alive prop
1  1         A   4     1    1
2  2         A   4     1    1
3  3         A   4     1    1
4  4         A   4     1    1
5  5         A   4     1    1
6  6         A   4     1    1


All together, there are 7 treatments (A:G) including a negative control (treatment A). Each replicate has an n=1 and prop is defined as the proportion of alive individuals / n. Therefore in prop, 1=alive, 0=dead.

I aim to look for differences in mean survival among treatments at day 4. Because my response variable is limited to values between 0-1, I believe I need to build a binomial glm. Here's what I have done:

df$$treatment <- as.factor(df$$treatment)
m <- glm(prop ~ treatment, data=df, family=binomial())
summary(m)

output:
Call:
glm(formula = prop ~ treatment, family = binomial(), data = df)

Deviance Residuals:
Min        1Q    Median        3Q       Max
-2.52113   0.00013   0.29175   0.51678   0.96954

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.135e+00  1.022e+00   3.069  0.00214 **
treatmentB  -2.625e+00  1.105e+00  -2.375  0.01755 *
treatmentC   1.543e+01  1.331e+03   0.012  0.99075
treatmentD  -3.644e-15  1.445e+00   0.000  1.00000
treatmentE  -1.190e+00  1.193e+00  -0.997  0.31890
treatmentF  -1.526e+00  1.159e+00  -1.317  0.18797
treatmentG  -7.376e-01  1.261e+00  -0.585  0.55845
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 122.65  on 167  degrees of freedom
Residual deviance: 101.86  on 161  degrees of freedom
AIC: 115.86

Number of Fisher Scoring iterations: 17


Side questions: #1: Not all of my coefficients are statistically significant. Is this an indication of how well the model fits the data for each coefficient?

#2: The deviance residuals at each quartile are not similar. Does this indicate that the model does not fit my data well? Is there a normality test equivalent for a binomial distribution that my data must pass before modeling?

To (attempt) to answer by question: 'Is there significant difference in mean survival among treatments?' I use anova using test='Chisq'.

a <- anova(m, test='Chisq')

output:

> a
Analysis of Deviance Table

Response: prop

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL                        167     122.65
treatment  6   20.785       161     101.86 0.002005 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Based upon the p-value of the anova, I see that there is significant difference in mean survival among treatments.

My next question is to use a post-hoc by looking at pair-wise comparisons mean survival between treatments to detect separation of statistically similar groups.

g <- summary(glht(m, mcp(treatment="Tukey")))

output:
> g

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glm(formula = prop ~ treatment, family = binomial(), data = df)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
B - A == 0 -2.625e+00  1.105e+00  -2.375    0.170
C - A == 0  1.543e+01  1.331e+03   0.012    1.000
D - A == 0 -3.644e-15  1.445e+00   0.000    1.000
E - A == 0 -1.190e+00  1.193e+00  -0.997    0.941
F - A == 0 -1.526e+00  1.159e+00  -1.317    0.809
G - A == 0 -7.376e-01  1.261e+00  -0.585    0.996
C - B == 0  1.806e+01  1.331e+03   0.014    1.000
D - B == 0  2.625e+00  1.105e+00   2.375    0.170
E - B == 0  1.435e+00  7.475e-01   1.920    0.409
F - B == 0  1.099e+00  6.912e-01   1.589    0.636
G - B == 0  1.887e+00  8.504e-01   2.219    0.238
D - C == 0 -1.543e+01  1.331e+03  -0.012    1.000
E - C == 0 -1.662e+01  1.331e+03  -0.012    1.000
F - C == 0 -1.696e+01  1.331e+03  -0.013    1.000
G - C == 0 -1.617e+01  1.331e+03  -0.012    1.000
E - D == 0 -1.190e+00  1.193e+00  -0.997    0.941
F - D == 0 -1.526e+00  1.159e+00  -1.317    0.809
G - D == 0 -7.376e-01  1.261e+00  -0.585    0.996
F - E == 0 -3.365e-01  8.252e-01  -0.408    1.000
G - E == 0  4.520e-01  9.625e-01   0.470    0.999
G - F == 0  7.885e-01  9.195e-01   0.857    0.972
(Adjusted p values reported -- single-step method)


And none of the pairwise comparisons have significant differences in mean survival. Which... surprised me.

So, I made a bar graph to visualize mean survival of each group.

t_means <- aggregate(df$$prop, by=list(df$$treatment), mean)
t_sd <- aggregate(df$$prop, by=list(df$$treatment), sd)
t_n <- aggregate(df$$prop, by=list(df$$treatment), length)
se <- t_sd$$x / sqrt(t_n$$x)
df_means <- cbind(t_means, se)
colnames(df_means) <- c('treatment', 'mean.survival', 'se')

ggplot(data=df_means, aes(x=treatment, y=mean.survival)) +
geom_bar(stat='identity') +
geom_errorbar(ymin=df_means$$mean.survival - df_means$$se, ymax=df_means$$mean.survival+df_means$$se)


And if I were to only look at this graph, it seems like at least B and C should separate, as in, have statistically different mean survival.

For lack of better words: What am I doing wrong?

My suspicions: Perhaps, the reason they are not separating is because all treatments have values 0 and 1, and there isn't 100% mortality and 100% survival.

Here is deput()

> dput(df)
structure(list(ID = 1:168, treatment = structure(c(1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L), levels = c("A", "B", "C", "D", "E", "F",
"G"), class = "factor"), day = c(4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L), alive = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L,
1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L,
1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L,
0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L), prop = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L,
1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L,
1L, 1L, 1L, 1L, 1L, 0L, 1L)), row.names = c(NA, -168L), class = "data.frame")


Here I make a contingency table:

alive<- aggregate(df$$alive, by=list(df$$treatment), sum)
dead <- 24 - alive$$x data <- cbind(alive$$x, dead)
row.names(data) <- c('A', 'B', 'C', 'D', 'E', 'F', 'G')
data

output:

A    23    1
B    15    9
C    24    0
D    23    1
E    21    3
F    20    4
G    22    2


From here I use Fisher's Exact Test for a significant relationship between treatment and outcome (alive/dead).

test <- fisher.test(data)
test

output:
> test

Fisher's Exact Test for Count Data

data:  data
p-value = 0.002345
alternative hypothesis: two.sided



From this I gather that there is a relationship between treatment and outcome...at day 4.

Nobody died in C. Great! It's easier to see what's going on with a cross tabs like here:

> table(dd$$alive, dd$$treatment)

A  B  C  D  E  F  G
0  1  9  0  1  3  4  2
1 23 15 24 23 21 20 22


The odds ratio relating death to treatment has as the comparison of C to A $$24 * 1 / (23 * 0) = \infty$$. The formula for an OR is one every statistician should know off hand, in addition the variance of the OR = AD/(BC) is given by:

$$\text{var}(\log \widehat{OR}) = 1/A + 1/B + 1/C + 1/D$$

so you tell me what happens when one of $$A, B, C, D$$ is 0!

R's fitter tells you this in weird ways. The actual OR estimate is $$\exp(15.4) \approx 4,876,801$$. In other words, this is some artificial value out in la-la-land when R decided to kill the GLM machine from spinning off into space. Unfortunately, just as the OR exploded, so did the standard error estimates:

> vcov(m)
(Intercept) treatmentB    treatmentC treatmentD treatmentE treatmentF treatmentG
(Intercept)    1.043478  -1.043478 -1.043478e+00  -1.043478  -1.043478  -1.043478  -1.043478
treatmentB    -1.043478   1.221256  1.043478e+00   1.043478   1.043478   1.043478   1.043478
treatmentC    -1.043478   1.043478  1.77e+06(wtf)   1.043478   1.043478   1.043478   1.043478
treatmentD    -1.043478   1.043478  1.043478e+00   2.086957   1.043478   1.043478   1.043478
treatmentE    -1.043478   1.043478  1.043478e+00   1.043478   1.424431   1.043478   1.043478
treatmentF    -1.043478   1.043478  1.043478e+00   1.043478   1.043478   1.343478   1.043478
treatmentG    -1.043478   1.043478  1.043478e+00   1.043478   1.043478   1.043478   1.588933


You can either collect more data, or use some methods for sparse categorical data, like Agresti Coulli correction, median unbiased estimation, or Fisher's Exact Test.

• Hi @Adam, thank you so much for your response. I have a few follow up questions. Please correct me where my logic is wrong. In my original post, I first built a glm which allowed me to detect significant differences in mean proportion alive at day 4 by comparing them to the reference (in this case, A, the control). After that, I wanted to test for significant differences in mean survival among treatments at day 4 among all treatments with an ANOVA. And finally, I wanted to detect separation among treatments using pairwise comparisons. I discovered an unexpected result when I failed to... Nov 11, 2022 at 15:27
• ... detect any separation among groups. Your insight suggested that because my data contains sparse categorical data (0 inflated data), that drove the OR into overdrive. I included what I tried next at the end of my question. Is the binomial glm still an appropriate model to use to compare to a reference group, given sparse categorical data? To me, it seems Fisher's exact test seems to answer the questions that ANOVA would. And lastly, is there a tool that I could use to handle pairwise comparisons of sparse categorical data? Sorry for the long response, I am doing my best to learn more stats. Nov 11, 2022 at 15:35
• @scott.pilgrim.vs.r your understanding is wrong on a few counts. Having lots of zeroes is not "zero inflation". ANOVAs do not factor in; whereas an ANOVA is just a linear regression, your binomial GLM suffices to estimate proportion differences by-group. The output from this GHLT is completely redundant; look closely that it's identical to your summary.glm output. This is what achieves precisely the inference you're after. FET is widely known and easily fit and understood in R and should handle the 0 failure cases easily. Nov 11, 2022 at 15:48
• thank you! That helps clarify a lot. Last question, if I wanted to ask "what treatments are statistically similar to say, treatment B (which is actually my positive control), how would I do that with this data. This was my reasoning for trying to use GHLT for multiple pairwise comparisons on top of the glm (should of included that, my apologies). Thanks! Nov 11, 2022 at 17:22