Cómo calcular el pseudo- $R^2$ de la regresión logística de R?

Question

Christopher Manning's escritura sobre la 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)

Algunos resultados:

> 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

A continuación, entra en detalles sobre cómo interpretar los coeficientes, comparar diferentes modelos, etc. Bastante útil.

Sin embargo, ¿cuánta varianza tiene en cuenta el modelo? A Página de Stata sobre regresión logística dice:

Técnicamente, $R^2$ no puede calcularse de la misma manera en la regresión logística que en la regresión OLS. La pseudo- $R^2$ en la regresión logística, se define como $1 - \frac{L1}{L0}$ , donde $L0$ representa la probabilidad logarítmica para el modelo "sólo constante" y $L1$ es la probabilidad logarítmica para el modelo completo con la constante y los predictores.

Lo entiendo a alto nivel. El modelo de sólo constante sería sin ninguno de los parámetros (sólo el término de intercepción). La probabilidad logarítmica es una medida de cómo los parámetros se ajustan a los datos. De hecho, Manning insinúa que la desviación podría ser $-2 \log L$ . Tal vez la desviación nula sea sólo constante y la desviación residual sea $-2 \log L$ del modelo? Sin embargo, no lo tengo muy claro.

¿Puede alguien verificar cómo se calcula realmente el pseudo- $R^2$ en R utilizando este ejemplo?

Answer 1

No olvides el rms paquete, por Frank Harrell. Encontrará todo lo que necesita para ajustar y validar los MLG.

He aquí 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, utilizando el lrm función,

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

Pronto se obtienen muchos índices de ajuste del modelo, incluyendo el de Nagelkerke $R^2$ con print(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í, $R^2=0.445$ y se calcula como $\left(1-\exp(-\text{LR}/n)\right)/\left(1-\exp(-(-2L_0)/n)\right)$ donde LR es el $\chi^2$ (comparando los dos modelos anidados que has descrito), mientras que el denominador es sólo el valor máximo de $R^2$ . Para un modelo perfecto, esperaríamos $\text{LR}=2L_0$ Es decir $R^2=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 habló del uso de $R^2$ con GLM, en su libro Modelos de predicción clínica (Springer, 2009, § 4.2.2 pp. 58-60). Básicamente, la relación entre el estadístico LR y el de Nagelkerke $R^2$ es aproximadamente lineal (será más lineal con baja incidencia). Ahora bien, como se discutió en el hilo anterior que enlacé en mi comentario, se pueden utilizar otras medidas como el $c$ que es equivalente a la estadística AUC (también hay una bonita ilustración en la referencia anterior, véase la Figura 4.6).

Answer 2

Para obtener fácilmente un pseudo McFadden $R^2$ para un modelo ajustado en R, utilice el paquete "pscl" de Simon Jackman y use el comando pR2. http://cran.r-project.org/web/packages/pscl/index.html

Answer 3

Tenga cuidado con el cálculo de Pseudo- $R^2$ :

McFadden's Pseudo- $R^2$ se calcula como $R^2_M=1- \frac{ln\hat{L}_{full}}{ln\hat{L}_{null}}$ , donde $ln\hat{L}_{full}$ es la log-verosimilitud del modelo completo, y $ln\hat{L}_{full}$ es la probabilidad logarítmica del modelo con sólo el intercepto.

Dos enfoques para calcular el Pseudo- $R^2$ :

1. Utilizar la desviación: ya que $deviance = -2*ln(L_{full})$ , $null.deviance = -2*ln(L_{null})$

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

Pero el planteamiento anterior no funciona para los Pseudo $R^2$

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

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

Esto puede modificarse ligeramente para calcular el Pseudo $R^2$

Ejemplo:

pseudo-R fuera de la muestra

Por lo general, el pseudomuestreo fuera de la muestra $R^2$ se calcula como $$R_p^2=1−\frac{L_{est.out}}{L_{null.out}},$$ donde $L_{est.out}$ es la probabilidad logarítmica para el período fuera de la muestra basada en los coeficientes estimados del período dentro de la muestra, mientras que y $L_{null.out}$ es la probabilidad logarítmica del modelo de sólo intercepción para el período fuera de la muestra.

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)

Answer 4

si la desviación fuera proporcional al logaritmo de la probabilidad, y se utiliza la definición (véase, por ejemplo, McFadden's aquí )

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

entonces el pseudo- $R^2$ arriba sería $1 - \frac{198.63}{958.66}$ = 0.7928

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

Answer 5

Si su fuera de la muestra Entonces creo que el $R^2$ debe calcularse con las probabilidades logarítmicas correspondientes como $R^2=1-\frac{ll_{full}}{ll_{constant}}$ , donde $ll_{full}$ es la log-verosimilitud de los datos de prueba con el modelo predictivo calibrado en el conjunto de entrenamiento, y $ll_{constant}$ es la log-verosimilitud de los datos de prueba con un modelo con sólo una constante ajustada en el conjunto de entrenamiento, y luego utilizar la constante ajustada para predecir en el conjunto de prueba calculando las probabilidades y por lo tanto obtener la log-verosimilitud.

Nótese que en una regresión lineal, es análogo, el fuera de muestra $R^2$ se calcula como $R^2=1-\frac{\sum_{i}(y_{i}-\hat{y}_i)^2}{\sum_{i}(y_{i}-\overline{y}_{train})^2}$ donde, en particular, si nos fijamos en el término del denominador $\sum_{i}(y_{i}-\overline{y}_{train})^2$ la predicción utiliza la media del conjunto de entrenamiento, $\overline{y}_{train}$ . Esto es como si ajustáramos un modelo en los datos de entrenamiento con sólo una constante, por lo que tenemos que minimizar $\sum_{i}(y_i-\beta_0)^2$ que se traduce en $\hat{\beta}_0=\overline{y}_{train}$ entonces, este modelo de predicción simple y constante es el que se utiliza como benchamrk (es decir, en el denominador del oos $R^2$ ) para el cálculo de la muestra fuera de la muestra $R^2$ .