Usando lm para la prueba de proporción de la muestra 2

Question

He estado usando modelos lineales para llevar a cabo 2 muestra la proporción de las pruebas por un tiempo, pero se han dado cuenta de que podría no ser del todo correcto. Parece que el uso de un modelo lineal generalizado con un binomio de la familia + de identidad enlace da exactamente el unpooled 2-proporción de la muestra los resultados de la prueba. Sin embargo, el uso de un modelo lineal (o glm con gaussiano de la familia) le da un poco diferente resultado. Yo soy la racionalización de que esto podría ser debido a la forma de R soluciona glm para binomial vs gaussiano familias, pero podría haber otra causa?

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

No tiene que ver con cómo resolver los problemas de optimización que corresponden a la colocación de los modelos, que tiene que ver con los reales problemas de optimización de los modelos de la pose.

Específicamente, en muestras grandes, se puede considerar como la comparación de dos mínimos cuadrados ponderados problemas

El modelo lineal (lm) se supone (cuando no ponderado) que la varianza de la proporción es constante. El glm se supone que la varianza de la proporción proviene de la binomial suposición $\text{Var}(\hat{p})=\text{Var}(X/n) = p(1-p)/n$. Esto de los pesos de los puntos de datos de manera diferente, y así llega a algo diferente de las estimaciones* y diferentes de la varianza de las diferencias.

* al menos en algunas situaciones, aunque no necesariamente en una recta de comparación de proporciones

Answer 2

En términos de cálculo, comparar el error estándar de la treatmentB coeficiente de lm vs binomio glm. Usted tiene la fórmula para el error estándar de la treatmentB coeficiente en el binomio glm (el denominador de z_unpooled). El error estándar de la treatmentB coeficiente en el estándar de la película es (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))

Ver este post para una derivación, con la única diferencia de que aquí el error de muestreo se encuentra en lugar de $\sigma^2$ (es decir, restar 2 a $n_A+n_B$ por la pérdida de grados de libertad). Sin ese $-2$, la lm y binomial glm errores estándar en realidad parecen coincidir al $n_A = n_B$.