Regresion lineal
INS FETP
4/02/2024
knitr::opts_chunk$set(echo = FALSE)
height<-survey$Height
hand<-survey$Wr.Hnd
plot(height~hand,pch=18,xlab="Writing handspan (cm)", ylab="Height (cm)", col="#ac70a9",
main="Writing handspan (cm) Vs Height (cm) Dispersion")
Relación entre la cuarta palmar y estatura
Como era de esperar, existe una asociación positiva entre la cuarta de un estudiante y su altura. Esa relación parece ser de naturaleza lineal. Para evaluar la fuerza de la relación lineal, puede encontrar el coeficiente de correlación estimado
Correlación
| Write.Hand | Hight | |
|---|---|---|
| Write.Hand | 1.0000000 | 0.6009909 |
| Hight | 0.6009909 | 1.0000000 |
Vemos que la correlación entre la cuarta de la mano de un estudiante y su altura es de 0.6009909 positivas, lo que indica que la relación es directa, es decir que cuando la longitud de la mano aumenta, la altura también lo hará; por otro lado su valor numérico nos indica que esta correlación es fuerte, es decir que los puntos tenderán a tener una distancia cercana a la recta que los pueda llegar a modelar.
Aunque hay 237 registros en el dataframe, la gráfica en realidad no muestra 237 puntos. Esto se debe a que faltan observaciones (codificadas como NA). De forma predeterminada, R elimina cualquier par “incompleto” al producir un gráfico como este. Para averiguar cuántas observaciones incompletas se han eliminado, puede utilizar el operador lógico abreviado | junto con is.na Luego use la función length para descubrir que faltan 29 pares de observaciones.
| Observaciones_NA |
|---|
| 29 |
Conceptos generales
El propósito de un modelo de regresión lineal es generar una función que estime la media de una variable dado un valor particular de otra variable. Estas variables se conocen como la variable de respuesta (la variable de “resultado” cuya media está tratando de encontrar) y la variable explicativa (la variable “predictora” cuyo valor ya tiene). En términos del ejemplo de la encuesta de estudiantes, podría preguntar algo como “¿Cuál es la altura esperada de un estudiante si su cuarta palmar mide 14,5 cm?” Aquí la variable de respuesta es la altura, y la variable explicativa es la cuarta.
Definición del Modelo
Suponga que está buscando determinar el valor de la variable de respuesta Y dado el valor de una variable explicativa X. El modelo de regresión lineal simple establece que el valor de una respuesta se expresa con la siguiente ecuación:
Para un modelo de regresión lineal con dos variables cuantitativas continuas, el número de posibles modelos se calcula usando la fórmula \(2^(n+1) - 1\), donde “n” es el número de predictores (variables independientes). En este caso, la medida de la cuarta de los estudiantes, por lo tanto “n” es igual a 1.
Entonces, el número de posibles modelos sería.
Así que, hay 3 posibles modelos diferentes.
- Modelo con la medida de la mano como predictor.
- Modelo con la constante \(\beta_0\) como predictor
- Modelo sin predictor, es decir decimos que la altura se predice por medio del promedio de esta.
Supuestos de los residuales
La validez de las conclusiones que puede extraer con base en el modelo depende de manera crítica de los supuestos hechos sobre \(\varepsilon_{i}\) , que se definen de la siguiente manera:
• Se supone que el valor de \(\varepsilon_{i}\) tiene una distribución normal tal que \(\varepsilon_{i}\) y \(\sim N(0, \sigma^{2})\).
• \(\varepsilon_{i}\) está centrada (es decir, tiene una media de) cero.
• La varianza de \(\varepsilon_{i}\) , \(\sigma^{2}\), es constante.
El término \(\varepsilon_{i}\) representa un error aleatorio. En otras palabras, se asume que cualquier valor de la respuesta se debe a un cambio lineal en un valor dadode X , más o menos alguna variación residual aleatoria o normalmente distribuida.
Parámetros
El valor denotado por \(\beta_{0}\) se llama intersecto, y el de \(\beta_{1}\) se llama pendiente. Juntos, también se conocen como coeficientes de regresión y se interpretan de la siguiente manera:
• El intercepto, \(\beta_{0}\), se interpreta como el valor esperado de la variable de respuesta cuando el predictor es cero.
• Generalmente, la pendiente, \(\beta_{1}\), es el foco de interés. Esto se interpreta como el cambio en la respuesta media por cada aumento de una unidad en el predictor. Cuando la pendiente es positiva, la línea de regresión aumenta de izquierda a derecha (la respuesta media es mayor cuanto mayor es el predictor); cuando la pendiente es negativa, la línea disminuye de izquierda a derecha (la respuesta media es más baja cuando el predictor es más alto). Cuando la pendiente es cero, esto implica que el predictor no tiene efecto sobre el valor de la respuesta. Cuanto más extremo es el valor de \(\beta_{1}\), (es decir, lejos de cero), más inclinada se vuelve la línea creciente para el caso positivo o decreciente en el caso negativo de \(\beta_{1}\).
Estimando el intercepto y pendiente
El objetivo es utilizar sus datos para estimar los parámetros de regresión, lo que arroja las estimaciones \(\beta_{0}\) y \(\beta_{1}\); esto se conoce como ajuste del modelo lineal. En este caso, los datos comprenden n pares de observaciones para cada individuo. El modelo ajustado de interés se refiere al valor medio de respuesta, denotado \(\hat{y_{i}}\), para un valor específico del predictor, x, y se escribe de la siguiente manera:
\(\hat{y_{i}} = \hat{\beta_{0}} + \hat{\beta_{1}}x_{i}\)
A veces, se usa una notación alternativa como \(E\left [ Y \right ]\) o \(E\left [ Y|X= x \right ]\) para enfatizar el hecho de que el modelo da la media (es decir, el valor esperado) de la respuesta. Mmuchos simplemente usan algo como \(\hat{y}\). Note \(x_{i}\) y \(y_{i}\) y sus n pares de datos observados para las variables predictoras y de respuesta, respectivamente; \(i\) = 1, . . . , n. Entonces, las estimaciones de los parámetros para la función de regresión lineal simple son
\(\hat{\beta}_{1}=\rho_{xy}\frac{S_{y}}{S_{x}}\) y \(\hat{\beta}_{0}=\bar{y} - \beta_{1}\bar{x}\)
Donde \(\bar{y}\) y \(\bar{x}\) son las medias de los \(x_{i}\) y \(y_{i}\)
\(S_{x}\) y \(S_{y}\) son las desviaciones estándard
\(\rho _{xy}\) es la correlación entre xy
Ajuste de modelos lineales con lm
En R, el comando lm realiza la estimación. Por ejemplo, en el siguiente código se crea una línea, un objeto de modelo lineal ajustado de la altura media de los estudiantes por la cuarta de la palma y se almacena en su entorno global como survfit:
##
## Call:
## lm(formula = Height ~ Wr.Hnd, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.7276 -5.0706 -0.8269 4.9473 25.8704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 113.9536 5.4416 20.94 <2e-16 ***
## Wr.Hnd 3.1166 0.2888 10.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.909 on 206 degrees of freedom
## (29 observations deleted due to missingness)
## Multiple R-squared: 0.3612, Adjusted R-squared: 0.3581
## F-statistic: 116.5 on 1 and 206 DF, p-value: < 2.2e-16El primer argumento es la respuesta ~ predictor, que especifica el modelo deseado. No tiene que usar el prefijo \(survey\) para extraer los vectores del dataframe porque le indica específicamente a lm que busque en el objeto proporcionado al argumento de datos.
El propio objeto de modelo lineal ajustado, survfit, tiene una clase especial en R, una de “lm”. Un objeto de la clase “lm” puede considerarse esencialmente como una lista que contiene varios componentes que describen el modelo. Si simplemente ingresa el nombre del objeto “lm”, proporcionará el resultado más básico: una repetición de su llamada y las estimaciones del intercepto \(\hat\beta {_{0}}\) y la pendiente \(\hat\beta{_{1}}\).
##
## Call:
## lm(formula = Height ~ Wr.Hnd, data = survey)
##
## Coefficients:
## (Intercept) Wr.Hnd
## 113.954 3.117Esto revela que el modelo lineal sería
\(\hat{y}\) = 113.954 + 3.117x
Si evalúa la función matemática para \(\hat{y}\) en un rango de diferentes valores para x, terminas con una línea recta cuando graficas los resultados. Teniendo en cuenta la definición de intersecto dada anteriormente como el valor esperado de la variable de respuesta cuando el predictor es cero, en el ejemplo actual, esto implicaría que la altura media de un estudiante con una envergadura de 0 cm es 113,954 cm (posiblemente menos). declaración que es útil ya que un valor de cero para la variable explicativa no tiene sentido; pero quiere decir que el modelo no parte de cero. La pendiente, es el cambio en la respuesta media por cada incremento de una unidad en estatura que es 3.117.
Esto establece que, en promedio, por cada 1 cm de aumento en la longitud de la mano, se estima que la altura de un estudiante aumenta en 3,117 cm.
Al ejecutar una vez más la línea para graficar los datos sin procesar y se muestra en la Figura, pero ahora agregue la línea de regresión ajustada usando abline. Cuando pasa un objeto de clase “lm” que representa un modelo lineal simple, como survfit, se agregará la línea de regresión ajustada.
## (Intercept) Wr.Hnd
## 113.95 3.12
Realizando una interpretación ajustada al ejemplo, no podemos decir que \(\beta_0\) sea el valor cuando la cuarta de la mano sea 0, lo que si podemos afirmar es que \(\beta_0=113.954\) es el valor estimado de la altura cuando no se tiene en cuenta la medida de la mano.
Los residuales
Cuando los parámetros se estiman, la línea ajustada es considerada como una implementación de la regresión de mínimos cuadrados porque es la línea que minimiza la diferencia cuadrática promedio entre los datos observados y ella misma. Este concepto es más fácil de entender dibujando las distancias entre las observaciones y la línea ajustada, formalmente llamados residuales, para un par de observaciones individuales. Primero, extraigamos dos registros específicos de los vectores de datos Wr.Hnd y Height y llamemos a los vectores resultantes obsA y obsB.
## [1] 15.00 170.18## [1] 21.50 172.72A continuación, inspeccione brevemente los nombres de los miembros del objeto survfit.
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "na.action" "xlevels" "call" "terms"
## [13] "model"Estos son los componentes que conforman automáticamente un objeto de modelo “lm” ajustado, mencionado brevemente anteriormente. Tenga en cuenta que hay un componente llamado “coefficients”. Este contiene un vector numérico de las estimaciones de la intersecto y la pendiente. Puede extraer este componente (y, de hecho, cualquiera de los otros enumerados aquí) de la misma manera que realizaría una referencia de miembro en una lista con nombre: ingresando survfit$coficients en el indicador. Sin embargo, siempre que sea posible, es técnicamente preferible para fines de programación extraer dichos componentes utilizando una función de “acceso directo”. Para el componente de coeficientes de un objeto “lm”, la función que usa es coef.
## (Intercept) Wr.Hnd
## 113.953623 3.116617Aquí, los coeficientes de regresión se extraen del objeto y luego se asignan por separado a los objetos beta0.hat y beta1.hat. Otras funciones comunes de acceso directo incluyen resid y fitted; estos dos pertenecen a los componentes “residuals” y “fitted.values”, respectivamente. uso segmentos para dibujar las líneas discontinuas verticales presentes en la Figura
Las líneas punteadas se encuentran con la línea ajustada en las ubicaciones del eje vertical pasadas a y0, que, con el uso de los coeficientes de regresión beta0.hat y beta1.hat, refleja la ecuación. Ahora, imagina una colección de líneas de regresión alternativas dibujadas a través de los datos (logrados al alterar el valor de la intersección y la pendiente). Luego, para cada una de las líneas de regresión alternativas, imagina que calculas los residuos (distancias verticales) entre el valor de respuesta de cada observación y el valor ajustado de esa línea. La línea de regresión lineal estimada es la línea que se encuentra “más cerca de todas las observaciones”. Con esto se quiere decir que el modelo de regresión ajustado está representado por la línea estimada que pasa por la coordenada proporcionada por las medias variables (\(\bar{x}\), \(\bar{y}\)), y es la línea que produce la medida total más pequeña del residuo cuadrático distancias Por esta razón, otro nombre para una ecuación de regresión estimada por mínimos cuadrados como esta es la línea de mejor ajuste.
Inferencia estadística
La estimación de una ecuación de regresión es relativamente sencilla, pero esto es solo el comienzo. En la regresión lineal simple, hay una pregunta natural que siempre se debe hacer: ¿Existe evidencia estadística que respalde la presencia de una relación entre el predictor y la respuesta? Para decirlo de otra manera, ¿hay evidencia de que un cambio en la variable explicativa afecte el resultado medio?
Resumen del modelo ajustado
R lleva a cabo automáticamente este tipo de inferencia basada en modelos cuando se procesan los objetos lm. El uso de la función de resumen en un objeto creado por lm le proporciona una salida mucho más detallada que simplemente imprimir el objeto en la consola. Por el momento, se centrará sólo en dos aspectos de la información presentada en resumen: las pruebas de significación asociadas con los coeficientes de regresión y la interpretación del llamado coeficiente de determinación (etiquetado como R-cuadrado). Use el resumen en el ajuste del objeto del modelo actual y verá lo siguiente:
##
## Call:
## lm(formula = Height ~ Wr.Hnd, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.7276 -5.0706 -0.8269 4.9473 25.8704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 113.9536 5.4416 20.94 <2e-16 ***
## Wr.Hnd 3.1166 0.2888 10.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.909 on 206 degrees of freedom
## (29 observations deleted due to missingness)
## Multiple R-squared: 0.3612, Adjusted R-squared: 0.3581
## F-statistic: 116.5 on 1 and 206 DF, p-value: < 2.2e-16Pruebas de significancia del coeficiente de regresión
La primer columna de la tabla de coeficientes contiene las estimaciones puntuales del intersecto y la pendiente; la tabla también incluye estimaciones de los errores estándar de estas estadísticas. Puede demostrarse que los coeficientes de regresión lineal simple, cuando se estiman utilizando mínimos cuadrados, siguen una distribución t con n-2 grados de libertad. El valor t estandarizado y un valor p se informan para cada parámetro. Estos representan los resultados de una prueba de hipótesis de dos colas definida formalmente como \(H_{0}=\beta _{j} = 0\)
\(H_{A}=\beta _{j} \neq 0\)
donde j = 0 para el intersecto y j = 1 para la pendiente.
Con el valor nulo de cero, la verdad de \(H_{0}\) implica que el predictor no tiene efecto sobre la respuesta. La afirmación aquí está interesada en si hay algún efecto de la covariable, no en la dirección de este efecto, por lo que HA tiene dos caras (a través de \(\neq 0\)). Al igual que con cualquier prueba de hipótesis, cuanto menor sea el valor p, mayor será la evidencia en contra de \(H_{0}\). Con un valor p pequeño (< 2 10−16) asociado a este estadístico de prueba en particular (T = (3.116 0)/0.2888 = 10.79), por lo tanto, concluiría que hay pruebas sólidas en contra de la afirmación de que el predictor no tiene efecto sobre el nivel medio de la respuesta.
La misma prueba se lleva a cabo para el intercepto, pero la prueba para el el parámetro de pendiente \(\beta _{1}\) suele ser más interesante (dado que el rechazo de la hipótesis nula para \(\beta _{0}\) simplemente indica evidencia de que la línea de regresión no choca con el eje vertical en cero), especialmente cuando los datos observados no incluyen x = 0,
A partir de esto, puede concluir que el modelo ajustado sugiere que hay evidencia de que un aumento en palma está asociado con un aumento en la altura entre la población que se estudia. Por cada centímetro adicional de la palma, el aumento promedio de altura es de aproximadamente 3,12 cm. También podría generar intervalos de confianza para sus estimaciones pero R proporciona una función conveniente para que un objeto de clase “lm”.
| 2.5 % | 97.5 % | |
|---|---|---|
| (Intercept) | 103.225178 | 124.682069 |
| Wr.Hnd | 2.547273 | 3.685961 |
A la función confint le pasa a su objeto modelo como primer argumento su nivel de confianza deseado como nivel. Esto indica que debe tener un 95 por ciento de confianza en que el verdadero valor de \(\beta_{1}\) se encuentra entre 2,55 y 3,69. Como de costumbre, la exclusión del valor nulo de cero refleja el resultado estadísticamente significativo de antes.
Coeficiente de Determinación
La salida de resumen también le proporciona los valores de R-cuadrado y R-cuadrado Ajustado, que son particularmente interesantes. Ambos se denominan coeficiente de determinación; describen la proporción de la variación en la respuesta que se puede atribuir al predictor. Para la regresión lineal simple, la primera medida (sin ajustar) se obtiene simplemente como el cuadrado del coeficiente de correlación estimado. Para el ejemplo de la altura del estudiante, primero almacene la correlación estimada entre Wr.Hnd y Height como rho.xy, y luego la eléva al cuadrado.
## [1] 0.3611901Obtiene el mismo resultado que el valorR-cuadrado. Esto le dice que alrededor del 36,1 por ciento de la variación en las alturas de los estudiantes se puede atribuir a la envergadura de la mano. La medida ajustada es una estimación alternativa que tiene en cuenta el número de parámetros que requieren estimación. La medida ajustada generalmente es importante solo si está utilizando el coeficiente de determinación para evaluar la “calidad” general del modelo ajustado en términos de un equilibrio entre la bondad del ajuste y la complejidad.
Otros resultados de resumen
El resumen del modelo le proporciona información aún más útil. El “error estándar residual” es el error estándar estimado del término \(\varepsilon_{i}\) (en otras palabras, la raíz cuadrada de la varianza estimada de \(\varepsilon_{i}\), es decir, \(\sigma^{2}\)). El resultado también proporciona un resumen de cinco números para las distancias residuales.
Como resultado final, se proporciona una determinada prueba de hipótesis realizada utilizando la distribución F. Esta es una prueba global del impacto de su(s) predictor(es) en la respuesta; esto se explorará junto con la regresión lineal múltiple.
Puede acceder a todos los resultados proporcionados por el resumen directamente, como objetos R individuales, en lugar de tener que leerlos fuera de la pantalla del resumen impreso completo. Así como los nombres (survfit) le brindan una indicación del contenido del objeto survfit independiente, el siguiente código le proporciona los nombres de todos los componentes accesibles después de usar summary para procesar survfit.
## [1] "call" "terms" "residuals" "coefficients"
## [5] "aliased" "sigma" "df" "r.squared"
## [9] "adj.r.squared" "fstatistic" "cov.unscaled" "na.action"se pueden extraer utilizando el operador pesos como de costumbre. El error estándar residual, por ejemplo, se puede recuperar directamente con esto:
## [1] 7.90878Predicción
La capacidad de ajustar un modelo estadístico significa que no solo puede comprender y cuantificar la naturaleza de las relaciones en sus datos (como el aumento estimado de 3,1166 cm en la altura media por cada 1 cm de aumento de la palma para el ejemplo del estudiante), sino que también puede predecir valores de el resultado de interés. Sin embargo, como con cualquier estadística, siempre existe la necesidad de acompañar cualquier estimación puntual o predicción con una medida de dispersión.
¿Intervalo de confianza o intervalo de predicción?
Con un modelo lineal simple ajustado, puede calcular una estimación puntual del valor de respuesta medio, condicionado al valor de una variable explicativa. Para hacer esto, simplemente conecta (en la ecuación del modelo ajustado) el valor de x que le interesa. Una estadística como esta siempre está sujeta a variaciones, así que al igual use el intervalo de confianza para la respuesta media (CI) y medir esta incertidumbre.
\(\hat{y}=t_{1-\alpha/2,n-2}S_{\varepsilon}\sqrt{\frac{1}{n}}+\frac{(x-\bar{x})^{2}}{(n-1)s_{x}^{2}}\)
donde se obtiene el límite inferior por resta, el límite superior por suma. Aquí, \(\hat{y}\) es el valor ajustado (de la línea de regresión) en x; t(1−α/2,n−2) es el valor crítico apropiado de una distribución t con n 2 grados de libertad (en otras palabras, resultando en un área de cola superior de exactamente α/2); sE es el error estándar residual estimado; y \(\hat{x}\) y \(s^{2}\) representan la media muestral y la varianza de las observaciones del predictor, respectivamente. Un intervalo de predicción (PI) para una respuesta observada es diferente del intervalo de confianza en términos de contexto. Cuando los IC se utilizan para describir la variabilidad de la respuesta media, se utiliza un IP para proporcionar el posible rango de valores que podría tomar una realización individual de la variable de respuesta, dado X. Esta distinción es sutil pero importante: el IC corresponde a una media y el PI corresponde a una observación individual. Se puede demostrar que 100(1 α) el intervalo de predicción porcentual para una respuesta individual dado un valor de x se calcula con la siguiente ecuación
\(\hat{y}=t_{1-\alpha/2,n-2}S_{\varepsilon}\sqrt{{1+\frac{1}{n}}+\frac{(x-\bar{x})^{2}}{(n-1)s_{x}^{2}}}\)
Interpretando Intervalos
Supongamos que desea determinar la altura media para estudiantes con cuarta de 14,5 cm y para estudiantes con un cuarta de 24 cm. Las estimaciones puntuales en sí mismas son fáciles: simplemente ingrese los valores x deseados en la ecuación de regresión
## [1] 159.1446## [1] 188.7524Según el modelo, se pueden esperar alturas medias de 159,14 y 188,75 cm para palmas de 14,5 y 24 cm, respectivamente. La función de coerción as.numeric se usa simplemente para eliminar el resultado de los nombres anotativos que de otro modo estarían presentes en los objetos beta0.hat y beta1.hat
Intervalos de confianza para alturas medias
Para encontrar intervalos de confianza para estas estimaciones, se pueden calcular manualmente usando (20.5), pero, R tiene un comando de predicción incorporado para hacerlo. Para usar predict, primero debe almacenar sus valores de x de una manera particular: como una columna en un nuevo dataframe. El nombre de la columna debe coincidir con el predictor utilizado originalente para crear el objeto de modelo ajustado. En este ejemplo, crearé un nuevo dataframe, xvals, con la columna denominada Wr.Hnd, que contiene solo dos valores de interés: los palmas de 14,5 y 24 cm.
| Wr.Hnd |
|---|
| 14.5 |
| 24.0 |
Ahora, cuando se llama a predict, el primer argumento debe ser el objeto de interés del modelo ajustado, survfit para este ejemplo. A continuación, en el argumento newdata, pasa el dataframe especialmente construido y el contiene los valores de predicción especificados. Para el argumento de intervalo, debe especificar “confianza” como un valor de cadena de caracteres. El nivel de confianza, aquí establecido en 95 por ciento, se pasa (en la escala de probabilidad) al nivel.
| fit | lwr | upr |
|---|---|---|
| 159.1446 | 156.4956 | 161.7936 |
| 188.7524 | 185.5726 | 191.9323 |
Con una matriz con tres columnas, cuyo número (y orden) de filas corresponde a los valores predictores que proporcionó en el dataframe newdata. La primera columna, con un encabezado de ajuste, es la estimación puntual en la línea de regresión; puede ver que estos números coinciden con los valores que calculó anteriormente. Las otras columnas proporcionan los límites de CI superior e inferior como las columnas lwr y upr, respectivamente. En este caso, interpretaría esto como un 95 por ciento de confianza de que la estatura media de un estudiante con un palma de 14,5 cm se encuentra entre 156,5 cm y 161,8 cm y se encuentra entre 185,6 cm y 191,9 cm para un palma de 24 cm (cuando se redondea a 1 d.p.).
Intervalos de predicción para observaciones individuales
La función de predicción también proporcionará sus intervalos de predicción. Para encontrar el intervalo de predicción para posibles observaciones individuales con una cierta probabilidad, simplemente necesita cambiar el argumento de intervalo a “predicción”.
| fit | lwr | upr |
|---|---|---|
| 159.1446 | 143.3286 | 174.9605 |
| 188.7524 | 172.8390 | 204.6659 |
Observe que los valores ajustados siguen siendo los mismos, como indican las Ecuaciones. Sin embargo, los anchos de los PI son significativamente mayores que los de los correspondientes IC; esto se debe a que las propias observaciones sin procesar, en un valor x específico, naturalmente serán más variables que su media.
La interpretación cambia en consecuencia. Los intervalos describen dónde se pronostica que estarán las alturas brutas de los estudiantes “95 por ciento del tiempo”. Para una palma de 14,5 cm, el modelo predice que las observaciones individuales estarán entre 143,3 cm y 175,0 cm con una probabilidad de 0,95; por una palma de 24 cm, el mismo PI se estima en 172,8 cm y 204,7 cm.
Trazado de intervalos
Tanto los CI como los PI se adaptan bien a la visualización de modelos de regresión lineal simple. Con el siguiente código, puede comenzar la Figura trazando los datos y la línea de regresión estimada como en la Figura anterior, pero esta vez usando xlim e ylim en el gráfico para ampliar los límites x o y, además acomoda el largo y ancho total del CI y el PI.
observo que en este grafico en la medida que aumenta la longitud de la palma aumenta la estatura
A esto, agregue las ubicaciones de los valores ajustados para x = 14,5 y x = 24, así como dos conjuntos de líneas verticales que muestran los IC y los PI.
Los puntos marcan los valores ajustados para estos dos valores particulares de x. El primer segmento establece los PI como líneas grises verticales gruesas, y la segunda establece los CI como líneas negras verticales más cortas. Las coordenadas de estos segmentos de línea trazados se toman directamente de los objetos mypred.pi y mypred.ci, respectivamente.
También puede producir “bandas” alrededor de la línea de regresión ajustada que marcan uno o ambos intervalos sobre todos los valores del predictor. Desde el punto de vista de la programación, esto no es técnicamente posible para un sistema continuo, pero puede lograrlo prácticamente definiendo una secuencia fina de valores a lo largo del eje x (usando seq con un valor de length alto) y evaluando el IC y el PI en cada punto de esta secuencia final. Luego, simplemente une los puntos resultantes como líneas al trazar la gráfica.
En R, esto requiere que vuelva a ejecutar el comando de predicción de la siguiente manera
La primera línea de este código crea la secuencia final de valores predictores y la almacena en el formato requerido por el argumento newdata. Las coordenadas del eje y para las bandas CI y PI se almacenan como la segunda y tercera columna de los objetos de matriz ci.band y pi.band. Finalmente, se utiliza líneas para sumar cada una de las cuatro líneas discontinuas correspondientes a los límites superior e inferior de los dos intervalos, y una leyenda añade el toque final.
Tenga en cuenta que las bandas de CI discontinuas negras se encuentran con las líneas negras verticales y las bandas PI discontinuas grises se encuentran con las líneas grises verticales para los dos valores x individuales de antes, tal como era de esperar
La curvatura “inclinada hacia adentro” de los intervalos es característica de este tipo de gráfica y es especialmente visible en el IC. Esta curva se produce porque naturalmente hay menos variación si está prediciendo dónde hay más datos. Para obtener más información sobre la predicción de objetos de modelo lineal, eche un vistazo a la archivo de ayuda ?predict.lm.
Interpolación frente a extrapolación
Antes de terminar esta introducción a la predicción, es importante aclarar las definiciones de dos términos clave: interpolation y extrapolation. Estos términos describen la naturaleza de una predicción dada. Una predicción se conoce como interpolación si el valor de x que especifica se encuentra dentro del rango de sus datos observados; la extrapolación es cuando el valor de x de interés se encuentra fuera de este rango. De las predicciones puntuales que acaba de hacer, puede ver que la ubicación x = 14.5 es un ejemplo de interpolación, y x = 24 es un ejemplo de extrapolación.
En general, la interpolación es preferible a la extrapolación: tiene más sentido usar un modelo ajustado para la predicción en la vecindad de los datos que ya se han observado. Sin embargo, la extrapolación que no está demasiado lejos de esa vecindad aún puede considerarse confiable. La extrapolación para el ejemplo de la altura del estudiante en x = 24 es un buen ejemplo. Esto está fuera del rango de los datos observados, pero no mucho en términos de escala, y los intervalos estimados para el valor esperado de \(\hat{y}\) = 188.75 cm parecen, al menos visualmente, razonables dada la distribución de las otras observaciones. Por el contrario, tendría menos sentido utilizar el modelo ajustado para predecir la altura de los estudiantes en una palma de, digamos, 50 cm:
| fit | lwr | upr |
|---|---|---|
| 269.7845 | 251.9583 | 287.6106 |
Una extrapolación tan extrema sugiere que la altura media de un individuo con una envergadura de 50 cm es casi 270 cm, siendo ambas medidas bastante poco realistas. Lo mismo es cierto en la otra dirección; el intercepto \(\hat{\beta_{0}}\) no tiene una interpretación práctica particularmente útil, indicando que la altura media de un estudiante con un palmo de 0 cm es de alrededor de 114 cm. El mensaje principal aquí es usar el sentido común al hacer cualquier predicción a partir del ajuste de un modelo lineal. En cuanto a la fiabilidad de los resultados, son preferibles las predicciones realizadas en valores dentro de una proximidad adecuada de los datos observados.
Comprender los predictores categóricos
Hasta ahora, ha visto modelos de regresión lineal simple que se basan en variables explicativas continuas, pero también es posible usar una variable explicativa discreta o categórica, formada por k grupos o niveles distintos, para modelar la respuesta media. Debe poder hacer los mismos supuestos que las observaciones son todas independientes entre sí y que los residuos se distribuyen normalmente con una varianza homogénea. Para empezar, el caso más simple en el que k=2 (un predictor de valor binario), que forma la base de la situación un poco más complicada en la que el predictor categórico tiene más de dos niveles o un predictor: k > 2).
Variables binarias: k = 2
El modelo de regresión se especifica como \(Y|X=\beta_{0}+\beta_{1}X_{i}+\varepsilon_{i}\) para una variable de respuesta Y y predictor X , y \(\varepsilon_{i}\sim N(0,\varepsilon_{i})\). Ahora, suponga que su variable predictora es categórica, con solo dos niveles posibles (binario; k = 2) y observaciones codificadas como 0 o 1. Para este caso, (20.1) aún se cumple, pero la interpretación de los parámetros del modelo, \(\beta_{0}\) y \(\beta_{1}\), ya no es realmente el intercepto o la pendiente. En cambio, es mejor pensar en ellos como algo así como dos interceptos, donde \(\beta_{0}\) proporciona el valor de referencia o de referencia de la respuesta cuando X = 0 y \(\beta_{1}\) representa el efecto aditivo en la respuesta media si X = 1. En otras palabras, si X = 0, entonces Y = \(\beta_{0}\) + \(\varepsilon_{i}\); si X = 1, entonces Y = \(\beta_{0}\) + \(\beta_{1}\) + \(\varepsilon_{i}\) . Como de costumbre, la estimación es en términos de encontrar la respuesta media \(\hat{y}=\mathbb{E}[Y|X=x]\) según la ecuación se convierte en \(\hat{y}\) = \(\beta_{0}\) + \(\beta_{1}x\).
Regrese al dataframe de la encuesta y observe que tiene una variable Sexo, donde los estudiantes registraron su género. Mire la documentación en la página de ayuda ?survey o ingrese algo como esto:
## [1] "factor"| Female | Male |
|---|---|
| 118 | 118 |
Verá que la columna de datos de sexo es un vector de factores con dos niveles, Femenino y Masculino, y que resulta que hay un número igual de los dos (uno de los 237 registros tiene un valor faltante para esta variable). Vas a determinar si existe evidencia estadística de que la altura de un estudiante se ve afectada por el sexo. Esto significa que nuevamente está interesado en modelar la altura como la variable de respuesta, pero esta vez, es con la variable sexo categórica como predictor. Para visualizar los datos, se realiza de la siguiente manera y obtendrá un par de diagramas de caja o boxplot.
Esto se debe a que la variable de respuesta especificada a la izquierda de ~ es numérica y la variable explicativa a la derecha es un factor, y el comportamiento predeterminado de R en esa situación es producir diagramas de caja uno al lado del otro.
Para enfatizar aún más la naturaleza categórica de la variable explicativa, puede superponer las observaciones brutas de altura y sexo sobre los diagramas de caja. Para hacer esto, simplemente convierta el vector factor a numérico así as.numeric; esto se puede hacer directamente.
Recuerde que los diagramas de caja marcan la mediana como la línea central en negrita, pero que las regresiones lineales de mínimos cuadrados se definen por el resultado medio, por lo que es útil mostrar también las alturas medias según el sexo.
## Female Male
## 165.6867 178.8260
Con tapply el argumento na.rm=TRUE se compara con los puntos suspensivos en la definición de tapply y se pasa a mean (lo necesita para asegurarse de que los valores faltantes presentes en los datos no terminen produciendo NA como resultados). Un código adicional a los puntos agrega esas coordenadas (como símbolos) a la imagen; La gráfica indica, en general, que los hombres tienden a ser más altos que las mujeres, pero ¿hay evidencia estadística de una diferencia que respalde esto?
Modelo de Regresión Lineal de Variables Binarias
Para responder esto con un modelo de regresión lineal simple, puede usar lm para producir estimaciones de mínimos cuadrados al igual que con cualquier otro modelo que haya ajustado hasta ahora.
##
## Call:
## lm(formula = Height ~ Sex, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.886 -5.667 1.174 4.358 21.174
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 165.687 0.730 226.98 <2e-16 ***
## SexMale 13.139 1.022 12.85 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.372 on 206 degrees of freedom
## (29 observations deleted due to missingness)
## Multiple R-squared: 0.4449, Adjusted R-squared: 0.4422
## F-statistic: 165.1 on 1 and 206 DF, p-value: < 2.2e-16Sin embargo, debido a que el predictor es un vector de factores en lugar de un vector numérico, el informe de los coeficientes es ligeramente diferente. La estimación de \(\beta_{0}\) se informa nuevamente como (intersecto); esta es la estimación de la altura media si un estudiante es mujer. La estimación de \(\beta_{1}\) se reporta como SexMale. El coeficiente de regresión correspondiente de 13.139 es la diferencia estimada que se imparte sobre la altura media de un estudiante si es hombre. Si nos fijamos en la ecuación de regresión correspondiente
\(\hat{y}\) = \(\beta_{0}\) + \(\beta_{1}x\) = 165.687 + 13.139\(X_{i}\)
puede ver que el modelo se ha ajustado suponiendo que la variable x se define como “el individuo es hombre”: 0 para no/falso, 1 para sí/verdadero. Es decir, se asume como referencia el nivel de “femenino” para la variable sexo, y ese es el efecto de “ser hombre” sobre la estatura media que se estima explícitamente. La prueba de hipótesis para \(\beta_{0}\) y \(\beta_{1}x\) se realiza con las mismas hipótesis definidas antes
\(H_{0}=\beta_{j} = 0\)
\(H_{A}=\beta_{j} \neq 0\)
Nuevamente, es la prueba para \(\beta_{1}x\) la que generalmente es de mayor interés ya que es este valor el que le dice si existe evidencia estadística de que la variable de respuesta media se ve afectada por la variable explicativa, es decir, si \(\beta_{1}x\) es significativamente diferente de cero.
Predicciones de una variable categórica binaria
Debido a que solo hay dos valores posibles para x, la predicción es sencilla aquí. Cuando evalúa la ecuación, la única decisión que se debe tomar es si se debe usar \(\beta_{1}x\) (en otras palabras, si un individuo es hombre) o no (si un individuo es mujer). Por ejemplo, puede ingresar el siguiente código para crear un factor de cinco observaciones adicionales con los mismos nombres de nivel que los datos originales y almacenar los nuevos datos en extra.obs:
## [1] Female Male Male Male Female
## Levels: Female MaleLuego, use predict para encontrar las alturas medias en esos valores adicionales del predictor. (Recuerde que cuando pasa nuevos datos para predict usando el argumento newdata, los predictores deben tener la misma forma que los datos que se usaron para ajustar el modelo en primer lugar).
| fit | lwr | upr |
|---|---|---|
| 165.6867 | 164.4806 | 166.8928 |
| 178.8260 | 177.6429 | 180.0092 |
| 178.8260 | 177.6429 | 180.0092 |
| 178.8260 | 177.6429 | 180.0092 |
| 165.6867 | 164.4806 | 166.8928 |
Puede ver en el resultado de las predicciones son diferentes solo entre los dos conjuntos de valores: las estimaciones puntuales de las dos instancias de Femenino son idénticas, simplemente \(\beta_{0}\) con IC del 90 por ciento. Las estimaciones puntuales y los IC para las instancias de Male también son iguales entre sí, según una estimación puntual de \(\beta_{0}\)+ \(\beta_{1}\) Por sí solo, es cierto, este ejemplo no es demasiado emocionante. Sin embargo, es fundamental entender cómo R presenta los resultados de la regresión cuando se usan predictores categóricos, especialmente cuando se considera la regresión múltiple.
Variables multinivel: k > 2
Es común trabajar con datos donde las variables predictoras categóricas que tienen más de dos niveles, de modo que (k > 2). También pueden denominarse variables categóricas multinivel. Para lidiar con esta situación más complicada mientras conserva la interpretabilidad de sus parámetros, primero debe codificar de forma ficticia su predictor en k − 1 variables binarias. #### Codificación ficticia de variables multinivel Para ver cómo se hace esto, suponga que desea encontrar el valor de la variable de respuesta Y cuando se le da el valor de una variable categórica X , donde X tiene k > 2 niveles (asumir también las condiciones de validez de la regresión lineal).
En el modelado de regresión, la codificación ficticia es el procedimiento utilizado para crear varias variables binarias a partir de una variable categórica como X. En lugar de la única variable categórica con posibles realizaciones
X = 1, 2, 3, . . . , k lo recodifica en varias variables sí/no, una para cada nivel, con posibles realizaciones:
X(1) = 0, 1; X(2) = 0, 1; X(3) = 0, 1; . . . ; X(k) = 0, 1
Como puede ver, X(i) representa una variable binaria para el i-ésimo nivel de la X original. Por ejemplo, si un individuo tiene X = 2 en la variable categórica original, entonces X(2) = 1 (sí) y todos los demás (X(1), X(3), . . . , X( k)) será cero (no).
Suponga que X es una variable que puede tomar cualquiera de los k = 4 valores 1, 2, 3 o 4, y ha realizado seis observaciones de esta variable: 1, 2, 2, 1, 4, 3. la Tabla 20 -1 muestra estas observaciones y sus equivalentes codificados ficticiamente X(1) , X(2) , X(3) y X(4) .
Tabla 20-1: Ejemplo ilustrativo de codificación ficticia para seis observaciones de una variable categórica con k = 4 grupos
X X(1) X(2) X(3) X(4) 1 1 0 0 0 2 0 1 0 0 2 0 1 0 0 1 1 0 0 0 4 0 0 0 1 3 0 0 1 0
Al ajustar el modelo, generalmente solo usa k - 1 de las variables binarias ficticias: una de las variables actúa como referencia o nivel de referencia, y se incorpora a la intersección general del modelo. En la práctica, terminaría con un modelo estimado como este,
\(\hat{y}\) = \(\beta_{0}\) + \(\beta_{1}X_{2}\)+ \(\beta_{1}X_{2}\)…++ \(\beta_{k-1}X_{k}\)
suponiendo que 1 es el nivel de referencia. Como puede ver, además del término de intersección general\(\beta_{0}\), tiene k-1 otros términos de intersección estimados que modifican el coeficiente de referencia\(\beta_{0}\) , dependiendo de las categorías originales asume una observación. si una observación tiene X(3) = 1 y, por lo tanto, todos los demás valores binarios son cero (de modo que la observación habría tenido un valor de X = 3 para la variable categórica original), el valor medio predicho de la respuesta sería \(\hat{y}\) = \(\beta_{0}\) + \(\beta_{1}X_{2}\). Por otro lado, debido a que el nivel de referencia se define como 1, si una observación tiene valores de cero para todas las variables binarias, implica que la observación originalmente tenía X = 1, y la predicción sería simplemente \(\hat{y}\) = \(\beta_{0}\).
La razón por la que es necesario crear un código ficticio para las variables categóricas de esta naturaleza es que, en general, las categorías no se pueden relacionar entre sí en el mismo sentido numérico que las variables continuas. A menudo no es apropiado, por ejemplo, pensar que una observación en la categoría 4 es “el doble” que una en la categoría 2, es lo que asumirían los métodos de estimación. Sin embargo, las variables binarias de presencia/ausencia son válidas y se pueden incorporar fácilmente en el marco de modelado. La elección del nivel de referencia generalmente tiene una importancia secundaria: los valores específicos de los coeficientes estimados cambiarán en consecuencia, pero cualquier interpretación general que haga en función del modelo ajustado será la misma independientemente.
Implementar este enfoque de codificación ficticia es técnicamente una forma de regresión múltiple, ya que ahora incluye varias variables binarias en el modelo. Sin embargo, es importante ser consciente de la naturaleza un tanto artificial de la codificación ficticia; aún debe pensar en los coeficientes múltiples como si representaran una única variable categórica, ya que las variables binarias X(1), . . . , X(k) no son independientes entre sí. Es por eso que he elegido definir estos modelos. #### Modelo de Regresión Lineal de Variables Múltiples categorias R hace que trabajar con predictores categóricos de esta manera sea bastante simple, ya que automáticamente genera códigos ficticios para cualquier variable explicativa de este tipo cuando llama a lm. Sin embargo, hay dos cosas que debe verificar antes de ajustar su modelo. 1. La variable categórica de interés debe almacenarse como un factor (formalmente desordenado). 2. Debe comprobar que está satisfecho con la categoría asignada como nivel de referencia (para fines interpretativos. Por supuesto, también debe estar satisfecho con la validez de los supuestos familiares de normalidad e independencia de \(\varepsilon_{i}\).
Para demostrar todas estas definiciones e ideas, regresemos a los datos de la encuesta de estudiantes del paquete MASS y mantengamos la “altura del estudiante” como la variable de respuesta de interés. Entre los datos se encuentra la variable Humo. Esta variable describe el tipo de fumador que cada estudiante reporta, definido por frecuencia y dividido en cuatro categorías: “intenso”, “nunca”, “ocasional” y “regular”.
## [1] TRUE| Heavy | Never | Occas | Regul |
|---|---|---|---|
| 11 | 189 | 19 | 17 |
Aquí, el resultado de is.factor(survey$Smoke) indica que efectivamente tiene un vector de factores, la llamada a la tabla arroja el número de estudiantes en cada una de las cuatro categorías y puede explicar - Solicitar oportunamente el atributo de niveles de cualquier factor R a través de niveles. Preguntémonos si existe evidencia estadística que respalde una diferencia en la estatura promedio de los estudiantes según la frecuencia con que fuman. Puede crear un conjunto de diagramas de caja de estos datos con las siguientes dos líneas; La figura muestra el resultado.
Tenga en cuenta de la salida anterior de R que, a menos que se defina explícitamente en la creación, los niveles de un factor aparecen en orden alfabético de forma predeterminada, como es el caso de Smoke, y R establecerá automáticamente el primero (como se muestra en la salida de una llamada a niveles). como el nivel de referencia cuando ese factor se utiliza como predictor en el ajuste posterior del modelo. Teniendo en cuenta el modelo lineal usando lm, puede ver que, de hecho, el primer nivel de Smoke, para “pesado”, se ha utilizado como referencia.
##
## Call:
## lm(formula = Height ~ Smoke, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.02 -6.82 -1.64 8.18 28.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 173.7720 3.1028 56.005 <2e-16 ***
## SmokeNever -1.9520 3.1933 -0.611 0.542
## SmokeOccas -0.7433 3.9553 -0.188 0.851
## SmokeRegul 3.6451 4.0625 0.897 0.371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.812 on 205 degrees of freedom
## (28 observations deleted due to missingness)
## Multiple R-squared: 0.02153, Adjusted R-squared: 0.007214
## F-statistic: 1.504 on 3 and 205 DF, p-value: 0.2147Como se describe en la ecuación, se obtienen estimaciones de los coeficientes correspondientes a las variables binarias ficticias para tres de las cuatro categorías posibles en este ejemplo: los tres niveles de no referencia. La observación en la categoría de referencia Pesado está representado únicamente por \(\beta_{0}\), designado primero intersección, con los otros coeficientes proporcionando los efectos asociados con una observación en una de las otras categorías.
Predicciones de una variable categórica multples categorias
Se encuentran estimaciones puntuales a través de la predicción, como de costumbre.
## [1] Heavy Never Occas Regul
## Levels: Heavy Never Occas Regul| fit | lwr | upr |
|---|---|---|
| 173.7720 | 167.6545 | 179.8895 |
| 171.8200 | 170.3319 | 173.3081 |
| 173.0288 | 168.1924 | 177.8651 |
| 177.4171 | 172.2469 | 182.5874 |
Aquí, se ha creado el objeto one.of.each con fines ilustrativos; representa una observación en cada una de las cuatro categorías, almacenadas como un objeto que coincide con la clase (y los niveles) de los datos originales de Smoke. Se prevé que un estudiante en la categoría Occas, por ejemplo, tenga una estatura media de 173,772 0,7433 = 173,0287.
Sin embargo, el resultado del resumen del modelo anterior muestra que ninguno de los coeficientes de la variable ficticia binaria se considera estadísticamente significativo desde cero (porque todos los valores p son demasiado grandes). Los resultados indican, como podría haber sospechado, que no hay evidencia de que la frecuencia de fumar (o más específicamente, tener una frecuencia de fumar que sea diferente desde el nivel de referencia) afecta la estatura media de los estudiantes con base en esta muestra de individuos. Como es común, el coeficiente de referencia \(\beta_{0}\) es estadísticamente muy significativo, pero eso solo sugiere que el intercepto general probablemente no es cero. (Debido a que la variable de respuesta es una medida de altura y claramente no estará centrada en 0 cm, ese resultado tiene sentido). Los intervalos de confianza proporcionados se calculan de la manera habitual basada en t. El pequeño valor de R-Squared refuerza esta conclusión, indicando que casi ninguna de las variaciones en la respuesta puede explicarse cambiando la categoría de frecuencia de tabaquismo. Además, el valor p general de la prueba F es bastante grande, alrededor de 0.215, lo que sugiere un efecto general no significativo del predictor de la respuesta.
Cambio del nivel de referencia
A veces, puede decidir cambiar el nivel de referencia seleccionado, en comparación con el cual se estiman los efectos de asumir cualquiera de los otros niveles. Cambiar la línea de base dará como resultado la estimación de diferentes coeficientes, lo que significa que los valores p individuales están sujetos a cambios, pero el resultado general (en términos de importancia global del factor) no se verá afectado.
Debido a esto, la alteración del nivel de referencia solo se realiza con fines interpretativos; a veces hay una línea de base intuitivamente natural del predictor (por ejemplo, “Placebo” versus “Medicamento A” y “Medicamento B” como una variable de tratamiento en el análisis de algún ensayo clínico) a partir del cual se quiere estimar la desviación en la respuesta media con respecto a las demás categorías posibles. La redefinición del nivel de referencia se puede lograr rápidamente utilizando la función de renivelación incorporada en R. Esta función le permite elegir qué nivel viene primero en la definición de un objeto de vector de factor dado y, por lo tanto, se designará como el nivel de referencia en el ajuste posterior del modelo. . En el ejemplo actual, supongamos que prefiere tener los no fumadores como nivel de referencia.
## [1] "Never" "Heavy" "Occas" "Regul"La función de renivelación ha movido la categoría Nunca a la primera posición en el nuevo vector de factores. Si sigue ajustando el modelo nuevamente utilizando SmokeReordered en lugar de la columna original de Smoke de la encuesta, proporcionará estimaciones de los coeficientes asociados con los tres niveles diferentes de fumadores. Vale la pena señalar las diferencias en el tratamiento de vectores de factores desordenados versus ordenados en aplicaciones de regresión. Podría parecer sensato ordenar formalmente la variable de tabaquismo, por ejemplo, aumentando la frecuencia de tabaquismo al crear un nuevo vector de factores. Sin embargo, cuando se proporciona un vector de factor ordenado en un código lm, R reacciona de una manera diferente: no realiza la codificación ficticia relativamente simple que se analiza aquí, donde un efecto se asocia con cada nivel opcional a la línea de base (técnicamente denominado como contrastes ortogonales). En cambio, el comportamiento predeterminado es ajustar el modelo basado en algo llamado contrastes polinómicos, donde el efecto de la variable categórica ordenada en la respuesta se define en una forma funcional más complicada. Basta con decir que este enfoque puede ser beneficioso cuando su interés radica en la naturaleza funcional específica de “ascender” a través de un conjunto ordenado de categorías.
Tratamiento de variables categóricas como numéricas
La forma en que lm decide definir los parámetros del modelo ajustado depende principalmente del tipo de datos que pasa a la función. Como se discutió, lm impone una codificación ficticia solo si la variable explicativa es un vector de factores desordenado. A veces, los datos categóricos que desea analizar no se han almacenado como un factor en su objeto de datos. Si la variable categórica es un vector de caracteres, lm implícitamente la convertirá en un factor. Sin embargo, si la variable categórica buscada es numérica, entonces lm realiza la regresión lineal exactamente como si fuera un predictor numérico continuo; estima un único coeficiente de regresión, que se interpreta como un “cambio por unidad” en la respuesta media. Esto puede parecer inapropiado si se supone que la variable explicativa original está compuesta por distintos grupos. En algunos entornos, sin embargo, especialmente cuando la variable puede tratarse naturalmente como numérica discreta, este tratamiento no solo es válido estadísticamente sino que también ayuda con la interpretación. Supongamos que está interesado en las variables kilometraje, mpg (continuo) y número de cilindros, cyl (discreto; el conjunto de datos contiene automóviles con 4, 6 u 8 cilindros). Ahora, es perfectamente sensato pensar automáticamente en cyl como una variable categórica. Tomando mpg como la variable de respuesta, los diagramas de caja son adecuados para reflejar la naturaleza agrupada de cyl como predictor; el resultado de la siguiente línea se da a la izquierda de la figura
Al ajustar el modelo de regresión asociado, debe tener en cuenta lo que le está indicando a R que haga. Dado que la columna cyl de mtcars es numérica y no un vector de factores per se, lm la tratará como continua si solo accede directamente al marco de datos.
## [1] "numeric"##
## Call:
## lm(formula = mpg ~ cyl, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9814 -2.1185 0.2217 1.0717 7.5186
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.8846 2.0738 18.27 < 2e-16 ***
## cyl -2.8758 0.3224 -8.92 6.11e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.206 on 30 degrees of freedom
## Multiple R-squared: 0.7262, Adjusted R-squared: 0.7171
## F-statistic: 79.56 on 1 and 30 DF, p-value: 6.113e-10Al igual que en las secciones anteriores, ha recibido una estimación de intercepto y la pendiente; el último es estadísticamente significativo, lo que indica que existe evidencia en contra de que el verdadero valor de la pendiente sea cero. Su línea de regresión ajustada es
\(\hat{y} = \beta_{0} + \beta_{1}X_{i}\) = 37.88 − 2.88\(X_{i}\)
donde \(\hat{y}\) es el kilometraje promedio y x es numérico: el número de cilindros. Por cada cilindro adicional, el modelo dice que su millaje disminuirá en 2.88 MPG, en promedio. Es importante reconocer el hecho de que ha ajustado una línea continua a lo que efectivamente son datos categóricos.
Diagrama de dispersión de los mismos datos con línea de regresión ajustada (tratando cyl como numérico-continuo) superpuesto.
Algunos investigadores ajustan predictores categóricos o discretos como variables continuas a propósito. Primero, permite la interpolación; por ejemplo, podría usar este modelo para evaluar el MPG promedio de un automóvil de 5 cilindros. Segundo, significa que hay menos parámetros que requieren estimación; en otras palabras, en lugar de k 1 intercepto para una variable categórica con k grupos, solo necesita un parámetro para la pendiente. Puede ser engañoso proceder de esta manera si las diferencias en la respuesta media de acuerdo con la categoría de predictor de una observación no están bien representadas linealmente; la detección de efectos significativos puede perderse por completo. Como mínimo, es importante reconocer esta distinción al ajustar modelos. Si acaba de reconocer que R ha ajustado la variable cyl como continua y desea ajustar realmente el modelo con cyl como categórico, tendría que convertirlo explícitamente en un factor vector de antemano
##
## Call:
## lm(formula = mpg ~ factor(cyl), data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.2636 -1.8357 0.0286 1.3893 7.2364
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.6636 0.9718 27.437 < 2e-16 ***
## factor(cyl)6 -6.9208 1.5583 -4.441 0.000119 ***
## factor(cyl)8 -11.5636 1.2986 -8.905 8.57e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.223 on 29 degrees of freedom
## Multiple R-squared: 0.7325, Adjusted R-squared: 0.714
## F-statistic: 39.7 on 2 and 29 DF, p-value: 4.979e-09Aquí, al convertir cyl en un factor para la fórmula para lm, puede ser que ha obtenido estimaciones de coeficientes de regresión para los niveles de cyl correspondientes a automóviles de 6 y 8 cilindros (con el nivel de referencia establecido automáticamente en 4 -coches cilíndricos).
Equivalencia con ANOVA unidireccional
Hay una observación final que hacer sobre los modelos de regresión con un solo predictor categórico nominal. Piense en el hecho de que estos modelos describen un valor de respuesta medio para los k grupos diferentes. aqui en realidad está haciendo lo mismo que en ANOVA unidireccional comparar más de dos medias y determinar si existe evidencia estadística de que al menos una media es diferente de las otras. De hecho, la regresión lineal simple con un solo predictor categórico, implementada mediante la estimación de mínimos cuadrados, es solo otra forma de realizar ANOVA unidireccional. O, quizás de manera más concisa, ANOVA es un caso especial de regresión de mínimos cuadrados. El resultado de una prueba ANOVA unidireccional es un único valor p que cuantifica un nivel de evidencia estadística contra la hipótesis nula que establece que las medias de los grupos son iguales. Cuando tiene un predictor categórico en una regresión, es exactamente ese valor p el que se informa al final de la regresión.
La prueba de significación global para la altura del estudiante modelada por el ejemplo del estado de fumador: tuvo un valor p de 0.2147. Esto provino de un estadístico de prueba F de 1.504 con df1 = 3 y df2 = 205. Ahora, suponga que le acaban de entregar los datos y le piden que realice un ANOVA de una vía de la altura sobre el tabaquismo.
## Df Sum Sq Mean Sq F value Pr(>F)
## Smoke 3 434 144.78 1.504 0.215
## Residuals 205 19736 96.27
## 28 observations deleted due to missingnessEsos mismos valores tambien estan aquí; puedes encontrar la raíz cuadrada del MSE:
## [1] 9.811728De hecho, este es el “error estándar residual” dado en el resumen de la salida. Las dos conclusiones que sacaría sobre el impacto del tabaquismo en la estatura (una para la producción de la salida del codigo y la otra para la prueba ANOVA) son, por supuesto, las mismas. La prueba global que proporciona lm no solo está ahí para el beneficio de confirmar los resultados de ANOVA. Como una generalización de ANOVA, los modelos de regresión de mínimos cuadrados proporcionan más que solo pruebas de coeficientes específicos.
Esa prueba mundial se conoce formalmente como la prueba F ómnibus, y si bien es equivalente a ANOVA unidireccional en la configuración de “predictor categórico único”, también es una prueba útil, general e independiente de la contribución estadística de varios predictores. al valor del resultado.
Taller para ser ejecutado
El conjunto de datos de la encuesta tiene una variable llamada Exer, un factor con k = 3 niveles que describen la cantidad de tiempo de ejercicio físico que hace cada estudiante: ninguno, algo o frecuente. Obtenga un recuento del número de estudiantes en cada categoría y produzca diagramas de caja uno al lado del otro de la altura de los estudiantes divididos por ejercicio.
Asumiendo la independencia de las observaciones y la normalidad como de costumbre, ajuste un modelo de regresión lineal con la altura como variable de respuesta y el ejercicio como variable explicativa (codificación ficticia). ¿Cuál es el nivel de referencia predeterminado del predictor? Producir un resumen del modelo
C. Saque una conclusión basada en el modelo ajustado de (b): ¿parece que la frecuencia del ejercicio tiene algún impacto en la altura media? ¿Cuál es la naturaleza del efecto estimado?
- Prediga las alturas medias de un individuo en cada una de las tres categorías de ejercicio, acompañadas de intervalos de predicción del 95 por ciento. ¿Llega al mismo resultado
e interprete el modelo de altura por ejercicio si construye una tabla ANOVA usando aov? F. ¿Hay algún cambio en el resultado de (e) si modifica el modelo para que el nivel de referencia de la variable de ejercicio sea “ninguno”? ¿Esperarías que lo hubiera?
Este es el “error estándar residual” dado. Las dos conclusiones que sacaría sobre el impacto del tabaquismo en la altura (una para la producción de película y la otra para la prueba ANOVA) son, por supuesto, las mismas. La prueba global que proporciona lm no solo está ahí para el beneficio de confirmar los resultados de ANOVA. Como una generalización de ANOVA, los modelos de regresión de mínimos cuadrados proporcionan más que solo pruebas de coeficientes específicos.
Es la prueba formalmente conocida como la prueba F ómnibus, y si bien es equivalente a ANOVA unidireccional en la configuración de “predictor categórico único”, también es una prueba útil, general e independiente de la contribución estadística de varios predictores al valor del resultado.
Usando el dataframe de la encuesta del paquete MASS
El conjunto de datos de la encuesta tiene una variable llamada Exer, un factor con k = 3 niveles que describen la cantidad de tiempo de ejercicio físico que hace cada estudiante: ninguno, algo o frecuente. Obtenga un recuento del número de estudiantes en cada categoría y produzca diagramas de caja uno al lado del otro de la altura de los estudiantes divididos por ejercicio.
Asumiendo la independencia de las observaciones y la normalidad como de costumbre, ajuste un modelo de regresión lineal con la altura como variable de respuesta y el ejercicio como variable explicativa (codificación ficticia).
¿Cuál es el nivel de referencia predeterminado del predictor?
Producir un resumen del modelo.
C. Saque una conclusión basada en el modelo ajustado de (b): ¿parece que la frecuencia del ejercicio tiene algún impacto en la altura media? ¿Cuál es la naturaleza del efecto estimado?
- Prediga las alturas medias de un individuo en cada una de las tres categorías de ejercicio, acompañadas de intervalos de predicción del 95 por ciento.
¿Llega al mismo resultado e interpretación para el modelo de altura por ejercicio si construye una tabla ANOVA usando aov? F.
¿Hay algún cambio en el resultado de (e) si modifica el modelo para que el nivel de referencia de la variable de ejercicio sea “ninguno”? ¿Esperarías que lo hubiera?
En el conjunto de datos mtcars. Una de las variables en este dataframe es qsec, descrita como el tiempo en segundos que se tarda en correr un cuarto de milla; otro es gear, el número de marchas hacia adelante (los automóviles en este conjunto de datos tienen 3, 4 o 5 marchas).
Utilizando los vectores directamente del dataframe, un modelo ajustado es un modelo de regresión lineal simple con qsec como variable de respuesta y engranaje como variable explicativa e interprete el resumen del modelo.
Regresión múltiple
La regresión lineal múltiple es una generalización directa del predictor único.
Permite modelar una variable de respuesta continua en términos de más de un predictor para que pueda medir el efecto conjunto de varias variables explicativas en la variable de respuesta. Su variable de respuesta de esta manera se usará R para ajustar el modelo usando mínimos cuadrados. También explorará otros aspectos estadísticos clave del modelado lineal en el entorno R, como la transformación de variables y la inclusión de efectos interactivos.
La regresión lineal múltiple representa una parte importante de la práctica de la estadística.
Le permite controlar o ajustar múltiples fuentes de influencia sobre el valor de la respuesta, en lugar de solo medir el efecto de una variable explicativa (en la mayoría de las situaciones, hay más de un contribuyente a las medidas de resultado).
En el corazón de esta clase de métodos está la intención de descubrir relaciones potencialmente causales entre su variable de respuesta y el efecto (conjunto) de cualquier variable explicativa. En realidad, la causalidad en sí misma es extremadamente difícil de establecer, pero puede fortalecer cualquier evidencia de causalidad mediante el uso de un estudio bien diseñado respaldado por una sólida recopilación de datos y ajustando modelos que puedan medir de manera realista las relaciones presentes en sus datos.
Algunos términos
• Una variable oculta influye en la respuesta, otro predictor o ambos, pero no se mide (o no se incluye) en un modelo predictivo. Por ejemplo, supongamos que un investigador establece un vínculo entre el volumen de basura arrojado por un hogar y si el hogar posee un trampolín.
La potencial variable al acecho aquí sería el número de niños en el hogar; es más probable que esta variable se asocie positivamente con un aumento en la basura y las posibilidades de tener un trampolín. Una interpretación que sugiera que poseer un trampolín es una causa de mayor desperdicio sería errónea.
• La presencia de una variable oculta puede conducir a conclusiones falsas sobre las relaciones causales entre la respuesta y los otros predictores, o puede enmascarar una verdadera asociación de causa y efecto; este tipo de error se conoce como confusión.
Para decirlo de otra manera, puede pensar en la confusión como el entrelazamiento de los efectos de uno o más predictores en la respuesta.
• Una variable molesta o extraña es un predictor de interés secundario o nulo que tiene el potencial de confundir las relaciones entre otras variables y, por lo tanto, afectar sus estimaciones de los otros coeficientes de regresión.
Las variables extrañas se incluyen en el modelo como una cuestión de necesidad, pero la naturaleza específica de su influencia en la respuesta no es el interés principal del análisis.
Estas definiciones se aclararán una vez que comience a ajustar e interpretar los modelos de regresión. El mensaje principal que quiero enfatizar aquí, una vez más, es que correlación no implica causalidad. Si un modelo ajustado encuentra una asociación estadísticamente significativa entre un predictor (o predictores) y una respuesta, es importante considerar la posibilidad de que las variables ocultas estén contribuyendo a los resultados e intentar controlar cualquier confusión antes de sacar conclusiones. Los modelos de regresión múltiple le permiten hacer esto.
Teoría
Antes de comenzar a usar R para ajustar modelos de regresión, se examinarán las definiciones técnicas de un modelo de regresión lineal con múltiples predictores. Se observarán cómo funcionan los modelos en un sentido matemático y obtendrá un vistazo a los cálculos que ocurren “detrás de escena” cuando se estiman los parámetros del modelo en R.
Ampliación del modelo simple a un modelo múltiple
En lugar de tener solo un predictor, desea determinar el valor de una variable de respuesta continua Y dados los valores de p > 1 variables explicativas independientes \(X_{1}, X_{2}, . . ., X_{P}\).
El modelo general se define como
\(Y_{i} = \beta_{0} + \beta_{1}X_{1} + ...+\beta_{p}X_{p} + \varepsilon_{i}\)
donde \(\beta_{0}\), . . . , \(\beta_{p}\) son los coeficientes de regresión y se asume residuos \(\varepsilon_{i}\) N(0, σ) distribuidos normalmente e independientes alrededor de la media.
En la práctica, tiene n registros de datos; cada registro proporciona valores para cada uno de los predictores X j ; j = 1, . . . .
El modelo a ajustar se da en términos de la respuesta media, condicionada a una realización particular del conjunto de variables explicativas
\(\mathbb{E}[Y|X] = \beta_{0} + \beta_{1}X_{1} + ...+\beta_{p}X_{p} + \varepsilon_{i}\)
En la regresión lineal simple, donde solo tiene una variable predictora, el objetivo es encontrar la “línea de mejor ajuste”. La idea de la estimación por mínimos cuadrados para modelos lineales con múltiples predictores independientes sigue la misma motivación. Ahora, sin embargo, en un sentido abstracto, puede pensar en la relación entre la respuesta y los predictores como un plano o superficie multidimensional. Desea encontrar la superficie que mejor se adapte a sus datos multivariados en términos de minimizar la distancia cuadrática general entre ella y los datos de respuesta sin procesar.
Más formalmente, para sus n registros de datos, los βˆj s se encuentran como los valores que minimizan la suma norte
\(\sum_{i}^{n}{y_{i}-(\beta_{0} + \beta_{1}X_{1} + ...+\beta_{p}X_{p} + \varepsilon_{i})}^{2}\)
Estimación en forma matricial
Los cálculos necesarios para minimizar esta distancia al cuadrado se facilitan mucho con una representación matricial de los datos. Cuando se trata de n observaciones multivariadas, puede escribir la ecuación de la siguiente manera,
\[\begin{bmatrix} \beta _{0}\\ \beta_{1} \\ .\\ .\\ \beta_{P}\\ \end{bmatrix}\] \[\begin{bmatrix} 1&x_{1,1}& . & . & x_{p,1}\\ 1& . & . & . & .\\ .& . & . & . & . \\ .& . & . & .& . \\ 1& x_{1,n}& . & . & x_{p,n} \\ \end{bmatrix}\] = \[\begin{bmatrix} \beta _{0}\\ \beta_{1} \\ .\\ .\\ \beta_{P}\\ \end{bmatrix}\]= \((X^{T}X)^{-1}X^{T}Y\)
Un ejemplo básico
| y | x1 | x2 |
|---|---|---|
| 1.55 | 1.13 | 1 |
| 0.42 | -0.73 | 0 |
| 1.29 | 0.12 | 1 |
| 0.73 | 0.52 | 1 |
| 0.76 | -0.54 | 0 |
| -1.09 | -1.15 | 1 |
| 1.41 | 0.20 | 0 |
| -0.32 | -1.09 | 1 |
| 1.55 |
| 0.42 |
| 1.29 |
| 0.73 |
| 0.76 |
| -1.09 |
| 1.41 |
| -0.32 |
| 1 | 1.13 | 1 |
| 1 | -0.73 | 0 |
| 1 | 0.12 | 1 |
| 1 | 0.52 | 1 |
| 1 | -0.54 | 0 |
| 1 | -1.15 | 1 |
| 1 | 0.20 | 0 |
| 1 | -1.09 | 1 |
| 1.2254572 |
| 1.0153004 |
| -0.6980189 |
Implementando en R e Interpretando
Siempre R construye automáticamente las matrices y lleva a cabo todos los cálculos necesarios cuando se le indica que se ajuste a un modelo de regresión lineal múltiple. Al igual que en los modelos de regresión simple, usa lm y solo incluye cualquier predictor adicional cuando se ñespecifica en la fórmula con el primer argumento. Por el momento se concentrará solo en los efectos principales y luego explorará relaciones más complejas.
Predictores adicionales
Cuando se trata de salida e interpretación, se puede trabajar con múltiples variables explicativas que siguen las mismas reglas de una regresión simple. Cualquier variable numérica continua (o una variable categórica que se trate como tal) tiene un coeficiente de pendiente que proporciona un cambio de unidad en la cantidad del desenlace. Cualquier variable categórica del grupo k (factores, formalmente desordenados) se codifica de forma ficticia y proporciona k − 1 intersectos. Primero confirmemos los cálculos matriciales manuales de hace un momento. Usando el dataframe demo.data, ajuste el modelo lineal múltiple y examine los coeficientes de ese objeto de la siguiente manera:
## (Intercept) x1 x2
## 1.2254572 1.0153004 -0.6980189Verá que obtiene exactamente las estimaciones puntuales almacenadas anteriormente en BETA.HAT, con la variable de respuesta a la izquierda como de costumbre, especifica los predictores en el lado derecho del símbolo ~; en conjunto, esto representa el argumento de la fórmula. Para ajustar un modelo con varios efectos principales, utilice + para separar las variables que desee incluir.
Para estudiar la interpretación de las estimaciones de los parámetros de un modelo de regresión lineal múltiple, regresemos al conjunto de datos de la encuesta en el paquete MASS. Se exploró varios modelos de regresión lineal simple basados en una variable de respuesta de la estatura del estudiante, así como predictores independientes de diámetro palmary sexo (categórico, k = 2) Descubrió que la envergadura era estadísticamente significativa, y el coeficiente estimado sugería un aumento promedio de aproximadamente 3,12 cm por cada 1 cm de aumento en el diámetro palmar. Cuando observó la misma prueba t usando el sexo como variable explicativa, el modelo también sugirió evidencia en contra de la hipótesis nula, con “ser hombre” agregando alrededor de 13.14 cm a la altura media en comparación con la media de las mujeres (la categoría utilizada). como nivel de referencia).
Lo que esos modelos no pueden decirle es el efecto conjunto del sexo y la envergadura en la predicción de la altura. Si incluye ambos predictores en un modelo lineal múltiple, puede (hasta cierto punto) reducir cualquier confusión que de otro modo podría ocurrir en los ajustes aislados del efecto de cualquiera de los predictores sobre la altura.
##
## Call:
## lm(formula = Height ~ Wr.Hnd + Sex, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.7479 -4.1830 0.7749 4.6665 21.9253
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 137.6870 5.7131 24.100 < 2e-16 ***
## Wr.Hnd 1.5944 0.3229 4.937 1.64e-06 ***
## SexMale 9.4898 1.2287 7.724 5.00e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.987 on 204 degrees of freedom
## (30 observations deleted due to missingness)
## Multiple R-squared: 0.5062, Adjusted R-squared: 0.5014
## F-statistic: 104.6 on 2 and 204 DF, p-value: < 2.2e-16El coeficiente para la envergadura de la palma ahora es solo alrededor de 1,59, casi la mitad de su valor correspondiente (3,12 cm) en la regresión lineal simple para la altura. A pesar de esto, sigue siendo estadísticamente muy significativo en presencia de sexo. El coeficiente de sexo también se ha reducido en magnitud en comparación con su modelo lineal simple y también sigue siendo significativo en presencia de palmo. Interpretarás estas nuevas cifras en un momento.
En cuanto al resto de la salida, el error estándar residual todavía le proporciona una estimación del error estándar del término de ruido aleatorio \(\varepsilon_{i}\), y también se le proporciona un valor R-cuadrado. Cuando se asocia con más de un predictor, este último se denomina formalmente coeficiente de determinación múltiple. El cálculo de este coeficiente, como en la configuración de predictor único, proviene de las correlaciones entre las variables del modelo, es importante tener en cuenta que R-cuadrado todavía representa la proporción de variabilidad en la respuesta que explica la regresión; en este ejemplo, se sitúa alrededor de 0,51.
Puede continuar agregando variables explicativas de la misma manera. Antes habia evaluado la frecuencia de tabaquismo como un predictor categórico independiente para la altura y encontró que esta variable explicativa no proporcionó evidencia estadística de un impacto en la respuesta media. Pero, ¿podría la variable fumar contribuir de manera estadísticamente significativa si se controla por el tamaño de la mano y el sexo?
##
## Call:
## lm(formula = Height ~ Wr.Hnd + Sex + Smoke, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.4869 -4.7617 0.7604 4.3691 22.1237
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 137.4056 6.5444 20.996 < 2e-16 ***
## Wr.Hnd 1.6042 0.3301 4.860 2.36e-06 ***
## SexMale 9.3979 1.2452 7.547 1.51e-12 ***
## SmokeNever -0.0442 2.3135 -0.019 0.985
## SmokeOccas 1.5267 2.8694 0.532 0.595
## SmokeRegul 0.9211 2.9290 0.314 0.753
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.023 on 201 degrees of freedom
## (30 observations deleted due to missingness)
## Multiple R-squared: 0.5085, Adjusted R-squared: 0.4962
## F-statistic: 41.59 on 5 and 201 DF, p-value: < 2.2e-16Dado que es una variable categórica con k > 2 niveles, Smoke tiene un código (con fumadores empedernidos como nivel de referencia predeterminado), lo que le brinda tres pendientes adicionales para los tres niveles que no son de referencia de la variable; el cuarto se incorpora como pendiente de comparación general. En el resumen del ajuste más reciente, puede ver que, si bien la envergadura de la palma y el sexo continúan arrojando valores de p muy pequeños, la frecuencia de fumar no sugiere tal evidencia en contra de las hipótesis de los coeficientes cero.
La variable tabaquismo ha tenido poco efecto sobre los valores de los demás coeficientes en comparación con el modelo anterior en superávit, y el coeficiente R cuadrado de determinación múltiple apenas ha aumentado. Una pregunta que podría hacerse ahora es, si la frecuencia con la que fuma no beneficia su capacidad para predecir la estatura media de manera sustancial, ¿debería eliminar esa variable del modelo por completo? Este es el objetivo principal de la selección de modelos: encontrar el “mejor” modelo para predecir el resultado, sin ajustar uno que sea innecesariamente complejo (al incluir más variables explicativas de las necesarias). Verá algunas formas comunes en que los investigadores intentan lograr esto.
Interpretación de los efectos marginales
En la regresión múltiple, la estimación de cada predictor tiene en cuenta el efecto de todos los demás predictores presentes en el modelo. Por lo tanto, un coeficiente para un predictor específico Z debe interpretarse como el cambio en la respuesta media para un aumento de una unidad en Z, mientras se mantienen constantes todos los demás predictores.
Como ha determinado que la frecuencia con la que fuma todavía parece no tener un impacto perceptible en la estatura media cuando se toma en consideración el sexo y la envergadura de la mano, vuelva a enfocarse en survmult, el modelo que incluye solo las variables explicativas de sexo y envergadura. Tenga en cuenta lo siguiente:
• Para estudiantes del mismo sexo (es decir, centrándose solo en hombres o mujeres), un aumento de 1 cm en la longitud de la mano conduce a un aumento estimado de 1,5944 cm en la altura media. • Para estudiantes de longitud similar, los hombres en promedio serán 9.4898 cm más altos que las mujeres. • La diferencia en los valores de los dos coeficientes predictores estimados al compararlos con sus respectivos ajustes del modelo lineal simple, más el hecho de que ambos continúan indicando evidencia en contra de la hipótesis nula de “ser cero” en el ajuste multivariado, sugiere que la confusión (en términos del efecto de la envergadura de la mano y el sexo sobre la variable de respuesta de la altura) está presente en los modelos de predictor único. El punto final destaca la utilidad general de la regresión múltiple. Muestra que, en este ejemplo, si usa solo modelos predictores únicos, la determinación del impacto “verdadero” que tiene cada variable explicativa en la predicción de la respuesta media es engañosa ya que parte del cambio en la altura está determinado por el sexo, pero algunos también se atribuye al diámetro de la palma. Vale la pena señalar que el coeficiente de determinación para el modelo de superávit es notablemente más alto que la misma cantidad en cualquiera de los modelos de una sola variable, por lo que en realidad está teniendo en cuenta una mayor parte de la variación en la respuesta. mediante el uso de regresión múltiple. El modelo ajustado en sí puede considerarse como “Altura media” = 137,687 + 1,594 × “palma” + 9,49 × “sexo” donde “handspan” es el handspan de escritura proporcionado en centímetros y “sexo” se proporciona como 1 (si es hombre) o 0 (si es mujer).
El intercepto en la línea de base (general) de alrededor de 137,687 cm representa la altura media de una mujer con una envergadura de 0 cm; nuevamente, esto claramente no se puede interpretar directamente en el contexto de la aplicación. Para este tipo de situación, algunos investigadores centran el predictor continuo (o predictores) en cero restando la media muestral de todas las observaciones en ese predictor de cada observación antes de ajustar el modelo. Los datos del predictor centrado se utilizan luego en lugar de los datos originales. El modelo ajustado resultante le permite usar el valor medio del predictor no traducido (en este caso, lapso de tiempo) en lugar de un valor cero para interpretar directamente la estimación del intercepto βˆ0. ### Visualización del modelo lineal múltiple Como se muestra aquí, “ser hombre” simplemente cambia el beta alrededor de 9,49 cm:
## (Intercept) Wr.Hnd SexMale
## 137.686951 1.594446 9.489814## [1] 147.1768Debido a esto, también podrías escribir dos ecuaciones. Aquí está la ecuación para estudiantes mujeres: “Altura media” = 137,687 + 1,594 × “palma” Esta es la ecuación para estudiantes varones: “Altura media” = (137,687 + 9,4898) + 1,594 × “palma” = 147.177 + 1.594 × “palma”
Esto es útil porque le permite visualizar el modelo multivariado de la misma manera que los modelos lineales simples. Este código produce la Figura siguiente
Se dibuja un diagrama de dispersión de las observaciones de altura y diámetro de la palma, divididas por sexo. Luego, abline agrega la línea correspondiente a las mujeres y agrega una segunda correspondiente a los hombres, con base en esas dos ecuaciones.
Aunque esta gráfica puede parecer dos ajustes de modelos lineales simples separados, uno para cada nivel de sexo, es importante reconocer que ese no es el caso. Efectivamente, está mirando una representación de un modelo multivariado en un diagrama bidimensional, donde las estadísticas que determinan el ajuste de las dos líneas visibles se han estimado “conjuntamente”, al considerar ambos predictores.
Búsqueda de intervalos de confianza
Se puede encontrar fácilmente los intervalos de confianza para cualquiera de los parámetros de regresión en modelos de regresión múltiple. Usando survmult2, el objeto del modelo ajustado para la altura del estudiante que incluye el predictor de frecuencia de fumar, el resultado de una llamada a confint se ve así:
| 2.5 % | 97.5 % | |
|---|---|---|
| (Intercept) | 124.5010442 | 150.310074 |
| Wr.Hnd | 0.9534078 | 2.255053 |
| SexMale | 6.9426040 | 11.853129 |
| SmokeNever | -4.6061148 | 4.517705 |
| SmokeOccas | -4.1312384 | 7.184710 |
| SmokeRegul | -4.8543683 | 6.696525 |
Tenga en cuenta que las variables Wr.Hnd y SexMale demostraron ser estadísticamente significativas al nivel del 5 por ciento en el resumen del modelo anterior y que sus niveles de confianza del 95 por ciento no incluyen el valor nulo de cero. Por otro lado, todos los coeficientes de las variables dummyes asociadas con el predictor de frecuencia de tabaquismo no son significativos y sus intervalos de confianza claramente incluyen cero. Esto refleja el hecho de que la variable tabaquismo no se considera, en su conjunto, estadísticamente significativa en este modelo en particular.
Prueba F ómnibus
En el contexto de predictores multinivel, puede pensar en la prueba F ómnibus de manera más general para modelos de regresión múltiple como una prueba con las siguientes hipótesis:
\(H_{0}:\beta_{1} = \beta_{2}= ...=\beta_{p}=0\)
\(H_{A}: Al menos un \beta_{j} \neq 0 (j= 1,2,...,p)\)
La prueba compara efectivamente la cantidad de error atribuida al modelo “nulo” (en otras palabras, uno con una pendiente) con la cantidad de error atribuida a los predictores cuando todos los predictores están presentes. En otras palabras, cuanto más capaces sean los predictores de modelar la respuesta, más error explicarán, brindándole una estadística F más extrema y, por lo tanto, un valor p más pequeño. El resultado único de la prueba es especialmente útil cuando tiene muchas variables explicativas. La prueba funciona igual con la combinación de predictores que tengan en un modelo dado: uno o más variables pueden ser continuas, discretas, binarias y/o categóricas con k > 2 niveles. Cuando se ajustan modelos de regresión múltiple, la salida por sí sola puede tomar tiempo para digerir e interpretar, y se debe tener cuidado para evitar errores de Tipo I (rechazo incorrecto de una hipótesis nula verdadera.
La prueba F ayuda a reducir todo eso, lo que permite concluir lo siguiente:
- Evidencia en contra de \(H_{0}\) si el valor p asociado es más pequeño que el nivel de significancia elegido α, lo que sugiere que su regresión, su combinación de las variables explicativas, hace un trabajo significativamente mejor para predecir la respuesta que si eliminara todas esas variables. predictores.
- No hay evidencia en contra de \(H_{0}\) si el valor p asociado es mayor que α, lo que sugiere que el uso de los predictores no tiene un beneficio tangible en comparación con tener un intercepto solo. La desventaja es que la prueba no le dice cuál de los predictores (o qué subconjunto de ellos) tiene un impacto beneficioso en el ajuste del modelo, ni le dice nada sobre sus coeficientes o errores estándar respectivos. Puede calcular la estadística de la prueba F utilizando el coeficiente de determinación, \(R^{2}\), del modelo de regresión ajustado. Sea p el número de parámetros de regresión que requieren estimación, excluyendo el intercepto \(\beta_{0}\).
\(\boldsymbol{F}=\frac{R^{2}(n-p-1)}{(1-R^{2})p}\)
donde n es el número de observaciones utilizadas para ajustar el modelo (después de que se hayan eliminado los registros con valores faltantes). Luego, bajo \(H_{0}\), sigue una distribución F con df1 = p, df2 = n- p - 1 grados de libertad. El valor p asociado se obtiene como el área de cola superior de esa distribución F. Como un ejercicio rápido para confirmar esto, vuelva su atención al modelo de regresión múltiple ajustado survmult2, que es el modelo para la estatura de los estudiantes por palma, sexo y condición de fumador de la encuesta. Puede extraer el coeficiente de determinación múltiple del informe resumido.
## [1] 0.508469Esto coincide con el valor múltiple de R-cuadrado. Luego, puede obtener n como el tamaño original del conjunto de datos en la encuesta menos los valores faltantes.
## [1] 207Obtiene p como el número de parámetros de regresión estimados (menos 1 para el intercepto).
## [1] 5A continuación, puede confirmar el valor de n-p-1, que coincide con el resumen salida (201 grados de libertad):
## [1] 201Finalmente, encuentre el estadístico de prueba y puede usar la función pf de la siguiente manera para obtener el valor p correspondiente para la prueba
## [1] 41.58529## [1] 0Puede ver que la prueba F para este ejemplo da un valor p que es tan pequeño que efectivamente es cero. Estos cálculos coinciden completamente con los resultados relevantes informados en la salida de summary(survmult2). Mirando hacia atrás en el ajuste de regresión múltiple de la altura del estudiante basado en palmo, sexo y tabaquismo en survmult2, no sorprende que con dos de los predictores que arrojan valores p pequeños, la prueba F sugiere una fuerte evidencia en contra de \(H_{0}\). Esto destaca la naturaleza “paraguas” de la prueba ómnibus: aunque la variable de frecuencia de tabaquismo en sí misma no parece aportar nada estadísticamente importante, la prueba F para ese modelo aún sugiere que se debe preferir survmult2 a un modelo “sin predictor”, porque tanto la envergadura de la palma como el sexo son importantes.
Predicción a partir de un modelo lineal múltiple
La predicción (o pronóstico) para la regresión múltiple sigue las mismas reglas que para la regresión simple. Es importante recordar que las predicciones puntuales encontradas para un perfil de covariable en particular, la colección de valores predictores para un individuo determinado, están asociadas con la media (o valor esperado) de la respuesta; que los intervalos de confianza proporcionan medidas para las respuestas medias; y que los intervalos de predicción proporcionan medidas para las observaciones sin procesar. También debe considerar el tema de la interpolación (predicciones basadas en valores de x que se encuentran dentro del rango de los datos de covariables observados originalmente) versus la extrapolación (predicción a partir de valores de x que se encuentran fuera del rango de dichos datos).
Como ejemplo, utilizando el modelo ajustado a la altura del estudiante como una función lineal de la longitud de la mano y el sexo (en superávit), puede estimar la altura media de un estudiante varón con una longitud de la mano de escritura de 16,5 cm, junto con un intervalo de confianza.
| fit | lwr | upr |
|---|---|---|
| 173.4851 | 170.9419 | 176.0283 |
De hecho, hay dos alumnas en el conjunto de datos con una palma de 13 cm. Con su conocimiento de subconjuntos de dataframe, puede inspeccionar estos dos registros y seleccionar las tres variables de interés.
| Sex | Wr.Hnd | Height | |
|---|---|---|---|
| 45 | Female | 13 | 180.34 |
| 152 | Female | 13 | 165.00 |
Ahora, la altura de la segunda mujer cae dentro del intervalo de predicción, pero la altura de la primera es significativamente más alta que el límite superior. Es importante darse cuenta de que, técnicamente, nada salió mal aquí en términos de ajuste e interpretación del modelo; todavía es posible que una observación pueda quedar fuera de un intervalo de predicción, incluso un amplio intervalo del 99 por ciento, aunque tal vez sea improbable. Puede haber varias razones para que esto ocurra. Primero, el modelo podría ser inadecuado. Por ejemplo, podrías estar excluyendo predictores importantes en el modelo ajustado y, por lo tanto, tienen menos poder predictivo. En segundo lugar, aunque la predicción está dentro del rango de los datos observados, se ha producido en un extremo del rango, donde es menos fiable porque los datos son relativamente escasos. Tercero, la observación en sí misma puede estar contaminada de alguna manera, tal vez el individuo registrado su longitud palmar puede ser incorrecta, en cuyo caso su observación inválida debe eliminarse antes del ajuste del modelo. Es con este ojo crítico que un buen estadístico evaluará datos y modelos.
Transformando Variables Numéricas
A veces, la función lineal definida estrictamente por la regresión estándar puede ser inadecuada cuando se trata de capturar las relaciones entre una respuesta y las covariables seleccionadas. Podría, por ejemplo, observar la curvatura en un diagrama de dispersión entre dos variables numéricas para las que una línea perfectamente recta no es necesariamente la más adecuada. Hasta cierto punto, el requisito de que sus datos muestren dicho comportamiento lineal para que un modelo de regresión lineal sea apropiado se puede relajar simplemente transformando (normalmente de forma no lineal) ciertas variables antes de que se realice cualquier estimación o ajuste del modelo. La transformación numérica se refiere a la aplicación de una función matemática a sus observaciones numéricas para cambiar su escala. Encontrar la raíz cuadrada de un número y convertir una temperatura de Fahrenheit a Celsius son ejemplos de una transformación numérica. En el contexto de la regresión, la transformación generalmente se aplica solo a variables continuas y se puede realizar de varias maneras. En esta sección, limitará su atención a los ejemplos que utilizan los dos enfoques más comunes: transformaciones polinómicas y logarítmicas. Sin embargo, tenga en cuenta que la idoneidad de los métodos utilizados para transformar las variables y cualquier beneficio de modelado que pueda ocurrir, solo pueden considerarse realmente caso por caso. La transformación en general no representa una solución universal para resolver problemas de no linealidad en las tendencias de sus datos, pero al menos puede mejorar la fidelidad con la que un modelo lineal puede representar esas tendencias.
Polinomio
Supongamos que observa una relación curva en sus datos, de modo que una línea recta no es una opción sensata para modelarlos. En un esfuerzo por ajustar mejor sus datos, se puede aplicar una transformación polinomial o de potencia a una variable predictora específica en su modelo de regresión. Esta es una técnica sencilla que, al permitir la curvatura polinomial en las relaciones, permite que los cambios en ese predictor influyan en la respuesta de formas más complejas de lo que sería posible de otro modo. Esto se logra al incluir términos adicionales en la definición del modelo que representan el impacto de potencias progresivamente más altas de la variable de interés en la respuesta. Para aclarar el concepto de curvatura polinomial, considere la siguiente secuencia entre −4 y 4, así como los vectores simples calculados a partir de ella:
Tomando el valor original de x y se calculan sus funciones específicas.
El vector y, como copia de x, es claramente lineal (en términos técnicos, se trata de un “polinomio de orden 1”). Asignas \(y^{2}\) para tomar un valor adicional al cuadrado de x, proporcionando un comportamiento cuadrático: un polinomio de orden 2. Por último, el vector \(y^{3}\) representa los resultados de una función cúbica de los valores de x, con la inclusión de x elevada a la potencia de un polinomio de orden 3.
Las siguientes tres líneas de código producen, por separado, los gráficos de izquierda a derecha de la figura.
Tal vez de manera un poco más general, digamos que tiene datos para un predictor continuo, X , que desea usar para modelar su respuesta, Y .
Siguiendo la estimación de la forma habitual, linealmente, el modelo simple es \(\hat{y} = \beta_{0} + \beta_{1}x_{i}\); a la tendencia cuadrática en X se puede modelar a través de la regresión múltiple \(\hat{y} = \beta_{0} + \beta_{1}x_{i} + \beta_{2}x^{2}\); una relación cúbica puede ser capturada por \(\hat{y} = \beta_{0} + \beta_{1}x_{i} + \beta_{2}x^{2} + \beta_{3}x^{3}\).
A partir de los gráficos de la figura, una buena manera de interpretar los efectos de incluir estos términos adicionales es la complejidad de las curvas que se pueden capturar. En el orden 1, la relación lineal no permite ninguna curvatura. En el orden 2, una función cuadrática de cualquier variable dada permite una “curva”. En el orden 3, el modelo puede hacer frente a dos curvas en la relación, y esto continúa si continúa agregando términos correspondientes a potencias crecientes de la covariable. Los coeficientes de regresión asociados con estos términos (todos implícitos como 1 en el código que produjo los gráficos anteriores) pueden controlar la apariencia específica (en otras palabras, la fuerza y la dirección) de la curvatura.
Ajuste de una transformación polinomial Vuelva a centrar su atención en el conjunto de datos integrado de mtcars. Considere la variable disp, que describe el volumen de desplazamiento del motor en pulgadas cúbicas, frente a una variable de respuesta de millas por galón. Si examina una gráfica de los datos en la Figura, puede ver que parece haber una curva leve pero notable en la relación entre el desplazamiento y el kilometraje.
¿La línea recta que proporcionaría un modelo de regresión lineal simple es realmente la mejor manera de representar esta relación? Para investigar esto, comience ajustando esa configuración lineal básica.
##
## Call:
## lm(formula = mpg ~ disp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8922 -2.2022 -0.9631 1.6272 7.2305
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.599855 1.229720 24.070 < 2e-16 ***
## disp -0.041215 0.004712 -8.747 9.38e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.251 on 30 degrees of freedom
## Multiple R-squared: 0.7183, Adjusted R-squared: 0.709
## F-statistic: 76.51 on 1 and 30 DF, p-value: 9.38e-10Es clara la evidencia estadística de un impacto lineal negativo del desplazamiento en el kilometraje: por cada pulgada cúbica adicional de desplazamiento, la respuesta media disminuye aproximadamente 0,041 millas por galón. Ahora, intente capturar la curva aparente en los datos agregando un término cuadrático en disp al modelo. Puede hacer esto de dos maneras. Primero, podría crear un nuevo vector en el espacio de trabajo simplemente elevando al cuadrado el vector mtcars$disp y luego proporcionando el resultado a la fórmula en lm. En segundo lugar, puede especificar \(disp^{2}\) directamente como un término aditivo en la fórmula. Si lo hace de esta manera, es esencial envolver esa expresión en particular en una llamada a I de la siguiente manera:
##
## Call:
## lm(formula = mpg ~ disp + I(disp^2), data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9112 -1.5269 -0.3124 1.3489 5.3946
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.583e+01 2.209e+00 16.221 4.39e-16 ***
## disp -1.053e-01 2.028e-02 -5.192 1.49e-05 ***
## I(disp^2) 1.255e-04 3.891e-05 3.226 0.0031 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.837 on 29 degrees of freedom
## Multiple R-squared: 0.7927, Adjusted R-squared: 0.7784
## F-statistic: 55.46 on 2 and 29 DF, p-value: 1.229e-10El uso de la función I alrededor de un término dado en la fórmula es necesario cuando dicho término requiere un cálculo aritmético, en este caso, \(disp^{2}\), antes de que el modelo en sí se ajuste. En cuanto al propio modelo de regresión múltiple ajustado, puede ver que la contribución del componente cuadrático es estadísticamente significativa: el resultado correspondiente a I(disp^2) muestra un valor p de 0,0031. Esto implica que incluso si se tiene en cuenta una tendencia lineal, el modelo que incluye un componente cuadrático (que introduce una curva) es un modelo de mejor ajuste. Esta conclusión está respaldada por un coeficiente de determinación notablemente más alto en comparación con el primer ajuste (0,7927 frente a 0,7183). Puede ver el ajuste de esta curva cuadrática en la figura 21-4 (cuyo código sigue a continuación).
Aquí podría preguntarse razonablemente si puede mejorar la capacidad del modelo para capturar aún más la relación agregando otro término de orden superior en la covariable de interés. Con ese fin:
##
## Call:
## lm(formula = mpg ~ disp + I(disp^2) + I(disp^3), data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0896 -1.5653 -0.3619 1.4368 4.7617
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.070e+01 3.809e+00 13.310 1.25e-13 ***
## disp -3.372e-01 5.526e-02 -6.102 1.39e-06 ***
## I(disp^2) 1.109e-03 2.265e-04 4.897 3.68e-05 ***
## I(disp^3) -1.217e-06 2.776e-07 -4.382 0.00015 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.224 on 28 degrees of freedom
## Multiple R-squared: 0.8771, Adjusted R-squared: 0.8639
## F-statistic: 66.58 on 3 and 28 DF, p-value: 7.347e-13El resultado muestra que un componente cúbico también ofrece una contribución estadísticamente significativa. Entonces, siendo \(\hat{y}\) las millas por galón y x el desplazamiento en pulgadas cúbicas, y expandiendo la notación e del resultado anterior, el modelo de regresión múltiple ajustado es
\(\hat{y} = 50.7 − 0.3372x + 0.0011x^{2} − 0.000001x^{3}\)
que es precisamente lo que refleja la línea de orden 3 en el panel izquierdo de la figura.
Trazar el ajuste polinomial
Para abordar el gráfico en sí, visualiza los datos y el primer modelo
(lineal simple) en car.order1 de la forma habitual. Para comenzar la
Figura, ejecute el siguiente código:
Es un poco más difícil agregar la línea correspondiente a cualquiera de los modelos denominados polinómicos ya que abline está equipado para manejar solo tendencias de línea recta. Una forma de hacer esto es hacer uso de predecir para cada valor en una secuencia que represente los valores deseados de la variable explicativa. (Estoy a favor de este enfoque porque también le permite calcular simultáneamente las bandas de confianza y predicción si lo desea). Para agregar la línea solo para el modelo de orden 2, primero cree la secuencia requerida sobre el rango observado de disp.
Aquí, la secuencia se ha ampliado un poco por menos y más 50 para predecir una pequeña cantidad a ambos lados del alcance de los datos de la covariable original, por lo que la curva se encuentra con los bordes del gráfico. Luego haces la predicción en sí y superpones la línea ajustada.
Usas la misma técnica, seguida de la suma final de la leyenda, para el polinomio de orden 3.
El resultado de todo esto está en el panel izquierdo de la Figura. Aunque haya utilizado datos sin procesar de una sola covariable, disp, el ejemplo ilustrado aquí se considera una regresión múltiple porque más de un parámetro (además del intercepto universal β0) requirió estimación en los modelos de orden 2 y 3. Los diferentes tipos de líneas de tendencia ajustadas a los datos de kilometraje y desplazamiento muestran claramente diferentes interpretaciones de la relación. Visualmente, podría argumentar razonablemente que el ajuste lineal simple es inadecuado para modelar la relación entre la respuesta y el predictor, pero es más difícil llegar a una conclusión clara al elegir entre las versiones de orden 2 y orden 3. El ajuste de orden 2 captura la curva que disminuye a medida que aumenta la disp; el ajuste de orden 3 permite además un bache (en términos técnicos, una silla de montar o una inflexión), seguido de una tendencia descendente más pronunciada en el mismo dominio. Entonces, ¿qué modelo es el “mejor”? En este caso, la significancia estadística de los parámetros sugiere que se debe preferir el modelo de orden 3. Habiendo dicho eso, hay otras cosas a considerar al elegir entre diferentes modelos.
Trampas de polinomios
Un inconveniente particular asociado con los términos polinómicos en los modelos de regresión lineal es la inestabilidad de la tendencia ajustada cuando se intenta realizar cualquier tipo de extrapolación. El diagrama de la derecha en la figura muestra los mismos tres modelos ajustados (MPG por desplazamiento), pero esta vez con una escala de desplazamiento mucho más amplia. Como puede ver, la validez de estos modelos es cuestionable. Aunque los modelos de orden 2 y 3 se ajustan aceptablemente a MPG dentro del rango de los datos observados, si se mueve aunque sea ligeramente fuera del umbral máximo de los valores de desplazamiento observados, las predicciones del kilometraje medio se desvían enormemente. El modelo de orden 2 en particular se vuelve completamente absurdo, lo que sugiere una rápida mejora en MPG una vez que la cilindrada del motor supera las 500 pulgadas cúbicas. Debe tener en cuenta este comportamiento matemático natural de las funciones polinómicas si está considerando usar términos de orden superior en sus modelos de regresión. Para crear este gráfico, se puede usar el mismo código que creó el gráfico de la izquierda; simplemente usa xlim para ampliar el rango del eje x y definir el objeto disp.seq en una secuencia correspondientemente más amplia (en este caso, solo configuro xlim=c(10,1000) con coincidencias desde y hasta límites en la creación de disp .seq). NOTA Los modelos como este todavía se conocen como modelos de regresión lineal, lo que puede parecer un poco confuso ya que las tendencias ajustadas para polinomios de orden superior son claramente no lineales. Esto se debe a que la regresión lineal se refiere al hecho de que la función que define la respuesta media es lineal en términos de los parámetros de regresión β0, β1, . . ., βp. Como tal, cualquier transformación aplicada a variables individuales no afecta la linealidad de la función con respecto a los coeficientes mismos.
Logarítmico
En situaciones de modelado estadístico donde tiene observaciones numéricas positivas, es común realizar una transformación logarítmica de los datos para reducir drásticamente el rango general de los datos y traer cambios extremos.
Las observaciones más cercanas a una medida de centralidad. En ese sentido, la transformación a una escala logarítmica puede ayudar a reducir la gravedad de los datos muy sesgados. En el contexto del modelado de regresión, las transformaciones logarítmicas se pueden usar para capturar tendencias donde las curvas aparentes se “aplanan”, sin el mismo tipo de inestabilidad fuera del rango de los datos observados que vio con algunos de los polinomios. Si necesita refrescar su memoria sobre logaritmos; aquí basta con señalar que el logaritmo es la potencia a la que debe elevar un valor base para obtener un valor de x. Por ejemplo, en 35 = 243, el logaritmo es 5 y 3 es la base, expresado como log3 243 = 5. Debido a la ubicuidad de la función exponencial en las distribuciones de probabilidad comunes, los estadísticos trabajan casi exclusivamente con el logaritmo natural (logaritmo en base e). A partir de aquí, asuma que todas las menciones de la transformación logarítmica se refieren al logaritmo natural. Para ilustrar brevemente el comportamiento típico de la transformación logarítmica, observe la Figura, lograda con lo siguiente:
Esto traza el logaritmo de los números enteros del 1 al 1000 frente a los valores sin procesar, así como el logaritmo negativo. Puede ver la forma en que los valores transformados logarítmicamente disminuyen y se aplanan a medida que aumentan los valores brutos.
Ajuste de la transformación de registro
Como se señaló, un uso de la transformación logarítmica en la regresión es permitir este tipo de curvatura en situaciones en las que una línea perfectamente recta no se ajusta a la relación observada. Como ilustración, regrese a los ejemplos de mtcars y considere el kilometraje como una función tanto de la potencia como del tipo de transmisión (variables hp y am, respectivamente). Cree un diagrama de dispersión de MPG contra caballos de fuerza, con diferentes colores que distingan entre autos automáticos y manuales.
Los puntos trazados que se muestran en la figura sugieren que las tendencias curvas en los caballos de fuerza pueden ser más apropiadas que las relaciones en línea recta. Tenga en cuenta que tiene que forzar explícitamente el vector numérico binario mtcars$am a un factor aquí para usarlo como selector para el vector de dos colores. Agregará las líneas después de ajustar el modelo lineal.
Hagámoslo usando la transformación logarítmica de caballos de fuerza para tratar de capturar la relación curva. Dado que, en este ejemplo, también desea tener en cuenta el potencial del tipo de transmisión para afectar la respuesta, esto se incluye como una variable de predicción adicional como de costumbre.
##
## Call:
## lm(formula = mpg ~ log(hp) + am, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9084 -1.7692 -0.1432 1.4032 6.3865
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.4842 5.2697 12.047 8.24e-13 ***
## log(hp) -9.2383 1.0439 -8.850 9.78e-10 ***
## am 4.2025 0.9942 4.227 0.000215 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.592 on 29 degrees of freedom
## Multiple R-squared: 0.827, Adjusted R-squared: 0.8151
## F-statistic: 69.31 on 2 and 29 DF, p-value: 8.949e-12El resultado indica conjuntamente los efectos estadísticamente significativos tanto del logaritmo de caballos de fuerza como del tipo de transmisión en el kilometraje. Manteniendo la transmisión constante, el promedio de MPG cae alrededor de 9.24 por cada unidad adicional de log-caballos de fuerza. Tener una transmisión manual aumenta el MPG medio en aproximadamente 4,2 (estimado en este orden debido a la codificación de am—0 para automático, 1 para manual; ver ?mtcars). El coeficiente de determinación muestra que el 82,7 por ciento de la variación en la respuesta se explica por esta regresión, lo que sugiere un ajuste satisfactorio.
Trazado del ajuste de transformación logarítmica
Para visualizar el modelo ajustado, primero debe calcular los valores ajustados para todos los valores predictores deseados. El siguiente código crea una secuencia de valores de caballos de fuerza (menos y más 20 caballos de fuerza) y realiza la predicción requerida para ambos tipos de transmisión.
En el código anterior, dado que desea trazar predicciones para ambos valores posibles de am, al usar newdata necesita replicar hp.seq dos veces. Luego, cuando proporciona valores para am a newdata, una serie de hp.seq se empareja con un valor am replicado apropiadamente de 0, el otro con 1. El resultado de esto es un vector de predicciones de longitud dos veces mayor que la de hp.seq , car.log.pred, correspondiendo los n primeros elementos a coches automáticos y los n últimos a manuales. Ahora puede agregar estas líneas a la Figura con lo siguiente:
Al examinar el diagrama de dispersión, puede ver que el modelo ajustado
parece hacer un buen trabajo al estimar la relación conjunta entre
potencia/transmisión y MPG. La significancia estadística del tipo de
transmisión en este modelo afecta directamente la diferencia entre las
dos líneas agregadas. Si am no fuera significativo, las líneas estarían
más juntas; en ese caso, el modelo estaría sugiriendo que una curva
sería suficiente para capturar la relación. Como de costumbre, la
extrapolación demasiado fuera del rango de los datos del predictor
observado no es una gran idea, aunque es menos inestable para las
tendencias transformadas logarítmicamente que para las funciones
polinómicas.
Otras transformaciones
La transformación puede involucrar más de una variable del conjunto de datos y tampoco se limita solo a las variables predictoras. En su investigación original sobre los datos de mtcars, Henderson y Velleman (1981) también notaron la presencia de las mismas relaciones curvas que ha descubierto entre la respuesta y variables como la potencia y el desplazamiento. Argumentaron que es preferible usar galones por milla (GPM) en lugar de MPG como variable de respuesta para mejorar la linealidad. Esto implicaría modelar una transformación de MPG, a saber, que GPM = 1/MPG. Los autores también comentaron sobre la influencia limitada que tanto la potencia como el desplazamiento tienen en GPM si el peso del automóvil se incluye en un modelo ajustado, debido a las correlaciones relativamente altas presentes entre estos tres predictores (conocidas como multicolinealidad). Para abordar esto, los autores crearon una nueva variable de predicción calculada como caballos de fuerza divididos por peso. Esto mide, en sus palabras, qué tan “superado” está un automóvil, y procedieron a usar ese nuevo predictor en lugar de solo la potencia o el desplazamiento. Esta es solo una parte de la experimentación que se llevó a cabo en la búsqueda de una forma adecuada de modelar estos datos. Con este fin, independientemente de cómo elija modelar sus propios datos, el objetivo de transformar variables numéricas siempre debe ser ajustar un modelo válido que represente los datos y las relaciones de manera más realista y precisa. Al alcanzar este objetivo, hay mucha libertad en cómo puede transformar las observaciones numéricas en aplicaciones de métodos de regresión.
Concepto y Motivación
Los diagramas como los que se encuentran en la figura se utilizan a menudo para ayudar a explicar el concepto de efectos interactivos. Estos diagramas muestran su valor medio de respuesta, \(\hat{y}\), en el eje vertical, como de costumbre, y un valor de predicción para la variable \(x_{1}\) en el eje horizontal. También muestran una variable categórica binaria \(x_{2}\), que puede ser cero o uno. Estas variables hipotéticas están etiquetadas como tales en las imágenes.
Al estimar modelos de regresión, siempre debe acompañar las interacciones con los efectos principales de los predictores relevantes, por razones de interpretabilidad. Dado que las interacciones en sí mismas se entienden mejor como un aumento de los efectos principales, no tiene sentido eliminar los últimos y dejar los primeros.
Para un buen ejemplo de una interacción, piense en la farmacología. Los efectos interactivos entre medicamentos son relativamente comunes, razón por la cual los profesionales de la salud a menudo preguntan sobre otros medicamentos que podría estar tomando.
Considere las estatinas, medicamentos comúnmente utilizados para reducir el colesterol. A los usuarios de estatinas se les recomienda evitar el jugo de toronja porque contiene compuestos químicos naturales que inhiben la eficacia de la enzima responsable de la correcta metabolización de la droga. Si una persona toma estatinas y no consume pomelo, se esperaría una relación negativa entre el nivel de colesterol y el uso de estatinas (piense en el “uso de estatinas” ya sea como una variable de dosis continua o categórica), a medida que aumenta el uso de estatinas o es afirmativo, el nivel de colesterol disminuye. Por otro lado, para una persona que toma estatinas y consume toronja, la naturaleza de la relación entre el nivel de colesterol y el uso de estatinas podría ser fácilmente diferente: negativo debilitado, neutral o incluso positivo. Si es así, dado que el efecto de las estatinas sobre el colesterol cambia según el valor de otra variable, ya sea que se consuma toronja o no, esto se consideraría una interacción entre esos dos predictores. Pueden ocurrir interacciones entre variables categóricas, variables numéricas capaces, o ambos. Es más común encontrar interacciones bidireccionales, interacciones entre exactamente dos predictores. Los efectos interactivos de tres vías y de orden superior son técnicamente posibles pero menos comunes, en parte porque son difíciles de interpretar en un contexto del mundo real.
Uno Categórico, Uno Continuo
En general, una interacción bidireccional entre un predictor categórico y uno continuo debe entenderse como un cambio en la pendiente del predictor continuo con respecto a los niveles de no referencia del predictor categórico. En presencia de un término para la variable continua, una variable categórica con k niveles tendrá k 1 términos de efecto principal, por lo que habrá más k 1 términos interactivos entre todos los niveles alternativos de la variable categórica y la variable continua.
Las diferentes pendientes de \(x_{1}\) por categoría de \(x_{2}\) para \(\hat{y}\). En tal situación, además de los efectos principales para x1 y x2, habría un término interactivo en el modelo ajustado correspondiente a x2 = 1. Esto define el término aditivo necesario para cambiar la pendiente en x1 para x2 = 0 a la nueva pendiente en x1 para x2 = 1. Por ejemplo, accedamos a un nuevo conjunto de datos. En este paquete, también encontrará el objeto de la diabetes: un conjunto de datos de enfermedades cardiovasculares que detalla las características de 403 afroamericanos (investigado e informado originalmente en Schorling et al., 1997; Willems et al., 1997). Instale faraway si aún no lo ha hecho y cárguelo con la biblioteca (“faraway”). Restrinja su atención al nivel de colesterol total (chol—continuo), la edad del individuo (edad: continuo) y tipo de marco corporal (marco: categórico con k = 3 niveles: “pequeño” como nivel de referencia, “mediano” y “grande”). Observará cómo modelar el colesterol total por edad y estructura corporal. Parece lógico esperar que el colesterol esté relacionado tanto con la edad como con el tipo de cuerpo, por lo que también tiene sentido considerar la posibilidad de que el efecto de la edad sobre el colesterol sea diferente para individuos con diferentes estructuras corporales. Para investigar, ajustemos la regresión lineal múltiple e incluyamos una interacción bidireccional entre las dos variables. En la llamada a lm, primero especifica los efectos principales, usando + como de costumbre, y luego especifica un efecto interactivo de dos predictores usando dos puntos (:) entre ellos.
##
## Call:
## lm(formula = chol ~ age + frame + age:frame, data = diabetes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -131.90 -26.24 -5.33 22.17 226.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 155.9636 12.0697 12.922 < 2e-16 ***
## age 0.9852 0.2687 3.667 0.00028 ***
## framemedium 28.6051 15.5503 1.840 0.06661 .
## framelarge 44.9474 18.9842 2.368 0.01840 *
## age:framemedium -0.3514 0.3370 -1.043 0.29768
## age:framelarge -0.8511 0.3779 -2.252 0.02490 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 42.34 on 384 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.07891, Adjusted R-squared: 0.06692
## F-statistic: 6.58 on 5 and 384 DF, p-value: 6.849e-06Al inspeccionar los parámetros del modelo estimado, puede verse un coeficiente de efecto principal para la edad, coeficientes de efecto principal para los dos niveles de la dummy (que no son el nivel de referencia) y dos términos más para el efecto interactivo de la edad. con esos mismos niveles de no referencia.
En realidad, hay un atajo para hacer esto en R: la notación de factor cruzado. El mismo modelo mostrado anteriormente podría haberse ajustado utilizando chol~ageframe en lm; el símbolo entre dos variables en una fórmula debe interpretarse como “incluye una intersección, todos los efectos principales y la interacción”. Usaré este atajo de ahora en adelante. El resultado muestra la importancia de la edad y alguna evidencia para apoyar la presencia de un efecto principal. También hay una ligera indicación de la importancia de la interacción, aunque es débil. Evaluar la significación en este caso, donde un predictor es categórico con k > 2 niveles, sigue la misma regla que se señaló en la discusión de las variables multinivel en la Sección 20.5.2: si al menos uno de los coeficientes es significativo, el efecto total debe ser considerado significativo. La ecuación general para el modelo ajustado se puede escribir directamente a partir de la salida. “Costerol total medio” = 155,9636 + 0,9852 × “edad” + 28.6051 × “marco medio” + 44,9474 × “marco grande” − 0,3514 × “edad: marco medio” − 0,8511 × “edad : complexión grande” (21,7) Se ha usado dos puntos (:) para indicar los términos interactivos para reflejar la salida R. Para el nivel de referencia del predictor categórico, tipo de cuerpo “pequeño”, el modelo ajustado se puede escribir directamente desde la salida. “Costerol total medio” = 155,9636 + 0,9852 × “edad” Para un modelo con solo los efectos principales, cambiar el tipo de cuerpo a “mediano” o “grande” afectaría solo a la intercepción; el efecto relevante simplemente se agrega al resultado. La presencia de la interacción, sin embargo, significa que además del cambio en el intercepto, la pendiente del efecto principal de la edad ahora también debe cambiarse de acuerdo con el término interactivo relevante. Para un individuo con un marco “medio”, el modelo es “Costerol total medio” = 155,9636 + 0,9852 × “edad” + 28,6051 − 0,3514 × “edad” = 184.5687 + (0.9852 − 0.3514) × “edad” = 184,5687 + 0,6338 × “edad” y para un individuo con un marco “grande”, el modelo es “Costerol total medio” = 155,9636 + 0,9852 × “edad” + 44,9474 − 0,8511 × “edad” = 200.911 + (0.9852 − 0.8511) × “edad” Puede calcularlos fácilmente en R accediendo a los coeficientes del objeto del modelo ajustado:
## (Intercept) age framemedium framelarge age:framemedium
## 155.9635868 0.9852028 28.6051035 44.9474105 -0.3513906
## age:framelarge
## -0.8510549A continuación, sumemos los componentes relevantes de este vector. Una vez que tenga las sumas, podrá trazar el modelo ajustado.
## (Intercept) age
## 155.9635868 0.9852028## (Intercept) age
## 184.5686904 0.6338122## (Intercept) age
## 200.9109973 0.1341479Las tres líneas se almacenan como vectores numéricos de longitud 2, con el intercepto y la pendiente. Esta es la forma requerida por el argumento coef opcional de abline, que le permite superponer estas líneas rectas en una gráfica de los datos sin procesar.
Si examina el modelo ajustado en la figura, está claro que la inclusión de una interaccion entre la edad y la estructura corporal ha permitido una mayor flexibilidad en la forma en que el colesterol total medio se relaciona con los dos predictores. La naturaleza no paralela de las tres líneas trazadas refleja el concepto. Revisé esto para ilustrar cómo funciona el concepto, pero en la práctica no es necesario seguir todos estos pasos para encontrar las estimaciones puntuales (y los intervalos de confianza asociados). Puede predecir a partir de un modelo lineal ajustado con interacciones de la misma manera que para los modelos de solo efecto principal mediante el uso de predict.
Dos categóricas
Descubrió evidencia de un efecto interactivo del tipo de lana y la tensión en el número medio de roturas de urdimbre en longitudes de hilo (basado en el marco de datos de roturas de urdimbre listas para usar). Implementemos el mismo modelo que el último ejemplo de warpbreaks.
##
## Call:
## lm(formula = breaks ~ wool * tension, data = warpbreaks)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.5556 -6.8889 -0.6667 7.1944 25.4444
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.556 3.647 12.218 2.43e-16 ***
## woolB -16.333 5.157 -3.167 0.002677 **
## tensionM -20.556 5.157 -3.986 0.000228 ***
## tensionH -20.000 5.157 -3.878 0.000320 ***
## woolB:tensionM 21.111 7.294 2.895 0.005698 **
## woolB:tensionH 10.556 7.294 1.447 0.154327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.94 on 48 degrees of freedom
## Multiple R-squared: 0.3778, Adjusted R-squared: 0.3129
## F-statistic: 5.828 on 5 and 48 DF, p-value: 0.0002772El símbolo de factor cruzado *, en lugar de lana + tensión + lana: tensión. Cuando ambos predictores en una interacción bidireccional son categóricos, habrá un término para cada nivel de no referencia del primer predictor combinado con todos los niveles de no referencia del segundo predictor. En este ejemplo, la lana es binaria con solo k = 2 niveles y la tensión tiene k = 3; por lo tanto, los únicos términos de interacción presentes son los niveles de tensión “medio” (M) y “alto” (H) (“bajo”, L, es el nivel de referencia) con lana tipo B (A es el nivel de referencia). Por lo tanto, en total en el modelo ajustado, hay términos para B, M, H, B:M y B:H. Existe evidencia estadística de un efecto interactivo entre el tipo de lana y la tensión en las roturas medias, además de los efectos principales contribuyentes de esos predictores. El modelo ajustado general puede entenderse como “Roturas de urdimbre medias” = 44.556 − 16.333 × “lana tipo B” − 20.556 × “media tensión” − 20.000 × “alta tensión” + 21.111 × “lana tipo B: tensión media” + 10.556 × “lana tipo B: alta tensión”
Los términos de interacción adicionales funcionan de la misma manera que los efectos principales: cuando solo están involucrados los predictores categóricos, el modelo puede verse como una serie de términos aditivos a la intercepción general. Exactamente cuáles usa en cualquier predicción determinada depende del perfil de covariable de un individuo determinado.
Una serie rápida de ejemplos: para la lana A a baja tensión, el número medio de roturas de urdimbre se predice simplemente como la intersección general; para la lana A a alta tensión, tiene el intercepto general y el término de efecto principal para alta tensión; para lana B a baja tensión, tiene la intercepción general y el efecto principal solo para lana tipo B; y para lana B en tensión media, tiene la intercepción general, el efecto principal para lana tipo B, el efecto principal para tensión media y el término interactivo para lana B con tensión media. Puede usar predict para estimar las roturas de deformación media para estos cuatro escenarios; están acompañados aquí con intervalos de confianza del 90 por ciento:
| fit | lwr | upr |
|---|---|---|
| 44.55556 | 38.43912 | 50.67199 |
| 24.55556 | 18.43912 | 30.67199 |
| 28.22222 | 22.10579 | 34.33866 |
| 28.77778 | 22.66134 | 34.89421 |
Dos continuos
Finalmente, la situación en la que los dos predictores son continuos. En este caso, un término de interacción opera como un modificador en el plano continuo que se ajusta usando solo los efectos principales. De manera similar a una interacción entre un predictor continuo y uno categórico, una interacción entre dos variables explicativas continuas permite que la pendiente asociada a una variable se vea afectada, pero esta vez, esa modificación se realiza de forma continua (es decir, , según el valor de la otra variable continua). Volviendo al dataframe mtcars, considere MPG una vez más como una función de la potencia y el peso. El modelo ajustado, que se muestra a continuación, incluye la interacción además de los efectos principales de los dos predictores continuos. Como puede ver, hay un solo término interactivo estimado y se considera significativamente diferente de cero.
##
## Call:
## lm(formula = mpg ~ hp * wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0632 -1.6491 -0.7362 1.4211 4.5513
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.80842 3.60516 13.816 5.01e-14 ***
## hp -0.12010 0.02470 -4.863 4.04e-05 ***
## wt -8.21662 1.26971 -6.471 5.20e-07 ***
## hp:wt 0.02785 0.00742 3.753 0.000811 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.153 on 28 degrees of freedom
## Multiple R-squared: 0.8848, Adjusted R-squared: 0.8724
## F-statistic: 71.66 on 3 and 28 DF, p-value: 2.981e-13El modelo se escribe como “MPG medio” = 49,80842 − 0,12010 × “caballos de fuerza” − 8,21662 × “peso” + 0.02785 × “caballos de fuerza: peso” = 49,80842 − 0,12010 × “caballos de fuerza” − 8,21662 × “peso” + 0,02785 × “caballos de fuerza” × “peso”
La segunda versión de la ecuación del modelo proporcionada aquí revela por primera vez una interacción expresada como el producto de los valores de los dos predictores, que es exactamente cómo se usa el modelo ajustado para predecir la respuesta. (Técnicamente, esto es lo mismo que cuando al menos uno de los predictores es categórico, pero la codificación dummy simplemente da como resultado ceros y unos para los términos respectivos, por lo que la multiplicación solo equivale a la presencia o ausencia de un término dado)
Puede interpretar una interacción entre dos predictores continuos considerando el signo (+ o ) del coeficiente. La negatividad sugiere que a medida que aumentan los valores de los predictores, la respuesta se reduce después de calcular el resultado de los efectos principales. La positividad, como es el caso aquí, sugiere que a medida que aumentan los valores de los predictores, el efecto es un aumento adicional, una amplificación, en la respuesta media. Contextualmente, los principales efectos negativos de hp y wt indican que el kilometraje se reduce naturalmente para automóviles más pesados y potentes. Sin embargo, la positividad del efecto interactivo sugiere que este impacto en la respuesta se “suaviza” a medida que aumenta la potencia o el peso. Para decirlo de otra manera, la relación negativa impartida por los efectos principales se vuelve “menos extrema” a medida que los valores de los predictores se hacen más y más grandes. La contrasta la versión del modelo de efectos principales solamente (obtenida usando lm con la fórmula mpg~hp+wt; no ajustada explícitamente en esta sección) con la versión de interacción del modelo ajustada justo arriba como el objeto car. adaptar.
Las superficies de respuesta muestran el promedio de MPG en el eje z vertical y las dos variables predictoras en los ejes horizontales como están marcadas. Puede interpretar el MPG medio pronosticado, basado en una potencia y un valor de peso dados, como un punto en la superficie. Tenga en cuenta que ambas superficies disminuyen en MPG (verticalmente a lo largo del eje z) a medida que avanza a valores más grandes de cualquiera de los predictores a lo largo de los respectivos ejes horizontales.
Interacciones de orden superior
Como se mencionó, las interacciones bidireccionales son el tipo más común de interacciones que encontrará en las aplicaciones de los métodos de regresión. Esto se debe a que para los términos de orden triple o superior, se necesitan muchos más datos para una estimación confiable de los efectos interactivos, y hay una serie de complejidades interpretativas que superar. Las interacciones de tres vías son mucho más raras que los efectos de dos vías, y las de cuatro vías y superiores son aún más raras.
En un conjunto de datos nucleares que se encuentra en el paquete de inicio que incluye datos sobre las construcciones de plantas de energía nuclear en los Estados Unidos. En los ejercicios, se concentró principalmente en los predictores de fecha y hora relacionados con los permisos de construcción para modelar el costo medio de construcción de las plantas de energía nuclear. Por el bien de este ejemplo, suponga que no tiene los datos sobre estos predictores. ¿Se puede modelar adecuadamente el costo de construcción usando solo las variables que describen las características de la planta en sí? Cargue el paquete de arranque y acceda a la página de ayuda de ?nuclear para encontrar detalles sobre las variables: cap (variable continua que describe la capacidad de la planta); cum.n (tratado como continuo, que describe el número de construcciones similares en las que los ingenieros habían trabajado previamente); ne (binario, que describe si la planta estaba en el noreste de los Estados Unidos); y ct (binario, que describe si la planta tenía una torre de enfriamiento). El siguiente modelo se ajusta con el costo final de construcción de la planta como respuesta; un efecto principal para la capacidad; y efectos principales de, y todas las interacciones de dos vías y la interacción de tres vías entre, cum.n, ne y ct:
##
## Call:
## lm(formula = cost ~ cap + cum.n * ne * ct, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -162.475 -50.368 -8.833 43.370 213.131
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 138.0336 99.9599 1.381 0.180585
## cap 0.5085 0.1127 4.513 0.000157 ***
## cum.n -24.2433 6.7874 -3.572 0.001618 **
## ne -260.1036 164.7650 -1.579 0.128076
## ct -187.4904 76.6316 -2.447 0.022480 *
## cum.n:ne 44.0196 12.2880 3.582 0.001577 **
## cum.n:ct 35.1687 8.0660 4.360 0.000229 ***
## ne:ct 524.1194 200.9567 2.608 0.015721 *
## cum.n:ne:ct -64.4444 18.0213 -3.576 0.001601 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 107.3 on 23 degrees of freedom
## Multiple R-squared: 0.705, Adjusted R-squared: 0.6024
## F-statistic: 6.872 on 8 and 23 DF, p-value: 0.0001264En este código, usted especifica las interacciones de orden superior ampliando el número de variables conectadas con un * (utilizando * en lugar de : ya que también desea incluir todos los efectos de orden inferior de esos tres predictores). En los resultados estimados, el principal efecto para cap es positivo, mostrando que una mayor capacidad de potencia está ligada a un mayor costo de construcción. Todos los demás efectos principales son negativos, lo que a primera vista parece implicar que un costo de construcción reducido está asociado con ingenieros más experimentados, plantas construidas en el noreste y plantas con una torre de enfriamiento. Sin embargo, esta no es una declaración precisa ya que aún no ha considerado los términos interactivos en esos predictores. Todos los efectos interactivos bidireccionales estimados son positivos: tener ingenieros más experimentados significa un costo de construcción más alto en el noreste independientemente de si hay una torre de enfriamiento, y tener ingenieros más experimentados también significa costos más altos para plantas con una torre de enfriamiento, independientemente de región. El costo también aumenta drásticamente para las plantas en el noreste con una torre de enfriamiento, independientemente de la experiencia del ingeniero. Habiendo dicho todo eso, la interacción negativa de tres vías sugiere que el mayor costo asociado con ingenieros más experimentados que trabajan en el noreste y en una planta con una torre de enfriamiento se reduce un poco después de que se calculan los efectos principales y los efectos interactivos de dos vías. Como mínimo, este ejemplo destaca las complejidades asociadas con la interpretación de los coeficientes del modelo para interacciones de orden superior. También es posible que surjan interacciones de alto orden estadísticamente significativas debido a variables ocultas que no se han tenido en cuenta, es decir, que las interacciones significativas son una manifestación espuria de patrones en los datos que los términos más simples que involucran esos predictores faltantes podrían explicar igual de bien. (si no mejor). En parte, esto motiva la importancia de la selección adecuada del modelo, que es lo siguiente en la discusión.
SELECCIÓN DEL MODELO LINEAL Y DIAGNÓSTICO
Ahora ha dedicado una buena cantidad de tiempo a muchos aspectos de los modelos de regresión lineal. La pregunta de cómo se pueden usar las herramientas y técnicas formales de R para investigar otros aspectos no menos importantes de la regresión: elegir un modelo apropiado para su análisis y evaluar la validez de las suposiciones que ha hecho.
Bondad de ajuste frente a complejidad
El objetivo general de ajustar cualquier modelo estadístico es representar fielmente los datos y las relaciones que se mantienen dentro de ellos. En general, ajustar modelos estadísticos se reduce a un acto de equilibrio entre dos cosas: bondad de ajuste y complejidad. La bondad de ajuste se refiere al objetivo de obtener un modelo que mejor represente las relaciones entre la respuesta y el predictor (o predictores). La complejidad describe cuán complicado es un modelo; esto siempre está ligado a la cantidad de términos en el modelo que requieren estimación; la inclusión de más predictores y funciones adicionales (como transformaciones polinómicas e interacciones) conduce a un modelo más complejo.
Principio de Parsimonia
Los estadísticos se refieren al acto de equilibrio entre la bondad de ajuste y la complejidad como el principio de parsimonia, donde el objetivo de la selección del modelo asociado es encontrar un modelo que sea lo más simple posible (en otras palabras, con una complejidad relativamente baja) , sin sacrificar demasiado la bondad de ajuste. Diríamos que un modelo que satisface esta noción es un ajuste parsimonioso. A menudo escuchará que los investigadores hablan de elegir el “mejor” modelo; en realidad, se refieren a la idea de parsimonia. Entonces, ¿cómo decide dónde trazar la línea en tal equilibrio? Naturalmente, la significación estadística juega un papel aquí, y la selección del modelo a menudo simplemente se reduce a evaluar la importancia del efecto de los predictores o las funciones de los predictores en la respuesta. En un esfuerzo por impartir cierta objetividad a dicho proceso, puede usar algoritmos de selección sistemática, como los que aprenderá en la Sección 22.2, para decidir entre múltiples variables explicativas y cualquier función asociada.
Directrices generales
Realizar cualquier tipo de selección de modelos o comparar varios modelos entre sí implica la toma de decisiones con respecto a la inclusión de variables predictoras disponibles. Sobre este tema, hay varias pautas que siempre debes seguir. • Primero, es importante recordar que no puede eliminar niveles individuales de un predictor categórico en un modelo dado; esto no tiene sentido. En otras palabras, si uno de los niveles que no son de referencia es estadísticamente significativo pero todos los demás no lo son, debe tratar la variable categórica, como un todo, como si hiciera una contribución estadísticamente significativa a la determinación de la respuesta media. Solo debe considerar la eliminación completa de ese predictor categórico si todos los coeficientes que no son de referencia están asociados con una falta de evidencia. Esto también es válido para los términos interactivos que implican predictores categóricos. • Si una interacción está presente en el modelo ajustado, todas las interacciones de orden inferior y los efectos principales de los predictores relevantes deben permanecer en el modelo. Cuando analicé la interpretación de los efectos interactivos como aumentos de los efectos de orden inferior. Como ejemplo, solo debe considerar eliminar el efecto principal de un predictor si no hay términos de interacción presentes en el modelo ajustado que involucre a ese predictor (incluso si ese efecto principal tiene un valor p alto). • En modelos en los que haya utilizado una transformación polinomial de una determinada variable explicativa, mantenga todos los términos polinómicos de orden inferior en el modelo si el mayor se considera significativo. Un modelo que contenga una transformación polinomial de orden 3 en un predictor, por ejemplo, también debe incluir las transformaciones de orden 1 y orden 2 de esa variable. Esto se debe al comportamiento matemático.
Funciones polinómicas: solo separando explícitamente los efectos lineales, cuadráticos y cúbicos (y así sucesivamente) como términos distintos puede evitar confundir dichos efectos entre sí.
Algoritmos de selección de modelo
El trabajo de un algoritmo de selección de modelos es filtrar a través de las variables explicativas disponibles de alguna manera sistemática para establecer cuáles pueden describir mejor la respuesta en conjunto, en lugar de ajustar modelos examinando combinaciones específicas de predictores de forma aislada. Los algoritmos de selección de modelos pueden ser controvertidos. Hay varios métodos diferentes, y ningún enfoque único es universalmente apropiado para cada modelo de regresión. Diferentes algoritmos de selección pueden dar como resultado diferentes modelos finales. En muchos casos, los investigadores tendrán información o conocimiento adicional sobre el problema que influye en la decisión, por ejemplo, que siempre se deben incluir ciertos predictores o que no tiene sentido incluirlos nunca. Esto debe considerarse al mismo tiempo que otras complicaciones, como la posibilidad de interacciones o variables ocultas no observadas que influyan en las relaciones significativas y la necesidad de garantizar que cualquier modelo ajustado sea estadísticamente válido. Es útil tener en cuenta esta famosa cita del célebre estadístico George Box (1919–2013): “Todos los modelos son incorrectos, pero algunos son útiles”. Nunca se puede asumir que cualquier modelo ajustado que produzca sea la verdad, pero un modelo que se ajuste y verifique cuidadosa y minuciosamente puede revelar características interesantes de los datos y, por lo tanto, tiene el potencial de revelar asociaciones y relaciones al proporcionar estimaciones cuantitativas de las mismas.
22.2.1 Comparaciones anidadas: la prueba F parcial La prueba F parcial es probablemente la forma más directa de comparar varios modelos diferentes. Examina dos o más modelos anidados, donde el modelo más pequeño y menos complejo es una versión reducida del modelo más grande y complejo. Formalmente, supongamos que ha ajustado dos modelos de regresión lineal de la siguiente manera:
\(\hat{y}redu = \beta_{0} + \beta_{1}x_{1} + ...+ \beta_{p}x_{p}\)
\(\hat{y}full = \beta_{0} + \beta_{1}x_{1} + ...+ \beta_{p}x_{p} + \beta_{q}x_{q}\)
Aquí, el modelo reducido, que predice \(\hat{y}redu\), tiene p predictores, más una intercepto. El modelo completo, que predice \(\hat{y}full\), tiene q términos predictores. La notación implica que q > p y que, junto con la inclusión estándar de un intercepto \(\beta_{0}\), el modelo completo incluye todos los p predictores del modelo reducido definido por \(\hat{y}redu\), así como q p términos adicionales. Esto enfatiza el hecho de que el modelo para \(\hat{y}redu\) está anidado dentro de \(\hat{y}full\).
Es importante tener en cuenta que aumentar la cantidad de predictores en un modelo de regresión siempre mejorará el \(R^{2}\) y cualquier otra medida de bondad de ajuste. La verdadera pregunta, sin embargo, es si esa mejora en la bondad de ajuste es lo suficientemente grande como para hacer que la complejidad adicional involucrada con la inclusión de cualquier término predictivo adicional “valga la pena”. Esta es precisamente la pregunta que intenta responder la prueba F parcial en el contexto de los modelos de regresión anidados. Su objetivo es probar si la inclusión de esos términos q p adicionales, que producen el modelo completo en lugar del modelo reducido, proporciona una mejora estadísticamente significativa en la bondad de ajuste. La prueba F parcial aborda estas hipótesis:
\(H_{0}:\beta_{p+1}=\beta_{p+2}=...=\beta_{q}=0\)
\(H_{A}\): Al menos unos de los \(\beta_{j}\neq 0\) (para j = p,…,q)
El cálculo de la prueba estadística para abordar estas hipótesis sigue las mismas ideas detrás de la prueba F ómnibus producida automáticamente por R al resumir un objeto de modelo lineal ajustado. Indique el coeficiente de determinación para los modelos reducido y completo con \(R^{2}\)redu y \(R^{2}\)full, respectivamente. Si n se refiere al tamaño de la muestra de los datos utilizados para ajustar ambos modelos, el estadístico de prueba dado por
\(F=\frac{R^{2}_{full}-R^{2}_{redu}(n-q-1)}{1-R^{2}_{full}(q-p)}\)
sigue una distribución F con df1 = q-p, df2 = n-q grados de libertad bajo el supuesto de \(H_{0}\). El valor p se encuentra como el área superior de la cola de la forma habitual; cuanto menor es, mayor es la evidencia en contra de la hipótesis nula, que establece que uno o más de los parámetros adicionales no tienen impacto en la variable de respuesta.
Tome los objetos modelo survmult y survmult2 como ejemplo. El modelo survmult tiene como objetivo predecir la estatura media de los estudiantes a partir de la palma de la mano y su longitud y el sexo en función del dataframe de la encuesta del paquete MASS; survmult2 agrega el estado de fumador a estos predictores. Si es necesario. La impresión de los objetos en la pantalla de la consola muestra una vista previa de los dos ajustes y facilita la confirmación de que el modelo más pequeño está realmente anidado dentro del modelo más grande en términos de sus variables explicativas:
##
## Call:
## lm(formula = Height ~ Wr.Hnd + Sex, data = survey)
##
## Coefficients:
## (Intercept) Wr.Hnd SexMale
## 137.687 1.594 9.490##
## Call:
## lm(formula = Height ~ Wr.Hnd + Sex + Smoke, data = survey)
##
## Coefficients:
## (Intercept) Wr.Hnd SexMale SmokeNever SmokeOccas SmokeRegul
## 137.4056 1.6042 9.3979 -0.0442 1.5267 0.9211Una vez que haya ajustado sus modelos anidados, R puede realizar pruebas F parciales utilizando la función anova (las pruebas F parciales se encuentran dentro del conjunto de metodologías de análisis de varianza). Para determinar si agregar fumar como predictor proporciona una mejora estadísticamente significativa en el ajuste, simplemente comience con el modelo reducido y proporcione los objetos del modelo como argumentos.
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 204 | 9959.182 | NA | NA | NA | NA |
| 201 | 9914.306 | 3 | 44.8756 | 0.3032653 | 0.823015 |
La salida proporciona las cantidades asociadas con el cálculo de R2full y el estadístico F, dado en la tabla resultante como F, que es de mayor interés. Usando los valores de p y q de la impresión survmult y survmult2, debería poder confirmar, por ejemplo, los valores de df1 y df2 que aparecen en la segunda fila de la tabla en las columnas Df y Res.Df, respectivamente. El resultado de esta prueba en particular, obtenido de un estadístico de prueba de = 0,3033 asociado con df1 = 3, df2 = 201, es un valor p alto de 0,823, lo que sugiere que no hay evidencia en contra de \(H_{0}\). Esto significa que agregar Smoke al modelo reducido, que incluye solo las variables explicativas Wr.Hnd y Sex, no ofrece una mejora tangible en el ajuste cuando se trata de modelar la altura de los estudiantes. Esa conclusión no es sorprendente, dados los valores p no significativos de todos los niveles de Fumar no de referencia. Así es como se utilizan las pruebas F parciales para la selección de modelos: en el actual Por ejemplo, el modelo reducido sería el ajuste más parsimonioso y se preferiría al modelo completo. Puede realizar comparaciones entre varios modelos anidados en una llamada determinada a anova, lo que puede ser útil para investigar cosas como la inclusión de términos interactivos o la inclusión de transformaciones polinómicas de predictores, ya que existe una jerarquía natural que requiere que conserve cualquier término de orden inferior. .
Por ejemplo, usemos el dataframe de diabetes, con el modelo ajustado para predecir el nivel de colesterol (chol) contra la edad (edad) y el marco corporal (marco) y la interacción entre esos dos predictores. Antes de usar las pruebas F parciales para comparar variantes anidadas, debe asegurarse de que está usando los mismos registros para cada modelo y que no faltan valores para ninguno de los predictores. Para hacer esto, solo necesita definir primero una versión de diabetes que elimine los registros con valores faltantes en los predictores que está utilizando. Cargue el paquete y use subconjuntos lógicos para identificar y eliminar a cualquier persona a la que le falte un valor para la edad o el marco. Defina esta nueva versión del objeto diabetes:
Ahora, ajusta los siguientes cuatro modelos usando tu nuevo objeto diab:
El primer modelo es solo una intersecto, el segundo agrega la edad como predictor, el tercero tiene la edad y el marco, y el cuarto incluye la interacción. El anidamiento es evidente y ahora puede comparar la importancia de las mejoras en la bondad de ajuste a medida que aumenta la complejidad del modelo en cada paso.
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 389 | 747264.8 | NA | NA | NA | NA |
| 388 | 712078.2 | 1 | 35186.579 | 19.630610 | 0.0000123 |
| 386 | 697527.5 | 2 | 14550.793 | 4.058948 | 0.0180134 |
| 384 | 688294.8 | 2 | 9232.677 | 2.575458 | 0.0774334 |
Si no hubiera eliminado los registros que contenían valores faltantes en esos predictores, habría recibido un error que le indicaría que los conjuntos de datos de los cuatro modelos no tenían el mismo tamaño. Los resultados en sí mismos sugieren que incluir la edad proporciona una mejora significativa para modelar chol; incluir un efecto principal para el marco proporciona una mejora leve adicional; y hay evidencia muy débil, si es que hay alguna, de que incluir un efecto interactivo es beneficioso para la bondad de ajuste. A partir de esto, es posible que prefiera usar dia.mod3, el modelo de solo efectos principales, como la representación más parsimoniosa del colesterol medio de estos cuatro modelos.
Selección Forward
Las pruebas F parciales son una forma natural de investigar modelos anidados, pero pueden ser difíciles de administrar si tiene muchos modelos diferentes para ajustar cuando, tiene muchas variables predictoras.
Aquí es donde entra en juego la selección Forward (también conocida como eliminación hacia adelante). La idea es comenzar con un modelo de solo intercepto y luego realizar una serie de pruebas independientes para determinar cuál de sus variables predictoras mejora significativamente la bondad de -adaptar.
Luego actualiza su modelo agregando ese término y ejecuta la serie de pruebas nuevamente para todos los términos restantes para determinar cuál de ellos mejoraría aún más el ajuste.
El proceso se repite hasta que no haya más términos que mejoren el ajuste de forma estadísticamente significativa.
Las funciones de R listas para usar add1 y update realizan la serie de pruebas y actualizan su modelo de regresión ajustado.
Utilizará dataframes nucleares. El objetivo es elegir el modelo más informativo para la predicción del costo de construcción. Cargue el arranque y acceda al archivo de ayuda ?nuclear para recordar las definiciones de las variables. Primero ajuste el modelo para el costo de construcción con un término de intersección general solamente.
##
## Call:
## lm(formula = cost ~ 1, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -254.05 -151.24 -13.46 150.40 419.68
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 461.56 30.07 15.35 4.95e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 170.1 on 31 degrees of freedomSabe por explicaciones anteriores que este modelo en particular es bastante inadecuado para la predicción confiable del costo. Entonces, considere la siguiente línea de código para iniciar la selección hacia adelante (he suprimido la salida, que mostraré por separado y discutiré en un momento):
| Df | Sum of Sq | RSS | AIC | F value | Pr(>F) | |
|---|---|---|---|---|---|---|
| NA | NA | 897172.3 | 329.7206 | NA | NA | |
| date | 1 | 334335.20129 | 562837.1 | 316.8003 | 17.8205309 | 0.0002071 |
| t1 | 1 | 186983.56088 | 710188.7 | 324.2416 | 7.8986141 | 0.0086296 |
| t2 | 1 | 27.38376 | 897144.9 | 331.7196 | 0.0009157 | 0.9760597 |
| cap | 1 | 199673.21054 | 697499.1 | 323.6647 | 8.5881062 | 0.0064137 |
| pr | 1 | 9036.67205 | 888135.6 | 331.3966 | 0.3052464 | 0.5847053 |
| ne | 1 | 128640.95163 | 768531.4 | 326.7680 | 5.0215629 | 0.0325885 |
| ct | 1 | 43042.23829 | 854130.1 | 330.1473 | 1.5117922 | 0.2284221 |
| bw | 1 | 16205.41516 | 880966.9 | 331.1373 | 0.5518510 | 0.4633402 |
| cum.n | 1 | 67938.09175 | 829234.2 | 329.2007 | 2.4578614 | 0.1274266 |
| pt | 1 | 305333.77985 | 591838.5 | 318.4081 | 15.4772171 | 0.0004575 |
El primer argumento para add1 es siempre el modelo que desea actualizar.
El segundo argumento, el alcance, es fundamental: debe proporcionar una fórmula que defina el modelo más “completo” y más complejo que consideraría ajustar.
Para esto, normalmente usaría el .~. notación, en la que los puntos se refieren a la definición del modelo en el primer argumento.
Específicamente, los puntos representan “lo que ya está ahí”. En otras palabras, a través del alcance le está diciendo a add1 que el modelo más completo que consideraría tiene un costo como respuesta, una intercepción y los efectos principales de todos los demás predictores en el marco de datos nucleares (restringiré el modelo completo a principal efectos sólo para facilitar la demostración).
No necesita proporcionar el dataframe como argumento, ya que esos datos están contenidos dentro del objeto modelo en el primer argumento. Por último, le indicas a add1 la prueba que debe realizar. Hay un puñado de variantes disponibles (ver ?add1), pero aquí te quedarás con test=“F” para pruebas F parciales. Ahora, concéntrese en la salida que se proporciona directamente después de la ejecución. de agregar1.
| Df | Sum of Sq | RSS | AIC | F value | Pr(>F) | |
|---|---|---|---|---|---|---|
| NA | NA | 897172.3 | 329.7206 | NA | NA | |
| date | 1 | 334335.20129 | 562837.1 | 316.8003 | 17.8205309 | 0.0002071 |
| t1 | 1 | 186983.56088 | 710188.7 | 324.2416 | 7.8986141 | 0.0086296 |
| t2 | 1 | 27.38376 | 897144.9 | 331.7196 | 0.0009157 | 0.9760597 |
| cap | 1 | 199673.21054 | 697499.1 | 323.6647 | 8.5881062 | 0.0064137 |
| pr | 1 | 9036.67205 | 888135.6 | 331.3966 | 0.3052464 | 0.5847053 |
| ne | 1 | 128640.95163 | 768531.4 | 326.7680 | 5.0215629 | 0.0325885 |
| ct | 1 | 43042.23829 | 854130.1 | 330.1473 | 1.5117922 | 0.2284221 |
| bw | 1 | 16205.41516 | 880966.9 | 331.1373 | 0.5518510 | 0.4633402 |
| cum.n | 1 | 67938.09175 | 829234.2 | 329.2007 | 2.4578614 | 0.1274266 |
| pt | 1 | 305333.77985 | 591838.5 | 318.4081 | 15.4772171 | 0.0004575 |
La salida consta de una serie de filas, comenzando con
##
## Call:
## lm(formula = cost ~ date, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -176.00 -105.27 -25.24 58.63 359.46
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6553.57 1661.96 -3.943 0.000446 ***
## date 102.29 24.23 4.221 0.000207 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 137 on 30 degrees of freedom
## Multiple R-squared: 0.3727, Adjusted R-squared: 0.3517
## F-statistic: 17.82 on 1 and 30 DF, p-value: 0.0002071La actualización, proporciona el modelo que desea actualizar como primer argumento, y el segundo argumento, la fórmula, le dice a la actualización cómo actualizar el modelo. Nuevamente usando el .~. notación, la instrucción es actualizar nuc.0 agregando la fecha como predictor, lo que da como resultado un objeto de modelo ajustado de la misma clase del primer argumento. Llame a un resumen del nuevo modelo, nuc.1, para ver esto. Entonces, ¡sigamos adelante! Vuelva a llamar a add1, pero ahora pase nuc.1 como su primer argumento.
| Df | Sum of Sq | RSS | AIC | F value | Pr(>F) | |
|---|---|---|---|---|---|---|
| NA | NA | 562837.1 | 316.8003 | NA | NA | |
| t1 | 1 | 15321.858 | 547515.2 | 317.9171 | 0.8115461 | 0.3750843 |
| t2 | 1 | 68161.195 | 494675.9 | 314.6695 | 3.9958983 | 0.0550606 |
| cap | 1 | 189731.747 | 373105.4 | 305.6442 | 14.7470963 | 0.0006163 |
| pr | 1 | 4027.160 | 558809.9 | 318.5705 | 0.2089935 | 0.6509638 |
| ne | 1 | 92256.403 | 470580.7 | 313.0716 | 5.6853918 | 0.0238671 |
| ct | 1 | 54793.671 | 508043.4 | 315.5228 | 3.1277177 | 0.0874906 |
| bw | 1 | 1240.053 | 561597.1 | 318.7297 | 0.0640344 | 0.8020147 |
| cum.n | 1 | 4658.212 | 558178.9 | 318.5344 | 0.2420159 | 0.6264574 |
| pt | 1 | 90586.683 | 472250.4 | 313.1849 | 5.5627558 | 0.0252997 |
Tenga en cuenta que ahora no hay una fila para agregar la fecha; ya está allí en nuc.1.
Parece que la siguiente adición más informativa sería cap. Actualice nuc.1 a tal efecto.
Ahora continúa, probando y actualizando.
Al llamar a add1 en nuc.2 (la salida no se muestra aquí), encontrará que la siguiente adición más importante es pt (por un pequeño margen). Actualice a un nuevo objeto llamado nuc.3, que incluye el siguiente término:
Luego pruebe nuevamente, usando add1 en nuc.3. Encontrará evidencia débil para incluir adicionalmente un efecto principal para ne, así que actualice con esa inclusión para crear nuc.4.
En este punto, puede estar razonablemente seguro de que no habrá más adiciones útiles, pero verifique con una última llamada a add1 en el ajuste más reciente para ser exhaustivo.
| Df | Sum of Sq | RSS | AIC | F value | Pr(>F) | |
|---|---|---|---|---|---|---|
| NA | NA | 222616.9 | 293.1191 | NA | NA | |
| t1 | 1 | 107.0452 | 222509.9 | 295.1037 | 0.0125081 | 0.9118095 |
| t2 | 1 | 19229.9190 | 203387.0 | 292.2282 | 2.4582587 | 0.1289998 |
| pr | 1 | 5230.8123 | 217386.1 | 294.3582 | 0.6256201 | 0.4361227 |
| ct | 1 | 15764.6566 | 206852.3 | 292.7688 | 1.9815159 | 0.1710750 |
| bw | 1 | 447.9715 | 222169.0 | 295.0546 | 0.0524252 | 0.8206872 |
| cum.n | 1 | 13819.9367 | 208797.0 | 293.0682 | 1.7208981 | 0.2010431 |
De hecho, parece que ninguna de las covariables restantes, si se incluyen en el modelo, produciría una mejora estadísticamente significativa en la bondad de ajuste, por lo que su modelo final se mantendrá en nuc.4.
##
## Call:
## lm(formula = cost ~ date + cap + pt + ne, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.894 -38.424 -2.493 35.363 267.445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.756e+03 1.286e+03 -3.699 0.000975 ***
## date 7.102e+01 1.867e+01 3.804 0.000741 ***
## cap 4.198e-01 8.616e-02 4.873 4.29e-05 ***
## pt -1.289e+02 4.950e+01 -2.605 0.014761 *
## ne 9.940e+01 3.864e+01 2.573 0.015908 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 90.8 on 27 degrees of freedom
## Multiple R-squared: 0.7519, Adjusted R-squared: 0.7151
## F-statistic: 20.45 on 4 and 27 DF, p-value: 7.507e-08Este método puede parecer un poco engorroso y, a veces, es difícil decidir cuál es el modelo más completo, pero es una excelente manera de mantenerse involucrado en cada etapa del proceso de selección para que pueda considerar cada adición cuidadosamente. Tenga en cuenta, sin embargo, que hay un elemento de subjetividad; es posible llegar a diferentes modelos finales eligiendo una adición sobre otra, como si hubiera agregado pt en lugar de fecha (tenían niveles similares de importancia en la primera llamada a add1).
Selección Backward
Después de aprender la selección Forward, comprender la selección Backward (o eliminación) no es una gran exageración. Como habrás adivinado, donde la selección Forward comienza con un modelo reducido y se abre paso hasta un modelo final agregando términos, la selección Backward comienza con tu modelo más completo y elimina términos sistemáticamente. Las funciones R para este proceso son drop1 para inspeccionar las pruebas F parciales y actualizar. La elección de la selección del modelo Forward o Backward generalmente se realiza caso por caso. Si no se conoce su modelo más completo o es difícil de definir y ajustar, generalmente se prefiere la selección Forward. Por otro lado, si tiene un modelo más completo natural y fácil de ajustar, entonces la selección Backward puede ser más conveniente de implementar. A veces, los investigadores realizarán ambos para ver si el modelo final al que llegan es diferente (una ocurrencia perfectamente posible). Revisa el ejemplo nuclear. Primero, defina el modelo más completo como aquel que predice el costo por los efectos principales de todas las covariables disponibles (como hizo en su uso del alcance en las selecciones Forward). Después de aprender la selección Forward, comprender la selección Backward (o eliminación) no es una gran exageración. Como habrás adivinado, donde la selección Forward comienza con un modelo reducido y se abre paso hasta un modelo final agregando términos, la selección Backward comienza con tu modelo más completo y elimina términos sistemáticamente. Las funciones R para este proceso son drop1 para inspeccionar las pruebas F parciales y actualizar. La elección de la selección del modelo Forward o Backward generalmente se realiza caso por caso. Si no se conoce su modelo más completo o es difícil de definir y ajustar, generalmente se prefiere la selección Forward. Por otro lado, si tiene un modelo más completo natural y fácil de ajustar, entonces la selección Backward puede ser más conveniente de implementar. A veces, los investigadores realizarán ambos para ver si el modelo final al que llegan es diferente (una ocurrencia perfectamente posible). Revisa el ejemplo nuclear. Primero, defina el modelo más completo como aquel que predice el costo por los efectos principales de todas las covariables disponibles (como hizo en su uso del alcance en las selecciones Forward.
##
## Call:
## lm(formula = cost ~ date + t1 + t2 + cap + pr + ne + ct + bw +
## cum.n + pt, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -128.608 -46.736 -2.668 39.782 180.365
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.135e+03 2.788e+03 -2.918 0.008222 **
## date 1.155e+02 4.226e+01 2.733 0.012470 *
## t1 5.928e+00 1.089e+01 0.545 0.591803
## t2 4.571e+00 2.243e+00 2.038 0.054390 .
## cap 4.217e-01 8.844e-02 4.768 0.000104 ***
## pr -8.112e+01 4.077e+01 -1.990 0.059794 .
## ne 1.375e+02 3.869e+01 3.553 0.001883 **
## ct 4.327e+01 3.431e+01 1.261 0.221008
## bw -8.238e+00 5.188e+01 -0.159 0.875354
## cum.n -6.989e+00 3.822e+00 -1.829 0.081698 .
## pt -1.925e+01 6.367e+01 -0.302 0.765401
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 82.83 on 21 degrees of freedom
## Multiple R-squared: 0.8394, Adjusted R-squared: 0.763
## F-statistic: 10.98 on 10 and 21 DF, p-value: 2.844e-06Claramente, hay varios predictores que parecen no contribuir significativamente a la respuesta, y estos mismos resultados son evidentes la primera vez que usa drop1 para examinar el impacto en la bondad de ajuste que ocurriría al eliminar cada variable.
| Df | Sum of Sq | RSS | AIC | F value | Pr(>F) | |
|---|---|---|---|---|---|---|
| NA | NA | 144064.9 | 291.1931 | NA | NA | |
| date | 1 | 51230.0415 | 195294.9 | 298.9290 | 7.4676825 | 0.0124702 |
| t1 | 1 | 2034.2579 | 146099.2 | 289.6418 | 0.2965290 | 0.5918028 |
| t2 | 1 | 28481.3581 | 172546.3 | 294.9659 | 4.1516605 | 0.0543902 |
| cap | 1 | 155942.7313 | 300007.6 | 312.6665 | 22.7314046 | 0.0001039 |
| pr | 1 | 27161.2608 | 171226.2 | 294.7201 | 3.9592330 | 0.0597943 |
| ne | 1 | 86580.9787 | 230645.9 | 304.2529 | 12.6207053 | 0.0018835 |
| ct | 1 | 10915.0306 | 154979.9 | 291.5301 | 1.5910583 | 0.2210075 |
| bw | 1 | 172.9693 | 144237.9 | 289.2315 | 0.0252133 | 0.8753538 |
| cum.n | 1 | 22939.4143 | 167004.3 | 293.9212 | 3.3438244 | 0.0816977 |
| pt | 1 | 626.9005 | 144691.8 | 289.3320 | 0.0913818 | 0.7654015 |
Una característica útil de drop1 es que su argumento de alcance es opcional. Si no incluye el alcance, el modelo predeterminado es solo de intercepto como el modelo “más reducido”, que suele ser una opción razonable.
Antes de sumergirse directamente en el proceso de eliminación, recuerde la interpretación de lo que está haciendo. Así como agregar cualquier término siempre mejorará la bondad de ajuste en la selección hacia adelante, eliminar cualquier término en la selección hacia atrás siempre empeorará la bondad de ajuste.
La verdadera pregunta es la importancia percibida de estos cambios en la calidad del ajuste. De la misma manera que antes, donde deseaba agregar solo aquellos términos que ofrecen una mejora estadísticamente significativa en la bondad de ajuste, al descartar términos, solo desea eliminar aquellos que no dan como resultado un detrimento estadísticamente significativo para la bondad. fuera de él. Como tal, la selección hacia atrás es el reverso completo de la selección hacia adelante en la forma en que se lleva a cabo. Entonces, de la salida de drop1, desea elegir el término para eliminar del modelo que tiene el efecto menos significativo de reducir la bondad del ajuste. En otras palabras, está buscando el término con el valor p no significativo más grande para su prueba F parcial, porque dejar caer un término con un valor p significativamente pequeño empeoraría significativamente la capacidad predictiva del modelo de regresión.
En el ejemplo actual, parece que el predictor bw tiene el efecto menos significativo en la reducción de la bondad de ajuste, así que comencemos la actualización eliminando ese término de nuc.0.
El uso de actualización en este algoritmo de selección es el mismo que antes; ahora, sin embargo, usas un - para significar la eliminación de un término siguiendo el estándar “lo que ya está allí” .~. notación. A continuación, se repite el proceso utilizando el último modelo nuc.1:
| Df | Sum of Sq | RSS | AIC | F value | Pr(>F) | |
|---|---|---|---|---|---|---|
| NA | NA | 144237.9 | 289.2315 | NA | NA | |
| date | 1 | 55941.7146 | 200179.6 | 297.7195 | 8.5325566 | 0.0079132 |
| t1 | 1 | 3123.8417 | 147361.7 | 287.9171 | 0.4764666 | 0.4972454 |
| t2 | 1 | 30717.4985 | 174955.4 | 293.4096 | 4.6852121 | 0.0415458 |
| cap | 1 | 159976.1541 | 304214.0 | 311.1120 | 24.4004962 | 0.0000610 |
| pr | 1 | 27139.5052 | 171377.4 | 292.7484 | 4.1394756 | 0.0541221 |
| ne | 1 | 86408.0127 | 230645.9 | 302.2529 | 13.1794541 | 0.0014792 |
| ct | 1 | 11814.7931 | 156052.7 | 289.7508 | 1.8020612 | 0.1931528 |
| cum.n | 1 | 24048.2963 | 168286.2 | 292.1659 | 3.6679864 | 0.0685573 |
| pt | 1 | 930.1263 | 145168.0 | 287.4372 | 0.1418683 | 0.7100391 |
Parecería que pt es el siguiente efecto principal más sensato a eliminar. Hazlo y nombra el objeto resultante nuc.2.
Ahora continúe, vuelva a verificar con una llamada a drop1 (no se muestra), y encontrará que el predictor t1 se revela como otra eliminación viable. Actualice su modelo con ese predictor eliminado; nombrar el objeto modelo nuc.3.
Vuelva a comprobar el nuevo nuc.3 con drop1. Ahora debería encontrar que el efecto de ct sigue siendo no significativo, así que elimínelo y actualícelo nuevamente, brindándole un nuevo nuc.4.
Realice otra verificación con drop1, esta vez en nuc.4. En este punto, es posible que dude en eliminar más predictores, ya que la importancia en diferentes intensidades está asociada con el efecto de su eliminación. Tenga en cuenta, sin embargo, que para al menos tres de los predictores restantes, t2, pr y cum.n, la significación estadística probablemente debería considerarse en el mejor de los casos en el límite: todos sus valores de p se encuentran entre los niveles de corte convencionales de α = 0.01 y α = 0,05.
Esto nuevamente enfatiza el papel activo que un investigador debe desempeñar en los algoritmos de selección de modelos, como la eliminación hacia adelante o hacia atrás; si debe eliminar más variables de aquí es una pregunta difícil de responder y se deja a su criterio. Quedémonos con nuc.4 como modelo final. En resumen, puede ver los parámetros de regresión estimados y las estadísticas posteriores al ajuste habituales.
##
## Call:
## lm(formula = cost ~ date + t2 + cap + pr + ne + cum.n, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -152.851 -53.929 -8.827 53.382 155.581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.702e+03 1.294e+03 -7.495 7.55e-08 ***
## date 1.396e+02 1.843e+01 7.574 6.27e-08 ***
## t2 4.905e+00 1.827e+00 2.685 0.012685 *
## cap 4.137e-01 8.425e-02 4.911 4.70e-05 ***
## pr -8.851e+01 3.479e+01 -2.544 0.017499 *
## ne 1.502e+02 3.400e+01 4.419 0.000168 ***
## cum.n -7.919e+00 2.871e+00 -2.758 0.010703 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 80.8 on 25 degrees of freedom
## Multiple R-squared: 0.8181, Adjusted R-squared: 0.7744
## F-statistic: 18.74 on 6 and 25 DF, p-value: 3.796e-08Inmediatamente, puede ver que su modelo final de la selección forward es diferente del modelo final seleccionado aquí, a pesar de que el modelo más completo es el mismo para ambos. ¿Cómo ha ocurrido eso? La respuesta, en pocas palabras, es que los predictores presentes en un modelo se afectan entre sí. Recuerde que los coeficientes estimados de los predictores presentes cambian fácilmente de valor a medida que controla las diferentes variables. A medida que aumenta el número de términos predictores, estas relaciones se vuelven más y más complejas, por lo que tanto el orden como la dirección del algoritmo de selección tienen el potencial de llevarlo por diferentes caminos a través del proceso de selección y llegar a diferentes destinos finales. que es exactamente lo que ha sucedido aquí. Como ejemplo perfecto de esto, considere el efecto principal de pt en los datos nucleares. En la selección forward, se agregó pt porque ofrecía la mejora “más significativa” al modelo costo~fecha+límite. En la selección hacia atrás, pt se eliminó antes, ya que ofrecía la menor reducción en la bondad de ajuste si se tomaba del modelo cost~date+t1+t2+cap+pr+ne+ct+cum.n+pt. Lo que esto significa es que para el último modelo, la contribución que pt podría hacer en términos de predicción del resultado ya está explicada por los otros términos predictores presentes. En el modelo más pequeño, ese efecto aún no se había explicado, por lo que pt era una adición atractiva. Todo esto sirve para resaltar la naturaleza voluble de la mayoría de los algoritmos de selección, a pesar de la forma sistemática en que se implementan. Es importante reconocer que el ajuste final del modelo probablemente variará entre los enfoques. y que debería ver estos métodos de selección más como guías útiles para encontrar el modelo más parsimonioso y no como una solución universal y definitiva.
Selección paso a paso de AIC
La aplicación de una serie de pruebas F parciales es el método de selección de modelos basado en pruebas más común, pero no es la única herramienta que un investigador tiene a su disposición. También puede localizar la parsimonia adoptando un enfoque basado en criterios. Una de las medidas de criterio más famosas se conoce como Criterio de información de Akaike (AIC). Habrá notado este valor como una de las columnas en la salida de add1 y drop1. Para un modelo lineal dado, AIC se calcula de la siguiente manera:
\(AIC = −2 × L + 2 × (p + 2)\)
Aquí, hay una medida de bondad de ajuste llamada log-verosimilitud y es L, p es el número de parámetros de regresión en el modelo, excluyendo el intercepto general. El valor es un resultado directo del procedimiento de estimación utilizado para ajustar el modelo, aunque su cálculo exacto está más allá del alcance de este texto. Lo que hay que saber es que toma valores más grandes para modelos que se ajustan mejor.
La ecuación produce un valor que recompensa la bondad de ajuste con el 2 pero simultáneamente penaliza la complejidad con el 2 (p+2). El signo negativo asociado con junto con el signo positivo de p + 2 significa que los valores más pequeños de AIC se refieren a modelos más parsimoniosos.
Para encontrar el AIC para un modelo lineal ajustado, utilice las funciones AIC o extractAIC en el objeto resultante de lm; eche un vistazo a los archivos de ayuda de estas funciones para ver las diferencias técnicas entre los dos.
El valor de (y por lo tanto también el AIC) no tiene una interpretación directa y es útil solo cuando lo compara con el AIC de otro modelo. Puede basar la selección del modelo en el AIC identificando el ajuste con el valor AIC más bajo. Esta es la razón por la que se informa directamente en la salida de add1 y drop1: puede decidir qué término agregar o quitar en función del cambio que resulte en un cambio al AIC más pequeño, en lugar de centrarse exclusivamente en la importancia del cambio a través de la prueba F.
Vayamos aún más lejos y combinemos las ideas de selección hacia adelante y hacia atrás. La selección del modelo paso a paso permite la opción de eliminar un término presente o agregar un término faltante y generalmente se implementa con respecto a la AIC. Es decir, se elige un término para agregarlo o eliminarlo en función del movimiento de todos los movimientos posibles que produce la mayor reducción individual en AIC. Esto le brinda más flexibilidad para explorar modelos candidatos en su camino hacia el ajuste del modelo final, determinado como el modelo del cual ninguna adición o eliminación reduciría aún más el valor de AIC. Es posible implementar la selección de AIC paso a paso usted mismo utilizando add1 o drop1 en cada etapa, pero afortunadamente R proporciona la función de paso incorporada para hacerlo por usted. Tome los datos de mtcars del paquete MASS de los últimos dos capítulos. Intentemos finalmente obtener un modelo para el kilometraje medio que ofrezca la oportunidad de incluir todos los predictores disponibles. Primero, eche un vistazo a la documentación en ?mtcars y una matriz de diagrama de dispersión de los datos nuevamente para recordar las variables y su formato en el objeto de marco de datos R. Luego defina el modelo inicial (a menudo llamado modelo nulo) como el modelo de solo intercepción.
Su modelo inicial puede ser cualquier cosa que desee, siempre que se encuentre dentro del dominio de los modelos descritos por su argumento de alcance que se proporcionará al paso. En este ejemplo, defina el alcance como el modelo más completo a considerar; configúrelo para que sea el modelo demasiado complejo con una interacción de cuatro vías entre wt, hp, cyl y disp (y todas las interacciones relevantes de orden inferior y los efectos principales, a través de el operador de factor cruzado), así como los efectos principales para am, gear, drat, vs, qsec y carb. Las dos variables categóricas multinivel, cyl y gear, se convierten explícitamente en factores para evitar que se traten como numéricos.
El potencial de interacciones en el modelo final servirá para resaltar una característica especialmente importante (y conveniente) de add1, drop1 y step. Todas estas funciones respetan la jerarquía impuesta por las interacciones y los efectos principales. Es decir, para add1 (y paso), no se proporcionará un término interactivo como opción para la adición a menos que todos los efectos relevantes de orden inferior ya estén presentes en el modelo ajustado actual; De manera similar, para drop1 (y paso), no se proporcionará un término interactivo o efecto principal como una opción para la eliminación a menos que todos los efectos relevantes de orden superior ya hayan desaparecido del modelo ajustado actual. La propia función de paso devuelve un objeto de modelo ajustado y, de forma predeterminada, proporciona un informe completo de cada etapa de selección. Llamémoslo ahora; por razones de impresión, se ha recortado parte de la salida, por lo que le animamos a que lo haga en su propia máquina.
## Start: AIC=115.94
## mpg ~ 1
##
## Df Sum of Sq RSS AIC
## + wt 1 847.73 278.32 73.217
## + disp 1 808.89 317.16 77.397
## + factor(cyl) 2 824.78 301.26 77.752
## + hp 1 678.37 447.67 88.427
## + drat 1 522.48 603.57 97.988
## + vs 1 496.53 629.52 99.335
## + factor(gear) 2 483.24 642.80 102.003
## + am 1 405.15 720.90 103.672
## + carb 1 341.78 784.27 106.369
## + qsec 1 197.39 928.66 111.776
## <none> 1126.05 115.943
##
## Step: AIC=73.22
## mpg ~ wt
##
## Df Sum of Sq RSS AIC
## + factor(cyl) 2 95.26 183.06 63.810
## + hp 1 83.27 195.05 63.840
## + qsec 1 82.86 195.46 63.908
## + vs 1 54.23 224.09 68.283
## + carb 1 44.60 233.72 69.628
## + disp 1 31.64 246.68 71.356
## + factor(gear) 2 40.37 237.95 72.202
## <none> 278.32 73.217
## + drat 1 9.08 269.24 74.156
## + am 1 0.00 278.32 75.217
## - wt 1 847.73 1126.05 115.943
##
## Step: AIC=63.81
## mpg ~ wt + factor(cyl)
##
## Df Sum of Sq RSS AIC
## + hp 1 22.281 160.78 61.657
## + wt:factor(cyl) 2 27.170 155.89 62.669
## <none> 183.06 63.810
## + qsec 1 10.949 172.11 63.837
## + carb 1 9.244 173.81 64.152
## + vs 1 1.842 181.22 65.487
## + disp 1 0.110 182.95 65.791
## + am 1 0.090 182.97 65.794
## + drat 1 0.073 182.99 65.798
## + factor(gear) 2 6.682 176.38 66.620
## - factor(cyl) 2 95.263 278.32 73.217
## - wt 1 118.204 301.26 77.752
##
## Step: AIC=61.66
## mpg ~ wt + factor(cyl) + hp
##
## Df Sum of Sq RSS AIC
## + wt:hp 1 34.623 126.16 55.897
## + hp:factor(cyl) 2 28.550 132.23 59.401
## + wt:factor(cyl) 2 25.385 135.39 60.158
## + am 1 9.752 151.03 61.655
## <none> 160.78 61.657
## + drat 1 2.438 158.34 63.168
## + carb 1 1.262 159.51 63.405
## + vs 1 0.655 160.12 63.527
## + disp 1 0.651 160.13 63.527
## + qsec 1 0.229 160.55 63.611
## - hp 1 22.281 183.06 63.810
## - factor(cyl) 2 34.270 195.05 63.840
## + factor(gear) 2 7.366 153.41 64.156
## - wt 1 116.390 277.17 77.084
##
## Step: AIC=55.9
## mpg ~ wt + factor(cyl) + hp + wt:hp
##
## Df Sum of Sq RSS AIC
## - factor(cyl) 2 3.606 129.76 52.799
## <none> 126.16 55.897
## + qsec 1 5.278 120.88 56.529
## + vs 1 2.345 123.81 57.297
## + hp:factor(cyl) 2 8.645 117.51 57.625
## + drat 1 0.229 125.93 57.839
## + am 1 0.082 126.07 57.876
## + carb 1 0.038 126.12 57.887
## + disp 1 0.032 126.12 57.889
## + factor(gear) 2 7.376 118.78 57.969
## + wt:factor(cyl) 2 1.070 125.08 59.624
## - wt:hp 1 34.623 160.78 61.657
##
## Step: AIC=52.8
## mpg ~ wt + hp + wt:hp
##
## Df Sum of Sq RSS AIC
## + qsec 1 8.720 121.04 52.573
## <none> 129.76 52.799
## + vs 1 4.324 125.44 53.714
## + disp 1 0.591 129.17 54.653
## + carb 1 0.176 129.59 54.755
## + am 1 0.042 129.72 54.788
## + drat 1 0.000 129.76 54.799
## + factor(gear) 2 7.219 122.54 54.967
## + factor(cyl) 2 3.606 126.16 55.897
## - wt:hp 1 65.286 195.05 63.840
##
## Step: AIC=52.57
## mpg ~ wt + hp + qsec + wt:hp
##
## Df Sum of Sq RSS AIC
## <none> 121.04 52.573
## - qsec 1 8.720 129.76 52.799
## + factor(gear) 2 9.482 111.56 53.962
## + am 1 1.939 119.10 54.056
## + carb 1 0.080 120.96 54.551
## + drat 1 0.012 121.03 54.570
## + vs 1 0.010 121.03 54.570
## + disp 1 0.008 121.03 54.571
## + factor(cyl) 2 0.164 120.88 56.529
## - wt:hp 1 65.018 186.06 64.331Cada bloque de salida muestra el ajuste del modelo actual, su valor
AIC y una tabla que muestra los posibles movimientos (ya sea agregando
+, eliminando - o sin hacer nada
A medida que avanza el algoritmo, verá que la fila
También tenga en cuenta que la adición del término de interacción bidireccional entre wt y factor(cyl) se considera solo en el tercer paso, después de que se hayan agregado los efectos principales de esos predictores.
Sin embargo, esa interacción bidireccional particular nunca termina
siendo incluida porque el efecto principal de hp es preferible en ese
tercer paso, y las interacciones posteriores que involucran hp ofrecen
una reducción mucho mejor en el valor AIC en el cuarto paso. De hecho,
en el quinto paso, se considera que eliminar el efecto principal para
factor(cyl) reduce al máximo el AIC, por lo que las tablas para los
pasos sexto y séptimo ya no incluyen ese término wt:factor(cyl) como una
opción. . El sexto paso sugiere que agregar el efecto principal para
qsec ofrece una reducción menor en el AIC, por lo que esto se hace. La
séptima mesa señala la final del algoritmo porque no hacer nada ofrece
el valor AIC más bajo y hacer cualquier otra cosa aumentaría el AIC (se
muestra a través de
##
## Call:
## lm(formula = mpg ~ wt + hp + qsec + wt:hp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8243 -1.3980 0.0303 1.1582 4.3650
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.310410 7.677887 5.250 1.56e-05 ***
## wt -8.681516 1.292525 -6.717 3.28e-07 ***
## hp -0.106181 0.026263 -4.043 0.000395 ***
## qsec 0.503163 0.360768 1.395 0.174476
## wt:hp 0.027791 0.007298 3.808 0.000733 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.117 on 27 degrees of freedom
## Multiple R-squared: 0.8925, Adjusted R-squared: 0.8766
## F-statistic: 56.05 on 4 and 27 DF, p-value: 1.094e-12A partir de esto, parece que valdría la pena investigar más a fondo este predictor, para establecer la validez del modelo ajustado y tal vez intentar transformaciones de los datos (como modelar GPM en lugar de MPG; para ver si este efecto persiste en ejecuciones posteriores del algoritmo paso a paso AIC.
La presencia de qsec en el modelo final ilustra el hecho de que la selección del modelo no se basó únicamente en la importancia de la contribución del predictor, sino en una medida basada en criterios que apuntaba a su propia definición de parsimonia. El AIC a veces es criticado por una tendencia a errar por el lado de una mayor complejidad y valores p más altos. Para equilibrar esto, puede aumentar el efecto de penalización de los predictores adicionales aumentando la contribución multiplicativa de (p + 2); aunque el factor multiplicativo estándar de 2 se usa en la mayoría de los casos (en el paso, puede usar el argumento opcional k para cambiar esto). Dicho esto, las medidas basadas en criterios son increíblemente útiles cuando tiene modelos que no están anidados (descartando la prueba F parcial) y desea compararlos para identificar rápidamente el que, potencialmente, proporciona la representación más parsimoniosa de los datos
Otros algoritmos de selección
Cualquier algoritmo de selección de modelos siempre tendrá como objetivo definir cuantitativamente la parsimonia y sugerir un modelo que optimice esa definición a la luz de los datos disponibles. Existen alternativas al AIC, como el AIC corregido (AICc) o el Criterio de Información Bayesiano (BIC), los cuales imponen sanciones más severas a la complejidad que el AIC predeterminado . A veces es tentador monitorear simplemente \(R^{2}\), el coeficiente de determinación, para una serie de modelos. Sin embargo, esto por sí solo es inadecuado para elegir entre modelos, ya que no penaliza la complejidad y, por lo general, siempre aumentará a medida que continúe agregando predictores, ya sea que tengan un impacto estadísticamente significativo o no. La estadística \(R^{2}\) ajustada, denotada como \(R^{2}\) y reportada como R-cuadrada ajustada en resumen, es una transformación simple de la \(R^{2}\) original que incorpora una penalización por complejidad relativa al tamaño de la muestra n; calculado como
\(R^{2} = 1 - \frac{(1-R^{2})(n-1)}{n-p-1}\)
donde p es el número de términos predictores (excluyendo la intercepto). Los algoritmos basados en pruebas y criterios siempre son preferibles, pero monitorear \(R^{2}\) puede ser útil como una verificación rápida entre modelos anidados: un valor más alto indica un modelo preferido. Independientemente del enfoque que emplee, recuerde siempre que cualquier modelo final alcanzado mediante el uso de estos algoritmos aún debe estar sujeto a escrutinio.
Diagnósticos residuales
Se examinó los aspectos prácticos de los modelos de regresión lineal múltiple, como el ajuste y la interpretación, la codificación, la transformación, etc., pero aún no ha analizado los métodos que son esenciales para determinar la validez de su modelo. modelo. La parte final es presentar el diagnóstico de modelos, cuyo objetivo principal es garantizar que su modelo de regresión sea válido y represente con precisión las relaciones en sus datos.
Como repaso, en general, cuando ajuste estos modelos, recuerde tener en cuenta estas cuatro cosas: Residuales Se supone que el término de error \(\varepsilon\), que define la desviación de cualquier observación del modelo de resultado medio ajustado, tiene una distribución normal con una media de cero y una varianza constante denotada con \(\sigma^{2}\). También se supone que el residual asociado con una observación dada es independiente del error de cualquier otra observación. Si un modelo ajustado sugiere una violación de cualquiera de estas suposiciones, deberá investigar más a fondo (lo que generalmente implica reajustar una variación del modelo). la linealidad Es fundamental poder suponer que la respuesta media como función es lineal en términos de los parámetros de regresión \(\beta_{0}, \beta_{1},...,\beta_{p}\). Aunque las transformaciones de variables individuales y la presencia de interacciones pueden relajar un poco la naturaleza específica de las tendencias estimadas, se debe investigar cualquier sugerencia de diagnóstico de que una relación no es lineal (y, por lo tanto, no está capturada por el modelo ajustado disponible).
Observaciones extremas o inusuales: Inspeccione siempre los puntos de datos extremos o los puntos de datos que influyen fuertemente en el modelo ajustado; por ejemplo, los puntos que se han registrado incorrectamente deben eliminarse del análisis. Los predictores de colinealidad altamente correlacionados entre sí pueden afectar negativamente a todo un modelo, lo que significa que puede ser fácil malinterpretar los efectos de cualquier predictor incluido. Esto debe evitarse en cualquier regresión.
Investigar los primeros tres después de ajustar el modelo usando herramientas de diagnóstico. Cualquier violación de estos supuestos disminuye la confiabilidad de su modelo, a veces severamente. La colinealidad y/o las observaciones extremas se pueden descubrir mediante exploraciones estadísticas básicas (por ejemplo, visualización de matrices de diagramas de dispersión) del preajuste de datos sin procesar, pero los efectos consecuentes se evalúan después del ajuste. Hay algunas pruebas estadísticas que puede realizar para diagnosticar un modelo estadístico, pero normalmente una inspección de diagnóstico se reduce a la interpretación de los resultados de las herramientas gráficas diseñadas para apuntar a supuestos específicos. Interpretar estas tramas puede ser bastante difícil y solo se vuelve más fácil con la experiencia.
Inspección e interpretación de residuos
Una buena demostración de la importancia de interpretar los resultados dados por \(\hat{y}\) como un valor de respuesta medio. Bajo el modelo asumido, cualquier desviación de las observaciones sin procesar de la línea ajustada se considera que es el resultado de los errores (normalmente distribuidos) definidos por el término \(\varepsilon\).
Por supuesto, en la práctica, no tiene los verdaderos valores de error porque no conoce el verdadero modelo de sus datos. Para la i-ésima observación de respuesta \(y_{i}\) y su valor ajustado del modelo, \(\hat{y}_{i}\), normalmente evaluaría las gráficas de diagnóstico usando los residuos estimados \(\hat{\varepsilon }_{i} = y_{i}-\hat{y_{i}}\). Una llamada al resumen incluso alienta un análisis posterior al ajuste de los residuales al proporcionarle un resumen de cinco números de la \(\hat{\varepsilon }_{i}\) arriba de la tabla de coeficientes estimados. Esto le permite echar un vistazo a sus valores y hacer una evaluación numérica preliminar de la simetría de su distribución (como lo exige el supuesto de normalidad.
Además de una inspección de diagnóstico de los residuos brutos \(\hat{\varepsilon }_{i}\), también se pueden realizar algunas comprobaciones de diagnóstico utilizando sus valores estandarizados (o estudentizados).
Los residuales estandarizados modifican la escala de los residuos sin procesar \(\hat{\varepsilon }_{i}\) para garantizar que todos tengan la misma varianza, lo cual es importante si necesita compararlos directamente entre sí.
Formalmente, esto se logra con el cálculo \(\frac{\varepsilon_{i}}{(\hat{\sigma({1-h_{ii})^{0.5}}})}\), donde \(\hat{\sigma}\) es la estimación del error estándar residual y \(h_{ii}\) es el apalancamiento de la i-ésima observación.
Podría decirse que la herramienta gráfica más común utilizada para un análisis posterior al ajuste de los residuos es un diagrama de dispersión simple de los residuos sin procesar “observados menos ajustados” en el eje vertical contra sus valores correspondientes del modelo ajustado de la regresión.
Si los supuestos relativas a \(\varepsilon\) son válidos, entonces el \(\hat{\varepsilon }_{i}\) debería aparecer disperso aleatoriamente alrededor de cero (ya que no se supone que los errores estén relacionados de ninguna manera con el valor de la respuesta).
Cualquier patrón sistemático en la gráfica sugiere que los residuos no concuerdan con los supuestos de los residuales; esto podría deberse a relaciones no lineales en sus datos o a la presencia de observaciones dependientes (en otras palabras, sus puntos de datos están correlacionados y por lo tanto no independientes entre sí).
La gráfica también se puede usar para detectar heteroscedasticidad, una variación no constante en los residuos, comúnmente vista como un “abanico” de los residuos alrededor de 0 a medida que aumentan los valores ajustados.
Nótese una vez más que es importante que estos supuestos teóricos sean válidos porque afectan la validez de las estimaciones de los coeficientes de regresión y la confiabilidad de sus errores estándar (y, por lo tanto, estadísticos), en otras palabras, la corrección de su interpretación de su impacto en la respuesta.
Es importante saber que incluso si su diagnóstico gráfico no proporciona un gráfico de buen comportamiento, esta no es una razón para abandonar inmediatamente el análisis.
Estos gráficos pueden formar una parte integral para encontrar un modelo apropiado para sus datos. A menudo, puede reducir la no linealidad al incluir predictores o interacciones adicionales, cambiar el tratamiento de una variable categórica o realizar transformaciones no lineales de ciertos predictores continuos.
La heterocedasticidad, donde la variabilidad es mayor para valores ajustados más altos, es común en algunos campos de investigación. Un primer paso para remediar este problema a menudo implica una simple transformación de registro de la respuesta seguida de una reinspección de los diagnósticos.
Es hora de un ejemplo. Se usó la selección paso a paso de AIC para elegir un modelo para MPG para los datos de mtcars, creando el objeto car.step. Ahora diagnostiquemos ese mismo ajuste para ver si hay algún problema con los supuestos del modelo.
Cuando aplica la función de gráfico directamente a un objeto de
película, puede producir convenientemente seis tipos de gráficos de
diagnóstico del ajuste. Por defecto, cuatro de estos gráficos se
producen en sucesión. Siga la señal al usuario en la consola, presione
En los ejemplos que siguen, sin embargo, seleccionará cada gráfico individualmente usando el argumento opcional which (especificado por los números enteros del 1 al 6; vea ?plot.lm para la documentación). La gráfica de residuos versus ajustada se da con which=1; la siguiente línea produce el gráfico de la izquierda
Como puede ver, R agrega una línea suavizada para ayudar al usuario a interpretar cualquier tendencia. De forma predeterminada, se anotan los tres puntos más extremos desde cero (de acuerdo con el atributo de nombres de fila del dataframe utilizado en la llamada que se ajustó al modelo).
La fórmula del modelo en sí se especifica debajo de la etiqueta del eje horizontal.
A partir de esto, se ve que la gráfica de residuos versus ajustada para el paso del automóvil ofrece pocos motivos de preocupación, si es que los hay. No hay mucha tendencia perceptible, y puede estar más tranquilo con el hecho de que los errores (ei) parecen homocedásticos en su distribución.
La gráfica de ubicación de escala es similar a la gráfica de residuos versus ajustada, aunque en lugar de la ei sin procesar en el eje vertical, la gráfica de ubicación de escala proporciona ei/(σˆ 1 hii 0.5) 0.5, es decir, la raíz cuadrada de la ecuación absoluta. valor (indicado por ; esto hace que todos los valores negativos sean positivos) de los residuos estandarizados. Estos se trazan contra los respectivos valores ajustados en el eje horizontal.
Al restringir la atención a la magnitud de cada residual de esta manera, la gráfica de ubicación de escala se usa para revelar tendencias en el tamaño de la desviación de cada punto de datos de su valor ajustado, a medida que aumentan los valores ajustados.
Esto significa que una gráfica de este tipo puede, por ejemplo, ser más útil que la gráfica de residuos sin procesar versus gráfica ajustada para detectar cosas como la heteroscedasticidad. Al igual que con la gráfica original de residuos versus gráfica ajustada, usted está buscando una gráfica sin un patrón perceptible como indicación de que no se han violado los supuestos de error.
El diagrama de la derecha de la figura 22-2 muestra el diagrama de ubicación de escala para car.step, seleccionado con which=3. Este gráfico también demuestra la capacidad de eliminar la línea de tendencia suavizada predeterminada con el argumento add.smooth y controlar cuántos puntos extremos se etiquetan con el argumento id.n.
Al igual que con la gráfica original de residuos versus gráfica ajustada, no parece haber mucho de qué preocuparse en la gráfica de escala-ubicación para este modelo mtcars. Regrese a los datos de Galileo sobre la bola rodante que se presentaron vez en su uso en el siguiente ejemplo, la variable de respuesta “distancia recorrida” se da como columna d, y la variable explicativa “altura” es la columna h, en el marco de datos gal.
Volveré a crear parte del ejercicio para brindarle un par de ejemplos sencillos de motivos de preocupación en las gráficas de diagnóstico residual. Ejecute el siguiente código para definir el marco de datos de las siete observaciones y ajuste dos modelos de regresión: el primero lineal en altura y el segundo cuadrático.
La primwe gráfica muestra los datos y proporciona la línea recta del modelo lineal simple. Aunque esto captura claramente la tendencia creciente, este gráfico sugiere que también hay cierta curvatura.
La segunda gráfica de residuos de diagnóstico versus gráfica ajustada muestra que el modelo de solo tendencia lineal es inadecuado: el patrón sistemático arroja una bandera roja con respecto a los supuestos que rodean los errores del modelo lineal. La tercer imagen in muestra la gráfica de residuos versus ajustada basada en la versión cuadrática del modelo en gal.mod2.
Incluir un término cuadrático en “altura” elimina esta curva prominente en los residuos. Sin embargo, estos últimos valores de \(\hat{\varepsilon }_{i}\) todavía parecen exhibir un comportamiento sistemático en forma de onda, lo que quizás sugiera que se debe probar un modelo cúbico, lo cual es difícil con un tamaño de muestra tan pequeño.