# When to treat a number as a count

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.

You're asking a question about methods here, but I would rather start an answer from your data and what you want to know.

Here's a version of your data that may be useful to people who don't routinely use R; the opening and closing lines are specifically for Stata, but users of other software should be able to adapt according to need. The periods are Stata's generic code for numeric missings and correspond to NA in R.

clear
input byte(id group baseline outcome coverage)
1 1 28 28 28
2 1 28  0 28
3 1 28 28 28
4 1 28  0 28
5 1 28  0 28
6 1 28  . 28
7 1 28  . 28
8 1 28 16 28
9 1 28 28 28
10 1 28 10 28
11 1 12 12 28
12 1 28  0 28
13 1 28 28 28
14 1 28 12 28
15 1 28  0 28
16 1 28  0 28
17 1 28 28 28
18 1 24  8 28
19 1 28  0 28
20 1 28 28 28
21 1 28 28 28
22 1 28  0 28
23 1 28  4 28
24 1 28  . 28
25 0 28  . 28
26 0 28  0 28
27 0 20  . 28
28 0 28 28 28
29 0 28  . 28
30 0 24 20 28
31 0 24  1 28
32 0 28  3 28
33 0 28 28 28
34 0 28 26 28
35 0 28  . 28
36 0 28  0 28
37 0 28 20 28
38 0 28 16 28
39 0 24 16 28
40 0 28  0 28
41 0 28  . 28
42 0 28  3 28
43 0 28  0 28
44 0 28  1 28
45 0 16 20 28
46 0 28  0 28
end

The core of the problem is comparing mean outcome for two values of group. A distraction is that baseline varies and it seems to be simplest at least at the outset to ignore cases that are not 28 days for baseline. It is not obvious to me that adding baseline as a predictor is the best way to adjust for varying baseline; an alternative is to scale outcome to fraction of baseline, but that is likely to be confusing too.

Comparing means can naturally be re-cast as a regression problem. Here is a graph with the regression line superimposed for the regression of outcome on group for baseline 28 days. With this simplification, the line just connects the two group means.

I am rotating your histograms and treating the data as what they are, data points in a regression problem comparing means. Stacking of identical outcomes is a graphical convention only and does not affect the regression results.

Although you refer to "Gaussian regression" the ideal condition of Gaussian or normal errors is the least important aspect of linear regression. The recent text by Gelman and friends

https://www.cambridge.org/core/books/regression-and-other-stories

even advises against normal quantile plots of residuals as a waste of time. I wouldn't go that far, but it is a serious point of view.

A glance at the graph and regression results points up a difference of 2.9 days; my lay guess is that a difference of that magnitude would be clinically or scientifically interesting, but the regression results show that the sample is too small to confirm it as significant at conventional levels. If you are worried that such an indication is over-dependent on the implicit assumption of normal errors, some bootstrapping of those regression results implies that a difference of zero is well inside just about any confidence interval for the difference of means. Retreat to Kruskal-Wallis seems to me a counsel of despair; why use 1950s technology when 1970s technology (bootstrap) is available and allows you to focus on the difference of means that is of prime interest?

In general, it is a really good idea to be sensitive to whether your data are counted or measured; to think about their conditional distributions; and to note whether an outcome is necessarily bounded. In this particular case, these plain regression results imply that it hardly matters what you assume or is what is assumed or ideal for the methods used. The difference between means looks interesting but is not conventionally significant and that indication is robust to whatever you do by way of analysis.

However, if I try to match your binomial regression, but focusing on baseline equal to 28, I find similarly that it is enough to flip the difference to conventionally significant. I didn't at first understand why there is such a big difference in indication.

But we should worry about what is assumed about distributions. I note that binomials can't be U-shaped. I first doubted whether that was the issue, but such thinking was visceral, not logical. If you repeat the analysis with robust (Eicker-Huber-White) standard errors, then the significance evaporates.

In short, in applying binomial regression rather than plain regression you're replacing a distributional assumption that doesn't bite -- even though it seems quite wrong -- with a distributional assumption that does bite! That's my diagnosis.

FWIW, the use of days here as an integer count is partly natural (people have daily rhythms they follow, sometimes rigidly and sometimes loosely) and partly a convention (in principle the data might be based on times of day too, yielding fractional days).

Finally, comparison of means is not the only game in town. I note that in group 0 just 2 out of 13 but in group 1 7 out of 19 people reported the full 28 days. Those differences naturally affected the means, but the detail may be important too.

Nitty-gritty follows as Stata output. R people expect that we're smart enough to decode R output if we're perverse enough not to use it (not to use it routinely, in my case) and I return the compliment. The minimalism of R output is admirable, except that not showing the sample size used in even the default summary is puzzling to me.

. set seed 2803

. quietly bootstrap diff=_b[1.group], reps(1000) : regress outcome i.group if baseline == 28
(running regress on estimation sample)

Linear regression                               Number of obs     =         32
Replications      =      1,000

command:  regress outcome i.group
diff:  _b[1.group]

------------------------------------------------------------------------------
|   Observed   Bootstrap                         Normal-based
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
diff |   2.910931   4.409327     0.66   0.509    -5.731191    11.55305
------------------------------------------------------------------------------

. estat bootstrap, percentile  normal bc

Linear regression                               Number of obs     =         32
Replications      =       1000

command:  regress outcome i.group
diff:  _b[1.group]

------------------------------------------------------------------------------
|    Observed               Bootstrap
|       Coef.       Bias    Std. Err.  [95% Conf. Interval]
-------------+----------------------------------------------------------------
diff |   2.9109312   .1026972   4.4093271   -5.731191   11.55305   (N)
|                                      -5.055556   11.84828   (P)
|                                      -5.582857   11.58442  (BC)
------------------------------------------------------------------------------
(N)    normal confidence interval
(P)    percentile confidence interval
(BC)   bias-corrected confidence interval

. glm outcome i.group baseline, f(binomial coverage)

Iteration 0:   log likelihood = -530.29406
Iteration 1:   log likelihood = -516.62802
Iteration 2:   log likelihood = -516.61552
Iteration 3:   log likelihood = -516.61552

Generalized linear models                         Number of obs   =         38
Optimization     : ML                             Residual df     =         35
Scale parameter =          1
Deviance         =  980.8498432                   (1/df) Deviance =   28.02428
Pearson          =  748.2307194                   (1/df) Pearson  =   21.37802

Variance function: V(u) = u*(1-u/coverage)        [Binomial]
Link function    : g(u) = ln(u/(coverage-u))      [Logit]

AIC             =   27.34819
Log likelihood   =  -516.615519                   BIC             =   853.5343

------------------------------------------------------------------------------
|                 OIM
outcome |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.group |   .2522059   .1263387     2.00   0.046     .0045866    .4998252
baseline |   -.038664   .0188569    -2.05   0.040    -.0756228   -.0017053
_cons |   .5471053   .5090758     1.07   0.283    -.4506649    1.544875
------------------------------------------------------------------------------

. glm outcome i.group if baseline == 28, f(binomial coverage)

Iteration 0:   log likelihood = -485.92872
Iteration 1:   log likelihood = -481.53913
Iteration 2:   log likelihood = -481.53724
Iteration 3:   log likelihood = -481.53724

Generalized linear models                         Number of obs   =         32
Optimization     : ML                             Residual df     =         30
Scale parameter =          1
Deviance         =  931.0323116                   (1/df) Deviance =   31.03441
Pearson          =  708.6313527                   (1/df) Pearson  =   23.62105

Variance function: V(u) = u*(1-u/coverage)        [Binomial]
Link function    : g(u) = ln(u/(coverage-u))      [Logit]

AIC             =   30.22108
Log likelihood   = -481.5372359                   BIC             =   827.0602

------------------------------------------------------------------------------
|                 OIM
outcome |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.group |   .4368407   .1406668     3.11   0.002     .1611389    .7125425
_cons |  -.6481498   .1103816    -5.87   0.000    -.8644938   -.4318058
------------------------------------------------------------------------------

. predict predicted
(option mu assumed; predicted mean outcome)

. tabdisp group, c(predicted)

--------------------------------
group |            predicted
----------+---------------------
0 |             9.615385
1 |             12.52632
--------------------------------

. glm outcome i.group if baseline == 28, f(binomial coverage) robust

Iteration 0:   log pseudolikelihood = -485.92872
Iteration 1:   log pseudolikelihood = -481.53913
Iteration 2:   log pseudolikelihood = -481.53724
Iteration 3:   log pseudolikelihood = -481.53724

Generalized linear models                         Number of obs   =         32
Optimization     : ML                             Residual df     =         30
Scale parameter =          1
Deviance         =  931.0323116                   (1/df) Deviance =   31.03441
Pearson          =  708.6313527                   (1/df) Pearson  =   23.62105

Variance function: V(u) = u*(1-u/coverage)        [Binomial]
Link function    : g(u) = ln(u/(coverage-u))      [Logit]

AIC             =   30.22108
Log pseudolikelihood = -481.5372359               BIC             =   827.0602

------------------------------------------------------------------------------
|               Robust
outcome |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.group |   .4368407   .6659552     0.66   0.512    -.8684075    1.742089
_cons |  -.6481498   .5129588    -1.26   0.206    -1.653531     .357231
------------------------------------------------------------------------------
• (+1) The distributional assumption bites in the logistic regression primarily because of the aggregation of drug-day counts over all the patients in a treatment group - as if it were immaterial whether there were 17 controls & 21 treated each studied over 28 days or a single control observed over 182 days & a single treated over 258 days. More appropriate would be to have 'patient' as a random effect nested in 'group' (though the assumption that a patient's decision to take a drug on any given day is independent of their decision on other days is still dubious). Commented Aug 25, 2020 at 12:17
• By the way, I tried a generalized linear mixed-effects model with the glmer` function from R's lme4 package, but the unconditional fits were way off - I imagine the assumption of normally distributed patient effects is inappropriate, as a Gaussian distribution can't be U-shaped either. Commented Aug 25, 2020 at 14:36
• Thank you so much for your answer Nick Cox. I have been grappling with this issue for some time as distributions of frequency of drug use data are bimodal as a rule. Thank you for reminding me about bootstrapping. I guess because I have not done it a lot I am not as comfortable with it, but it actually makes perfect sense. I must confess I am not sure what you mean by 'bite'. Is 'biting' a good or bad thing, or do you just mean 'returns a singificant p-value'. So can I clarfiy you are recommending bootstrapping a Guassian regression? To me that seems very sensible... Commented Aug 26, 2020 at 1:41
• ...I completely distrusted the results of the binomial regression. Even if a mean difference of ~2 was clinically significant (which I doubt) as you say the sample size is so small it would be absolute folly to make any inferences about group differences. This is all just part of my ongoing quest to find a way to deal with this sort of strange data that doesn't fit comfortably into any one category. I imagine I will be working with it for a long time to come so I need to get to grips with it somehow. Thank you again. A great, thoughtful answer. Commented Aug 26, 2020 at 1:52
• Thanks for your comments. "bite" is bad; I mean that an assumption (I would rather say "ideal condition") for a technique you used isn't even roughly satisfied, and that matters, because your results are seriously affected by that failure. As before, I wouldn't say bootstrapping Gaussian regression; you're bootstrapping linear regression. Note that you could bootstrap the difference between means directly. Commented Aug 26, 2020 at 7:28