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)