¿Cómo calcular pseudo- partir de la regresión logística de R?

46

El artículo de Christopher Manning sobre regresión logística en R muestra una regresión logística en R de la siguiente manera:

ced.logr <- glm(ced.del ~ cat + follows + factor(class), 
  family=binomial)

Alguna salida:

> summary(ced.logr)
Call:
glm(formula = ced.del ~ cat + follows + factor(class),
    family = binomial("logit"))
Deviance Residuals:
Min            1Q    Median       3Q      Max
-3.24384 -1.34325   0.04954  1.01488  6.40094

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.31827    0.12221 -10.787 < 2e-16
catd          -0.16931    0.10032  -1.688 0.091459
catm           0.17858    0.08952   1.995 0.046053
catn           0.66672    0.09651   6.908 4.91e-12
catv          -0.76754    0.21844  -3.514 0.000442
followsP       0.95255    0.07400  12.872 < 2e-16
followsV       0.53408    0.05660   9.436 < 2e-16
factor(class)2 1.27045    0.10320  12.310 < 2e-16
factor(class)3 1.04805    0.10355  10.122 < 2e-16
factor(class)4 1.37425    0.10155  13.532 < 2e-16
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 958.66 on 51 degrees of freedom
Residual deviance: 198.63 on 42 degrees of freedom
AIC: 446.10
Number of Fisher Scoring iterations: 4

Luego entra en algunos detalles sobre cómo interpretar los coeficientes, comparar diferentes modelos, etc. Bastante útil.

Sin embargo, ¿cuánta varianza representa el modelo? Una página de Stata sobre regresión logística dice:

Técnicamente, no se puede calcular de la misma manera en la regresión logística que en la regresión OLS. El pseudo- , en regresión logística, se define como , donde representa la probabilidad logarítmica para el modelo "solo constante" y es la probabilidad logarítmica para el modelo completo con constante y predictores.R 2 1 - L 1R2R2 L0L11L1L0L0L1

Entiendo esto en el alto nivel. El modelo de solo constante no tendría ninguno de los parámetros (solo el término de intercepción). La probabilidad de registro es una medida de qué tan cerca se ajustan los parámetros a los datos. De hecho, Manning tipo de indicios de que la desviación podría ser . ¿Quizás la desviación nula es constante y la desviación residual es del modelo? Sin embargo, no lo tengo claro.- 2 log L2logL2logL

¿Alguien puede verificar cómo se calcula realmente el pseudo- en R usando este ejemplo?R2

dfrankow
fuente
55
Las páginas de computación estadística UCLA, por lo general excelentes, han cometido un error raro aquí: no debe haber paréntesis en la expresión para pseudo- , es decir, debe ser . (Lo siento por no responder a sus consultas, estoy a punto de la cabeza de la cama - Estoy seguro de que alguien más va a haber contestado antes Estoy lo suficientemente despierto como para hacerlo.) 1 - L 1 / L 0R21L1/L0
onestop
3
Esta página discute varios pseudo-R ^ 2s.
dfrankow
2
Nota: a la pregunta relacionada no le gustan los pseudo-R ^ 2s, pero prefiere la validación cruzada o la predicción de prueba de resistencia.
dfrankow

Respuestas:

49

No olvides el paquete rms , de Frank Harrell. Encontrarás todo lo que necesitas para ajustar y validar GLM.

Aquí hay un ejemplo de juguete (con un solo predictor):

set.seed(101)
n <- 200
x <- rnorm(n)
a <- 1
b <- -2
p <- exp(a+b*x)/(1+exp(a+b*x))
y <- factor(ifelse(runif(n)<p, 1, 0), levels=0:1)
mod1 <- glm(y ~ x, family=binomial)
summary(mod1)

Esto produce:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.8959     0.1969    4.55 5.36e-06 ***
x            -1.8720     0.2807   -6.67 2.56e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 258.98  on 199  degrees of freedom
Residual deviance: 181.02  on 198  degrees of freedom
AIC: 185.02

Ahora, usando la lrmfunción,

require(rms)
mod1b <- lrm(y ~ x)

Pronto obtendrá muchos índices de ajuste del modelo, incluido Nagelkerke , con :R2print(mod1b)

Logistic Regression Model

lrm(formula = y ~ x)

                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       

Obs           200    LR chi2      77.96    R2       0.445    C       0.852    
 0             70    d.f.             1    g        2.054    Dxy     0.705    
 1            130    Pr(> chi2) <0.0001    gr       7.801    gamma   0.705    
max |deriv| 2e-08                          gp       0.319    tau-a   0.322    
                                           Brier    0.150                     


          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept  0.8959 0.1969  4.55  <0.0001 
x         -1.8720 0.2807 -6.67  <0.0001 

Aquí, y se calcula como , donde LR es la estadística (comparando los dos modelos anidados que describió), mientras que el denominador es solo el valor máximo para . Para un modelo perfecto, esperaríamos , es decir .( 1 - exp ( - LR / n ) ) / ( 1 - exp ( - ( - 2 L 0 ) / n ) )R2=0.445(1exp(LR/n))/(1exp((2L0)/n))R 2 LR = 2 L 0 R 2 = 1χ2R2LR=2L0R2=1

A mano,

> mod0 <- update(mod1, .~.-x)
> lr.stat <- lrtest(mod0, mod1)
> (1-exp(-as.numeric(lr.stat$stats[1])/n))/(1-exp(2*as.numeric(logLik(mod0)/n)))
[1] 0.4445742
> mod1b$stats["R2"]
       R2 
0.4445742 

Ewout W. Steyerberg discutió el uso de con GLM, en su libro Clinical Prediction Models (Springer, 2009, § 4.2.2 pp. 58-60). Básicamente, la relación entre el estadístico LR y el Nagelkerke es aproximadamente lineal (será más lineal con baja incidencia). Ahora, como se discutió en el hilo anterior al que me vinculé en mi comentario, puede usar otras medidas como el estadístico que es equivalente al estadístico AUC (también hay una buena ilustración en la referencia anterior, consulte la Figura 4.6).R2R2c

chl
fuente
¿Puede explicar cómo obtuvo .445? Usé 1-exp (-77.96 / 200) pero obtuve .323. ¿Que estoy haciendo mal? Gracias.
2
¿Cuál es Nagelkerke R2?
JetLag
1
@JetLag Bajo los índices de discriminación, el Nagelkerke se abrevia como R2 (es decir, 0.445). Puede verificar esto usando la función NagelkerkeR2 () del paquete fmsb.
Chernoff
7

Tenga cuidado con el cálculo de Pseudo-R2 :

El Pseudo- McFadden se calcula como , donde es la probabilidad logarítmica del modelo completo, y es la probabilidad logarítmica del modelo con solo intercepción.R2RM2=1lnL^fulllnL^nulllnL^fulllnL^full

Dos enfoques para calcular Pseudo- :R2

  1. Usar desviación: dado que ,deviance=2ln(Lfull)null.deviance=2ln(Lnull)

    pR2 = 1 - mod$deviance / mod$null.deviance # works for glm

Pero el enfoque anterior no funciona para Pseudo fuera de muestraR2

  1. Utilice la función "logLik" en R y definición (también funciona para muestras)

    mod_null <- glm(y~1, family = binomial, data = insample) 1- logLik(mod)/logLik(mod_null)

Esto puede modificarse ligeramente para calcular Pseudo fuera de muestraR2

Ejemplo:

pseudo-R fuera de muestra

Por lo general, el pseudo- fuera de muestra se calcula como donde es el probabilidad de registro para el período fuera de la muestra basado en los coeficientes estimados del período dentro de la muestra, mientras que es la probabilidad de registro para el modelo de solo intercepción para el período fuera de la muestra.R2Lest. outLnull. out

Rp2=1Lest.outLnull.out,
Lest.outLnull.out

Códigos:

pred.out.link <- predict(mod, outSample, type = "link") mod.out.null <- gam(Default~1, family = binomial, data = outSample) pR2.out <- 1 - sum(outSample$y * pred.out.link - log(1 + exp(pred.out.link))) / logLik(mod.out.null)

Xiaorui Zhu
fuente
deviance=2ln(Lfull) no se cumple para binomio, solo ver model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp, data = esoph, family = binomial)y llamar model1$deviancey -2*logLik(model1).
Curioso
6

si la desviación era proporcional a la probabilidad de registro, y uno usa la definición (ver, por ejemplo, McFadden aquí )

pseudo R^2 = 1 - L(model) / L(intercept)

entonces el pseudo- anterior sería = 0.7928 1 - 198,63R21198.63958.66

La pregunta es: ¿se informa la desviación proporcional a la probabilidad logarítmica?

dfrankow
fuente
3
Este pseudo-R ^ 2 no concuerda en absoluto con la respuesta Nagelkerke R ^ 2 de @ chl.
dfrankow
La desviación se definió como -2 * LL cuando estaba en la escuela.
DWin
@dfrankow no está de acuerdo, porque Nagelkerke es una normalización de Cox y Snell R2, que es diferente de McFaddens R2.
Colin
0

Si está fuera de la muestra , entonces creo que debe calcularse con las probabilidades de registro correspondientes como , donde es la probabilidad logarítmica de los datos de prueba con el modelo predictivo calibrado en el conjunto de entrenamiento, y es la probabilidad logarítmica de los datos de prueba con un modelo con solo una constante ajustada en el conjunto de entrenamiento, y luego usar el constante para predecir en el conjunto de pruebas calculando las probabilidades y, por lo tanto, obtener la probabilidad logarítmica.R2R2=1llfullllconstantllfullllconstant

Tenga en cuenta que en una regresión lineal, es análogo, el resultado de la muestra se calcula como , donde en particular si observamos el término denominador , la predicción usa el promedio sobre el conjunto de entrenamiento, . Esto es como si ajustamos un modelo en los datos de entrenamiento con solo una constante, por lo que tenemos que minimizar , lo que resulta en , entonces, este modelo predictivo constante simple es el que se usa como benchamrk (es decir, en el denominador de los oosR2R2=1i(yiy^i)2i(yiy¯train)2i(yiy¯train)2y¯traini(yiβ0)2 β 0= ¯ y trainR2Rβ^0=y¯trainR2plazo) para el cálculo de la muestra .R2

cthraves
fuente