Para un conjunto de \(p-1\) variables regresoras, se pueden construir \(2^{p-1}\) modelos. En la mayoría de circunstancias es imposible para un analista hacer una revisión detallada de todos los modelos posibles. Por ejemplo, si se tienen \(10\) variables regresoras, se podrían construir \(2^{10}=1024\) modelos. En este sentido la utilización de software estadístico toma importancia.
El procedimiento de selección de modelos, también conocido como selección de variables, se han desarrollado para identificar un grupo de modelos de regresión que son buenos según un criterio o métrica
Muchos criterios han sido desarrollados para comparar modelos de regresión. Algunos de ellos serán vistos en esta sección del curso.
Notación.
\(P-1\): número total de variables regresoras \(X\).
Se asume que los modelos a comparar poseen término \(\beta_0\). Por lo tanto, la función de regresión contiene a todas las variables regresoras potencias \(X\) y los \(P\) parámetros.
El número de variables regresoras \(X\) en un submodelo será denotado por \(p-1\). De modo que hay \(p\) parámetros en la función de regresión para un submodelo con un subconjunto de variables \(X\). Esto es:
\[\begin{align}
1 \leq p \leq P
\end{align} \tag{1}\]
Se asume que el número de observaciones excede el máximo número de parámetros \(\beta\) potenciales.
\[\begin{align}
n>p
\end{align} \tag{2}\]
De hecho, es deseable que \(n\) sea mucho mayor que \(p\) para obtener resultados sólidos.
Criterio \(R_p^2\)
El criterio \(R_p^2\) usa el coeficiente de determinación múltiple, que se muestra en la Ecuación 3:
\(R_p^2\) indica que hay \(p\) parámetros o \(p-1\) variables \(X\).
El criterio \(R_p^2\) no pretende identificar los subconjuntos que maximizan este criterio. \(R_p^2\) nunca puede disminuir a medida que se incluyen variables \(X\) adicionales en el modelo. Por lo tanto, \(R_p^2\) será máximo cuando todas las \(P - 1\) variables \(X\) potenciales se incluyan en el modelo de regresión. La intención al utilizar el criterio \(R_p^2\) es encontrar el punto en el que añadir más variables \(X\) no merece la pena porque conduce a un aumento muy pequeño de \(R_p^2\). A menudo, este punto se alcanza cuando sólo se incluye un número limitado de variables \(X\) en el modelo de regresión. Evidentemente, la determinación del punto en el que aparecen los rendimientos decrecientes es una cuestión de juicio.
Call:
lm(formula = datos$ln_tiempo_supervivencia ~ datos$coagulacion_score +
datos$indice_pronostico + datos$test_enzimas)
Residuals:
Min 1Q Median 3Q Max
-0.46994 -0.17938 -0.03116 0.17959 0.59105
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.766441 0.226757 16.610 < 2e-16 ***
datos$coagulacion_score 0.095475 0.021692 4.401 5.66e-05 ***
datos$indice_pronostico 0.013344 0.002035 6.558 2.95e-08 ***
datos$test_enzimas 0.016444 0.001630 10.089 1.19e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2493 on 50 degrees of freedom
Multiple R-squared: 0.7572, Adjusted R-squared: 0.7427
F-statistic: 51.99 on 3 and 50 DF, p-value: 2.137e-15
Criterio \(R^2_p~ajustado\)
Como el \(R_p^2\) no tiene en cuenta el número de parámetros del modelo de regresión y el \(max(R_p^2)\) no decrece cuando la cantidad de parámetros \(p\) incrementa, el coeficiente de determinación ajustado, \(R^2_{ap}\) (Ecuación 5), resulta de mucha utilidad
Donde \(SSE_p\) es la suma de cuadrados del error para el submodelo ajustado con \(p\) parámetros.
Valores de \(C_p\) cercanos a \(p\) son interpretados como que no muestran sesgo.
En el contexto de la regresión lineal el sesgo En el contexto de modelos de regresión y otros modelos predictivos, el sesgo es el error que proviene de la simplificación excesiva del modelo, como cuando no se incluyen variables importantes. Por ejemplo, si se usa un modelo lineal para datos que en realidad tienen una relación no lineal, el modelo tendrá un sesgo alto.
Si se usa el criterio \(C_p\) para selección de variables y modelos se debe observar:
El valor de \(C_p\) es bajo.
\(C_p\) es cercano a \(p\).
Ejemplo de cálculo de \(C_p\) de Mallow.
#Datoslibrary(readxl)datos <-read_excel("data_cirugia.xlsx")datos$genero <-as.factor(datos$genero)datos$uso_alcohol_moderado <-as.factor(datos$uso_alcohol_moderado)datos$uso_alcohol_exceso <-as.factor(datos$uso_alcohol_exceso)#Modelo completo con X1, X2, X3, X4, p=5 modelo_completo <-lm(datos$ln_tiempo_supervivencia~datos$coagulacion_score+datos$indice_pronostico+datos$test_enzimas+datos$test_higado)#Modelo que solo tiene X4modelo <-lm(datos$ln_tiempo_supervivencia~datos$test_higado)# MSE del modelo completoanova_completo <-aov(modelo_completo)resumen_completo <-summary(anova_completo)MSE <- resumen_completo[[1]]$`Mean Sq`[5]# SSE modelo con X4anova <-aov(modelo)resumen <-summary(anova)SSE <- resumen[[1]]$`Sum Sq`[2]# Cpn <-nrow(datos)p<-length(coef(modelo))c_p <- (SSE/MSE)-(n-(2*p))c_p
[1] 67.69589
#Con sentencias de Rlibrary(olsrr)resultados <-ols_mallows_cp(modelo,modelo_completo)# Cp con todos los posibles modeloslibrary(olsrr)resultados_todos <-ols_step_all_possible(modelo_completo)cp_todos<-c(resultados_todos$result[3],resultados_todos$result[8])print(cp_todos)
Se observa que \(C_p = 3.387945< p=4\), para ese modelo con \(X1, X2~y~ X3\), la relación indica un sesgo uy bajo (o no sesgo).
Criterio \(AIC_p\)
El criterio Akaike’s Information Criterion\(AIC_p\), que se calcula como se muestra en la Ecuación 7, penaliza a los modelos con muchas variables regresoras (como el \(R^2_{ap}\) y \(C_p\)).
Note que el primer término \(n~ln(SSE_p)\)decrece cuando \(p\) crece, el segundo término \(n~ln(n)\) es fijo para un valor de \(n\), el tercer término \(2p\), aumenta con el número de parámetros \(p\).
En este sentido, los modelos con \(SSE_p\) bajos tienen buen desempeño en este criterio, siempre que la penalización \(2p\) no sea tan grande.
Se escogen los modelos con menor \(AIC_p\).
Ejemplo de cálculo de \(AIC_p\).
#Datoslibrary(readxl)datos <-read_excel("data_cirugia.xlsx")datos$genero <-as.factor(datos$genero)datos$uso_alcohol_moderado <-as.factor(datos$uso_alcohol_moderado)datos$uso_alcohol_exceso <-as.factor(datos$uso_alcohol_exceso)#Modelo que solo tiene X4modelo <-lm(datos$ln_tiempo_supervivencia~datos$test_higado)# SSE modelo con X4anova <-aov(modelo)resumen <-summary(anova)SSE <- resumen[[1]]$`Sum Sq`[2]# AIC_pn <-nrow(datos)p<-length(coef(modelo))AIC_p <- (n*log(SSE))-(n*log(n))+(2*p)AIC_p
[1] -103.2679
#Con sentencias de Rlibrary(olsrr)resultados <-ols_aic(modelo, method="SAS", corrected =FALSE)resultados
[1] -103.2679
# AIC con todos los posibles modelos#Modelo completo con X1, X2, X3, X4, p=5 modelo_completo <-lm(datos$ln_tiempo_supervivencia~datos$coagulacion_score+datos$indice_pronostico+datos$test_enzimas+datos$test_higado)library(olsrr)resultados_todos <-ols_step_all_possible(modelo_completo)aic_todos <-c(resultados_todos$result[3],resultados_todos$result[9])print(aic_todos)
El criterio Bayesian Information Criterion\(BIC_p\), que se calcula como se muestra en la Ecuación 8, penaliza a los modelos con muchas variables regresoras (como el \(R^2_{ap}\) y \(C_p\)).
Note que el primer término \(n~ln(SSE_p)\)decrece cuando \(p\) crece, el segundo término \(n~ln(n)\) es fijo para un valor de \(n\), el tercer término \(ln(n)~p\), aumenta con el número de parámetros \(p\).
En este sentido, los modelos con \(SSE_p\) bajos tienen buen desempeño en este criterio, siempre que la penalización \(ln(n)~p\) no sea tan grande.
Se escogen los modelos con menor \(AIC_p\).
Ejemplo de cálculo de \(BIC_p\).
#Datoslibrary(readxl)datos <-read_excel("data_cirugia.xlsx")datos$genero <-as.factor(datos$genero)datos$uso_alcohol_moderado <-as.factor(datos$uso_alcohol_moderado)datos$uso_alcohol_exceso <-as.factor(datos$uso_alcohol_exceso)#Modelo que solo tiene X4modelo <-lm(datos$ln_tiempo_supervivencia~datos$test_higado)# SSE modelo con X4anova <-aov(modelo)resumen <-summary(anova)SSE <- resumen[[1]]$`Sum Sq`[2]# BIC_pn <-nrow(datos)p<-length(coef(modelo))BIC_p <- (n*log(SSE))-(n*log(n))+(log(n)*p)BIC_p
[1] -99.28994
#Con sentencias de Rlibrary(olsrr)resultados <-ols_sbc(modelo, method="SAS")resultados
[1] -99.28994
# BIC con todos los posibles modelos#Modelo completo con X1, X2, X3, X4, p=5 modelo_completo <-lm(datos$ln_tiempo_supervivencia~datos$coagulacion_score+datos$indice_pronostico+datos$test_enzimas+datos$test_higado)library(olsrr)resultados_todos <-ols_step_all_possible(modelo_completo)bic_todos <-c(resultados_todos$result[3],resultados_todos$result[10])print(bic_todos)
El método Stepwise funciona de la siguiente forma:
Se empieza con las \(p-1\) variables regresoras.
Se elimina la variable regresora menos explicativa. Usando el siguiente estadístico, (Ecuación 9):
\[\begin{align}
F = \frac{\frac{SSE_0-SSE_1}{DFE_0-DFE_1}}{\frac{SSE_1}{DF_1}}
\end{align} \tag{9}\]
Donde:
\(SSE_1\): suma de cuadrados del error del modelo ampliado.
\(SSE_0\): suma de cuadrados del error del modelo reducido.
\(DFE_1\): grados de libertad del error del modelo ampliado.
\(DFE_0\): suma de cuadrados del error del modelo reducido.
Se ajusta el modelo con la con los \(p-2\) variables restantes.
Se repite numeral 2 hasta que todas las variables del modelo tengan un \(F\) predeterminado (olsrr cuando todos los \(p-values>0.3\), se puede modificar con el argumento penter).
Eliminación progresiva en R con olsrr
#Datoslibrary(readxl)datos <-read_excel("data_cirugia.xlsx")datos$genero <-as.factor(datos$genero)datos$uso_alcohol_moderado <-as.factor(datos$uso_alcohol_moderado)datos$uso_alcohol_exceso <-as.factor(datos$uso_alcohol_exceso)#Modelo completomodelo <-lm(datos$ln_tiempo_supervivencia~datos$coagulacion_score+datos$indice_pronostico+datos$test_enzimas+datos$test_higado+datos$edad+datos$genero+datos$uso_alcohol_moderado+datos$uso_alcohol_exceso)#Backwardlibrary (olsrr)ols_step_backward_p(modelo, details=TRUE) #details es para mostrar el paso a paso
Stepwise Summary
----------------------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
----------------------------------------------------------------------------------------
0 Full Model -5.532 14.358 -155.257 0.84611 0.81875
1 datos$test_higado -7.497 10.403 -157.634 0.84601 0.82258
2 datos$uso_alcohol_moderado -8.613 7.299 -159.406 0.84347 0.82348
----------------------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
--------------------------------------------------------------
R 0.918 RMSE 0.193
R-Squared 0.843 MSE 0.043
Adj. R-Squared 0.823 Coef. Var 3.211
Pred R-Squared 0.784 AIC -8.613
MAE 0.162 SBC 7.299
--------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
-------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
-------------------------------------------------------------------
Regression 10.800 6 1.800 42.209 0.0000
Residual 2.004 47 0.043
Total 12.805 53
-------------------------------------------------------------------
Parameter Estimates
-----------------------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
-----------------------------------------------------------------------------------------------------
(Intercept) 4.054 0.235 17.272 0.000 3.582 4.527
datos$coagulacion_score 0.072 0.019 0.233 3.839 0.000 0.034 0.109
datos$indice_pronostico 0.014 0.002 0.473 8.051 0.000 0.010 0.017
datos$test_enzimas 0.015 0.001 0.653 10.909 0.000 0.012 0.018
datos$edad -0.003 0.003 -0.078 -1.343 0.186 -0.009 0.002
datos$genero1 0.087 0.058 0.089 1.512 0.137 -0.029 0.203
datos$uso_alcohol_exceso1 0.351 0.076 0.280 4.597 0.000 0.197 0.505
-----------------------------------------------------------------------------------------------------
Introducción progresiva (Forwards Selection = FS)
El método Stepwise funciona de la siguiente forma:
Se empieza con \(0\) variables regresoras.
Se introduce la variable con mayor estadístico \(F\), (Ecuación 9). Esto equivale a introducir la variable con menor \(p-value\) para los contrastes individuales \(\beta_k=0\)
Se repite numeral 2 hasta que el valor de \(F\) de todas las variables que no están incluidas en el modelo es menor que un nivel predeterminado. (olsrr cuando los \(p-values>0.3\), se puede modificar con el argumento penter)
Introducción progresiva en R con olsrr
#Datoslibrary(readxl)datos <-read_excel("data_cirugia.xlsx")datos$genero <-as.factor(datos$genero)datos$uso_alcohol_moderado <-as.factor(datos$uso_alcohol_moderado)datos$uso_alcohol_exceso <-as.factor(datos$uso_alcohol_exceso)#Modelo completomodelo <-lm(datos$ln_tiempo_supervivencia~datos$coagulacion_score+datos$indice_pronostico+datos$test_enzimas+datos$test_higado+datos$edad+datos$genero+datos$uso_alcohol_moderado+datos$uso_alcohol_exceso)#Forwardlibrary (olsrr)ols_step_forward_p(modelo, details=TRUE) #details es para mostrar el paso a paso
Stepwise Summary
--------------------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
--------------------------------------------------------------------------------------
0 Base Model 79.529 83.507 -76.676 0.00000 0.00000
1 datos$test_enzimas 51.434 57.401 -105.440 0.42725 0.41624
2 datos$indice_pronostico 24.767 32.723 -131.597 0.66318 0.64997
3 datos$uso_alcohol_exceso 4.243 14.188 -150.402 0.77805 0.76473
4 datos$coagulacion_score -8.131 3.803 -160.533 0.82992 0.81603
5 datos$genero -8.580 5.343 -160.229 0.83746 0.82053
6 datos$edad -8.613 7.299 -159.406 0.84347 0.82348
--------------------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
--------------------------------------------------------------
R 0.918 RMSE 0.193
R-Squared 0.843 MSE 0.043
Adj. R-Squared 0.823 Coef. Var 3.211
Pred R-Squared 0.784 AIC -8.613
MAE 0.162 SBC 7.299
--------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
-------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
-------------------------------------------------------------------
Regression 10.800 6 1.800 42.209 0.0000
Residual 2.004 47 0.043
Total 12.805 53
-------------------------------------------------------------------
Parameter Estimates
-----------------------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
-----------------------------------------------------------------------------------------------------
(Intercept) 4.054 0.235 17.272 0.000 3.582 4.527
datos$test_enzimas 0.015 0.001 0.653 10.909 0.000 0.012 0.018
datos$indice_pronostico 0.014 0.002 0.473 8.051 0.000 0.010 0.017
datos$uso_alcohol_exceso1 0.351 0.076 0.280 4.597 0.000 0.197 0.505
datos$coagulacion_score 0.072 0.019 0.233 3.839 0.000 0.034 0.109
datos$genero1 0.087 0.058 0.089 1.512 0.137 -0.029 0.203
datos$edad -0.003 0.003 -0.078 -1.343 0.186 -0.009 0.002
-----------------------------------------------------------------------------------------------------
Método exhaustivo o análisis de los mejores subconjuntos (Exhaustive or All-Possible-Subsets Method = APS)
Consiste en examinar todos los posibles subconjuntos de regresores con un cardinal \(r≤k\) predeterminado. Esto significa que hay que examinar \(2^{r−1}\) subconjuntos de regresores (ignoramos el conjunto vacío).
Los posibles submodelos se pueden comparar entre sí mediante un criterio de selección de modelos COMO \(C_p, AIC~o~BIC\)
Método exhaustivo en R con leaps
#Datoslibrary(readxl)datos <-read_excel("data_cirugia.xlsx")datos$genero <-as.factor(datos$genero)datos$uso_alcohol_moderado <-as.factor(datos$uso_alcohol_moderado)datos$uso_alcohol_exceso <-as.factor(datos$uso_alcohol_exceso)#Exhaustivolibrary(leaps)exhaustivo<-regsubsets(datos$ln_tiempo_supervivencia~., data = datos, nvmax =3, force.in =NULL, force.out =NULL, method ="exhaustive") #nbest = 1, un único modelo óptimo para cada número de predictores. nvmax = NULL para que no haya límite en el número de regresorea. # force.in: regresores que deberían estar en todos los modelos. force.out: variables que no deberían estar en ningún modelosummary(exhaustivo)