06.-REGRESION

OBJETIVO

Lleva a cabo análisis de regresión múltiple por mínimos cuadrados ordinadorios (OLS) o robusta (M-estimator) para una variable criterio en p variables predictoras de forma simultánea o mediante el procedimiento de paso a paso o jerárquica.

1. Regresión simultánea y paso a paso

2. Regresión y casos extremos

3. Regresión robusta

4. Estimación Bootstrap

5. Regresión Jerárquica

6. Análisis de Dominancia

7. Regresión con Interacción

8. Regresión con modelo introducido por el usuario

9. Análisis de mediación

1. Regresión simultánea y paso a paso

El análisis de regresión se lleva a cabo mediante la función regresion.multiple.fnc. Para ejemplificar su uso partiremos de la base de datos zprostate de la librería instalada en el toolbox bestglm.

head(zprostate)

lcavol lweight age lbph svi lcp

1 -1.6373556 -2.0062118 -1.8624260 -1.024706 -0.5229409 -0.8631712

2 -1.9889805 -0.7220088 -0.7878962 -1.024706 -0.5229409 -0.8631712

3 -1.5788189 -2.1887840 1.3611634 -1.024706 -0.5229409 -0.8631712

4 -2.1669171 -0.8079939 -0.7878962 -1.024706 -0.5229409 -0.8631712

5 -0.5078745 -0.4588340 -0.2506313 -1.024706 -0.5229409 -0.8631712

6 -2.0361285 -0.9339546 -1.8624260 -1.024706 -0.5229409 -0.8631712

gleason pgg45 lpsa train

1 -1.0421573 -0.8644665 -0.4307829 TRUE

2 -1.0421573 -0.8644665 -0.1625189 TRUE

3 0.3426271 -0.1553481 -0.1625189 TRUE

4 -1.0421573 -0.8644665 -0.1625189 TRUE

5 -1.0421573 -0.8644665 0.3715636 TRUE

6 -1.0421573 -0.8644665 0.7654678 TRUE

Vemos que disponemos de 9 variables cuantitativas y una de caracter lógico (train). Se trata de un estudio realizado a 97 hombres con cancer de próstata a los que se ha medido el antígeno prostático PSA así como otras variables de carácter clínico (?zprostate te informará sobre esta base de datos y sus variables). Deseamos encontrar la mejor combinación lineal de variables que predigan la variable lpsa (logaritmo natural de la psa).

Podemos en primer lugar valorar el cumplimiento de relación lineal para dichas variables mediante la función grafica.dispersion.fnc.

grafica.dispersion.fnc(zprostate, variables=1:9)

Si observamos la última fila vemos que las variables con una relación lineal mas clara son: lcavol, lweigth y svi.

Seguidamente solicitaremos un modelo de regresión para la variable lpsa como criterio y el resto de las variables cuantitativas como predictoras.

regre1=regresion.multiple.fnc(zprostate, variables=1:9,

vd='lpsa', grafica=T)

Por defecto la función realizará una estimación paso a paso con el argumento paso.a.paso actuando como verdadero por defecto. Si deseas una estimación de todas las variables incluidas en el argumento variables, deberás incluir el argumento paso.a.paso=F

Al finalizar la estimación del modelo lineal obtendrás cuatro gráficas desde las cuales valorar el cumplimiento de los supuestos relativos a los residuales o errores del modelo . Así mismo podemos ver los casos que contribuyen al incremento del error en el modelo estimado. En este sentido podemos ver que los casos 39, 95 y 96 presentan los mayores valores de residual.

Por defecto la estimación de un modelo de regresión múltiple generá la siguiente lista de resultados:

1. Determinación de la transformación de potencia necesaria para promover la normalidad de la variable criterio.

2. Determinación de la transformación de potencia necesaria para promover varianza constante en la variable criterio.

3. Matriz de correlación de pearson de las variables predictoras y criterio.

4. Resultado de la estimación paso a paso del modelo.

5. Test de autocorrelación de los errores de Durbin-Watson.

6. Coeficientes no estandarizados y estandarizados del modelo.

7. Diagnóstico de colinealidad por tolerancia y factor de inflación de varianza (VIF).

8. ïndices de correlación parcial y semiparcial de las variables predictoras seleccionadas en el modelo.

9. Indice de determinación puntual y ajustado.

10. Anova de la varianza explicada por el modelo estimado.

2. Regresión y casos extremos

Tal y como hemos podido ver en el ejemplo anterior tenemos una serie de casos que contribuyen negativamente al ajuste del modelo. Por ese motivo vamos a re-estimarlo, eliminando los casos que en las gráficas de residuales anteriores hemos detectado como extremos incluyendo el argumento elimina.extremos con el número de los casos que queremos eliminar del análisis (pero no de la base de datos)

regresion.multiple.fnc(zprostate, variables=1:8,

vd='lpsa', elimina.extremos=c(81,95,96,39))

3. Regresión robusta

Los mínimos cuadrados ordinarios son muy sensibles a la presencia de casos extremos multivariados en los datos. En el apartado anterior vimos una posible solución consistente en eliminar dichos casos. Otro acercamiento al problema se denomina regresión robusta donde se emplea un criterio de ajuste no tan vulnerable a los casos extremos. El procedimiento mas común es el introcido por Huber denominado estimador M. Lo activaremos en la función regresion.multiple.fnc incluyendo el argumento robusta=T.

regre3=regresion.multiple.fnc(zprostate, variables=1:8,

vd='lpsa', paso.a.paso=T, robusta=T, grafica=T)

Esto da lugar a la estimación de pesos donde los casos con peso inferior a la unidad son presentados al usuari@

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

# ANALISIS DE REGRESION - EST. ROBUSTA (estimador M de Huber)

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

*** Casos con ponderacion (w) inferior a 1 ***

lpsa lcavol lweight svi w indice

5 0.3715636 -0.5078745 -0.45883405 -0.5229409 0.5459589 5

47 2.5687881 1.1690261 0.85549191 1.8925480 0.6009827 47

38 2.1916535 -0.7573103 -2.92717970 -0.5229409 0.6049781 38

97 5.5829322 1.8003666 0.80776440 1.8925480 0.6173111 97

57 2.7880929 -0.3185491 -1.78307338 -0.5229409 0.6321376 57

18 1.4929041 0.7962471 0.04765594 -0.5229409 0.6341435 18

7 0.7654678 -0.5199666 -0.36279315 -0.5229409 0.7236293 7

1 -0.4307829 -1.6373556 -2.00621178 -0.5229409 0.7350848 1

85 3.5876769 0.1801563 0.18898744 -0.5229409 0.7369823 85

2 -0.1625189 -1.9889805 -0.72200876 -0.5229409 0.7407559 2

8 0.8544153 -0.5573125 -0.20875657 -0.5229409 0.7608354 8

72 3.0373539 -0.1611952 -0.67190036 -0.5229409 0.7731158 72

4 -0.1625189 -2.1669171 -0.80799390 -0.5229409 0.8406158 4

87 3.6800909 0.5720085 0.23985445 -0.5229409 0.8491387 87

22 1.6582281 0.6017430 -0.29854413 -0.5229409 0.9186895 22

90 3.9936030 0.1801563 0.15444819 1.8925480 0.9317534 90

69 2.9626924 -1.5240614 1.81975701 -0.5229409 0.9496610 69

4.- 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, 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))

regre1=regresion.multiple.fnc(zprostate, variables=1:9,

vd='lpsa', boot=T)

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

# Extrayendo 1000 muestras bootstrap

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

$Bootstrap

Original Sesgo ET.boot t.boot IC.BCa.inf IC.BCa.sup

(Intercept) 2.47839 -0.002697595 0.074035 33.4758 2.340560 2.62906

lcavol 0.61978 0.000045504 0.086693 7.1491 0.453712 0.79317

lweight 0.28351 0.001545922 0.078142 3.6281 0.130208 0.43789

svi 0.27558 0.000376978 0.092892 2.9667 0.085199 0.45219

En la tabla de resultados Bootstrap podemos ver en primer lugar los valores originales de la estimación, seguido del sesgo (diferencia entre el valor estimado y el promedio de las 1ooo 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).

5. Regresión Jerárquica

El usuario puede solicitar un modelo de regresión jerárquica donde se define la entrada según el orden definido por el usuario en el nuevo argumento jerarquica. Esta metodología permitirá valorar el efecto que la primera variable predictora tiene sobre la criterio obviando cualquier covarianza con las restantes. De esta forma el usuario obtener una aproximación al Anova con suma de cuadrados tipo I, modificando el orden de entrada de las variables en el modelo de regresión.

regresion.multiple.fnc(zprostate, vd='lpsa',

elimina.extremos=c(5,38,39,95,96),

jerarquica=c('lcavol','lweight','svi'))

El procedimiento de regresión jerárquica ha sido sustituida por el análisis de regresión de dominancia, que puedes ver en el siguiente apartado, dado que aporta mucha más información y de mayor calidad acerca de la importancia que cada variable tiene en términos de varianza explicada.

6. Análisis de Dominancia

Si se desea estimar la importancia relativa de los regresores del modelo lineal estimado Groemping, U. 2006, utilizaremos el argumento dominancia=T. Esto dará lugar a una descomposición de la varianza total explicada a partir de las variables incluídas en el modelo finalmente estimado. Si se incluyese el argumento ajustada.por se entiende que se desea descomponer dicha varianza descontando previamente lo que estas variables incluídas en ese argumento explican del índice de determinación.

regresion.multiple.fnc(zprostate, variables=1:9, vd='lpsa', dominancia=T)

lmg puede verse como el promedio ponderado de la varianza explicada por cada variable utilizando para ello tanto los efectos directos así como los ajustados para todas las variables del modelo.

pmvd a diferencia del estadístico anterior elimina la posibilidad de que una variable con peso igual a cero pueda presentar una varianza promediada positiva si esta presenta una correlación importante con las otras variables predictoras pero no con la criterio (supresión cooperativa).

last te indica el valor de la varianza cuando esa variable es la última en entrar al modelo de regresión. De igual modo first nos informa del valor de varianza explicada cuando esa es la primera variable en entrar al modelo.

Si incluímos el argumento ajustada.por con el nombre de las variables del modelo estimado sobre las que deseamos ajustar la varianza residual (varianza que queda por explicar una vez controlado el efecto de las variables incluídas en ajustada.por).

regresion.multiple.fnc(zprostate, variables=1:9,

vd='lpsa', dominancia=T,

ajustada.por=c('age','lweight'))

$Dominancia

Response variable: lpsa

Total response variance: 1.332476

Analysis based on 97 observations

3 Regressors:

Proportion of variance explained: 63.59%

One Regressor always included in model:

lweight

18.78 % of variance explained by this regressor

Relative importance of 2 regressors assessed:

lcavol svi

44.82 % of variance decomposed among these

Metrics are not normalized (rela=FALSE).

Relative importance metrics:

lmg pmvd

lcavol 0.3004382 0.37057959

svi 0.1477460 0.07760463

Podemos ver que la variable lweight (age no entró en el modelo estimado paso a paso) explica el 18.8% de la varianza de la criterio. Las restantes variables explican en 30 y 15% respectivamente según lmg y el 37 y 7 según pmvd una vez controlado su efecto por la variable lweight.

lmg puede verse como el promedio ponderado de la varianza explicada por cada variable utilizando para ello tanto los efectos directos así como los ajustados para todas las variables del modelo.

pmvd a diferencia del estadístico anterior elimina la posibilidad de que una variable con peso igual a cero pueda presentar una varianza promediada positiva si esta presenta una correlación importante con las otras variables predictoras pero no con la criterio (supresión cooperativa).

last te indica el valor de la varianza cuando esa variable es la última en entrar al modelo de regresión. De igual modo first nos informa del valor de varianza explicada cuando esa es la primera variable en entrar al modelo.

Si incluímos el argumento ajustada.por con el nombre de las variables del modelo estimado sobre las que deseamos ajustar la varianza residual (varianza que queda por explicar una vez controlado el efecto de las variables incluídas en ajustada.por).

regresion.multiple.fnc(zprostate, variables=1:9,

vd='lpsa', dominancia=T,

ajustada.por=c('age','lweight'))

$Dominancia

Response variable: lpsa

Total response variance: 1.332476

Analysis based on 97 observations

3 Regressors:

Proportion of variance explained: 63.59%

One Regressor always included in model:

lweight

18.78 % of variance explained by this regressor

Relative importance of 2 regressors assessed:

lcavol svi

44.82 % of variance decomposed among these

Metrics are not normalized (rela=FALSE).

Relative importance metrics:

lmg pmvd

lcavol 0.3004382 0.37057959

svi 0.1477460 0.07760463

Podemos ver que la variable lweight (age no entró en el modelo estimado paso a paso) explica el 18.8% de la varianza de la criterio. Las restantes variables explican en 30 y 15% respectivamente según lmg y el 37 y 7 según pmvd una vez controlado su efecto por la variable lweight.

7. 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.

Para el análisis de la interacción la función hace uso de la librería sjPlot.

Utilizaremos la base de datos swiss que contiene 6 variables medidas en 47 observaciones relativas a la tasa de fertilidad y su relación con 5 variables características de 47 poblaciones (agricultura, educacion, catolicismo, etc).

head(swiss)

Fertility Agriculture Examination Education Catholic Infant.Mortality

Courtelary 80.2 17.0 15 12 9.96 22.2

Delemont 83.1 45.1 6 9 84.84 22.2

Franches-Mnt 92.5 39.7 5 5 93.40 20.2

Moutier 85.8 36.5 12 7 33.77 20.3

Neuveville 76.9 43.5 17 15 5.16 20.6

Porrentruy 76.1 35.3 9 7 90.57 26.6

Debemos ahora centrar las variables que van a entrar como predictoras en el modelo (todas). Para ello utilizaremos la función centra.variable.fnc conjuntamente con un bucle for (para hacerlo todo en un solo paso).

for(i in 2:6) swiss = centra.variable.fnc(swiss, variable=i)

nombres.var(swiss)

nombre numero tipo

1 Fertility 1 numeric

2 Agriculture 2 numeric

3 Examination 3 integer

4 Education 4 integer

5 Catholic 5 numeric

6 Infant.Mortality 6 numeric

7 c.Agriculture 7 numeric

8 c.Examination 8 numeric

9 c.Education 9 numeric

10 c.Catholic 10 numeric

11 c.Infant.Mortality 11 numeric

regresion.multiple.fnc(swiss, vd='Fertility', variables=7:11, interaccion=T)

8. Modelo definido por el usuario

En ocasiones el modelo lineal por defecto no es el adecuado. Pensemos en modelos no lineales aunque intrínsecamente lineales como los polinomiales, logarítmicos, exponenciales, inversos, etc. Asimismo podemos desear estimar un modelo con interacción entre dos o mas variables cualitativas, cuantitativas o mixtas. Para cubrir ese rango de modelos mas inusuales existe el argumento modelo donde el usuario definirá el vínculo que existe entre la variable criterio y las predictoras. Para demostrar su uso simularemos algunas matrices de datos que nos permitirán exponer el modo correcto de utilizarlo.

Crearemos una relación no lineal cuadrática entre una variable con distribución normal x y otra y que resulta de aplicar el modelo que puedes observar en las líneas siguientes.

x=rnorm(100)

y=20 + 0.8*x^2 + runif(100)

datos=data.frame(y,x)

head(datos)

y x

1 20.759 0.53510

2 22.265 1.27190

3 20.518 -0.28641

4 20.126 0.23650

5 20.689 -0.30591

6 20.950 0.68272

Podemos estimar inicialmente ese modelo de forma clásica, aunque ya vemos en el diagrama de dispersión que se solicita a continuación que no existe relación lineal entre ambas variables.

grafica.dispersion.fnc(datos, variables=1:2)

regresion.multiple.fnc(datos, vd='y', grafica=T)

$Coeficientes

Beta Estimate Std. Error t value Pr(>|t|) ic.inf ic.sup

(Intercept) 0.000 21.2059 0.1051 201.6865 0.0000 20.9973 21.4145

x 0.237 0.2607 0.1081 2.4126 0.0177 0.0462 0.4752

$Indices

R R2 R2.adj

0.236776 0.056063 0.046431

Tal y como podemos ver en la figura de la derecha, los residuales no tienen la distribución aleatoria esperada. Por el contrario se observa una clara relación no lineal y cuadrática entre la variable predictora y la criterio tal y como el diagrama de dispersión nos mostraba.

Seguidamente estimaremos un modelo tentativo para esa relación haciendo uso del argumento modelo. Observa como en la definición de ese modelo utilizamos la funcion I( ) la cual puede traducirse de forma literal como "tal cual". Esta función es obligatoria cuando queremos indicar que la variable debe ser tratada tal y como se indica dentro de la función I( ).

Todos los modelos deber tener a la variable criterio y las predictoras separados por el caracter ~, el cual puede obtenerse en general mediante la combinación de teclas AltGr + 4.

modelo = 'y ~ x + I(x^2)'

regresion.multiple.fnc(datos, modelo=modelo, grafica=T)

$Coeficientes

Beta Estimate Std. Error t value Pr(>|t|) ic.inf ic.sup

(Intercept) 0.0000 20.5341 0.0366 561.6473 0.0000 20.4615 20.6067

x 0.0054 0.0059 0.0321 0.1844 0.8541 -0.0578 0.0696

I(x^2) 0.9591 0.7793 0.0237 32.8860 0.0000 0.7323 0.8263

$Indices

R R2 R2.adj

0.96037 0.92231 0.9207

Podemos ver ahora en la figura de la derecha, que nuestro modelo ha mejorado considerablemente. En primer lugar hemos incrementado la varianza explicada desde el exiguo 5% al 92%. Por otra parte los residuales que en el modelo lineal anterior no solo no tenían un patrón aleatorio sino que daba lugar a valores muy alejados del 0. Ahora puedes ver como estos han perdido el patrón no lineal así como se mantienen en valores muy bajos como corresponde a la clara mejoría de este modelo evidenciada en la varianza explicada en el índice de determinación.

Ahora plantearemos una serie de modelos no lineales sobre la base de datos disponible por defecto attitude que el usuario puede probar a estimar. Para ello crearemos artificialmente dos variables cualitativas en esa base de datos.

attitude$sexo=c('Hombre','Mujer')

attitude$zona=rep(c('A','B'), each=15)

head(attitude)

rating complaints privileges learning raises critical advance sexo zona

1 43 51 30 39 61 92 45 Hombre A

2 63 64 51 54 63 73 47 Mujer A

3 71 70 68 69 76 86 48 Hombre A

4 61 63 45 47 54 84 35 Mujer A

5 81 78 56 66 71 83 47 Hombre A

6 43 55 49 44 54 49 34 Mujer A

A continuación presentaremos múltiples modelos alternativos, creados exclusivamente por motivos didácticos y no guiadas por la naturaleza real de los datos. El usuario puede a partir de estos modelos estimarlos y comprobar el uso del argumento modelo.

modelo='rating ~ sexo*zona + complaints + I(1/privileges)'

modelo='rating ~ complaints + I(complaints^2) + privileges'

modelo='rating ~ complaints + log(critical) + sqrt(privileges)'

regresion.multiple.fnc(datos=attitude, modelo=modelo,

grafica=T, robusta=T)

Subir->