Uso de glm() como sustituto de la prueba simple de chi cuadrado

Question

Estoy interesado en cambiar las hipótesis nulas utilizando glm() en R.

Por ejemplo:

x = rbinom(100, 1, .7)  
summary(glm(x ~ 1, family = "binomial"))

pone a prueba la hipótesis de que $p = 0.5$ . ¿Y si quiero cambiar el nulo por $p$ = algún valor arbitrario, dentro de glm() ?

Sé que esto se puede hacer también con prop.test() y chisq.test() pero me gustaría explorar la idea de utilizar glm() para comprobar todas las hipótesis relativas a los datos categóricos.

Answer 1

Puede utilizar un offset : glm con family="binomial" estima los parámetros en la escala log-odds o logit, por lo que $\beta_0=0$ corresponde a un logaritmo de 0 o a una probabilidad de 0,5. Si se quiere comparar con una probabilidad de $p$ se desea que el valor de referencia sea $q = \textrm{logit}(p)=\log(p/(1-p))$ . El modelo estadístico es ahora

\begin{split} Y & \sim \textrm{Binom}(\mu) \\ \mu & =1/(1+\exp(-\eta)) \\ \eta & = \beta_0 + q \end{split}

donde sólo la última línea ha cambiado de la configuración estándar. En el código R:

• utilice offset(q) en la fórmula
• la función logit/log-odds es qlogis(p)

• Un poco molesto, usted tiene que proporcionar un valor de desplazamiento para cada elemento en la variable de respuesta - R no replicará automáticamente un valor constante para usted. Esto se hace a continuación mediante la creación de un marco de datos, pero usted podría simplemente utilizar rep(q,100) .

  x = rbinom(100, 1, .7) dd <- data.frame(x, q = qlogis(0.7)) summary(glm(x ~ 1 + offset(q), data=dd, family = "binomial"))

Answer 2

Mire el intervalo de confianza para los parámetros de su GLM:

> set.seed(1)
> x = rbinom(100, 1, .7)
> model<-glm(x ~ 1, family = "binomial")
> confint(model)
Waiting for profiling to be done...
    2.5 %    97.5 % 
0.3426412 1.1862042 

Se trata de un intervalo de confianza para las probabilidades logarítmicas.

Para $p=0.5$ tenemos $\log(odds) = \log \frac{p}{1-p} = \log 1 = 0$ . Así que la prueba de la hipótesis de que $p=0.5$ equivale a comprobar si el intervalo de confianza contiene 0. Éste no lo contiene, por lo que se rechaza la hipótesis.

Ahora, para cualquier $p$ se puede calcular el logaritmo de las probabilidades y comprobar si está dentro del intervalo de confianza.

Answer 3

No es (totalmente) correcto/exacto utilizar los valores p basados en los valores z/t de la función glm.summary como prueba de hipótesis.

1. Es un lenguaje confuso. Los valores indicados se denominan valores z. Pero en este caso utilizan el estimado error estándar en lugar de la desviación real. Por lo tanto, en realidad están más cerca de los valores t . Compara las tres salidas siguientes:
   1) resumen.glm
   2) Prueba t
   3) Prueba z

   > set.seed(1)
   > x = rbinom(100, 1, .7)
   
   > coef1 <- summary(glm(x ~ 1, offset=rep(qlogis(0.7),length(x)), family = "binomial"))$coefficients
   > coef2 <- summary(glm(x ~ 1, family = "binomial"))$coefficients
   
   > coef1[4]  # output from summary.glm
   [1] 0.6626359
   > 2*pt(-abs((qlogis(0.7)-coef2[1])/coef2[2]),99,ncp=0) # manual t-test
   [1] 0.6635858
   > 2*pnorm(-abs((qlogis(0.7)-coef2[1])/coef2[2]),0,1) # manual z-test
   [1] 0.6626359

2. No son valores p exactos. Un cálculo exacto del valor p utilizando la distribución binomial funcionaría mejor (con la potencia de cálculo actual, esto no es un problema). La distribución t, asumiendo una distribución gaussiana del error, no es exacta (sobreestima p, superando el nivel alfa ocurre con menos frecuencia en la "realidad"). Véase la siguiente comparación:

   # trying all 100 possible outcomes if the true value is p=0.7
   px <- dbinom(0:100,100,0.7)
   p_model = rep(0,101)
   for (i in 0:100) {
     xi = c(rep(1,i),rep(0,100-i))
     model = glm(xi ~ 1, offset=rep(qlogis(0.7),100), family="binomial")
     p_model[i+1] = 1-summary(model)$coefficients[4]
   }
   
   # plotting cumulative distribution of outcomes
   outcomes <- p_model[order(p_model)]
   cdf <- cumsum(px[order(p_model)])
   plot(1-outcomes,1-cdf, 
        ylab="cumulative probability", 
        xlab= "calculated glm p-value",
        xlim=c(10^-4,1),ylim=c(10^-4,1),col=2,cex=0.5,log="xy")
   lines(c(0.00001,1),c(0.00001,1))
   for (i in 1:100) {
     lines(1-c(outcomes[i],outcomes[i+1]),1-c(cdf[i+1],cdf[i+1]),col=2)
   #  lines(1-c(outcomes[i],outcomes[i]),1-c(cdf[i],cdf[i+1]),col=2)
   }
   
   title("probability for rejection as function of set alpha level")

   

   La curva negra representa la igualdad. La curva roja está por debajo de ella. Esto significa que para un valor p calculado por la función de resumen glm, encontramos esta situación (o una diferencia mayor) con menos frecuencia en la realidad de lo que indica el valor p.