¿Cómo puedo calcular Pearson $\chi^2$ estadística de prueba de falta de ajuste en un modelo de regresión logística en R?

Question

El cociente de probabilidad (un.k.una. desviación) $G^2$ estadística y la falta de ajuste (o de bondad de ajuste) de la prueba es bastante sencillo para obtener un modelo de regresión logística (ajuste utilizando el glm(..., family = binomial) función) en R. sin Embargo, puede ser fácil de tener algunos de los recuentos de células terminar lo suficientemente bajo para que la prueba no es fiable. Una manera de verificar la fiabilidad de la prueba de razón de verosimilitud para la falta de ajuste para comparar el estadístico de prueba y el P-valor para aquellos de Pearson chi cuadrado (o $\chi^2$) la falta de prueba de ajuste.

Ni el glm objeto ni su summary() método de informe de la prueba estadística de Pearson chi-cuadrado de la prueba de falta de ajuste. En mi búsqueda, lo único que se me ocurrió es la chisq.test() de la función (en stats paquete): su documentación dice "chisq.test realiza chi-cuadrado de la tabla de contingencia de las pruebas y de bondad de ajuste de las pruebas." Sin embargo, la documentación es escasa sobre cómo realizar estas pruebas:

Si x es una matriz con una fila o columna, o si x es un vector y y no se da, entonces, una bondad de ajuste prueba se realiza (x es tratado como un one-dimensional de la tabla de contingencia). Las entradas de x deben ser enteros no negativos. En este caso, la hipótesis planteada es si la población probabilidades iguales a las que en p, o son todos iguales si p no se da.

Me imagino que usted podría utilizar el y componente de la glm objeto de la x argumento de chisq.test. Sin embargo, no puede usar la fitted.values componente de la glm objeto de la p argumento de chisq.test, porque obtendrá un error: "probabilities must sum to 1."

¿Cómo puedo (en R) al menos calcular la prueba de Pearson $\chi^2$ estadística de prueba de falta de ajuste, sin tener que ejecutar a través de los pasos de forma manual?

Answer 1

La suma de los cuadrados de los residuos de Pearson es exactamente igual a la de Pearson $\chi^2$ estadística de prueba de falta de ajuste. Así que si tu modelo ajustado (es decir, el glm objeto) se llama logistic.fit, el siguiente código de retorno de la estadística de prueba:

sum(residuals(logistic.fit, type = "pearson")^2)

Consulte la documentación en residuals.glm para obtener más información, incluyendo lo que otros residuos están disponibles. Por ejemplo, el código de

sum(residuals(logistic.fit, type = "deviance")^2)

va a llegar la $G^2$ estadística de prueba, de la misma manera deviance(logistic.fit) proporciona.

Answer 2

La prueba de Pearson estadístico tiene una distribución degenerada por lo que no se recomienda en general para el modelo logístico de bondad de ajuste. Yo prefiero estructurado pruebas (linealidad, la aditividad). Si quieres un ómnibus de la prueba de ver el único grado de libertad que le Cessie - van Houwelingen - Copas - Hosmer no ponderado suma de los cuadrados de la prueba, como el implementado en el R rms paquete residuals.lrm función.

Answer 3

También puede utilizar c_hat(mod) que va a dar el mismo resultado como sum(residuals(mod, type = "pearson")^2).

Answer 4

Gracias, no me di cuenta de que era tan simple como: suma(residuos(f1, type="pearson")^2) Sin embargo, por favor tenga en cuenta que Pearsons residual varía dependiendo de si se calcula por la covariable grupo o individuales. Un ejemplo sencillo:

el m1 es una matriz (esta es la cabeza de una gran matriz):

m1[1:4,1:8]

    x1 x2 x3 obs    pi   lev  yhat y
obs  1  1 44   5 0.359 0.131 1.795 2
obs  0  1 43  27 0.176 0.053 4.752 4
obs  0  1 53  15 0.219 0.062 3.285 1
obs  0  1 33  22 0.140 0.069 3.080 3

Donde x1-3 son predictores, obs es no. observaciones en cada grupo, pi es la probabilidad de pertenencia al grupo (predecir a partir de la ecuación de regresión), lev es el apalancamiento, la diagonal de la matriz hat, yhat la predicción no. (y=1) en el grupo e y el real no.

Esto le dará a usted de Pearson por el grupo. Nota cómo es diferente si y==0: '$ fun1 <- function(j){ si (m1[j,"y"] ==0){ # y=0 para esta covariable patrón Pr1 <- sqrt( m1[i,"pi"] / (1-m1[i,"pi"])) Pr2 <- -sqrt (m1[i,"obs"]) res <- round( Pr1 * Pr2, 3) return(res) } else { Pr1 <- m1[j,"y"] - m1[j,"yhat"] Pr2 <- sqrt( m1[j,"yhat"] * ( 1-(m1[j,"pi"]) ) ) res <- round( Pr1/Pr2, 3) return(res) } } $'

Así

nr <-nrow(m1)
nr1 <- seq(1,nr)
Pr <- sapply(1:nrow[m1], FUN=fun1)
PrSj <- sum(Pr^2)

Si hay un gran número de sujetos con y=0 covariable patrones, entonces Pearons residual será mucho mayor cuando se calcula utilizando el 'grupo' en lugar de la 'individual' método.

Ver, por ejemplo, de Hosmer & Lemeshow ", que se Aplica de Regresión Logística", Wiley, 200.