Using lm for 2-sample proportion test

Question

I have been using linear models to perform 2-sample proportion tests for a while, but have realized that might not be completely correct. It appears that using a generalized linear model with a binomial family + identity link gives exactly the unpooled 2-sample proportion test results. However, using a linear model (or glm with gaussian family) gives a slightly different result. I'm rationalizing that this might be due to how R solves glm for binomial vs. gaussian families, but could there be another cause?

## prop.test gives pooled 2-sample proportion result
## glm w/ binomial family gives unpooled 2-sample proportion result
## lm and glm w/ gaussian family give unknown result

library(dplyr)
library(broom)
set.seed(12345)

## set up dataframe -------------------------
n_A <- 5000
n_B <- 5000

outcome <- rbinom(
  n = n_A + n_B,
  size = 1,
  prob = 0.5
)
treatment <- c(
  rep("A", n_A),
  rep("B", n_B)
)

df <- tbl_df(data.frame(outcome = outcome, treatment = treatment))


## by hand, 2-sample prop tests ---------------------------------------------
p_A <- sum(df$outcome[df$treatment == "A"])/n_A
p_B <- sum(df$outcome[df$treatment == "B"])/n_B

p_pooled <- sum(df$outcome)/(n_A + n_B)
z_pooled <- (p_B - p_A) / sqrt( p_pooled * (1 - p_pooled) * (1/n_A + 1/n_B) )
pvalue_pooled <- 2*(1-pnorm(abs(z_pooled)))

z_unpooled <- (p_B - p_A) / sqrt( (p_A * (1 - p_A))/n_A + (p_B * (1 - p_B))/n_B )
pvalue_unpooled <- 2*(1-pnorm(abs(z_unpooled)))


## using prop.test --------------------------------------
res_prop_test <- tidy(prop.test(
  x = c(sum(df$outcome[df$treatment == "A"]), 
        sum(df$outcome[df$treatment == "B"])),
  n = c(n_A, n_B),
  correct = FALSE
))
res_prop_test # same as pvalue_pooled
all.equal(res_prop_test$p.value, pvalue_pooled)
# [1] TRUE


# using glm with identity link -----------------------------------
res_glm_binomial <- df %>%
  do(tidy(glm(outcome ~ treatment, family = binomial(link = "identity")))) %>%
  filter(term == "treatmentB")
res_glm_binomial # same as p_unpooled
all.equal(res_glm_binomial$p.value, pvalue_unpooled)
# [1] TRUE


## glm and lm gaussian --------------------------------

res_glm <- df %>%
  do(tidy(glm(outcome ~ treatment))) %>%
  filter(term == "treatmentB")
res_glm 
all.equal(res_glm$p.value, pvalue_unpooled)
all.equal(res_glm$p.value, pvalue_pooled)

res_lm <- df %>%
  do(tidy(lm(outcome ~ treatment))) %>% 
  filter(term == "treatmentB")
res_lm
all.equal(res_lm$p.value, pvalue_unpooled)
all.equal(res_lm$p.value, pvalue_pooled)

all.equal(res_lm$p.value, res_glm$p.value)
# [1] TRUE

Answer 1

It's not to do with how they solve the optimization problems that correspond to fitting the models, it's to do with the actual optimization problems the models pose.

Specifically, in large samples, you can effectively consider it as comparing two weighted least squares problems

The linear model (lm) one assumes (when unweighted) that the variance of the proportions is constant. The glm assumes that the variance of the proportions comes from the binomial assumption $\text{Var}(\hat{p})=\text{Var}(X/n) = p(1-p)/n$. This weights the data points differently, and so comes to somewhat different estimates* and different variance of differences.

* at least in some situations, though not necessarily in a straight comparison of proportions

Answer 2

In terms of calculation, compare the standard error of the treatmentB coefficient for lm vs. binomial glm. You have the formula for the standard error of the treatmentB coefficient in the binomial glm (the denominator of z_unpooled). The standard error of the treatmentB coefficient in the standard lm is (SE_lm):

    test = lm(outcome ~ treatment, data = df)
    treat_B =  as.numeric(df$treatment == "B")
    SE_lm = sqrt( sum(test$residuals^2)/(n_A+n_B-2) / 
              sum((treat_B - mean(treat_B))^2))

See this post for a derivation, the only difference being that here the sample error is found instead of $\sigma^2$ (i.e. subtract 2 from $n_A+n_B$ for lost degrees of freedom). Without that $-2$, the lm and binomial glm standard errors actually seem to match when $n_A = n_B$.