In using the cbind() function in R for a logistic regression on a $2 \times 2$ table, what is the explicit functional form of the regression equation?

Question

Suppose I have a $2 \times 2$ table that looks like:

            Disease       No Disease
Treatment         55                67
Control           42                34

I would like to do a logistic regression in R on this table. I understand that the standard way is to use the glm function with a cbind function in the response. In other words, the code looks like:

glm(formula = cbind(c(55,67),c(42,34)) ~ as.factor(c(1, 0)), family = binomial())

I am wondering why R requires us to use the cbind function and why simply using proportions is not sufficient. Is there a way to write this out explicitly as a formula? What would it look in the form of:

$$ log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1X $$

where $X = 1$ if we have treatment and $X=0$ for control?

Right now it seems like I am regressing on a matrix for the dependent value.

Answer 1

First I show how you can specify a formula using aggregated data with proportions and weights. Then I show how you could specify a formula after dis-aggregating your data to individual observations.

Documentation inglm indicates that:

"For a binomial GLM prior weights are used to give the number of trials when the response is the proportion of successes"

I create new columns total and proportion_disease in df for the 'number of trials' and 'proportion of successes' respectively.

library(dplyr)
df <- tibble(treatment_status = c("treatment", "no_treatment"),
       disease = c(55, 42),
       no_disease = c(67,34)) %>% 
  mutate(total = no_disease + disease,
         proportion_disease = disease / total) 

model_weighted <- glm(proportion_disease ~ treatment_status, data = df, family = binomial("logit"), weights = total)

The above weighted approach takes in aggregated data and will give the same solution as the cbind method but allows you to specify a formula. (Below is equivalent to Original Poster's method but cbind(c(55,42), c(67,34)) rather than cbind(c(55,67),c(42,34)) so that 'Disease' rather than 'Treatment' is the response variable.)

model_cbinded <- glm(cbind(disease, no_disease) ~ treatment_status, data = df, family = binomial("logit"))  

You could also just dis-aggregate your data into individual observations and pass these into glm (allowing you to specify a formula as well).

df_expanded <- tibble(disease_status = c(1, 1, 0, 0), 
                      treatment_status = rep(c("treatment", "control"), 2)) %>%
                        .[c(rep(1, 55), rep(2, 42), rep(3, 67), rep(4, 34)), ]

model_expanded <- glm(disease_status ~ treatment_status, data = df_expanded, family = binomial("logit"))

Let's compare these now by passing each model into summary. model_weighted and model_cbinded both produce the exact same results. model_expanded produces the same coefficients and standard errors, though outputs different degrees of freedom, deviance, AIC, etc. (corresponding with the number of rows/observations).

    > lapply(list(model_weighted, model_cbinded, model_expanded), summary)
[[1]]

Call:
glm(formula = proportion_disease ~ treatment_status, family = binomial("logit"), 
    data = df, weights = total)

Deviance Residuals: 
[1]  0  0

Coefficients:
                          Estimate Std. Error z value Pr(>|z|)
(Intercept)                 0.2113     0.2307   0.916    0.360
treatment_statustreatment  -0.4087     0.2938  -1.391    0.164

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.9451e+00  on 1  degrees of freedom
Residual deviance: 1.0658e-14  on 0  degrees of freedom
AIC: 14.028

Number of Fisher Scoring iterations: 2


[[2]]

Call:
glm(formula = cbind(disease, no_disease) ~ treatment_status, 
    family = binomial("logit"), data = df)

Deviance Residuals: 
[1]  0  0

Coefficients:
                          Estimate Std. Error z value Pr(>|z|)
(Intercept)                 0.2113     0.2307   0.916    0.360
treatment_statustreatment  -0.4087     0.2938  -1.391    0.164

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.9451e+00  on 1  degrees of freedom
Residual deviance: 1.0658e-14  on 0  degrees of freedom
AIC: 14.028

Number of Fisher Scoring iterations: 2


[[3]]

Call:
glm(formula = disease_status ~ treatment_status, family = binomial("logit"), 
    data = df_expanded)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.268  -1.095  -1.095   1.262   1.262  

Coefficients:
                          Estimate Std. Error z value Pr(>|z|)
(Intercept)                 0.2113     0.2307   0.916    0.360
treatment_statustreatment  -0.4087     0.2938  -1.391    0.164

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 274.41  on 197  degrees of freedom
Residual deviance: 272.46  on 196  degrees of freedom
AIC: 276.46

Number of Fisher Scoring iterations: 3

(See R bloggers for conversation on weights parameter in glm in the regression context.)

Answer 2

Binomial distribution versus Bernoulli distribution

The binomial distribution can be viewed as multiple Bernoulli distributions with the same probability parameter $p$. Their probability distributions differ up to a constant.

For observations $x_i \in \lbrace 0,1 \rbrace$, where $k$ is the number of observations that are to $1$

• $1$ Bernoulli experiment PMF $$f(x_1) = p^{x_1} (1-p)^{1-x_1}$$
• $n$ Bernoulli experiments PMF $$f(x_1,x_2,\dots,x_n) = \prod p^{x_i} (1-p)^{1-x_i} = p^{k} (1-p)^{n-k}$$
• Binomial distribution PMF $$f(x_1,x_2,\dots,x_n) = f(k,n) = {n \choose k} p^{k} (1-p)^{n-k}$$

The difference is that the Binomial distribution sort of 'ignores' the order of the observations $x_i$ and only looks at the total number of $x_i=1$ and $x_i=0$.

Influence of sample size on likelihood distribution and error

why simply using proportions is not sufficient

The maximum likelihood estimate is maximimizing $p^k(1-p)^{n-k}$ and this happens at $p_{max} = k/n$ the proportion of observations with $x_i=1$. So indeed for finding the maximum likelihood the proportion is sufficient.

However, when we wish to estimate some expression for the error of the estimate (for instance standard error or confidence intervals) then the number of observations becomes important.

An intuitive view to see how the sample size matters might be the following: Let's look at the sample distribution of the outcome of the experiment as a function of the sample size. In the image below we plot the binomial distribution for $p=0.8$ with different size parameters $n$. On the left you see that the curves shift to the right for larger $n$, which is obvious as you would expect more counts when the sample size $n$ is increased. But the distance from the expected rate $p$ times $n$ is becoming relatively shorter. You see this on the right where we divide the outcome (the counts) by the total sample size $n$, which gives a relative rate. You see that the observed relative rate is more likely to be close to the true rate (0.8 in this graph) when the size $n$ is larger. (A related question is: How to estimate a probability of an event to occur based on its count?).

Input for the GLM function

I am wondering why R requires us to use the cbind function

You are not required to use the cbind function. There are 2/3 methods (two or three depending on how you look at it: two ways to present the data and the last one splits up in two ways to feed it to the glm function).

• We can also use GLM as if it was doing a Bernoulli variable. In that case the data would consist of $55+67+42+34 = 198$ rows. One row for each single count (each count can be considered a Bernoulli distributed variable)

  So the data would look like a long table

  $$\begin{array}{rrr} \text{ID} & \text{Observation of disease} & \text{Treatment}\\ 1 & 1 & 1\\ 2 & 1 & 1\\ 3 & 1 & 1\\ \vdots & \vdots & \vdots \\ 53 & 1 & 1\\ 54 & 1 & 1\\ 55 & 1 & 1\\ 56 & 0 & 1\\ 57 & 0 & 1\\ 58 & 0 & 1\\ \vdots & \vdots & \vdots \\ 120 & 0 & 1\\ 121 & 0 & 1\\ 122 & 0 & 1\\ 123 & 1 & 0\\ 124 & 1 & 0\\ 125 & 1 & 0\\ \vdots & \vdots & \vdots \\ 162 & 1 & 0\\ 163 & 1 & 0\\ 164 & 1 & 0\\ 165 & 0 & 0\\ 166 & 0 & 0\\ 167 & 0 & 0\\ \vdots & \vdots & \vdots \\ 196 & 0 & 0\\ 197 & 0 & 0\\ 198 & 0 & 0\\ \end{array}$$

   data <- data.frame(ID = 1:198, 
                      observation = c(rep(1,55),rep(0,67),rep(1,42),rep(0,34)),
                      treatment = c(rep(1,55+67), rep(0,42+34)))
   glm1 <- glm(observation ~ treatment, family = binomial, data = data) 
  

• The previous method is often used when the independent variable is continuous. In the case of a categorical variable, we can group all the 1s and 0s for a category. In this case, the data is not 0s and 1s but the counts of 0s and 1s.

   data <- data.frame(disease = c(55,42),
                      healthy = c(67,34),
                      treatment = c(1,0))
  
   ### the 'cbind method' gives glm the numbers in each categories
   glm2 <- glm(cbind(disease,healthy) ~ treatment, family = binomial, data = data)
  
   ### this method gives glm the fraction and uses 'weight' to describe the number of data
   glm3 <- glm(disease/(disease+healthy) ~ treatment, family = binomial, data = data, weight = disease+healthy)
  

The summaries of the above three different calls to glm give the same results. There is a difference in statistics like deviance and AIC, that is because the likelihood function is different by a factor $n \choose k$ but that is not important because the relevant measures in comparing models are the ratios of likelihood.

> summary(glm1)

Call:
glm(formula = observation ~ treatment, family = binomial, data = data)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.268  -1.095  -1.095   1.262   1.262  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.2113     0.2307   0.916    0.360
treatment    -0.4087     0.2938  -1.391    0.164

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 274.41  on 197  degrees of freedom
Residual deviance: 272.46  on 196  degrees of freedom
AIC: 276.46

Number of Fisher Scoring iterations: 3

> summary(glm2)

Call:
glm(formula = cbind(disease, healthy) ~ treatment, family = binomial, 
    data = data)

Deviance Residuals: 
[1]  0  0

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.2113     0.2307   0.916    0.360
treatment    -0.4087     0.2938  -1.391    0.164

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.9451e+00  on 1  degrees of freedom
Residual deviance: 1.0658e-14  on 0  degrees of freedom
AIC: 14.028

Number of Fisher Scoring iterations: 2

> summary(glm3)

Call:
glm(formula = disease/(disease + healthy) ~ treatment, family = binomial, 
    data = data, weights = disease + healthy)

Deviance Residuals: 
[1]  0  0

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.2113     0.2307   0.916    0.360
treatment    -0.4087     0.2938  -1.391    0.164

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.9451e+00  on 1  degrees of freedom
Residual deviance: 1.0658e-14  on 0  degrees of freedom
AIC: 14.028

Number of Fisher Scoring iterations: 2

>