# Regression model for count data with restricted upper bound

What is the recommended model for count data with a known/restricted upper bound on the number of counts? I know that Poisson is used for count data and the negative binomial in case of overdispersion, but if I am not mistaken, they both assume a non-zero count that has no upper limit (0 to inf).

More details: the dependent variable $$y$$ is the number of events attended per customer (per year). There is a maximum of 25 events per year (which is the upper bound on thecount). Customers attend 1 to 25 events so $$y$$ is an integer belonging to [1,25].

The independent variable $$x$$ is a categorical variable representing the location of the customer and that can take three values (local, non-local and foreigner). I want to use a glm(y ~ x) to predict the number of events attended based on the value of $$x$$ (I am using R).

Also note that I have deducted the number of attended events by 1 so that my $$y$$ is now between 0 and 24 to avoid the need for a zero-truncated model (in case I need to use a Poisson or negative binomial). I then intend to use the inverse link function to get the prediction of number of events, and add one to it.

• Binomial distribution may be applicable. – Nick Cox Sep 6 '17 at 11:13
• Can you explain more about what the variable is/why there's an upper limit? – Glen_b Sep 6 '17 at 11:57
• @Glen_b: The dependent variable "y" is the number of events attended per customer (per year). There is a maximum of 25 events per year (which is the upper bound on count). Customers attend 1 to 25 events so "y" is an integer belonging to [1,25]. The independent variable "x" is a categorical variable representing the location of the customer and that can take three values (local, non-local and foreigner). I want to use a glm(y ~ x) to predict the number of events attended based on the value of x (I am using R). Also note that I have deducted the number of attended events by 1 so that my "y" is – NadTeX Sep 7 '17 at 1:52
• @NickCox applicable for use as what? – Tilefish Poele Sep 8 '17 at 6:36
• Please join your accounts if NedTex and Nader are one person. Also if you have a new question ask a new question. – Ferdi Sep 8 '17 at 8:06

It seems like a binomial GLM might work reasonably well (as Nick also mentioned).

Failing that a quasi binomial glm (to pick up the different dispersion due to a mix of $p$'s across people), or a beta-binomial model, or a binomial glmm (generalized linear mixed model).

There's questions relating to each of these on site.

## poisson models for proportions

The trick to this problem is that it is more a proportion than it is a count. You may have been told that a Poisson model is "used for count data". In fact, you could use Poisson data to model proportions (or binomial) outcomes. The only issue is that the Poisson model tends to over estimate the variance, (Binomial: var = p(1-p), Poisson: var = p) leading to wider confidence intervals and larger p-values, on average. It can also predict beyond the range of normal values. Is that a bad thing? Just truncate the over predictions, and wide confidence intervals just mean reduced power.

## independent data models for over/under-dispersed binomial outcomes

Diving into the problem more, binomial variables come from independent and identically distributed Bernoulli events. It's important to verify those assumptions. One way is to inspect the dispersion of Bernoulli outcomes. For instance, let's simulate 3 meetings and their attendance, but 1st not so popular, and 3rd is very popular.

set.seed(123)
p <- c(0.3, 0.5, 0.7)
x <- matrix(rbinom(100*3, 1, p), 100, 3, byrow = T)
y <- rowSums(x)

summary(glm(cbind(y, 3) ~ 1, family=quasibinomial))


gives

Call:
glm(formula = cbind(y, 3) ~ 1, family = quasibinomial)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.5555  -0.3563   0.3188   0.3188   0.8486

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.69984    0.05662  -12.36   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.3191032)

Null deviance: 41.543  on 99  degrees of freedom
Residual deviance: 41.543  on 99  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4


A dispersion parameter of 0.32 which is highly under-dispersed. Another way to check in this univariate case is simply:

> var(y/3)
[1] 0.07958474
> mean(y/3)*(1-mean(y/3))
[1] 0.2499


Note: 0.08 (actual variance) / 0.25 (expected variance under binomial probability model) = 0.32. I consider the glm approach superior because it allows for adjusting for covariates.

One way of conceptualizing the problem of underdispersion (in this case) is that if you simulated outcomes from a binomial GLM, their distributions would not match the target distributions. The calibration of the probability model for the response is:

> # expected
> dbinom(x = 0:3, size = 3, prob = mean(y/3))
[1] 0.1275167 0.3774832 0.3724834 0.1225166

> # observed
> prop.table(table(y))
y
0    1    2    3
0.11 0.37 0.44 0.08


(see figure at bottom of the page for a calibration plot)

Here, using a quasibinomial GLM will naturally widen confidence intervals to compensate for the miscalibration noted above. That way we can get correct p-values and confidence intervals. Predicting and simulating outcomes remains a problem, however.

## proportional odds models for over/under-dispersed binomial outcomes.

Another approach is to use a proportional odds model. These "cumulative link" models allow for adjacent counts to have a unique latent log-odds difference determining the likelihood of achieving at least one higher outcome. This is why this type of model is preferred for ordinal data.

A "novariate" (intercept only) model for the response gives the following log-odds parameters for endorsing an at-least-one-higher response

polr(factor(y) ~ 1) Call: polr(formula = factor(y) ~ 1)

No coefficients

Intercepts:
0|1         1|2         2|3
-2.09049882 -0.07990206  2.44238850

Residual Deviance: 234.7927
AIC: 240.7927


Importantly, this leads to a probability model with a much higher agreement with the empirical frequencies (calibration)

> predict(polr(factor(y) ~ 1), newdata = data.frame('1'=1), type = 'probs')
0          1          2          3
0.11002372 0.37001138 0.43996795 0.07999695


So I find that this method is superior for handling the case when the Bernoulli outcomes contributing the binomial event are possibly non-independent and possible non-identical in their probability.

This seems like an ideal case for an Ordinal Regression. The number of events attended is an ordered categorical variable with a defined range, so for this reason I would shy away from a binomial GLM. I would suggest starting with the Ordered Logit model.

While I'm not very familiar with R here are two resources which should be helpful for you:

• I don't follow the reasoning behind the implication that "The number of events attended is an ordered categorical variable with a defined range" means a Binomial GLM is contraindicated, because (under some additional assumptions explored in comments to the question) this is a principal part of the characterization of a binomial response. – whuber Dec 21 '18 at 16:57
• See my answer for a rationale of using proportion odds models (so called ordinal regression) in the case of non-independent or non-identical Bernoulli trials. – AdamO Dec 21 '18 at 19:49

You essentially just have categorical dependent data, just with 25 categories, it doesn't really matter that it's count data or that it has an upper limit in this case. Thus a multinomial logistic regression should work.

Theoretically Poisson regression would be technically more interpretable, as essentially the errors are from a count distribution, however the implication of the limit at 25 could give you erroneous predictions, whereas you can't in the multinomial case.