Using offset in binomial model to account for increased numbers of patients

Question

Two related questions from me. I have a data frame which contains numbers of patients in one column (range 10 - 17 patients) and 0s and 1s showing whether an incident happened that day. I'm using a binomial model to regress probability of incident on number of patients. However, I would like to adjust for the fact that when there are more patients, there will inevitably be more incidents because the total amount of patient time on the ward is higher on that day.

So I'm using an offset binomial model like this (R-code):

glm(Incident~Numbers, offset=Numbers, family=binomial, data=threatdata)

My questions are:

1. Is it okay to have exactly the same variables predicting and in the offset? I want to partial out the tonic increase in incident probability and see if there's anything left, essentially. It makes sense to me but I'm a little cautious in case I'm wrong.

2. Is the offset specified correctly? I know that in poisson models it would read

   offset=log(Numbers)
   

I don't know if there's an equivalent here and I can't seem to find any binomial offsets with Google (major problem being that I keep getting negative binomial which of course is no good).

Answer 1

If you are interested in the probability of an incident given N days of patients on ward then you want a model either like:

mod1 <- glm(incident ~ 1, offset=patients.on.ward, family=binomial)

the offset represents trials, incident is either 0 or 1, and the probability of an incident is constant (no heterogeneity in tendency to generate incidents) and patients do not interact to cause incidents (no contagion). Alternatively, if the chance of an incident is small, which it is for you (or you've thresholded the incident counts without mentioning it to us) then you might prefer the Poisson formulation

log.patients.on.ward <- log(patients.on.ward)
mod2 <- glm(incident ~ 1, offset=log.patients.on.ward, family=poisson)

where the same assumptions apply. The offset is logged because the number of patients on ward has a proportional/multiplicative effect.

Expanding on the second model, maybe you think there are more incidents than would be otherwise expected simply due to increased patient numbers. That is, perhaps patients do interact or are heterogenous. So you try

mod3 <- glm(incident ~ 1 + log.patients.on.ward, family=poisson)

If the coefficient on log.patients.on.ward is significantly different from 1, where it was fixed in mod2, then something may indeed be wrong with your assumptions of no heterogeneity and no contagion. And while you cannot of course distinguish these two (nor either one from other missing variables), you do now have an estimate of how much increasing the number of patients on ward increases the rate / probability of an incidents over and above what you'd expect from chance. In the space of parameters it's 1-coef(mod3)[2] with interval derivable from confint.

Alternatively you can just work with the log quantity and its coefficient directly. If you just want to predict the probability of incident using the number of patients on ward, then this model would be a simple way to do it.

The Questions

1. Is it ok to have dependent variables in your offset? It sounds like a very bad idea to me, but I don't see that you have to.

2. The offset in Poisson regression models for exposure is indeed log(exposure). Perhaps confusingly the use of offset in R's Binomial regression models is basically way to indicate the number of trials. It can always be replaced by a dependent variable defined as cbind(incidents, patients.on.ward-incidents) and no offset. Think of it like this: in the Poisson model it enters on the right hand side behind the log link function, and in the Binomial model it enters on the left hand side in front of the logit link function.

Answer 2

Offsets in Poisson regressions

Lets start by looking at why we use an offset in a Poisson regression. Often we want to due this to control for exposure. Let $\lambda$ be the baseline rate per unit of exposure and $t$ be the exposure time in the same units. The expected number of events will be $\lambda \times t$.

In a GLM model we are modelling the expected value using a link function $g$, that is

$$g(\lambda t_i) = \log(\lambda t_i) = \beta_0 + \beta_1x_{1,i} + \dots $$

where $t_i$ is the exposure duration for individual $i$ and $x_i$ is the covariate value for individual $i$. The ellipsis simply indicates additional regression terms we may want to add.

We can simplify simplifying the above expression

$$\log(\lambda) = \log(t_i) + \beta_0 +\beta_1x_{1,i} + \dots$$

The $\log(t_i)$ is simply an "offset" added to the Poisson regression as it is not a product of any of the model parameters which we will be estimating.

Binomial Regression

In a binomial regression, which typically use a logit link, that is:

$$g(p_i) = \textrm{logit}(p_i) = log\left(\frac{p_i}{1-p_i}\right) = \beta_0 +\beta_1x_{1,i}+\dots $$

You can see it will be difficult to derive a model for $p_i$ that will produce a constant offset.

For example, if $p_i$ is the probability that one any patient on day $i$ has an incident. It will be a function of the the individual patients available on that day. As jboman stated it is easier to derive the compliment of no incidence, rather than directly determine probability for at least one incident.

Let $p_{i,j}^*$ be the probability of a patient $j$ having an incident on day $i$. The probability of no patients having an incident on day $i$ will be $\prod_{j=1}^{N_i}(1-p^*_{i,j})$, where $N_i$ is the number of patients on day $i$. By the compliment, the probability of at least one patient having an incident will be, $$p_i = 1-\prod_{j=1}^{N_i}(1-p^*_{i,j}).$$

If we are willing to assume the probability of any patient having an incident on any day is the same we can simplify this to $$p_i = 1-(q^*)^{N_i},$$ where $q^*= 1-p^*$ and $p^*$ is the shared incidence probability.

If we substitute this new definition of $p_i$ back into our logit link function $g(p_i)$, the best we can do in terms of simplification and rearranging is $\log\left((q^*)^{-N} -1 \right)$. This still does not leave us with a constant term that can be factored out.

As a result we cannot use an offset in this case.

The best you can do is discretize the problem (as suggested by jboman) you can create bins for the number of patients and estimate a separate value for $p$ for each of these bins. Otherwise you will need to derive a more complicated model.

Answer 3

This answer comes in two parts, the first a direct answer to the question and the second a commentary on the model you're proposing.

The first part relates to the use of Numbers as an offset along with having it on the r.h.s. of the equation. The effect of doing this will simply be to subtract 1 from the estimated coefficient of Numbers, thereby reversing out the effect of the offset, and will not otherwise change the results. The following example, with a few lines of irrelevant output removed, demonstrates this:

library(MASS)
Numbers <- rpois(100,12)
p <- 1 / (1 + exp(0.25*Numbers))
y <- rbinom(100, Numbers, p)
Incident <- pmin(y, 1) 

> summary(glm(Incident~Numbers, family="binomial"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3121  -1.0246  -0.8731   1.2512   1.7465  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  0.99299    0.80624   1.232   0.2181  
Numbers     -0.11364    0.06585  -1.726   0.0844 . <= COEFFICIENT WITH NO OFFSET TERM
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 135.37  on 99  degrees of freedom
Residual deviance: 132.24  on 98  degrees of freedom
AIC: 136.24

> summary(glm(Incident~Numbers, offset=Numbers, family="binomial"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3121  -1.0246  -0.8731   1.2512   1.7465  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.99299    0.80624   1.232    0.218    
Numbers     -1.11364    0.06585 -16.911   <2e-16 *** <= COEFFICIENT WITH OFFSET TERM
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 342.48  on 99  degrees of freedom
Residual deviance: 132.24  on 98  degrees of freedom
AIC: 136.24

Note how everything is the same except the coefficient of Numbers and the null deviance (and the t-statistic, because it's still testing against 0 instead of -1.)

The second part relates to the model you are building. Since the incidents are recorded not as the number of incidents in a day but whether there were any incidents in a day, the probability of observing a 1 on day $t$ is $1-(1-p_t)^{N_t}$, where $N_t$ is the number of patients on day $t$ and $p_t$ is the per-patient probability of an incident on day $t$. The usual link function, the logit, would parameterize this as $\log(1-(1-p_t)^{N_t})/N_t\log(1-p_t)$. This indicates that the relationship between the probability of observing a 1 on day $t$ and $N_t$ may not be well-modeled by a linear function on the logit scale. (This may be the case anyway, as one might expect some rough "threshold" below which the quality of patient care is OK but above which the quality of patient care drops rapidly.) Reversing the definition of the probabilities so as to move the $N_t$ in the denominator instead of the numerator still leaves you with that awkward exponential inside the log.

One might also suspect that the per-patient probability varies from patient to patient, which would lead to a more complex, hierarchical model, but I won't go into that here.

In any case, given this and the limited range of the number of patients you observe, rather than use a model that is linear on the logit scale, it might be better to be nonparametric about the relationship and group the number of patients into three or four groups, for example, 10-11, 12-13, 14-15, and 16-17, construct dummy variables for those groups, then run the logistic regression with the dummy variables on the right hand side. This will better enable the capture of nonlinear relationships such as "the system is overloaded around 16 patients and incidents start to ramp up significantly." If you had a much wider range of patients, I'd suggest a generalized additive model, e.g., 'gam' from the 'mgcv' package.

Answer 4

Seems simplest to specify a log-link and keep the offset as for a Poisson model.