regresion.logistica.fnc
Objetivo
Lleva a cabo la estimación de modelos de regresión logística binaria y multinomial (J categorías) a partir de las variables covariantes definidas por el usuario.
1.- Regresión Logística (V. cuantitativas)
Partiremos una base de datos simulada que se encuentra en el archivo binario de R logistica.Rdata. En este ejemplo tenemos una variable criterio: rendimiento con dos niveles alto y bajo y cuatro variables predictoras (covariantes) continuas: X2 a X5.
datos=lee.archivo.fnc('logistica.Rdata')
*** Se ha leido correctamente el archivo logistica.Rdata
*** Esta es la cabecera:
El objeto es una matriz que pertenece a la clase data.frame
y tiene 100 filas y 5 columnas (variables)
rendimiento X2 X3 X4 X5
1 low -1.2431936 -0.05121179 -0.9686906 -2.10458829
2 low -2.4816524 -0.68465902 -0.3790959 -1.25686392
3 low 1.1117752 1.74769268 2.2902208 1.10005473
4 low 0.5235463 1.18290975 0.6442378 -0.43535722
5 low 1.0170960 0.50613817 0.2914288 -0.02599073
6 low -0.9165170 -1.64128257 -0.8047344 -2.10915093
Solicitamos la distribución de frecuencias de la variable rendimiento con el objetivo de comprobar tanto el número de niveles como el orden que presentan en la base de datos leída.
frecuencias.fnc(datos, variables='rendimiento')
#------------------------------------------------------------------
# TABLA DE FRECUENCIAS
#------------------------------------------------------------------
$rendimiento
$rendimiento$n.total
[1] 100
$rendimiento$tabla
low hig
21 79
Observa que el primer nivel de rendimiento es low (bajo) y no hig (alto) como esperamos (por defecto los niveles se ordenan alfabéticamente) ya que deseamos modelar el riesgo de bajo rendimiento y no el de alto. Para solucionarlo reordenamos los niveles de la variable mediante la función reordena.factor.fnc.
datos=reordena.factor.fnc(datos, que.factor='rendimiento', niveles=c('hig','low'))
A partir de este momento la variable rendimiento tendrá como primer nivel de la misma el nivel hig y como segundo low. Esto hará que el modelo que estimamos lo sea sobre el riesgo de estar en el nivel low (bajo) de dicha variable.
Estimaremos un modelo de regresión logística binaria utilizando como variable de agrupación (argumento grupo) la variable rendimiento y el resto de las variables como posibles covariantes. Por defecto la función realizará una estimación paso a paso hacia adelante condicional (paso.a.paso=T). Si incluimos el argumento grafica=T obtendremos en primer lugar las gráficas de residuos que nos permitirá evaluar la presencia de casos extremos o de palanca (leverages).
regresion.logistica.fnc(datos, grupo='rendimiento', variables=2:5, grafica=T)
#------------------------------------------------------------------
# REGRESION LOGISTICA
#------------------------------------------------------------------
Degrees of Freedom: 99 Total (i.e. Null); 96 Residual
Null Deviance: 102.8
Residual Deviance: 55.5 AIC: 63.5
$variables.seleccionadas
variables
1 X2
2 X5
3 X4
$ajuste
Chi.cuadrado gl p
47.2929 3.0000 0.0000
$coeficientes
B ET z p exp(B) ic.inf ic.sup
(Intercept) 3.3450 0.6856 4.8786 0.0000 28.36 7.40 108.72
X2 -2.0760 0.5088 -4.0799 0.0000 0.13 0.05 0.34
X4 1.1308 0.4429 2.5532 0.0107 3.10 1.30 7.38
X5 -2.4602 0.6939 -3.5452 0.0004 0.09 0.02 0.33
$cada.predictora
Single term deletions
Model:
rendimiento ~ X2 + X5 + X4
Df Deviance AIC LRT Pr(>Chi)
<none> 55.498 63.498
X2 1 84.843 90.843 29.3449 0.00000006058 ***
X4 1 62.805 68.805 7.3064 0.006871 **
X5 1 75.512 81.512 20.0139 0.00000768826 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
$clasificacion
Predicho
Observado hig low
hig 15 6
low 4 75
$prop.clasificacion
Predicho
Observado hig low
hig 0.71428571 0.28571429
low 0.05063291 0.94936709
La primera columna de coeficientes tienen explicación inicialmente poco clara porque en realidad predicen el logaritmo de la odd ratio (P/1-P). Precisamente, una forma de eliminar el logaritmo consiste en exponenciar ambos lados de la ecuación con lo que la última columna puede ser interpretada como el efecto multiplicativo que la unidad de cambio de la variable asociada tiene sobre la odd ratio (riesgo) y no sobre su logaritmo. En este sentido, vemos que las variables X2 y X5 tienen un efecto protector sobre el riesgo (exponente inferior a 1) de 0.13 y 0.09 respectivamente sobre el riesgo de tener un bajo rendimiento. La variable X4 por el contrario, tiene un efecto multiplicativo e incremental del riesgo de 3.1 .
Sin embargo, tal y como podemos ver en la primera figura hay claramente 3 residuales tipificados superiores a 2: 83, 60 y 28. Repetiremos el análisis eliminando estos 3 casos del análisis.
regresion.logistica.fnc(datos[-c(28,15,85),], grupo='rendimiento',
variables=2:5, grafica=T)
#------------------------------------------------------------------
# REGRESION LOGISTICA
#------------------------------------------------------------------
Degrees of Freedom: 96 Total (i.e. Null); 93 Residual
Null Deviance: 98.72
Residual Deviance: 51.23 AIC: 59.23
$variables.seleccionadas
variables
1 X2
2 X5
3 X4
$ajuste
Chi.cuadrado gl p
47.4836 3.0000 0.0000
$coeficientes
B ET z p exp(B)
(Intercept) 3.4630 0.7270 4.7632 0.0000 31.91 7.68 132.68
X2 -2.0542 0.5158 -3.9826 0.0001 0.13 0.05 0.35
X4 1.1362 0.4516 2.5160 0.0119 3.11 1.29 7.55
X5 -2.6338 0.7393 -3.5624 0.0004 0.07 0.02 0.31
Podemos ver que al eliminar esos 3 casos la deviance ha sido ligeramente disminuida aunque con escaso impacto sobre los parámetros estimados.
Las 3 figuras presentadas evidencian la característica distribución sigmoidad de la estimación logística. Podemos ver como para la variable X2 los valores inferiores a 2 sitúan a los sujetos en zona de probabilidad de bajo rendimiento, mientras que para X5 los valores deberán ser inferiores a 1. Como puedes ver los efectos de X2 y X5 son protectores dado que a mayor valor esas variables menor es el riesgo de rendimiento bajo. No ocurre así con la variable X4 que como puedes ver tiene un efecto multiplicativo del riesgo.
2.- Estimación Bootstrap
Si incluimos el argumento bootstrap=T, la función llevará a cabo la extracción de 1000 muestras bootstrap para calcular el sesgo en la estimación de los coeficientes de la regresión logística, así como el error típico bootstrap el cual se utilizará para estimar el estadístico "t" así como el intervalo de confianza con corrección de sesgo (adjusted bootstrap percentile BCa method) Efron (1987))
regresion.logistica.fnc(datos, grupo='rendimiento', variables=2:5, bootstrap=T)
#------------------------------------------------------------------
# Extrayendo 1000 muestras bootstrap
#------------------------------------------------------------------
$Bootstrap
Original Sesgo ET.boot t.boot IC.BCa.inf IC.BCa.sup
(Intercept) 3.344973 0.3094383 1.2612623 2.652084 2.1433701 5.042513
X2 -2.076020 -0.2384329 0.6592583 -3.149023 -3.0507864 -1.020164
X4 1.130758 0.1310681 0.6246963 1.810092 0.1626953 2.143294
X5 -2.460152 -0.2800133 1.3109993 -1.876547 -4.0415771 -1.290363
En la tabla de resultados Bootstrap podemos ver en primer lugar los valores originales de la estimación de los coeficientes de la regresión logística, seguido del sesgo (diferencia entre el valor estimado y el promedio de las 1000 replicaciones bootstrap), el error típico boostrap (desviación típica a partir de las 1000 muestras), el estadístico "t" (valor del parámetro dividido por el error típico bootstrap) y el intervalo de confianza corregido por el sesgo de asimetría (Efron 1987).
3.- Regresión Logística (V. cualitativas)
Para demostrar el uso de la regresión logística con variables cualitativas crearemos una versión discretizada percentílica (3 niveles) de la variable con efecto protector sobre el riesgo de bajo rendimiento X4.
datos=discretiza.variable.fnc(datos, variable='X4', ntiles=3)
#------------------------------------------------------------------
# DISCRETIZADO DE VARIABLE CONTINUA
#------------------------------------------------------------------
$resumen
Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.4980 -0.3903 0.3279 0.3032 1.0080 2.5900
$metodo
[1] "Cortes Percentilicos"
$tabla.n
< P33 P33-P66 > P66
Freq. 33 33 34
$tabla.medias
< P33 P33-P66 > P66
Media -0.776708 0.260619 1.392602
-----------------------------------------------------------------------------
*** Se ha creado la varible X4.r discretizado de la variable X4 ***
----------------------------------------------------------------------------
Ahora la variable X4.r tiene tres niveles <P33 (menor o igual al percentil 33), P33-P66 y >P66.
Dado que vamos a utilizar la variable discretizada, obviamente no incluiremos en el modelo la variable original X4.
head(datos)
rendimiento X2 X3 X4 X5 X4.r
1 low -1.2431936 -0.05121179 -0.9686906 -2.10458829 < P33
2 low -2.4816524 -0.68465902 -0.3790959 -1.25686392 < P33
3 low 1.1117752 1.74769268 2.2902208 1.10005473 > P66
4 low 0.5235463 1.18290975 0.6442378 -0.43535722 P33-P66
5 low 1.0170960 0.50613817 0.2914288 -0.02599073 P33-P66
6 low -0.9165170 -1.64128257 -0.8047344 -2.10915093 < P33
regresion.logistica.fnc(datos, grupo='rendimiento', paso.a.paso=F,
variables=c(2,5,6), grafica=T)
$coeficientes
B ET z p exp(B) ic.inf ic.sup
(Intercept) 2.2850 0.7281 3.1383 0.0017 9.83 2.36 40.94
X2 -2.0138 0.5053 -3.9856 0.0001 0.13 0.05 0.36
X5 -2.5048 0.7274 -3.4434 0.0006 0.08 0.02 0.34
X4.rP33-P66 1.2906 0.9765 1.3216 0.1863 3.63 0.54 24.64
X4.r> P66 2.9108 1.1832 2.4602 0.0139 18.37 1.81 186.77
Vemos que se han creado automáticamente dos variables dummys : X4.rP33-P66 y X4.r> P66 que miden el efecto del cambio del nivel de referencia (<P33) al nivel definido en cada dummy. En este sentido podemos ver que el auténtico efecto multiplicador del riesgo, es el que se produce al cambiar desde el percentil mas bajo al mas alto (> P66 p < 0.05) y no desde el mas bajo al medio (P33-P66 p > 0.05)
4.- Análisis de Interacción
Mediante la inclusión del argumento interaccion=T, la función lleva a cabo una estimación paso a paso de todas las combinaciones dos a dos de las variables incluidas en el modelo. Es muy importante (y así te lo indicará la propia función) que las variables cuantitativas entren centradas en el modelo (para dar significado al término Bo). En caso contrario se producirían estimaciones de la variable criterio fuera de los límites de los valores reales de la variable.
for(i in 2:5) datos=centra.variable.fnc(datos, variable=i)
nombres.var(datos)
nombre numero tipo
1 rendimiento 1 factor
2 X2 2 numeric
3 X3 3 numeric
4 X4 4 numeric
5 X5 5 numeric
6 X4.r 6 factor
7 c.X2 7 numeric
8 c.X3 8 numeric
9 c.X4 9 numeric
10 c.X5 10 numeric
regresion.logistica.fnc(datos, grupo='rendimiento',
variables=7:10, interaccion=T)
$Anova
Analysis of Deviance Table (Type II tests)
Response: rendimiento
LR Chisq Df Pr(>Chisq)
c.X2 18.9314 1 0.0000135505 ***
c.X3 1.3947 1 0.237608
c.X4 8.1846 1 0.004225 **
c.X5 25.1984 1 0.0000005173 ***
c.X2:c.X3 2.9161 1 0.087700 .
c.X2:c.X4 0.1031 1 0.748198
c.X2:c.X5 0.6214 1 0.430520
c.X3:c.X4 6.6576 1 0.009873 **
c.X3:c.X5 3.0991 1 0.078337 .
c.X4:c.X5 1.9977 1 0.157539
---
$Graficas.Interaccion
"*** Se ha creado el archivo graf_interac_reg_log.pdf con las gráficas"
"*** de las interacciones que han resultado significativas."
