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 postthis 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$.