Unexpected post-hoc result using binomial GLM

Question

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

Model: binomial, link: logit

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.

Thanks in advance!

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")

Edit for @AdamO's response.

Here I make a contingency table:

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


output: 

  alive dead
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.

Answer 1

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.