I am testing group differences in number of days participants used a drug in the previous 28 days. These are the data, but I am having trouble deciding on which approach to use: standard Gaussian regression or aggregated binomial regression. I have asked similar questions before on CV (e.g. here) but am still a bit unsure.
I have provided R code for replicability, but of course anyone who wants to weigh in - R user or otherwise - is more than welcome.
df <- data.frame(group = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
baseline = as.integer(c(28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 12, 28, 28, 28, 28, 28, 28, 24, 28, 28, 28, 28, 28, 28, 28, 28, 20, 28, 28, 24, 24, 28, 28, 28, 28, 28, 28, 28, 24, 28, 28, 28, 28, 28, 16, 28)),
outcome = as.integer(c(28, 0, 28, 0, 0, NA, NA, 16, 28, 10, 12, 0, 28, 12, 0, 0, 28, 8, 0, 28, 28, 0, 4, NA, NA, 0, NA, 28, NA, 20, 1, 3, 28, 26, NA, 0, 20, 16, 16, 0, NA, 3, 0, 1, 20, 0)),
coverage = 28)
group is the treatment participants received; baseline the number of days they used in the 28 days prior to commencing the study; outcome the number of days they used during the 28 day study (coverage is the number of days in the trial).
Here are the summary statistics:
library(tidyverse)
df %>%
group_by(group) %>%
drop_na(outcome) %>%
summarise(mean = mean(outcome, na.rm = T),
sd = sd(outcome, na.rm = T),
median = median(outcome, na.rm = T),
firstQuartile = quantile(outcome, probs = 0.25, na.rm = T),
thirdQuartile = quantile(outcome, probs = 0.75, na.rm = T),
tally = n())
# output
# group mean sd median firstQuartile thirdQuartile tally
# <dbl> <dbl> <dbl> <int> <dbl> <dbl> <int>
# 0 10.7 11.3 3 0 20 17
# 1 12.3 12.3 10 0 28 21
And the distribution of the outcomes in each group
ggplot(df, aes(x = outcome, group = group)) + geom_histogram() + facet_wrap(~group) + scale_x_continuous(breaks = seq(0,28,7))
As is typical for substance use data the outcomes are distributed quite bimodally.
When I analyse the outcome, regressing days used, treated as a continuous variable, on treatment group and baseline days used...
summary(contMod <- lm(formula = outcome ~ group + baseline,
data = df,
na.action = na.exclude))
# output
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 17.7783 16.0047 1.111 0.274
# group 1.7020 3.9248 0.434 0.667
# baseline -0.2660 0.5919 -0.449 0.656
The model indicates no significant difference between groups in mean days used when controlling for baseline days used. However, the model residuals are very non-normal:
hist(resid(contMod))
The quantile-quantile plot tells the same story
plot(contMod,2)
So to me it looks like the standard Gaussian regression may not be appropriate to model these data.
Given that the data are basically integer counts of occurrences of a binary event (used on day x vs did not use on day x) within a known number of 'trials' (28 days). I wondered if an aggregated binomial regression might be a more appropriate way to model the data?
summary(contMod <- glm(formula = cbind(outcome, coverage-outcome) ~ group + baseline,
data = df,
family = binomial,
na.action = na.exclude))
# output
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 0.54711 0.50908 1.075 0.2825
# group 0.25221 0.12634 1.996 0.0459 *
# baseline -0.03866 0.01886 -2.050 0.0403 *
Now the group difference is significant when controlling for baseline.
Such a dramatic difference in results from two different models of the same is quite surprising to me. Of course I was aware it was possible but had never encountered it in the wild before.
So I have several questions for the clever CV users
1. Is aggregated binomial regression a better way to model these data given the extreme non-normality of both the outcome and the model residuals? And if it is appropriate did I do it correctly? And, even if I did do it correctly is there another even better way (nonparametric for example? Kruskal-Wallis test kruskal.test(outcome ~ group, data = df) yielded similar results to the Gaussian, $\chi^2 = 0.07, p = 0.80$, but doesn't control for baseline).
2. How do I interpret the output from the aggregated logistic regression? If the outcome was a Bernoulli process I would use simple binary logistic regression and interpreting the results would be straightforward, exponentiate the group coefficient and that represents how much greater the odds are of using on the single day in question in the 1 group than the 0 group. But I am less sure of how one would report the outcome from the aggregated binomial.
Help much appreciated as always.
