logistic.regression.fnc

Objetive

Lleva a cabo la estimación de modelos de regresión logística binaria a partir de las variables covariantes definidas por el usuario.

Logistic Regression (numeric variables)

Partiremos una base de datos simulada que se encuentra en el archivo binario de R logistic.Rdata. En este ejemplo tenemos una variable criterio: perfomance con dos niveles high y low y cuatro variables predictoras (covariantes) continuas: X2 a X5.

dat=read.file.fnc('logistic.Rdata')

*** The external file logistic.Rdata has been correctly read. ***

*** This is it head:

The object is a matrix that belong to class: data.frame

and has 100 rows and 5 columns (variables)

perfomance 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 perfomance con el objetivo de comprobar tanto el número de niveles como el orden que presentan en la base de datos leída.

frequency.fnc(dat, 'perfomance')

#------------------------------------------------------------------

# TABLE OF FRECUENCY

#------------------------------------------------------------------

$perfomance

$perfomance$n.total

[1] 100

$perfomance$table

high low

21 79

Observa que el primer nivel de rendimiento es alto y no bajo 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 reorder.factor.fnc.

dat=reorder.factor.fnc(dat, which.factor='perfomance',

new.levels=c('low','high'))

Estimaremos un modelo de regresión logística binaria utilizando como variable de agrupación la variable perfomance y el resto de las variables como posibles covariantes. Por defecto la función realizará una estimación paso a paso hacia adelante condicional (stepwise=T). Si incluimos el argumento graphic=T obtendremos en primer lugar las gráficas de residuos que nos permitirá evaluar la presencia de casos extremos o de palanca (leverages).

logistic.regression.fnc(dat, group='perfomance',

variables=2:5, graphic=T)

#------------------------------------------------------------------

# LOGISTIC REGRESSION

#------------------------------------------------------------------

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)

(Intercept) -3.3450 0.6856 -4.8786 0.0000 0.0353

X2 2.0760 0.5088 4.0799 0.0000 7.9725

X5 2.4602 0.6939 3.5452 0.0004 11.7072

X4 -1.1308 0.4429 -2.5532 0.0107 0.3228

$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 6.058e-08 ***

X5 1 75.512 81.512 20.0139 7.688e-06 ***

X4 1 62.805 68.805 7.3064 0.006871 **

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$clasificacion

Predicho

Observado bajo alto

bajo 75 4

alto 6 15

$prop.clasificacion

Predicho

Observado bajo alto

bajo 0.94936709 0.05063291

alto 0.28571429 0.71428571

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 multiplicativo de 8 y 11 aproximadamente sobre el riesgo de tener un bajo rendimiento. La variable X4 tiene un efecto protector sobre el riesgo (divide por 3 dicho riesgo (1/0.3).

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.

logistic.regression.fnc(dat[-c(28,60,83),], group='perfomance',

variables=2:5, graphic=T)

#------------------------------------------------------------------

# LOGISTIC REGRESSION

#------------------------------------------------------------------

Degrees of Freedom: 96 Total (i.e. Null); 93 Residual

Null Deviance: 93.07

Residual Deviance: 37.68 AIC: 45.68

$variables.seleccionadas

variables

1 X2

2 X5

3 X4

$ajuste

Chi.cuadrado gl p

55.3854 3.0000 0.0000

$coeficientes

B ET z p exp(B)

(Intercept) -4.9781 1.1866 -4.1951 0e+00 0.0069

X2 2.8567 0.7426 3.8469 1e-04 17.4040

X5 3.6254 1.0565 3.4315 6e-04 37.5397

X4 -1.6117 0.5866 -2.7476 6e-03 0.1995

Podemos ver que al eliminar esos 3 casos la deviance ha sido claramente disminuida. El impacto sobre los efectos de las variables predictoras es mas que evidente. Vemos que la unidad de cambio de X2 multiplica por 17 (mas del doble que antes) el riesgo de tener un bajo rendimiento, mientras que el efecto de X5 afecta multiplicando dicho riesgo por 37. Si observas ahora el efecto protector de X4 divide por 5 el riesgo.

Las 3 gráficas presentadas evidencian la característica distribución sigmoidad de la estimación logística. Podemos ver como para la variable X2 los valores cercanos a 2 sitúan a los sujetos en zona de probabilidad de bajo rendimiento, mientras que para X5 los valores deberán ser superiores a 1.5. Si observas la última gráfica verás el efecto protector de X4 que consiste precisamente en una disminución de la probabilidad de bajo rendimiento a medida que aumenta de valor dicha variables.

Logistic Regression (factor variables)

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.

dat=cut.variable.fnc(dat, variable='X4', ntiles=3)

#------------------------------------------------------------------

# CUT OF A NUMERIC VARIABLE

#------------------------------------------------------------------

$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

-----------------------------------------------------------------------------

*** The variable X4.r discretized by categories of the variable X4

*** has been created.

-----------------------------------------------------------------------------

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(dat)

perfomance 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

logistic.regression.fnc(dat[-c(28,60,83),], group='perfomance',

variables=c(2,3,5,6), graphic=T)

$coeficientes

B ET z p exp(B)

(Intercept) -3.3988 1.0834 -3.1372 0.0017 0.0334

X2 2.7784 0.7443 3.7326 0.0002 16.0933

X5 3.7600 1.1488 3.2729 0.0011 42.9484

X4.rP33-P66 -1.9707 1.2807 -1.5387 0.1239 0.1394

X4.r> P66 -4.3002 1.6397 -2.6225 0.0087 0.0136

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 protector 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)

Subir ->