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.