Menggunakan lm untuk uji proporsi 2 sampel

Saya telah menggunakan model linier untuk melakukan tes proporsi 2 sampel untuk sementara waktu, tetapi telah menyadari bahwa mungkin tidak sepenuhnya benar. Tampaknya menggunakan model linier umum dengan tautan keluarga + identitas binomial memberikan hasil tes proporsi 2 sampel yang tidak dikumpulkan. Namun, menggunakan model linier (atau glm dengan keluarga gaussian) memberikan hasil yang sedikit berbeda. Saya merasionalisasi bahwa ini mungkin karena bagaimana R memecahkan glm untuk keluarga binomial vs gaussian, tetapi mungkinkah ada penyebab lain?

## 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

Ini tidak ada hubungannya dengan bagaimana mereka menyelesaikan masalah optimisasi yang sesuai dengan pemasangan model, itu berkaitan dengan masalah optimisasi aktual yang ditimbulkan oleh model.

Khususnya, dalam sampel besar, Anda dapat secara efektif menganggapnya sebagai membandingkan dua masalah kuadrat terkecil

Model linier ( lm) kita mengasumsikan (ketika tidak tertimbang) bahwa varians dari proporsi adalah konstan. Glm mengasumsikan bahwa varian proporsi berasal dari asumsi binomial . Ini menimbang titik data secara berbeda, sehingga muncul estimasi yang agak berbeda * dan perbedaan ragam yang berbeda.$\text{Var}(\hat{p})=\text{Var}(X/n) = p(1-p)/n$

* setidaknya dalam beberapa situasi, meskipun tidak harus dalam perbandingan proporsi yang lurus

Dalam hal perhitungan, bandingkan kesalahan standar dari koefisien treatmentB untuk lm vs binomial glm. Anda memiliki rumus untuk kesalahan standar dari koefisien treatmentB dalam glomer binomial (penyebut z_unpooled). Kesalahan standar dari koefisien treatmentB dalam lm standar adalah (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))

Lihat posting ini untuk derivasi, satu-satunya perbedaan adalah bahwa di sini kesalahan sampel ditemukan alih-alih (yaitu kurangi 2 dari untuk kehilangan derajat kebebasan). Tanpa itu , kesalahan standar lm dan binomial glm benar-benar cocok ketika .$\sigma^2$ - 2 n A = n $n_A+n_B$$-2$$n_A = n_B$