Logistic regression is most commonly used with a Bernoulli response. If \(Y \sim Bernoulli(p)\), then \(P(Y = 1) = p = 1 - P(Y = 0)\). Logistic regression can be used to estimate \(p|\textbf{x}\) where \(\textbf{x}\) is a collection of predictor variables. This is accomplished by assuming the relationship
\[log(\frac{p}{1-p}) = \textbf{x}\beta\]
and then choosing \(\beta\) to maximize the joint likelihood:
\[\Pi_i p_i^{y_i}(1-p_i)^{(1-y_i)}\]
Note that the function linking the predictor variables to \(p\) is called the logit function, which is what gives logistic regression its name. In \(\textbf{R}\), this model can be fit using glm by specifying the options family = "binomial".
The binomial distribution is a generalization of the Bernoulli, typically conceptualized as describing the number of “successes” (here denoted \(y\)) out of some number of attempts or trials (here denoted \(n\)). The probability mass function for a binomial random variable, \(Y\), is
\[P(Y = y|n, p) = {x \choose n} p^x(1 - p)^{n - x}\]
From this, we can see that by setting \(n = 1\), we recover the Bernoulli PMF.
Since the binomial distribution is a generalization of the Bernoulli distribution, it stands to reason that it may be possible to generalize logistic regression for Bernoulli variables to binomial variables. Indeed, this is the case, and it turns out it is relatively easy. Using the logit link function again, maximize the joint likelihood:
\[\Pi_i {y_i \choose n_i} p_i^{y_i}(1-p_i)^{(1-y_i)}\]
R is also capable of fitting binomial regression models using glm with family = "binomial". However, there are a few syntactic quirks. First, in the model statement, the response variable must be given as a proportion of successes for a given run. This can be done either by providing the proportions as a column in the data matrix or by specifying a ratio of the number of successes to the number of attempts. An example here might be instructive.
df <- data.frame(
x = runif(30, -1, 1),
n = rbinom(30, 6, .8)) %>%
mutate(y = rbinom(30, n, (exp(x)/(1 + exp(x)))))Consider the data set, df, given below.
## x n y
## 1 -0.01132681 5 1
## 2 -0.26414193 4 2
## 3 -0.81380816 4 0
## 4 0.55575400 4 3
## 5 0.65899013 6 3
## 6 0.89954037 4 3Here, \(x\) is some predictor variable related to the number of successes, \(y\), for the given number of attempts, \(n\). We can fit a logistic regression to this data to estimate this relationship:
##
## Call: glm(formula = y/n ~ x, family = "binomial", data = df, weights = n)
##
## Coefficients:
## (Intercept) x
## 0.1041 0.9450
##
## Degrees of Freedom: 29 Total (i.e. Null); 28 Residual
## Null Deviance: 51.94
## Residual Deviance: 42.44 AIC: 93.45We get the same results using different syntax:
## x n y prob
## 1 -0.01132681 5 1 0.20
## 2 -0.26414193 4 2 0.50
## 3 -0.81380816 4 0 0.00
## 4 0.55575400 4 3 0.75
## 5 0.65899013 6 3 0.50
## 6 0.89954037 4 3 0.75##
## Call: glm(formula = prob ~ x, family = "binomial", data = dfProb, weights = n)
##
## Coefficients:
## (Intercept) x
## 0.1041 0.9450
##
## Degrees of Freedom: 29 Total (i.e. Null); 28 Residual
## Null Deviance: 51.94
## Residual Deviance: 42.44 AIC: 93.45Notably, these fitted values are equivalent to what we would get if we were to transform each of our binomial responses into equivalent sets of Bernoulli trials. For example, the first observation in df is a single success out of 5 attempts with a covariate of -0.011. Another way to think of this is five Bernoulli trials, one successful, four unsuccessful, all with a covariate value of -0.011.
For example:
## x y
## 1 -0.01132681 1
## 2 -0.01132681 0
## 3 -0.01132681 0
## 4 -0.01132681 0
## 5 -0.01132681 0
## 6 -0.26414193 1##
## Call: glm(formula = y ~ x, family = "binomial", data = dfTall)
##
## Coefficients:
## (Intercept) x
## 0.1041 0.9450
##
## Degrees of Freedom: 139 Total (i.e. Null); 138 Residual
## Null Deviance: 194.1
## Residual Deviance: 184.6 AIC: 188.6Note that while the degrees of freedom, information criteria, etc. are different, the coefficient estimates are the same, as are the standard errors for these coefficients.
##Using ciTools for binomial regression
ciTools supports logistic regression with both Bernoulli and binomial response variables. For both types, add_ci has intuitive and expected functionality. It produces a point estimate and confidence intervals around the estimated probability of success. However, the other functions of ciTools (add_pi, add_quantile, and add_probs) produce different behaviors, because prediction intervals, quantiles, and probability values are not very useful for Bernoulli variables.
For example, consider prediction intervals. These are used to quantify the observation-to-observation variability and estimate its range. For a Bernoulli variable, only two values are possible, so this is not particularly meaningful – the prediction interval will always be \([0,1]\). Quantiles are similarly problematic. And probability estimates (that is, \(P(Y>y|x)\)) are equivalent to estimates of \(p\) or \(1-p\) . Thus, asking for any of these will produce either an error (in the case of add_pi.glm(family = "binomial) or add_quantiles(family = "binomial)) or a warning (in the case of add_probs.glm(family = "binomial")).
For a binomial response, on the other hand, each of these functions produce meaningful output. Binomial response variables can take more than two values, making it sensible to consider prediction intervals, quantiles, or probabilities. For example, it may be interesting to consider the 90 percent quantile for a binomial regression:
## Warning in add_quantile.glm(df, fit, p = 0.9): Treating weights as indicating
## the number of trials for a binomial regression where the response is the
## proportion of successes## Warning in add_quantile.glm(df, fit, p = 0.9): The response variable is not
## continuous so Prediction Intervals are approximate## Warning in sim_quantile_other(df, fit, p, name, yhatName, nSims): For binomial
## models, add_quantile's column of fitted values reflect E(Y|X) rather than
## typical default for logistic regression, pHat## x n y pred quantile0.9
## 1 -0.01132681 5 1 2.616636 4
## 2 -0.26414193 4 2 1.854723 3
## 3 -0.81380816 4 0 1.358501 3
## 4 0.55575400 4 3 2.609291 4
## 5 0.65899013 6 3 4.044646 6
## 6 0.89954037 4 3 2.887789 4There is a lot to unpack here, including three warning messages. Let’s consider the output first. Unlike when we use add_ci.glm, the pred column includes the predicted values \(E(Y|\boldsymbol{x})\) rather than the estimated probability of success \(\hat{p}|\boldsymbol{x}\).
Next, there are three warnings. None of these indicate that anything has gone wrong. Rather, they are included to ensure that users understand the output they’re given. The first warning makes it clear to the user that the weights argument is assumed to indicate that they’re doing binomial rather than weighted Bernoulli regression. The second warning informs users that the estimated interval is approximate. This is because the current method used by ciTools for GLMs is based on simulation, and because ciTools forces quantile estimates to lie in the support set of the response distribution. The final warning points out that the pred column refers to the fitted values rather than estimated probability of success, which was mentioned in the previous paragraph.
Logistic regression can be used for both Bernoulli and Binomial response variables, and glm supports both. ciTools differentiates between the two based on whether the user has included a column in the weights argument. Aside from add_ci, none of the functions in ciTools produce useful information for Bernoulli response variables. For Binomial response variables, all of these functions produce useful information, though error messages are included to ensure that users understand the output presented.