MODELOS LINEALES I TAREA EXAMEN 1: REGRESIÓN LINEAL SIMPLE Y MÚLTIPLE
MODELOS LINEALES I TAREA EXAMEN 1: REGRESIÓN LINEAL SIMPLE Y MÚLTIPLE
- Pregunta 1: Gauss-Markov
- Pregunta 2: Significancia debida a la Regresión
- Pregunta 3: Mínimos Cuadrados y Máxima Verosimilitud.
- Pregunta 4: Intervalos de Predicción y Respuesta Media.
- Pregunta 5: Algunas Propiedades de los Estimadores
- Pregunta 6: Algunas Propiedades de los Estimadoes II
- Pregunta 7: Algunas Propiedades de los Estimadores III
- Pregunta 8: Modelo de Regresión Lineal sin Intercepto
- a). Obtenga el estimador de mínimos cuadrados de \(\beta\).
- b). Demuestre que es insesgado y obtenga su varianza.
- c). Encuentre una expresión para \(\hat{\sigma^2}\). ¿Cuántos grados de libertad tiene y por qué?
- d). Si se tienen dos modelos del mismo conjunto de datos, uno con y otro sin intercepto, investigue cuál es la estadística que permite comparar el mejor de los modelos.
- Pregunta 9: Algunas Propiedades del Modelo de Regresión Múltiple
- Pregunta 10: Algunas Propiedades del Modelo de Regresión Múltiple II
- Pregunta 11: Algunas Propiedades del Modelo de Regresión Múltiple III
- Pregunta 12: Algunas Propiedades del Modelo de Regresión Múltiple IV
- Pregunta 13: Ejemplo I
- Pregunta 14: Ejemplo II
Pregunta 1: Gauss-Markov
Enuncia el Teorema de Gauss-Markov.
\(\underline{\text{Solución:}}\)
Dado un modelo de Regresión para una serie de datos \(\left(X_i,Y_i\right)_{i=1}^n\):
\[Y=\beta_0+\beta_1X+\epsilon,~\text{i.e.}~\forall_{i=1,\cdots, n} ~Y_i=\beta_0+\beta_1X_i+\epsilon_i,\] bajo los supuestos de que \(\sum_{i=1}^n\left[\epsilon_i\right]=0,\) que \(\forall_{i=1,\cdots,n}~\text{Var}(\epsilon_i)=0\) y que los errores \(\epsilon\) no son correlacionados, entoces los estimadores de mínimos cuadrados para \(\beta_0\) y \(\beta_1\):
\[ \begin{aligned} \hat{\beta_0}&=\bar{Y}-\hat{\beta_1}\bar{X}:\bar{Y}=\frac{1}{n}\sum_{i=1}^n(Y_i),\\ \hat{\beta_1}&=\frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{\sum_{i=1}^n(X_i-\bar{X})^2}=\frac{S_{XY}}{S_{XX}}. \end{aligned} \] son insesgagos y de varianza mínima al compararlos con otros estimadores insesgados que son combinaciones lineales de las \(\left(Y_i\right)_{i=1}^n\). Es decir, \(\hat{\beta_0}\) y \(\hat{\beta_1}\) son considerados los Mejores Estimadores Lineales Insesgados en el sentido de varianza mínima.
\(_\blacksquare\)
Pregunta 2: Significancia debida a la Regresión
¿Qué significa la hipótesis de significancia debida a la regresión? ¿En qué supuestos está basada?
\(\underline{\text{Solución:}}\)
Dado un modelo de Regresión para una serie de datos \(\left(X_i,Y_i\right)_{i=1}^n\):
\[Y=\beta_0+\beta_1X+\epsilon,~\text{i.e.}~\forall_{i=1,\cdots, n} ~Y_i=\beta_0+\beta_1X_i+\epsilon_i,\]
La Hipótesis de Significancia debida a la Regresión se define como la prueba de hipótesis:
\[\text{H}_0:\beta_1=0\text{ vs. }\text{H}_1:\beta_1\ne 0.\] Básicamente se la prueba con la se pone en duda la exitencia de un modelo lineal para el conjunto de datos \(\left(X_i,Y_i\right)_{i=1}^n\), puesto que si no se rechaza la hipótesis nula \(\text{H}_0\), entonces se tendría evidencia para decir que el modelo de regresión lineal es una constante: \(\hat{Y}=\hat{\beta_1}X+\hat{\beta_0}=\hat{\beta_0}\), lo cual implica una regresión lineal trival ó, para fines prácticos, inútil.
Esta prueba de hipótesis se basa en el supuesto de que existe una relación entre las variables aleatorias \(X\) y \(Y\). Además, para ocupar un estadístico con el cual realizar la prueba son necesarios también los supuestos de independecia de los errores \(\epsilon_i=Y_i-\hat{Y_i}\), los errores tienen todo varinaza \(\sigma^2\) constante y desconocida, tienen media cero y cada uno tiene distribución \(\epsilon_i\sim\mathcal{N}(0,\sigma^2)\).Esto último se debe a que con esos supuestos uno puede obtener las distribuciones de los estimadores \(\hat{\beta_0},\hat{\beta_0},\hat{\sigma^2}\) al saber que las variables \(Y_i\sim\mathcal{N}(\beta_0+\beta_1X_i,\sigma^2)\) y que los estimadores son todos combinaciones lineales de las \(Y_i\).
\(_\blacksquare\)
Pregunta 3: Mínimos Cuadrados y Máxima Verosimilitud.
¿Qué diferencia existe entre la estimación de por mínimos cuadrados y los de máxima verosimilitud?
\(\underline{\text{Solución:}}\)
Hay dos diferencias principales:
1). Supuestos
El método de Mínimos Cuadrados [MC de aquí en adelante], no requiere suponer la distribución de ningún componente de la regresión.
Por otro lado, el método de Máxima Verosimilitud [MV de aquí en adelante], require suponer la distribución de los errores de la regresión \(\epsilon_i=Y_i-\hat{Y_i}\) y (por ende) la distribución de la variable de respuesta \(Y_i\). Usualmente suponemos que \(\epsilon_i\sim\mathcal{N}(0,\sigma^2)\).
2). Obtención
En MC por medio de cálculo es que buscamos los estimadores de \(\beta_0\) y \(\beta_1\) que minimicen la diferencia de cuadrados:
\[\inf_{\beta_0,\beta_1\in\mathbb{R}}\Biggr\{\sum_{i=1}^n\left(Y_i-(\beta_0+\beta_1X_i)\right)^2\Biggr\}\]
En MV ocupamos la función de verosimilitud para obtener los estimadores de MV patra \(\beta_0,\beta_1,\sigma^2\):
\[L(\beta_0,\beta_1,\sigma^2|X,Y)=\prod_{i=1}^{n}e^{\left[-\frac{1}{2\sigma^2}\left(Y_i-\beta_0-\beta_1X_i\right)^2\right]}\] Esas son las únicas diferencias, ya que al minimizar la suma de cuadrados de las diferencias y al maximizar la función de verosimilitud, obtenedremos que los estimadores por mínimos cuadrados y los de máxima verosimilitud de \(\beta_0,\beta_1\) son exactamente iguales:
\[ \begin{aligned} \hat{\beta_0}&=\bar{Y}-\hat{\beta_1}\bar{X}:\bar{Y}=\frac{1}{n}\sum_{i=1}^n(Y_i),\\ \hat{\beta_1}&=\frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{\sum_{i=1}^n(X_i-\bar{X})^2}=\frac{S_{XY}}{S_{XX}}. \end{aligned} \] Además tendremos que el estimador de MV para \(\sigma^2\) es:
\[\hat{\sigma^2}=\frac{1}{n}\sum_{i=1}^{n}\left(Y_i-\hat{\beta_0}-\hat{\beta_1}X_i\right)^2.\] \(_\blacksquare\)
Pregunta 4: Intervalos de Predicción y Respuesta Media.
¿Cuál es la diferencia entre el intervalo de predicción y el de la respuesta media?
\(\underline{\text{Solución:}}\)
Sea \(X_0\) el valor del regresor, independiente para el cual se pretende calcular el valor de la respuesta media de la variable de respuesta:
\[\mathbb{E}\left[Y|X_0\right]=\hat{Y_0}=\hat{\beta_0}+\hat{\beta_1}X_0\implies\hat{Y_0}\sim\mathcal{N}\left(\beta_0+\beta_1X_0,\sigma^2\left(\frac{1}{n}+\frac{\left(X_0-\bar{X}\right)^2}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right)\right)\] Por lo que, (después de estadarizar y obtener una catidad pivotal), tendremos que el intevalo de confianza de \(1-\alpha\) para la respuesta media es:
\[\mathbb{E}\left[Y|X_0\right]\in\left(\hat{Y_0}\pm t_{(n-2)}^{1-\frac{\alpha}{2}}\sqrt{\text{CME}\left(\frac{1}{n}+\frac{\left(X_0-\bar{X}\right)^2}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right)}\right)\]
Donde \(\text{CME}=\frac{1}{n}\sum_{i=1}^{n}\left(Y_i-\hat{\beta_0}-\hat{\beta_1}X_i\right)^2\) es el Cuadrado Medio del Error.
Mientras que para la predicción \(Y_0^*=\hat{\beta_0}+\hat{\beta_1}X_0+\epsilon\sim\mathcal{N}\left(Y_0,\sigma^2\left(1+\frac{1}{n}+\frac{1}{S_{XX}}(X_0-\bar{X})^2\right)\right)\)
Por lo que, (después de estadarizar y obtener una catidad pivotal), tendremos que el intevalo de confianza de \(1-\alpha\) para la predicción es:
\[Y_0\in\left(\hat{Y_0}\pm t_{(n-2)}^{1-\frac{\alpha}{2}}\sqrt{\text{CME}\left(1+\frac{1}{n}+\frac{\left(X_0-\bar{X}\right)^2}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right)}\right)\] La diferencia radica en que si comparamos el intervalo de \(\mathbb{E}\left[Y|X_0\right]\) con el de predicción, el primero es menor porque depende del error debido al ajuste del modelo y el segundo al del error asociado con las futuras observaciones.
\(_\blacksquare\)
Pregunta 5: Algunas Propiedades de los Estimadores
En un modelo de regresión simple, el estimador puntual de \(\hat{\beta_1}=\frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{\sum_{i=1}^n(X_i-\bar{X})^2}=\frac{S_{XY}}{S_{XX}}.\) Además, se sabe que \(\mathbb{E}\left[\hat{\beta_1}\right]=\beta_1,\text{Var}(\hat{\beta_1})=\frac{\sigma^2}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\). Para demostrar esto, necesitamos reconocer que \(\hat{\beta_1}\) es una combinación lineal de las observaciones \(Y_i\). Entonces, dado que \(\hat{\beta_1}=\sum_{i=1}^{n}k_iY_i,\) donde \(k_i=\frac{X_i-\bar{X}}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\), demuestre lo siguiente:
i). \(\sum_{i=1}^n k_i=0\)
\(\underline{\text{Demostración:}}\)
\[ \begin{aligned} \sum_{i=1}^nk_i&=\sum_{i=1}^n\left[\frac{X_i-\bar{X}}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right]=\frac{1}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\sum_{i=1}^n\left[X_i-\bar{X}\right]\\ &=\frac{1}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\left[\sum_{i=1}^n(X_i)-n\bar{X}\right]=\frac{1}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\left[n\bar{X}-n\bar{X}\right]=0. \end{aligned} \]
\(_\blacksquare\)
ii). \(\sum_{i=1}^n k_iX_i=1\)
\(\underline{\text{Demostración:}}\)
\[ \begin{aligned} \sum_{i=1}^nk_iX_i&=\sum_{i=1}^n\left[\frac{(X_i-\bar{X})X_i}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right]=\frac{1}{{\sum_{i=1}^{n}(X_i-\bar{X})^2}}\sum_{i=1}^n\left[X_i^2-\bar{X}X_i\right]\\ &=\frac{\sum_{i=1}^n\left[X_i^2-\bar{X}X_i\right]}{\sum_{i=1}^{n}\left[X_i^2-2X_i\bar{X}+\bar{X}^2\right]}=\frac{\sum_{i=1}^n\left[X_i^2-\bar{X}X_i\right]}{\sum_{i=1}^{n}\left[X_i^2\right]+-2n\bar{X}^2+n\bar{X}^2}\\ &=\frac{\sum_{i=1}^n\left[X_i^2-\bar{X}X_i\right]}{\sum_{i=1}^{n}\left[X_i^2\right]+-n\bar{X}^2}=\frac{\sum_{i=1}^n\left[X_i^2\right]-n\bar{X}^2}{\sum_{i=1}^{n}\left[X_i^2\right]-n\bar{X}^2}=1. \end{aligned} \]
\(_\blacksquare\)
iii). \(\sum_{i=1}^n k_i^2=\frac{1}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\)
\(\underline{\text{Demostración:}}\)
\[ \begin{aligned} \sum_{i=1}^nk_i^2&=\sum_{i=1}^n\left[\frac{X_i-\bar{X}}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right]^2=\sum_{i=1}^n\frac{\left[X_i-\bar{X}\right]^2}{\left[\sum_{i=1}^{n}(X_i-\bar{X})^2\right]^2}\\ &=\frac{1}{\left[\sum_{i=1}^{n}(X_i-\bar{X})^2\right]^2}\left[\sum_{i=1}^n \left[X_i-\bar{X}\right]^2\right]=\frac{1}{\sum_{i=1}^{n}(X_i-\bar{X})^2}. \end{aligned} \]
\(_\blacksquare\)
Pregunta 6: Algunas Propiedades de los Estimadoes II
Del modelo de regresión lineal simple \(Y_i=\beta_0+\beta_1X_i+\epsilon_i\) con \(\mathbb{E}[\epsilon_i]=0\) y \(\text{Var}[\epsilon_i]=\sigma^2\) y los \(\{\epsilon_i\}_{i=1}^n\) son non correlacionados.
i). Demuestre que \(\text{Cov}\left[\hat{\beta_0},\hat{\beta_1}\right]=\frac{-\bar{X}\sigma^2}{S_{XX}}\)
\(\underline{\text{Demostración:}}\)
Debido a que \(\hat{\beta_0},\hat{\beta_1}\) pueden expresarse como combinaciones lineales de las \(Y_i\), entonces existen \(\{a_i,b_i\}_{i=1}^{n}\subset\mathbb{R}\) tales que:
\[ \begin{aligned} \text{Cov}\left[\hat{\beta_0},\hat{\beta_1}\right]&=\mathbb{E}\left[\left(\hat{\beta_0}-\mathbb{E}\left[\hat{\beta_0}\right]\right)\left(\hat{\beta_1}-\mathbb{E}\left[\hat{\beta_1}\right] \right)\right]\\ &=\mathbb{E}\left[\left(\sum_{i=1}^na_iY_i-\mathbb{E}\left[\sum_{i=1}^na_iY_i\right]\right)\left(\sum_{i=1}^nb_iY_i-\mathbb{E}\left[\sum_{i=1}^nb_iY_i\right] \right)\right]\\ &=\sum_{i=1}^na_ib_i\mathbb{E}\left[(Y_i-\mathbb{E}[Y_i])(Y_i-\mathbb{E}[Y_i])\right]=\sum_{i=1}^na_ib_i\text{Var}[Y_i]=\sum_{i=1}^na_ib_i\sigma^2 \end{aligned} \] Sabemos que para toda \(i=1,\cdots,n\) tenemos que \(a_i=\frac{X_i-\bar{X}}{S_{XX}}\) y que \(b_i=\frac{1}{n}-\left(\frac{X_i-\bar{X}}{S_{XX}}\right)\bar{X}\), (), por lo tanto :
\[ \begin{aligned} \text{Cov}\left[\hat{\beta_0},\hat{\beta_1}\right]&=\sum_{i=1}^na_ib_i\sigma^2=\sum_{i=1}^n\left[\frac{X_i-\bar{X}}{S_{XX}}\right]\left[\frac{1}{n}-\left(\frac{X_i-\bar{X}}{S_{XX}}\right)\bar{X}\right]\sigma^2\\ &=\sum_{i=1}^n\left[\frac{X_i-\bar{X}}{nS_{XX}}-\left(\frac{X_i-\bar{X}}{S_{XX}}\right)^2\bar{X}\right]\sigma^2\\ &=\left[\frac{n\bar{X}-n\bar{X}}{nS_{XX}}-\frac{\bar{X}S_{XX}}{S_{XX}^2}\right]\sigma^2=\frac{-\bar{X}\sigma^2}{S_{XX}}. \end{aligned} \]
\(_\blacksquare\)
ii). Demuestre que \(\text{Cov}\left[\bar{Y},\hat{\beta_1}\right]=0\)
\(\underline{\text{Demostración:}}\)
\[ \begin{aligned} \text{Cov}\left[\bar{Y},\hat{\beta_1}\right]&=\text{Cov}\left[\hat{\beta_0}+\hat{\beta_1}\bar{X},\hat{\beta_1}\right]\\ &=\text{Cov}\left[\hat{\beta_0},\hat{\beta_1}\right]+\bar{X}\text{Cov}\left[\hat{\beta_1},\hat{\beta_1}\right]\\ &=\frac{-\bar{X}\sigma^2}{S_{XX}}+\bar{X}\text{Var}\left[\hat{\beta_1}\right]=\frac{-\bar{X}\sigma^2}{S_{XX}}+\frac{\bar{X}\sigma^2}{S_{XX}}=0. \end{aligned} \] \(_\blacksquare\)
Pregunta 7: Algunas Propiedades de los Estimadores III
Demuestre que el estadístico \(\left(t_0\right)^2\) para la prueba parcial de \(\beta_1\) es equivalente al estadístico \(F_0\) para la prueba de significancia de la regresión
\(\underline{\text{Demostración:}}\)
Sabemos que para la prueba de hipótesis parcial de \(\beta_1\), es decir:
\[\text{H}_0:\beta_1=\beta_{10}\text{ vs. }\text{H}_1:\beta_1\ne\beta_{10},\] la cantidad pivotal ocupada requiere del estimador \(\hat{\beta_1}\) y esa es:
\[t_0=\frac{\hat{\beta_1}-\beta_{10}}{\sqrt{\frac{\frac{1}{n-2}\sum_{i=1}^n\left(Y_i-\hat{Y_i}\right)^2}{\sum_{i=1}^n\left(X_i-\bar{X}\right)^2}}}=\frac{\hat{\beta_1}-\beta_{10}}{\sqrt{\frac{\text{CME}}{\sum_{i=1}^n\left(X_i-\bar{X}\right)^2}}}\sim t_{(n-2)}.\] De forma similar sabemos que para la prueba de hipótesis de significancia del modelo:
\[\text{H}_0:\beta_1=0\text{ vs. }\text{H}_1:\beta_1\ne 0,\] la cantidad pivotal ocupada es:
\[ F_0=\frac{\sum_{i=1}^{n}\left(\hat{Y_i}-\bar{Y_i}\right)^2}{\frac{1}{n-2}\sum_{i=1}^{n}\left(\hat{Y_i}-Y_i\right)^2}=\frac{\text{SCR}}{\frac{1}{n-2}\text{SCE}}\sim F_{(1,n-2)}. \] Por lo que si elevamos al cuadrado la primer cantidad pivotal y suponemos \(\beta_{10}=0\) tendremos que:
\[ \begin{aligned} (t_0)^2&=\left[\frac{\hat{\beta_1}}{\sqrt{\frac{\frac{1}{n-2}\sum_{i=1}^n\left(Y_i-\hat{Y_i}\right)^2}{\sum_{i=1}^n\left(X_i-\bar{X}\right)^2}}} \right]^2=\frac{\left[\frac{S_{XY}}{S_{XX}}\right]^2}{\frac{\frac{1}{n-2}\sum_{i=1}^n\left(Y_i-\hat{Y_i}\right)^2}{\sum_{i=1}^n\left(X_i-\bar{X}\right)^2}}=\frac{\left[\frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{\sum_{i=1}^n(X_i-\bar{X})^2}\right]^2}{\frac{\frac{1}{n-2}\sum_{i=1}^n\left(Y_i-\hat{Y_i}\right)^2}{\sum_{i=1}^n\left(X_i-\bar{X}\right)^2}}\\ &=\frac{\sum_{i=1}^n\left(X_i-\bar{X}\right)^2\left[\frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{\sum_{i=1}^n(X_i-\bar{X})^2}\right]^2}{\frac{1}{n-2}\text{SCE}}=\frac{\frac{\sum_{i=1}^n\left(X_i-\bar{X}\right)^2}{\sum_{i=1}^n(X_i-\bar{X})^2}\left(\frac{\left(\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})\right)^2}{\sum_{i=1}^n(X_i-\bar{X})^2}\right)}{\frac{1}{n-2}\text{SCE}}\\ &=\frac{\frac{(S_{XY})^2}{S_{XX}}}{\frac{1}{n-2}\text{SCE}}=\frac{\hat{\beta_1}S_{XY}}{\frac{1}{n-2}\text{SCE}}. \end{aligned} \] Nótese que:
\[ \text{SCR}=\sum_{i=1}^{n}\left(\hat{Y_i}-\bar{Y_i}\right)^2=\hat{\beta_1}S_{XY}. \] Por lo tanto:
\[ (t_0)^2=\frac{\hat{\beta_1}S_{XY}}{\frac{1}{n-2}\text{SCE}}=\frac{\text{SCR}}{\frac{1}{n-2}\text{SCE}}=F_0. \]
\(_\blacksquare\)
Pregunta 8: Modelo de Regresión Lineal sin Intercepto
Ocasionalmente, un modelo de regresión en el cual el intercepto tiene valor (a priori) cero puede ser ajustado por el siguiente modelo:
\[Y_i=\beta X_i+\epsilon_i.\] Resuelva lo siguiente
a). Obtenga el estimador de mínimos cuadrados de \(\beta\).
\(\underline{\text{Solución:}}\)
Se busca minimizar \(\sum_{i=1}^n(\epsilon_i)^2\), (en este caso \(\hat{Y_i}=\beta X_i\)) ,por ello:
\[S(\beta)=\sum_{i=1}^n(\epsilon_i)^2=\sum_{i=1}^n\left(Y_i-\hat{\beta}X_i\right)^2.\] Derivando con respecto a \(\hat{\beta}\) e igualando a 0:
\[\frac{\partial{S(\beta)}}{\partial\hat{\beta}}=-2\sum_{i=1}^n\left(Y_i-\hat{\beta}X_i\right)X_i=0\implies \sum_{i=1}^nX_iY_i-\hat{\beta}X_i^2\implies \hat{\beta}=\frac{\sum_{i=1}^nX_iY_i}{\sum_{i=1}^nX_i^2}.\] Volviendo a derivar para saber si es un máximo ó mínimo:
\[\frac{\partial^2{S(\beta)}}{\partial\hat{\beta}^2}=2\sum_{i=1}^n(Y_i-\hat{\beta}X_i)^{-1}X_i^2>0.\] Por lo tanto es un mínimo, por consiguiente, el estimador de \(\beta\) que minimiza la suma de cuadrados de los residuos es:
\[\hat{\beta}=\frac{\sum_{i=1}^nX_iY_i}{\sum_{i=1}^nX_i^2}.\]
\(_\blacksquare\)
b). Demuestre que es insesgado y obtenga su varianza.
\(\underline{\text{Solución:}}\)
\[ \begin{aligned} \mathbb{E}\left[\hat{\beta}\right]&= \mathbb{E}\left[\frac{\sum_{i=1}^nX_iY_i}{\sum_{i=1}^nX_i^2}\right]=\frac{\sum_{i=1}^nX_i}{\sum_{i=1}^nX_i^2}\mathbb{E}[Y_i]\\ &=\frac{\sum_{i=1}^nX_i}{\sum_{i=1}^nX_i^2}\beta X_i=\frac{\sum_{i=1}^nX_i^2}{\sum_{i=1}^nX_i^2}\beta\\ \implies \mathbb{E}\left[\hat{\beta}\right]=\beta \end{aligned} \] Por lo que \(\hat{\beta}\) sí es insesgado.
Ahora, nótese que:
\[ \begin{aligned} \hat{\beta}&=\frac{\sum_{i=1}^nX_iY_i}{\sum_{i=1}^nX_i^2}\\ \implies& \hat{\beta}= \left(\frac{X_1}{\sum_{i=1}^nX_i^2}\right)Y_1+\cdots+\left(\frac{X_n}{\sum_{i=1}^nX_i^2}\right)Y_n. \end{aligned} \] Por lo que:
\[ \begin{aligned} \text{Var}\left(\hat{\beta}\right)&=\text{Var}\left(\left(\frac{X_1}{\sum_{i=1}^nX_i^2}\right)Y_1+\cdots+\left(\frac{X_n}{\sum_{i=1}^nX_i^2}\right)Y_n\right)\\ &=\left(\frac{X_1}{\sum_{i=1}^nX_i^2}\right)^2\text{Var}\left(Y_1\right)+\cdots+\left(\frac{X_n}{\sum_{i=1}^nX_i^2}\right)^2\text{Var}\left(Y_n\right)\\ &=\left(\frac{X_1}{\sum_{i=1}^nX_i^2}\right)^2\sigma^2+\cdots+\left(\frac{X_n}{\sum_{i=1}^nX_i^2}\right)^2\sigma^2.\\ \implies \text{Var}\left(\hat{\beta}\right)&=\left(\frac{\sum_{i=1}^nX_i^2}{\left(\sum_{i=1}^nX_i^2\right)^2}\right)\sigma^2=\left(\frac{1}{\sum_{i=1}^nX_i^2}\right)\sigma^2. \end{aligned} \]
\(_\blacksquare\)
c). Encuentre una expresión para \(\hat{\sigma^2}\). ¿Cuántos grados de libertad tiene y por qué?
\(\underline{\text{Solución:}}\)
De forma totalmente análoga al modelo de regresión lineal con intercepto, (y ocupando el supuesto de normalidad de los errores), tendremos que un estimador para la varianza del modelo es:
\[\hat{\sigma^2}=\frac{1}{n-1}\sum_{i=1}^n\left(Y_i-\hat{\beta}X_i \right)^2=\frac{1}{n-1}\sum_{i=1}^n\left(Y_i-\hat{Y_i} \right)^2.\] Este estimador posee \(n-1\) grados de libertad y eso se debe a que sólo es un parámetro el que debe estimarse en esta ocasión \(\beta\), (a diferencia del modelo de regresión lineal con intercepto, en el cual eran dos lo parámetros a estimar: \(\beta_1,\beta_2\)) y son \(n\) datos los que componen nuestra muestra aleatoria \(\{(X_i,Y_i)\}_{i=1}^n\)
\(_\blacksquare\)
d). Si se tienen dos modelos del mismo conjunto de datos, uno con y otro sin intercepto, investigue cuál es la estadística que permite comparar el mejor de los modelos.
\(\underline{\text{Solución:}}\)
Sea \(D\) el conjunto de Datos, \(M_0\) el modelo de regresión lineal sin intercepto y \(M_1\) el modelo de regresión lineal con intercepto que ajustan a \(D\), entonces un posible estadístico que nos permita comparar la eficacia de los modelos y determinar cúal es el mejor es el Cuadrado Medio del Error \(CME_{M_1}=\frac{SCE}{n-2}=\frac{1}{n-2}\sum_{i=1}^{n}\left(Y_i-\hat{Y_i}\right)^2\sim\chi^2_{(n-2)}\) y \(CME_{M_0}=\frac{SCE}{n-2}=\frac{1}{n-2}\sum_{i=1}^{n}\left(Y_i-\hat{Y_i}\right)^2\sim\chi^2_{(n-1)}\) (eso debido a la cantidad de parámetro que se estiman en cada caso), puesto que el ‘mejor modelo’ debería ser aquel que (entre otras cosas) reduzca la suma de errores cuadrados y, por lo tanto, minimizar la esperanza de la suma de cuadrados de los errores, (eso es \(CME\)). La ventaja de ocupar este estadístico es el hecho de que su distribución \(\chi^2\) es estrictamente positiva, susu cuartiles están totalmente determinados y que la determinación del mejor modelo es simple:
\(M_0\) será mejor que \(M_1\) si\(CME_{M_0}<CME_{M_1}\) y viceversa.
\(_\blacksquare\)
Pregunta 9: Algunas Propiedades del Modelo de Regresión Múltiple
Se tiene que el modelo de regresión lineal múltiple \(Y=X\mathbb{B}+\epsilon\). Demostrar que el esimador de mínimos cuadrados se puede escribir de la forma:
\[\hat{\mathbb{B}}=\mathbb{B}+R\epsilon:R=(X^tX)^{-1}X^t\]
\(\underline{\text{Demostración:}}\)
Tenemos que
\[ Y_i=\beta_0+\beta_1X_{1i}+\beta_2X_{2i}+\cdots+\beta_rX_{ri}+\epsilon_i=\beta_0+\sum_{j=1}^{r}\beta_jX_{ji}+\epsilon_{i}~,~i=1,\cdots,n \]
La función de mínimo cuadrados es:
\[S(\mathbb{B})=S(\beta_0,\beta_1,\beta_2,\cdots,\beta_r)=\sum_{i=1}^n\left(Y_i-\hat{Y_i}\right)^2=\sum_{i=1}^n\left(Y_i-\left(\beta_0+\sum_{j=1}^{r}\beta_jX_{ji}\right)\right)^2\] Las ecuaciones normales para cada parámetro son:
\[ \begin{aligned} \frac{\partial S(\beta_0,\beta_1,\beta_2,\cdots,\beta_r)}{\partial \beta_0}&=-2\sum_{i=1}^{n}\left(Y_i-\beta_0-\sum_{j=1}^{r}\beta_jX_{ji}\right)=0\\ \frac{\partial S(\beta_0,\beta_1,\beta_2,\cdots,\beta_r)}{\partial \beta_i}&=-2\sum_{i=1}^{n}X_i\left(Y_i-\sum_{j=1}^{r}\beta_jX_{ji}\right)=0 \end{aligned} \]
Ocupando la notación de matrices:
\[ \begin{aligned} Y_{1\times n}^{t}&=[Y_1,\cdots,Y_n],\\ \beta_{1\times r}^t&=[\beta_0,\beta_1,\cdots,\beta_r],\\ \epsilon_{1\times n}&=[\epsilon_1,\cdots,\epsilon_n],\\ X_{n\times r}&=\begin{bmatrix} 1 & X_{11} & \cdots & X_{1r}\\ \vdots & \vdots & \cdots & \vdots\\ 1 & X_{n1} & \cdots & X_{nr} \end{bmatrix},\\ S(\beta_0,\beta_1,\beta_2,\cdots,\beta_r)&=S(\mathbb{B})=\epsilon^t\epsilon=(Y-X\mathbb{B})^t(Y-X\mathbb{B}),\\ \frac{\partial S(\beta_0,\beta_1,\beta_2,\cdots,\beta_r)}{\partial \mathbb{B}}\Bigg |_{\mathbb{B}}&=-2X^tY+2X^tX\mathbb{B}=0\\ \therefore \hat{\mathbb{B}}&=(X^tX)^{-1}X^tY \end{aligned} \] Ahora, sustituyendo la anterior expresión en \(Y\), (sea \(D=\left[(X^tX)^{-1}X^t\right]^{-1}\)) ,tendremos lo siguiente:
\[ \begin{aligned} Y&=X\mathbb{B}+\epsilon=D\hat{\mathbb{B}}\\ &\implies \hat{\mathbb{B}}=D^{-1}\left[X\mathbb{B}+\epsilon\right]=(X^tX)^{-1}X^tX\mathbb{B}+(X^tX)^{-1}X^t\epsilon=\mathbb{B}+R\epsilon.\\ \therefore& ~ \hat{\mathbb{B}}=\mathbb{B}+R\epsilon:R=(X^tX)^{-1}X^t. \end{aligned} \]
\(_\blacksquare\)
Pregunta 10: Algunas Propiedades del Modelo de Regresión Múltiple II
Demostrar que los residuos de un modelo de regresión lineal se pueden expresar de la forma \(e=(1-H)\epsilon\)
\(\underline{\text{Demostración:}}\)
Dado que \(e=Y-\hat{Y}\) y \(\hat{Y}=X\hat{\mathbb{B}}=(X^tX)^{-1}X^tY=RY\) entonces tenemos que:
\[ \begin{aligned} e&=Y-X\hat{\mathbb{B}}=Y-X\left[\mathbb{B}+R\epsilon\right]=Y-X\mathbb{B}-XR\epsilon=\epsilon-XR\epsilon=\left[I-XR\right]\epsilon\\ &=[I-X(X^tX)^{-1}X^t]\epsilon=[I-H]\epsilon.\\ \therefore~e&=[I-H]\epsilon. \end{aligned} \]
\(_\blacksquare\)
Pregunta 11: Algunas Propiedades del Modelo de Regresión Múltiple III
Considere un modelo de regresión con \(k\) variables y \(n\) observaciones. Demuestre que:
\[1-\text{R}^{2}_{\text{adj}}=\left(1-\text{R}^2\right)\cdot\frac{n-1}{n-k-1}.\]
\(\underline{\text{Demostración:}}\)
Sabemos que, (sea \(p-1=k+1-1\) el número de grados de libertad del modelo (eso debido a que tenemos \(k\) regresores ó variables)):
\[ \begin{aligned} & \text{R}^{2}=1-\frac{SCE}{SCT}\implies (1-\text{R}^{2})=\frac{SCE}{SCT},\\ & \implies (1-\text{R}^{2})\cdot\frac{n-1}{n-k-1}=\frac{\frac{1}{n-k-1}}{\frac{1}{n-1}}\frac{SCE}{SCT}=\frac{\frac{1}{n-k-1}}{\frac{1}{n-1}}\left(1-\frac{SCR}{SCT}\right)\\ & =\frac{\frac{1}{n-p}}{\frac{1}{n-1}}\left(1-\frac{SCR}{SCT}\right) \end{aligned} \] Por otro lado
\[ 1-\text{R}^{2}_{\text{adj}}=\frac{\frac{1}{n-p}}{\frac{1}{n-1}}\frac{SCE}{SCT}=\frac{\frac{1}{n-p}}{\frac{1}{n-1}}\left(1-\frac{SCR}{SCT}\right) \]
Por lo tanto:
\[ 1-\text{R}^{2}_{\text{adj}}=\left(1-\text{R}^2\right)\cdot\frac{n-1}{n-p}=\left(1-\text{R}^2\right)\cdot\frac{n-1}{n-k-1} \]
\(_\blacksquare\)
Pregunta 12: Algunas Propiedades del Modelo de Regresión Múltiple IV
Con las mismas consideraciones anteriores, demuestre que una forma equivalente de hacer a prueba de significancia de la regresión es (re)calculando el estadísitico \(F_0\) de la siguiente forma:
\[F_0=\frac{\text{R}^2}{1-\text{R}^2}\left(\frac{n-k-1}{k}\right)\]
\(\underline{\text{Demostración:}}\)
Desarrollando
\[ \begin{aligned} &\frac{1-\text{R}^2}{\text{R}^2}\left(\frac{n-k-1}{k}\right)=\frac{1-\text{R}^2}{\text{R}^2}\left(\frac{n-k-1}{k}\right)=\frac{\frac{SCE}{SCT}}{\frac{SCR}{SCT}}\left(\frac{n-k-1}{k}\right)=\frac{SCE}{SCR}\left(\frac{n-k-1}{k}\right)\\ &\therefore \frac{\text{R}^2}{1-\text{R}^2}\left(\frac{n-k-1}{k}\right)=\frac{SCR}{SCE}\frac{\frac{1}{k}}{\frac{1}{n-k-1}}=\frac{\frac{SCR}{k}}{\frac{SCE}{n-k-1}}=\frac{\frac{SCR}{p-1}}{\frac{SCE}{n-p}}=F_{0}\sim F_{(p-1,n-p)} \end{aligned} \]
Pregunta 13: Ejemplo I
La siguiente información sirve para determinar si existe una relación lineal entre el Gasto nacional de Salud como el Producto Interno Bruto (PIB) y el tiempo (t) siendo el modelo siguiente:
\[\text{PIB}=\beta_0+\beta_1t+\epsilon_i.\]
El PIB en términos porcentuales y el tiempo en años.
Con los siguientes datos, complete la siguiente tabla y estima los siguientes valores :
| Estimación | Valores |
|---|---|
| n | 10 |
| Tiempo Promedio | 1990.19 |
| PIB Promedio | 11.645 |
| SXX | 1983.636 |
| SYY | 102.027 |
| SXY | 445.709 |
a). Estimadores de \(\beta_0,\beta_1 \text{ y }\hat{\sigma^2}\).
\(\underline{\text{Solución:}}\)
\(\hat{\beta_1}\)
\[\hat{\beta_1}=\frac{S_{XY}}{S_{XX}}=\frac{445.709}{1983.636}\]
## [1] "beta_1_estimado= 0.224692937615571"\(\hat{\beta_0}\)
\[\hat{\beta_0}=\bar{Y}-\beta_1\bar{X}=11.645-\frac{445.709}{1983.636}1990.19\]
## [1] "beta_0_estimado= -435.536637513133"\(\hat{\sigma^2}\)
\[ \begin{aligned} \hat{\sigma^2}&=\frac{1}{n-2}\sum_{i=1}^{n}\left(Y_i-\hat{Y_i}\right)^2=\frac{1}{n-2}\left[\frac{S_{YY}S_{XX}-\left(S_{XY}\right)^2}{S_{XX}}\right]\\ &=\frac{1}{8}\left[\frac{102.027\cdot 1983.636-\left(445.709\right)^2}{1983.636}\right] \end{aligned} \]
## [1] "sigma_estimado= 0.234916933537706"\(_\blacksquare\)
b). Los coeficientes de correlación y determinación.
\(\underline{\text{Solución:}}\)
\(R^2\)
\[R^2=1-\frac{SCE}{SCT}=1-\frac{S_{YY}-\hat{\beta_1}S_{XY}}{\hat{\beta_1}S_{XY}}=1-\frac{102.027-0.2246929\cdot445.709}{0.2246929\cdot445.709}\]
## [1] "R^2= 0.981234355518028"\(\rho_{Y\hat{Y}}\)
\[\rho_{Y\hat{Y}}=\sqrt{R^2}=\sqrt{1-\frac{SCE}{SCT}}=\sqrt{0.9812344}\]
## [1] "rho= 0.990572741154342"\(_\blacksquare\)
c). El intervalo de confianza para \(\hat{\beta_1}\) con \(t_{n-2}^{0.975}=2.7515236\).
\(\underline{\text{Solución:}}\)
Sea \(S^2=\frac{S_{YY}-\hat{\beta_1}S_{XY}}{n-2}\), entonces tendremos que:
## [1] "S^2= 0.234916933537706"\[ \begin{aligned} \beta_1&\in\left(\hat{\beta_1}-t_{n-2}^{0.975}\sqrt{\frac{S^2}{S_{XX}}},\hat{\beta_1}+t_{n-2}^{0.975}\sqrt{\frac{S^2}{S_{XX}}}\right)\\ &=\left(0.2246929-2.7515236\cdot0.01088244,0.2246929+2.7515236\cdot0.01088244\right)\\ \therefore \beta_1&\in (0.1947496,0.2546362) \end{aligned} \]
\(_\blacksquare\)
d). Interprete todos los valores
| Fuente de Variación | Grados de Libertad | Suma de Cuadrados | Cuadrado Medio | F0 | p-value |
|---|---|---|---|---|---|
| Regresión | 1 | 100.1477 | 100.1477 | 426 | 3.172e-08 |
| Error | 8 | 1.879335 | 0.2349169 | ||
| Total | 9 | 102.027 | 11.33633 |
Se rechaza la hipótesis de pendiente nula del modelo sii \(F_0=426>F_{(1,8)}^{0.975}=7.570882\). Además tenemos un \(\text{p-value}=3.172e-08<\alpha=0.025\), por lo que tenemos mayor evidencia para garantizar que sí existe un MRL no trivial para los datos que poseemos sobre el Tiempo y el PIB.
\(\underline{\text{Solución:}}\)
\(\hat{\beta_1}\)
Donde \(\hat{\beta_1}=0.2246929\), nos quiere decir que estimador insesgado de varianza mínima de ‘el coeficiente de pendiente’ del MRL es positivo, por lo que se espera que la tasa del cambio inmediato del PIB (la derivada) con respecto al tiempo sea del \(0.2246929\).
\(\hat{\beta_0}\)
Donde \(\hat{\beta_0}=-435.5366\), nos quiere decir que estimador insesgado de varianza mínima de ‘la ordenada al origen’ del MRL es negativa, por lo que, (a menos que tengamos garantía de que uno de los datos del modelo en el PIB tiene asociada una coordenada t=0, en cuyo caso \(\hat{\beta_0}\) sería el valor esperado de la variable \(Y=PIB\)) no tiene interpretación sobre el modelo, simplemente es la ordenada al origen que mejor ajusta a los datos.
\(\hat{\sigma^2}\)
Donde \(\hat{\sigma^2}=0.2349169\), nos quiere decir que estimador insesgado de varianza mínima de ‘la varianza de los errores del modelo’ del MR es de \(0.2349169\), por lo que se espera que la dispersión de los residuos no es mucha, lo que da una ‘buena señal’ del ajuste.
\(R^2\)
Donde \(R^2=0.9812344\), nos quiere decir que ‘el coeficiente de variabilidad explicada por el modelo’ (ó el coefiente de determinación del modelo) es muy cercano a \(1\), por lo que se puede decir que el MRL explica un \(98.12344\%\) de la varibilidad de los datos. Esto nos indica que el modelo es adecuado y que está funcionando para ‘predecir’ ó ‘justificar’ el comportamiento de los datos.
\(\rho_{Y\hat{Y}}\)
Donde \(\rho_{Y\hat{Y}}=0.9905727\), nos quiere decir que ‘el coefiente de correlación el modelo’ (ó bien, de las variables \(X,Y\) del modelo) es muy cercano a \(1\), por lo que se puede decir que existe una ‘fuerte’ dependencia lineal entre el \(PIB\) y \(T\), lo cual da mayor peso a nuestra suposición de la existencia de un MRL.
\(\beta_1\in(0.1947496,0.2546362)\)
Este es el intervalo, con un \(97.5\%\) de probabilidad, podemos asegurar que el valor real de \(\beta_1\) estará- Nótese que este es un intervalo que no incluye al cero y tampoco valores negativos, lo cual nos dice dos cosas: 1). Con un \(97.5\%\) de probabilidad la pendiente del modelo lineal no es cero, lo cual implica que el modelo no sea una constante (es decir un caso trivial y descartable) y 2). Con un \(97.5\%\) de probabilidad la pendiente del modelo lineal no es negativa, lo cual implica una correlación positiva (eso y lo que obtuvimos de \(\rho_{Y\hat{Y}}\)), entre la varible \(PIB\) y \(T\), por lo que si aumenta el tiempo se espera que aumente el \(PIB\) de cierta nación.
\(_\blacksquare\)
Pregunta 14: Ejemplo II
El archivo BGSall.csv es un estudio de cohorte que se realizó entre Junio de 1928 a Junio de 1929 en Berkeley, California en niños que nacieron en ese periodo de tiempo. La base contiene mediciones periódicas hasta los 18 años. Los datos incluye información del Sexo (SEX,0=Hombre, 1=Mujer), Peso (WT2, WT9, WT18) en kilogramos, Talla HT2, HT9, HT18 en centímetros, Circunferencia de Pierna (LG2, LG9, LG18) en centímetros, Fuerza (ST2, ST9, ST18) en kilogramos. Otras dos variables han sido incluidas y son IMC a las 18 años (BMI18) y SOMA, escala que de 1: Muy delgado a 7: Obeso. Con base en dichas medidas, realizar lo siguiente:
I.
Efectuar un análisis exploratorio de datos de todas las variables que contiene usando las tablas y las gráficas que únicamente considere necesarias y útiles con base en lo que requiera visualizar y/o analizar. Notas: 1) Es importante recordar que tienen una variable implícita dentro de la base que es el tiempo, tienen mediciones a los 2, 9 y 18 años. 2) El análisis no es sólo poner gráficos ó estadísticos, sino explicarlos e interpretarlos. No utilice algo que no sepa interpretar.
\(\underline{\text{Solución:}}\)
DATA<-read.csv(file="BGSall.csv", header=TRUE, sep=",",
stringsAsFactors=FALSE)#Obtención de Datos como DF
head(DATA)#Primeras seis onservaciones ## Sujeto Sexo WT2 HT2 WT9 HT9 LG9 ST9 WT18 HT18 LG18 ST18 BMI18
## 1 1 0 13.6 90.2 41.5 139.4 31.6 74 110.2 179.0 44.1 226 34.4
## 2 2 0 12.7 91.4 31.0 144.3 26.0 73 79.4 195.1 36.1 252 20.9
## 3 3 0 12.6 86.4 30.1 136.5 26.6 64 76.3 183.7 36.9 216 22.6
## 4 4 0 14.8 87.6 34.1 135.4 28.2 75 74.5 178.7 37.3 220 23.3
## 5 5 0 12.7 86.7 24.5 128.9 24.2 63 55.7 171.5 31.0 200 18.9
## 6 6 0 11.9 88.1 29.8 136.0 26.7 77 68.2 181.8 37.0 215 20.6
## Soma
## 1 7.0
## 2 4.0
## 3 6.0
## 4 2.0
## 5 1.5
## 6 3.0## 'data.frame': 136 obs. of 14 variables:
## $ Sujeto: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Sexo : int 0 0 0 0 0 0 0 0 0 0 ...
## $ WT2 : num 13.6 12.7 12.6 14.8 12.7 11.9 11.5 15 13.2 16.9 ...
## $ HT2 : num 90.2 91.4 86.4 87.6 86.7 88.1 82.2 91.7 83.8 91 ...
## $ WT9 : num 41.5 31 30.1 34.1 24.5 29.8 26 32.1 30.1 37.9 ...
## $ HT9 : num 139 144 136 135 129 ...
## $ LG9 : num 31.6 26 26.6 28.2 24.2 26.7 26.5 26.3 27.6 29 ...
## $ ST9 : int 74 73 64 75 63 77 45 61 70 61 ...
## $ WT18 : num 110.2 79.4 76.3 74.5 55.7 ...
## $ HT18 : num 179 195 184 179 172 ...
## $ LG18 : num 44.1 36.1 36.9 37.3 31 37 39.1 36.4 37.3 33.9 ...
## $ ST18 : int 226 252 216 220 200 215 152 166 189 183 ...
## $ BMI18 : num 34.4 20.9 22.6 23.3 18.9 20.6 26.3 20.6 21.8 19.4 ...
## $ Soma : num 7 4 6 2 1.5 3 6 4 4 3 ...Debido a que Sujeto ,Sexo y Soma son categóricas, (Sujero más bien sería una varible de identificación de ) las convertiremos a factores
## Sujeto Sexo WT2 HT2 WT9 HT9 LG9 ST9 WT18 HT18 LG18 ST18 BMI18
## 1 1 0 13.6 90.2 41.5 139.4 31.6 74 110.2 179.0 44.1 226 34.4
## 2 2 0 12.7 91.4 31.0 144.3 26.0 73 79.4 195.1 36.1 252 20.9
## 3 3 0 12.6 86.4 30.1 136.5 26.6 64 76.3 183.7 36.9 216 22.6
## 4 4 0 14.8 87.6 34.1 135.4 28.2 75 74.5 178.7 37.3 220 23.3
## 5 5 0 12.7 86.7 24.5 128.9 24.2 63 55.7 171.5 31.0 200 18.9
## 6 6 0 11.9 88.1 29.8 136.0 26.7 77 68.2 181.8 37.0 215 20.6
## Soma
## 1 7
## 2 4
## 3 6
## 4 2
## 5 1.5
## 6 3Primero, vamos a buscar posibles ‘Missing Values’ en la base de datos.
## [1] 0lo cual quiere decir que no tenermos datos faltantes ó perdidos.
Ahora, desplegaremos la información estadística básica de cada varible eso es: valor mínimo, primer cuartil (), mediana (), media (ó valor esperado), tercer cuartil (), valor máximo, la cantidad de varones y mujeres. De esta forma tendremos una idea del comportamiento de cada variable del conjunto.
## Sujeto Sexo WT2 HT2 WT9
## 1 : 1 0:66 Min. :10.10 Min. :80.90 Min. :19.90
## 2 : 1 1:70 1st Qu.:12.20 1st Qu.:85.97 1st Qu.:27.88
## 3 : 1 Median :13.20 Median :87.70 Median :30.90
## 4 : 1 Mean :13.21 Mean :87.80 Mean :31.63
## 5 : 1 3rd Qu.:14.10 3rd Qu.:89.70 3rd Qu.:34.12
## 6 : 1 Max. :18.60 Max. :98.20 Max. :66.80
## (Other):130
## HT9 LG9 ST9 WT18
## Min. :121.4 Min. :21.80 Min. : 22.00 Min. : 42.90
## 1st Qu.:132.5 1st Qu.:26.30 1st Qu.: 57.00 1st Qu.: 56.88
## Median :135.7 Median :27.30 Median : 64.00 Median : 65.10
## Mean :135.5 Mean :27.68 Mean : 64.57 Mean : 64.87
## 3rd Qu.:139.4 3rd Qu.:28.65 3rd Qu.: 73.00 3rd Qu.: 70.85
## Max. :152.5 Max. :40.40 Max. :121.00 Max. :110.20
##
## HT18 LG18 ST18 BMI18
## Min. :153.6 Min. :30.00 Min. : 77.0 Min. :16.60
## 1st Qu.:166.2 1st Qu.:34.17 1st Qu.:124.0 1st Qu.:19.98
## Median :172.5 Median :35.75 Median :150.5 Median :21.30
## Mean :172.6 Mean :35.84 Mean :167.1 Mean :21.72
## 3rd Qu.:179.2 3rd Qu.:37.23 3rd Qu.:214.0 3rd Qu.:22.70
## Max. :195.1 Max. :44.10 Max. :260.0 Max. :36.90
##
## Soma
## 4 :29
## 5 :21
## 3 :18
## 2 :13
## 4.5 :12
## 5.5 :11
## (Other):32Otra forma de visualizar mejor lo anteriormente mencionado es a través de las densidades de cada varible en un solo plot:
DATA[,-c(1,2)]%>%
keep(is.numeric) %>%
gather() %>%
ggplot(aes(value)) +
facet_wrap(~ key, scales = "free") +
geom_density()+
theme(axis.text.x=element_blank(),
axis.title.x=element_blank(),
axis.title.y = element_blank())
Como puede observarse, cada variable ocupada tiene una densidad particular, algunas asemejan distribuciones ‘lognormal’ (BMI18,WT2,WT9,WT18) y otras normal (HT9,HT18,LG18,ST9)
Para analizar los outliers de las variables recurrimos a graficar los ‘boxplot’ de la base de datos.
DATA[,-c(1,2)]%>%
keep(is.numeric) %>%
gather() %>%
ggplot()+
facet_wrap(~ key, scales = "free")+
geom_boxplot(aes(y=value))+
theme(axis.text.x=element_blank(),
axis.title.x=element_blank(),
axis.title.y = element_blank())
Todas las variables aquí presentes (a excepción de HT18 y ST18) tienen valores atípicos, algunas por una sola cola (como lo son BMI18, HT2, LG18, WT18, WT2 y WT9) lo cual puede explicar el por qué sus densidades no son simétricas, (en términos de las medidas de Peso, Estatura (a los 2 años), Circunferencia de Pierna (a los 18 años) e Índice de Masa Corporal (a los 18 años) quiere decir que encontraremos más individuos que se encuentresn por debajo de la mediana de cada medición), aquí también coincide con el hecho de que la barra de la media en cada boxplot de estas medidas se encuentre ‘más hacia arriba’ lo que confirma lo anteriormente mencionado. Por otro lado tenemos variables que presentan valores atípicos por las dos colas (como LG9 y ST9), lo cual dota de cierta ‘simetría’ a las densidades de estas medidas y viendo el área acumulada por la ‘caja’ uno podría suponer que las densidades están ‘centradas’ (es decir, que es más probable encontrar individuos con LG9 y ST9 cercanos a sus medianas en la base de datos), no obstante si uno se fija en la cola superior de ST9 encontrará valores atípicos más alejados que para el caso de la cola inferior, por lo que no tenemos razón para creer que ST9 tenga densidad simétrica (es más probable encontrar individuos con ST9 menor a la mediana en la base de datos). Para el caso de ST18 un tiene que observar la posición en la que se encuentra la barra de mediana de la caja, como está más cerca de la ‘zona inferior’ de la caja y la región restante de la caja acumula un área mayor nos daría la impresión de que existe un primer valor de la variable en el cual encontraremos un 50% de las observaciones antes de este y otro valor distinto y mayor al anterior donde la proporción de observaciones con lecturas mayores a esta sea a lo mucho de 50%, es decir pareciera que tenemos ‘dos medias’.
\(_\blacksquare\)
II.
Ajustar un modelo de regresión utilizando como variables respuesta a BMI18 y como predictoras: SEX, HT2, WT2, HT9, WT9 y ST9.
##
## Call:
## lm(formula = BMI18 ~ Sexo + HT2 + WT2 + HT9 + WT9 + ST9, data = DATA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5304 -1.1631 -0.0903 1.0576 11.3718
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.955227 6.664366 5.245 6.22e-07 ***
## Sexo1 -1.047965 0.432792 -2.421 0.0169 *
## HT2 0.008348 0.095188 0.088 0.9303
## WT2 -0.395456 0.188740 -2.095 0.0381 *
## HT9 -0.147330 0.064834 -2.272 0.0247 *
## WT9 0.428558 0.051848 8.266 1.46e-13 ***
## ST9 -0.027774 0.016643 -1.669 0.0976 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.305 on 129 degrees of freedom
## Multiple R-squared: 0.3717, Adjusted R-squared: 0.3425
## F-statistic: 12.72 on 6 and 129 DF, p-value: 3.025e-11a). Obtener la ecuación de regresión e interpretarlo, así como los coeficientes estandarizados.
\(\underline{\text{Solución:}}\)
Ocupando la información del modelo que previamente se imprimió, tenemos que la ecuación de regresión es:
\[ \begin{aligned} \hat{\text{BMI18}}=&34.955227+-1.047965\cdot\text{Sexo}+ 0.008348\cdot\text{HT2} +-0.395456\cdot\text{WT2}+\\&-0.147330\cdot\text{HT9}+0.428558\cdot\text{WT9}+-0.027774\cdot\text{ST9} \end{aligned} \] A menos que tengamos una observación en la base de datos donde todos los valores de las variables regresoras fueran cero (en tal caso \(\hat{\beta_0}=34.955227\) es el valor esperado de \(\hat{\text{BMI18}}\) ), el intercepto del hiper-plano \(\hat{\beta_0}=34.955227\) no tiene interpretación alguna. Por otro lado, para los demás coeficientes podemos decir que son los cambios esperados en sus respectivos regresores cuando los demás coeficientes son constantes. Por ejemplo, cuando los demás coeficientes son constantes, tenemos que se espera un cambio inmediato de \(0.008348\) centímetros en la variable de estaura a los dos años (HT2). Nótese que para lavariable Sexo su coeficiente no tiene una interpretación adecuada ya que hablamos de una variable dicotómica. Ahora, no todos los coeficientes y regresores tienen el mismo impacto en el modelo, eso puede visualizarse en la columna del p-value, de tal forma que si ocupamos un nivel de significancia del \(\alpha=0.05\), podríamos decir que con una probabilidad del \(95\%\) ‘las pendientes’ de altura a los dos años (HT2) y fuerza a los 9 años (ST9) son insignificantes y, por ende, descartables (en el sentido de que con una probabilidad del \(95\%\) ‘las pendientes’ son iguales a cero).
Además, los coeficientes estandarizados, (los estimadores de los coeficientes del modelo donde la base de datos ha sido estandarizada en el sentido de que las varianzas de las variables dependientes e independientes son iguales a 1, es decir los coeficientes estandarizados no tienen unidad y se refieren a cuántas desviaciones estándar cambiará una variable dependiente, por aumento de la desviación estándar en la variable predictora. Coeficiente estandarizado) son:
modelostd<-lm(scale(BMI18)~scale(as.numeric(Sexo))+scale(HT2)+scale(WT2)+
scale(HT9)+scale(WT9)+scale(ST9),data=DATA)
summary(modelostd)##
## Call:
## lm(formula = scale(BMI18) ~ scale(as.numeric(Sexo)) + scale(HT2) +
## scale(WT2) + scale(HT9) + scale(WT9) + scale(ST9), data = DATA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6486 -0.4091 -0.0318 0.3720 3.9998
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.929e-16 6.953e-02 0.000 1.0000
## scale(as.numeric(Sexo)) -1.849e-01 7.636e-02 -2.421 0.0169 *
## scale(HT2) 9.867e-03 1.125e-01 0.088 0.9303
## scale(WT2) -2.235e-01 1.067e-01 -2.095 0.0381 *
## scale(HT9) -2.848e-01 1.253e-01 -2.272 0.0247 *
## scale(WT9) 8.997e-01 1.089e-01 8.266 1.46e-13 ***
## scale(ST9) -1.509e-01 9.043e-02 -1.669 0.0976 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8108 on 129 degrees of freedom
## Multiple R-squared: 0.3717, Adjusted R-squared: 0.3425
## F-statistic: 12.72 on 6 and 129 DF, p-value: 3.025e-11Dado que esta estandarización se ocupa para responder a la pregunta de cuál de las variables independientes tiene un mayor efecto sobre la variable dependiente Coeficiente estandarizado, se puede decir que, (de forma similar al caso sin estandarización) las variables de Sexo, Peso a los 9 años (WT9), Estatura a los a los 9 años (HT9) y Estatura a los a los 12 años (HT12) son significativas para el modelo en términos de la magintud de sus betas estandarizadas y sus p-values, (esas betas toman valores en el intervalo \((-2.848e-01 ,8.997e-01)\)).
\(_\blacksquare\)
b). Realizar la prueba ANOVA y obtener \(R2\) y \(R_{aj}^2\) e interpretarlos.
\(\underline{\text{Solución:}}\)
Nótese que al momento de ajustar un modelo de regresión múltiple en R, se estima automáticamente el valor de la estadística de Fisher en la prueba de hipótesis de la significancia del modelo a través del análisis de varianzas (ANOVA), es decir:
##
## Call:
## lm(formula = BMI18 ~ Sexo + HT2 + WT2 + HT9 + WT9 + ST9, data = DATA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5304 -1.1631 -0.0903 1.0576 11.3718
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.955227 6.664366 5.245 6.22e-07 ***
## Sexo1 -1.047965 0.432792 -2.421 0.0169 *
## HT2 0.008348 0.095188 0.088 0.9303
## WT2 -0.395456 0.188740 -2.095 0.0381 *
## HT9 -0.147330 0.064834 -2.272 0.0247 *
## WT9 0.428558 0.051848 8.266 1.46e-13 ***
## ST9 -0.027774 0.016643 -1.669 0.0976 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.305 on 129 degrees of freedom
## Multiple R-squared: 0.3717, Adjusted R-squared: 0.3425
## F-statistic: 12.72 on 6 and 129 DF, p-value: 3.025e-11Al llamar a la fución summary() en el modelo, parte de la respuesta que arroja es la siguiente información:
Multiple R-squared: 0.3717, Adjusted R-squared: 0.3425, F-statistic: 12.72 on 6 and 129 DF, p-value: 3.025e-11
Es decir, el valor de la estadística de Fisher valuada en el modelo es de 12.72 para una \(F_{(6,129)}\) y con un p-value de \(3.025e-11<\alpha=0.05\), por lo que con un bivel de confianza del \(95\%\) podemos asegurar la existencia de un modelo de regresión lineal múltiple del Índice de Masa Corporal a la edad de 18 años en términos del Sexo, Estatura a los 2 y 9 años, Peso a los 2 y 9 años y Fuerza a los 9 años (ya que con esa probabilidad es que no se anulan esos coeficientes de la regresión). Al mismo tiempo, obtenemos un coeficiente de variabilidad explicada por el modelo de \(0.3717\), por lo que tenemos un ajuste del \(37.17\%\) del modelo a la variable de Índice de Masa Corporal a la edad de 18 años, no es un resultado adecuado, esto nos indica que no es un modelo fiable. El ajuste del coeficiente de determinación nos dá in resultado totalmente análogo a lo anteriormente mencionado al tener un valor de \(0.3425\).
Ahora, ocuparemos la función nativa de R para realizar el análisis de varianzas individual, (es decir de cada una de las variables regresoras del modelo):
## Df Sum Sq Mean Sq F value Pr(>F)
## Sexo 1 5.2 5.2 0.987 0.3224
## HT2 1 4.5 4.5 0.842 0.3606
## WT2 1 17.8 17.8 3.352 0.0694 .
## HT9 1 13.7 13.7 2.569 0.1114
## WT9 1 349.7 349.7 65.796 3.39e-13 ***
## ST9 1 14.8 14.8 2.785 0.0976 .
## Residuals 129 685.6 5.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Siendo así que en las pruebas individuales de ANOVA, la varible aleatoria de Peso a los 9 años es la de mayor significancia al tener un valor del estadísitico de Fisher de \(F_{ind}=65.796\), con lo cual se obtiene un p-value de \(3.39e-13\), por lo que WT9 es la varible aleatoria regresora que más aporta a nuestra conclusión de existencia del MRLM, (es decir con una probabilidad del \(95\%\) podemos decir que el coeficiente del regresor de WT9 no es cero, desafortunadamente esto no es algo que podamos concluír sobre los demás regresores al ver que sus p-values son mayores a \(\alpha=0.05\))
\(_\blacksquare\)
c). Calcule las correlaciones parciales y concluya.
\(\underline{\text{Solución:}}\)
Ocuparemos la fución cor() en la base de datos DATA a la cual extraemos las variables Sexo, HT2, WT2, HT9, WT9 y ST9, ajustaremos a que el método sea el de Pearson:
D<-DATA[,c('Sexo','HT2','WT2','HT9','WT9','ST9')]
#D%>%mutate(Sexo=as.numeric(Sexo),ST9=as.numeric(ST9))
D<-lapply(D, as.numeric)
DE<-cbind(D$Sexo,D$HT2,D$WT2,D$HT9,D$WT9,D$ST9)
colnames(DE)<-c('Sexo','HT2','WT2','HT9','WT9','ST9')
CP<-pcor(DE,method = "pearson")
#Mostramos un 'HeatMap' de las correlaciones recién obtenidas para evidenciar
#de una manera práctica cúales son las variables aleatorias con mayor
#correlación.
P1<-ggcorrplot(CP$estimate,
hc.order = TRUE,
type = "full",
lab = TRUE,
lab_size = 3,
method="square",
#colors = c("#A40000", "#04004C", "#04894C"),
title = "Correlograma de v.a's",
ggtheme = theme_bw)
#Mostramos un 'HeatMap' de los p-values las correlaciones recién obtenidas para
#evidenciar de una manera práctica cúalespares de las variables aleatorias
#tienen el menor p-value en la pruea de correlación de Pearson,
#(mientras menor sea el p-value menor razón tendremos para suponer no
#correlación entre los pares de variables).
P2<-ggcorrplot(CP$p.value,
hc.order = TRUE,
type = "full",
lab = TRUE,
lab_size = 3,
method="square",
#colors = c("#A40000", "#04004C", "#04894C"),
title = "Pares de p-values",
ggtheme = theme_bw)Para visualizar las correlaciones, primero veámos el gráfico de dispersión de las variables:
Visualmente, puede notarse que existen correlaciones positivas (APARENTEMENTE) entre todos los pares generados por las variables HT2,WT2,HT9,WT9,ST9 eso tiene sentido ya que se espera que al paso del tiempo la estatura y el peso de un individuo promedio aumente, no sólo eso sino que también (en promedio) el peso de un individuoen el presente depende del peso que tuvo en el pasado, (lo mismo aplica para la estatura). De igual forma la fuerza del individuo (casi siempre) es influenciada por su peso y (por lo anteriormente mencionado) estatura. En resumen, existe una tendencia de crecimiento positivo en cada una de las variables anteriormente mencionada cuando hay un incremento encualquiera de las otras variables. Mientras que para el caso de Sexo, al ser dicotómica no podemos garantizar una correlación fuerte negativa/posiva ya que hablamos de correlaciones lineales, ningún modelo lineal ajustaría correctamete a las variables Sexo vs HT2,WT2,HT9,WT9,ST9, por lo que la única conclusión lógica que podemos llevar a cabo es el decir que las correlaciones entre Sexo y HT2,WT2,HT9,WT9,ST9 son muy cercanas al cero ó son nulas, tal como veremos a a continuación:
Aquí tenemos las correlaciones parciales recién obtenidas:
Aquí tenemos los p-values de las pruebas de correlaciones parciales ocuapndo el método de Pearson:
\(_\blacksquare\)
Como puede verse en el HeatMap de correlaciones parciales, las variables con mayor correlación (en cuanto a la mangnitud) son HT9 y HT2 con un valor de 0.56, mientras que las que tienen menor correlación parcial son las de Sexo y HT2 con -0.05, seguidas por Sexo y HT9 con -0.05. Por ende, podemos concluir que exsiten variables regresoras en el modelo que no son independientes entre sí, (tales como HT9 y HT2), por lo que estamos ocupando variables de más en el modelo. Por lo que sería recomendable descartar algunas de esas variables para ocupar sólo variables regresoras independientes (ó con correlación cercana a cero).
d). Obtener los intervalos de confianza de los parámetros estimados y explicarlos.
\(\underline{\text{Solución:}}\)
Por medio de la función confidint() es que vamos a desplegar los intervalos de confianza del modelo del \(95\%\):
## 2.5 % 97.5 %
## (Intercept) 21.76961576 48.140838994
## Sexo1 -1.90425458 -0.191676147
## HT2 -0.17998467 0.196680243
## WT2 -0.76888284 -0.022029036
## HT9 -0.27560591 -0.019054185
## WT9 0.32597546 0.531141191
## ST9 -0.06070244 0.005155202Como puede observarse los coeficientes de las variables ST9 y HT2 contienen al cero, por lo que es muy probable que el valor real del coefiente de esos estimadores sea cero, por lo cual no son significativos esos regresores. Por lo tanto, podemos descartar estos regresores con el objeto de mejorar nuestro modelo.
\(_\blacksquare\)
e). Graficar los distintos tipos de residuos y observar su comportamiento. ¿Siguen los supuestos? Menciónelos.
\(\underline{\text{Solución:}}\)
Como se demostró en los ejrcicios anteriores, sí existe una relación lineal entre la variable de Índice de Masa Corporal a los 18 años y las variables SEX, HT2, WT2, HT9, WT9 y ST9, (eso se evidencía en la prueba de la F). Además, (a pesar de existir correlación positiva y negativa no nula entre los regresores), no tenemos cantidades significativas para suponer una correlación fuerte entre las variables, (eso es, ninguna de las correlaciones tiene maginitud mayor a \(0.70\)).
Empezamos con el primer tipo de residuo:
\[e_i=Y_i-\hat{Y_i}\]
CME=mean(modelolineal$residuals^2)
Residuo_e=resid(modelolineal)
#Residuo_d=resid(modelo)/(CME)
#Residuo_t=resid(modelo)/sqrt(CME*(1-(1/136+)))
DATA<-DATA%>%mutate(Residuo_E=Residuo_e,Fitted=modelolineal$fitted.value)
head(DATA)## Sujeto Sexo WT2 HT2 WT9 HT9 LG9 ST9 WT18 HT18 LG18 ST18 BMI18
## 1 1 0 13.6 90.2 41.5 139.4 31.6 74 110.2 179.0 44.1 226 34.4
## 2 2 0 12.7 91.4 31.0 144.3 26.0 73 79.4 195.1 36.1 252 20.9
## 3 3 0 12.6 86.4 30.1 136.5 26.6 64 76.3 183.7 36.9 216 22.6
## 4 4 0 14.8 87.6 34.1 135.4 28.2 75 74.5 178.7 37.3 220 23.3
## 5 5 0 12.7 86.7 24.5 128.9 24.2 63 55.7 171.5 31.0 200 18.9
## 6 6 0 11.9 88.1 29.8 136.0 26.7 77 68.2 181.8 37.0 215 20.6
## Soma Residuo_E Fitted
## 1 7 8.8778885 25.52211
## 2 4 0.2059669 20.69403
## 3 6 0.8947258 21.70527
## 4 2 0.8839250 22.41608
## 5 1.5 -1.5157883 20.41579
## 6 3 -0.9803251 21.58033Revisaremos la ‘normalidad’ de los residuos:
Para tener un resultado mas exacto, ocuparemos el test de Shapiro-Wilk de normalidad de una muestra:
##
## Shapiro-Wilk normality test
##
## data: DATA$Residuo_E
## W = 0.89993, p-value = 4.454e-08Dado a que p-value = \(4.454e-08<\alpha=0.05\), rechazamos la hipótesis de normalidad de los residuos, por lo que aquí se está violando el principio de normalidad de los residuos.
Revisaremos la variabilidad constante de los residuos (homocedasticidad):
Al representar los residuos frente a los valores ajustados por el modelo (modelolineal$fitted.values), los primeros se tienen que distribuir de forma aleatoria en torno a cero, manteniendo aproximadamente la misma variabilidad a lo largo del eje X. Si se observa algún patrón específico, por ejemplo forma cónica o mayor dispersión en los extremos, significa que la variabilidad es dependiente del valor ajustado y por lo tanto no hay homocedasticidad.
ggplot(DATA, aes(Fitted,Residuo_E)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Nótese cómo la línea roja (el comportamiento de la varianza) empieza a trazar cierto patrón conforme aumenta el valor de los valores ajustados, ésto nos indica una variabilidad no constante de los residuos, es decir otro principio violado. Para tener una lectura más precisa ocuparemos el test studentizado de Breusch-Pagan en nuestro modelos para checar qué sucede con la homocedasticidad:
##
## studentized Breusch-Pagan test
##
## data: modelolineal
## BP = 30.666, df = 6, p-value = 2.935e-05Dado a que p-value = \(2.935e-05<\alpha=0.05\), rechazamos la hipótesis de homocedasticidad de los residuos.
Por otro lado, el supuesto de la suma de residuos igual a 0 sí se cumple:
## [1] 6.161738e-15Ahora, continuamos con el residuo estandarizado:
\[d_i=\frac{Y_i-\hat{Y_i}}{\sqrt{\text{CME}}}\]
n=dim(DATA)[1]
CME=mean(modelolineal$residuals^2)/(n-2)
#Residuo_e=resid(modelo)
#Residuo_d=resid(modelo)/(sqrt(CME))
#Residuo_t=resid(modelo)/sqrt(CME*(1-(1/136+)))
Residuo_d=MASS::stdres(modelolineal)
DATA<-DATA%>%mutate(Residuo_D=Residuo_d)Revisaremos la ‘normalidad’ de los residuos std:
Para tener un resultado mas exacto, ocuparemos el test de Shapiro-Wilk de normalidad de una muestra:
##
## Shapiro-Wilk normality test
##
## data: DATA$Residuo_D
## W = 0.89197, p-value = 1.692e-08Dado a que p-value = \(1.692e-08<\alpha=0.05\), rechazamos la hipótesis de normalidad de los residuos std, por lo que aquí se está violando el principio de normalidad de los residuos.
Revisaremos la variabilidad constante de los residuos (homocedasticidad):
Al representar los residuos frente a los valores ajustados por el modelo (modelo$fitted.values), los primeros se tienen que distribuir de forma aleatoria en torno a cero, manteniendo aproximadamente la misma variabilidad a lo largo del eje X. Si se observa algún patrón específico, por ejemplo forma cónica o mayor dispersión en los extremos, significa que la variabilidad es dependiente del valor ajustado y por lo tanto no hay homocedasticidad.
ggplot(DATA, aes(Fitted,Residuo_D)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Nuevamente, la línea roja (el comportamiento de la varianza) empieza a trazar cierto patrón conforme aumenta el valor de los valores ajustados, ésto nos indica una variabilidad no constante de los residuos, es decir otro principio violado.
Tampoco se satisface el supuesto de suma cero de los residuos.
## [1] -0.6752011Ahora, continuamos con el residuo studentizado:
\[r_i=\frac{Y_i-\hat{Y_i}}{\sqrt{\text{CME}\left(1-\left(\frac{1}{n}+\frac{(Y_i-\bar{Y})^2}{\sum_{i=1}^n(Y_i-\bar{Y})^2}\right)\right)}}\]
# yu<-(DATA$BMI18-mean(DATA$BMI18))^2
# ye<-sqrt(CME*(1-(1/n+(yu)/(sum(yu)))))
# Residuo_t=resid(modelo)/ye
Residuo_t=MASS::studres(modelolineal)
DATA<-DATA%>%mutate(Residuo_T=Residuo_t)Revisaremos la ‘normalidad’ de los residuos studentizados:
Para tener un resultado mas exacto, ocuparemos el test de Shapiro-Wilk de normalidad de una muestra:
##
## Shapiro-Wilk normality test
##
## data: DATA$Residuo_T
## W = 0.86727, p-value = 1.085e-09Dado a que p-value = \(1.085e-09<\alpha=0.05\), rechazamos la hipótesis de normalidad de los residuos studentizados, por lo que aquí se está violando el principio de normalidad de los residuos.
Revisaremos la variabilidad constante de los residuos (homocedasticidad):
Al representar los residuos frente a los valores ajustados por el modelo (modelolineal$fitted.values), los primeros se tienen que distribuir de forma aleatoria en torno a cero, manteniendo aproximadamente la misma variabilidad a lo largo del eje X. Si se observa algún patrón específico, por ejemplo forma cónica o mayor dispersión en los extremos, significa que la variabilidad es dependiente del valor ajustado y por lo tanto no hay homocedasticidad.
ggplot(DATA, aes(Fitted,Residuo_T)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Nuevamente, la línea roja (el comportamiento de la varianza) empieza a trazar cierto patrón conforme aumenta el valor de los valores ajustados, ésto nos indica una variabilidad no constante de los residuos, es decir otro principio violado.
Tampoco se satisface el principio de la suma cero de los residuos:
## [1] -0.02141444En conclusión, ninguno de los tres tipos de residuo satisface todos los principios, como puede observarse en los plots de ‘Fitted’ (valor estimado de Y) vs ‘Residual’, la variabilidad de los residuos no es constante, de hecho pareciera que el error aumenta conforme aumenta el valor de los valores estimados en todos los casos de distintos residuos.
\(_\blacksquare\)
f). ¿Hay observaciones atípicas? ¿A qué individuos corresponden? ¿Tienen sentido dichas observaciones?
\(\underline{\text{Solución:}}\)
Vamos a identificar ahora posibles valores atípicos ó influyentes. Para ello ocuparemos los residuos studentizados, ya que dentro de una mestra de distribución t de Student se considera valor atípico aque mayor ó igual a 3.
ggplot(DATA,aes(x=predict(modelolineal),y=abs(Residuo_T)))+
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
# se identifican en rojo observaciones con residuos estandarizados
#absolutos > 3
geom_point(aes(color = ifelse(abs(Residuo_T) >= 3, 'blue','red') ) )+
labs(title = "Distribución de los residuos studentized",
x = "predicción modelo")+
theme_bw() + theme(plot.title = element_text(hjust = 0.5),
legend.title = element_blank())
Se pueden observar 4 observaciones atípicas, siendo más precisos, hablamos de:
## Sujeto Sexo WT2 HT2 WT9 HT9 LG9 ST9 WT18 HT18 LG18 ST18 BMI18
## 1 1 0 13.6 90.2 41.5 139.4 31.6 74 110.2 179.0 44.1 226 34.4
## 19 19 0 13.3 90.0 31.6 130.3 27.5 68 88.6 172.9 40.4 230 29.6
## 60 60 0 18.6 95.2 66.8 147.5 40.4 121 86.1 188.0 40.0 233 24.4
## 134 134 1 13.3 86.4 41.4 138.2 32.5 44 97.7 162.8 42.5 125 36.9
## Soma Residuo_E Fitted Residuo_D Residuo_T
## 1 7 8.877889 25.52211 3.961330 4.210329
## 19 7 6.696304 22.90370 2.993805 3.091510
## 60 4 -7.530363 31.93036 -4.047793 -4.315434
## 134 7 11.371790 25.52821 5.141866 5.744260Por lo que estas observaciones atípicas corresponden a individuos con sobrepeso u obesidad (Sujetos 1, 19 y 134) su Índice de Masa Corporal a los 18 años es mayor a 25.0 y en los casos de los sujetos 1 y 134 el BMI18 es mayor a 30.0, al mismo tiempo el modelo está prediciendo sus de BMI18 con valores bastante ‘alejados’ del dato real: al Sujeto 1 le predice 25.5221 cuando en realidad tiene 34.4, al Sujeto 19 le predice 22.90370 cuando en realidad tiene 29.6, al Sujeto 134 le predice 25.52821 cuando en realidad tiene 36.9 .Mientras que en el caso del individuo 60 tenemos que sus registros en las variables regresoras son casi todos outliers de sus respectivas varibles, pero al mismo tiempo el modelo predice erroneamente que su BMI a los 18 será de 31.93036 cuando en realidad es de 24.4. Por lo tanto el modelo no está predicendo correctamente el BMI18 de estos individuos con observaciones atípicas en algunas de sus varibles regresoras y en el valor real del BMI18.
\(_\blacksquare\)
g). Al realizar la prueba F parcial, ¿la variable edad es significativa?
\(\underline{\text{Solución:}}\)
Dado que todas las varibles regresoras dependen de forma implícicta del tiempo (2,9 ó 18 años) y, sabiendo que en la prueba del ANOVA para cada regresor al menos la variable WT9 sí es significativa para un valor \(\alpha=0.05\), entonces implicitamente pareciera que sí es significativa la variable aleatoria de la Edad a pesar de no ser como tal un regresor de nuestro modelo. Por otro lado para probar lo anterior dicho de froma explícita:
Para realizar este inciso, vamos a trasponer la base de datos de forma que la edad sí sea una variable aleatoria:
DATA<-read.csv(file="BGSall.csv", header=TRUE, sep=",",stringsAsFactors=FALSE)
#DATA<-DATA%>%mutate(Sujeto=as.factor(Sujeto),
#Sexo=as.factor(Sexo),Soma=as.factor(Soma))
T<-gather(DATA,'Variable','Datos',3:14)
T<-T%>%mutate(Edad=case_when(
Variable=='WT2'|Variable=='HT2'|Variable=='ST2'|Variable=='LG2'~2,
Variable=='WT9'|Variable=='HT9'|Variable=='ST9'|Variable=='LG9'~9, Variable=='WT18'|Variable=='HT18'|Variable=='ST18'|Variable=='LG18'|
Variable=='BMI18'|Variable=='Soma'~18))
F<-T%>%spread(Variable,Datos)
F[is.na(F)] <- 0
F<-F%>%mutate(WT=WT2+WT9+WT18,HT=HT2+HT9+HT18,ST=ST9+ST18,LG=LG9+LG18,
BMI=BMI18)%>%select(Sujeto,Sexo,Edad,HT,WT,ST,LG,Soma,BMI)
F<-F%>%mutate(BMI=case_when(BMI!=0~BMI,BMI==0~WT/(HT/100)^2))
#Base de datos transformada de forma de que la Edad sea variable explícita.
head(F)## Sujeto Sexo Edad HT WT ST LG Soma BMI
## 1 1 0 2 90.2 13.6 0 0.0 0 16.71575
## 2 1 0 9 139.4 41.5 74 31.6 0 21.35613
## 3 1 0 18 179.0 110.2 226 44.1 7 34.40000
## 4 2 0 2 91.4 12.7 0 0.0 0 15.20237
## 5 2 0 9 144.3 31.0 73 26.0 0 14.88775
## 6 2 0 18 195.1 79.4 252 36.1 4 20.90000D<-F%>%select(Sexo,Edad,HT,WT,ST)
CP11<-pcor(D,method = "pearson")
#Mostramos un 'HeatMap' de las correlaciones recién obtenidas para evidenciar
#de una manera práctica cúales son las variables aleatorias con mayor correlación.
P111<-ggcorrplot(CP11$estimate,
hc.order = TRUE,
type = "full",
lab = TRUE,
lab_size = 3,
method="square",
#colors = c("#A40000", "#04004C", "#04894C"),
title = "Correlaciones de v.a's",
ggtheme = theme_bw)
plot(P111)
## Df Sum Sq Mean Sq F value Pr(>F)
## Sexo 1 7.9 7.9 7.959 0.00502 **
## HT 1 1343.5 1343.5 1357.739 < 2e-16 ***
## ST 1 318.7 318.7 322.041 < 2e-16 ***
## WT 1 2039.3 2039.3 2060.840 < 2e-16 ***
## Edad 1 0.3 0.3 0.317 0.57390
## Residuals 402 397.8 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Tal parece que la edad no es un regresor significativo puesto que la prueba de las F parciales nos indica que para toda significancia menor ó igual a 0.1815 no rechazamos la hipótesis de que el coeficiente del regresor Edad es igual a cero, por lo que puede ser descartado. Esto puede deberse a que esta variable aleatoria sólo tiene los valores de 2,9 y 18, tiene una alta correlación (mayor a 0.73) con la estatura (HT) por lo que (sabiendo que el Índice de Masa Corporal se calcula diréctamente en términos de la estatura), por lo que la variable Edad sale sobrando.
\(_\blacksquare\)
III.
Introducir las demás variables y utilizar los distintos criterios y algoritmos de selección de variables.
a). ¿Cuál o cuáles modelos son los mejores candidatos y por qué?
\(\underline{\text{Solución:}}\)
Para lorgrar nuestro cometido de obtener el mejor modelo en términos de: 1)Reducir el \(CMe()\), 2). Encontrar el Sesgo mínimo \(Cp\) y 3). Estabilizar \(R_p^2\), es que vamos a ocupar la función stepAIC() que determinará el mejor modelo por medio del Criterio de Información de Akaike. Una vetaja de esta función es que permite ajustar el tipo de selección, es decir Stepwise, Backward ó Forward Selection.
Stepwise regression
DATA<-read.csv(file="BGSall.csv", header=TRUE, sep=",",stringsAsFactors=FALSE)
DATA<-DATA[,-c(1)]#Sujeto es un ID, no una variable
modelo_completo<-lm(BMI18~.,data=DATA) #El modelo que estima BMI18 en términos
#de las demás variables de la base de datos como regresores lineales.
step.model<-stepAIC(modelo_completo,direction = "both",
trace = FALSE)
summary(step.model)##
## Call:
## lm(formula = BMI18 ~ Sexo + ST9 + WT18 + HT18, data = DATA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.09151 -0.08358 0.02979 0.12202 1.51000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.597508 0.668284 62.245 <2e-16 ***
## Sexo 0.134079 0.063291 2.118 0.0360 *
## ST9 -0.003063 0.001649 -1.857 0.0655 .
## WT18 0.340215 0.002700 126.025 <2e-16 ***
## HT18 -0.242296 0.004209 -57.562 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2591 on 131 degrees of freedom
## Multiple R-squared: 0.9919, Adjusted R-squared: 0.9917
## F-statistic: 4032 on 4 and 131 DF, p-value: < 2.2e-16Backward regression
DATA<-read.csv(file="BGSall.csv", header=TRUE, sep=",",stringsAsFactors=FALSE)
DATA<-DATA[,-c(1)]#Sujeto es un ID, no una variable
modelo_completo<-lm(BMI18~.,data=DATA) #El modelo que estima BMI18 en
#términos de las demás variables de la base de datos como regresores lineales.
bwd.model<-stepAIC(modelo_completo,direction = "backward",
trace = FALSE)
summary(bwd.model)##
## Call:
## lm(formula = BMI18 ~ Sexo + ST9 + WT18 + HT18, data = DATA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.09151 -0.08358 0.02979 0.12202 1.51000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.597508 0.668284 62.245 <2e-16 ***
## Sexo 0.134079 0.063291 2.118 0.0360 *
## ST9 -0.003063 0.001649 -1.857 0.0655 .
## WT18 0.340215 0.002700 126.025 <2e-16 ***
## HT18 -0.242296 0.004209 -57.562 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2591 on 131 degrees of freedom
## Multiple R-squared: 0.9919, Adjusted R-squared: 0.9917
## F-statistic: 4032 on 4 and 131 DF, p-value: < 2.2e-16Forward regression
DATA<-read.csv(file="BGSall.csv", header=TRUE, sep=",",stringsAsFactors=FALSE)
DATA<-DATA[,-c(1)]#Sujeto es un ID, no una variable
modelo_completo<-lm(BMI18~.,data=DATA) #El modelo que estima BMI18 en
#términos de las demás variables de la base de datos como regresores lineales.
fwd.model<-stepAIC(modelo_completo,direction = "forward",
trace = FALSE)
summary(fwd.model)##
## Call:
## lm(formula = BMI18 ~ Sexo + WT2 + HT2 + WT9 + HT9 + LG9 + ST9 +
## WT18 + HT18 + LG18 + ST18 + Soma, data = DATA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.09523 -0.10080 0.02605 0.12194 1.37734
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.7999596 1.0908432 38.319 <2e-16 ***
## Sexo 0.0906496 0.1312795 0.691 0.491
## WT2 0.0142571 0.0227610 0.626 0.532
## HT2 -0.0090977 0.0110073 -0.827 0.410
## WT9 -0.0122575 0.0145195 -0.844 0.400
## HT9 0.0145655 0.0116577 1.249 0.214
## LG9 0.0146574 0.0346915 0.423 0.673
## ST9 -0.0035000 0.0023493 -1.490 0.139
## WT18 0.3452658 0.0063672 54.226 <2e-16 ***
## HT18 -0.2508167 0.0075999 -33.003 <2e-16 ***
## LG18 -0.0052012 0.0228263 -0.228 0.820
## ST18 -0.0004834 0.0013409 -0.360 0.719
## Soma -0.0310598 0.0313917 -0.989 0.324
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2636 on 123 degrees of freedom
## Multiple R-squared: 0.9922, Adjusted R-squared: 0.9914
## F-statistic: 1299 on 12 and 123 DF, p-value: < 2.2e-16#Estadística de Mallow
Cpstep<-ols_mallows_cp(step.model,modelo_completo)
Cpbwd<-ols_mallows_cp(bwd.model,modelo_completo)
Cpfwd<-ols_mallows_cp(fwd.model,modelo_completo)
#Información de Akaike
AICstep<-ols_aic(step.model)
AICbwd<-ols_aic(bwd.model)
AICfwd<-ols_aic(fwd.model)
paste('Estadística de Mallow Step= ',Cpstep)## [1] "Estadística de Mallow Step= 0.570700675558115"## [1] "Estadística de Mallow Backward= 0.570700675558115"## [1] "Estadística de Mallow Forward= 13"## [1] "Información de Akaike Step= 25.4785145923758"## [1] "Información de Akaike Backward= 25.4785145923758"## [1] "Información de Akaike Forward= 37.5866440328671"Como puede observarse, los algoritmos de Backward y de Stepwise nos arrojan el mismo resultado en términos de variables seleccionadas, ocupando sólamente Sexo, ST9, WT18, HT18 (eso tiene todo el sentido al saber que BMI18 se calcula únicamente como el cociente de WT18 sobre (HT18)^2, como este es un modelo lineal, no podemos ocupar funciones cuadráticas ó racionales; por otro lado, la inclusión de Sexo tiene sentido al considerar la diferencia biológica que en promedio existe entre Varones y Mujeres y, finalmente, ST9 puede jugar cierto papel reelevante al considerar que la fuerza de un individuo depene de su Peso, pero para fines prácticos no tiene tanto sentido que ocupen la fuerza a los 9 años, personalmente creo que tendría más sentido ocupar la fuerza a los 18 años) arrojando los dos valores en los estadísticos iguales: \(F_c=4032, \text{p-value}< 2.2e-16\),\(CMe=0.2591\),\(R^2_{adj}=0.9917\),\(C_p=0.5707007\) y \(AIC=25.47851\). En este caso se ocupó el algoritmo que minimiza la entropía esperada por el modelo \(AIC\) (misma interpretación que \(R^2\) y \(C_p\) (Mallow)).NOTA la función ocupada (stepAIC()) no sólo se basa en minimizar \(AIC\) puede ajustarse para que se encuentre la \(R^2\) y \(C_p\) mínima a la hora de seleccionar los modelos.
Por otro lado el Forward nos arroja un pésimo resultado al ocupar las 12 variables de la base de datos como regresores, no sólo se ignora el hecho de que hay alta correlación de por medio, sino que se reduce el coefiente de determinaciónal, se reduce el valor del estadístico \(F_c\),aumenta el Cuadrado Medio del error, aumenta el valor de AIC y mejora el resultado de la \(C_p\):
\(F_c=1299, \text{p-value}< 2.2e-16\),\(CMe=0.2636\),\(R^2_{adj}=0.9914\),\(C_p=13\) y \(AIC=37.58664\)
\(_\blacksquare\)
b). ¿Se obtienen los mismos resultados con todos los algoritmos?
\(\underline{\text{Solución:}}\)
Sólo en el Backward y el Stepwise se obtienen el mismo resultado, en el Forward se tiene un resultado totalmente distinto.
\(_\blacksquare\)
c). ¿Es mejor o peor al modelo obtenido en el inciso II?
\(\underline{\text{Solución:}}\)
Por supuesto que el modelo obtenido por Stepwise Selection es mejor que el que se construyó en el inciso II, el lm() que tiene como regresoras a las variables Sexo, ST9, WT18, HT18 no sólo tiene una estadísitica de Fisher \(F_c\) mucho mayor (\(4032 >12.72\), y por lo tanto un p-value mucho menor),tiene un coeficiente de determinación mucho mayor (\(0.9914>0.3425\)) y su Cuadrado Medio del error es menor (\(0.2591<2.305\)), además de que ocupa menos variables regresoras, por lo que es más eficiente y por lo tanto se ajusta mucho mejor a la muestra de datos.
Ahora, veamos los resultados del primer modelo y comparemos lo que resta:
#Estadística de Mallow
Cplin<-ols_mallows_cp(modelolineal,modelo_completo)
#Información de Akaike
AIClin<-ols_aic(modelolineal)
paste('Estadística de Mallow Backward= ',Cpbwd)## [1] "Estadística de Mallow Backward= 0.570700675558115"## [1] "Información de Akaike Step= 25.4785145923758"## [1] "Estadística de Mallow Lineal= 9747.26656290146"## [1] "Información de Akaike Step= 621.946163212698"Puede notarse que \(C_p\) para el modelo lineal está mas alejado del número de variables que ocupa que el modelo que se obtuvo por Backwards Selection. De manera similar el criterio de Akaike es más favorable para el modelo Backwards. En conclusión, sí es mucho mejor el modelo que sólo ocupa las cuatro variables regresoras Sexo, ST9, WT18 y HT18 que aquel que ocupa las seis variables regresoras Sexo, HT2, WT2, HT9, WT9 y ST9.
\(_\blacksquare\)
IV. Ajuste un solo modelo de regresión lineal “final”. Justifique e interprete dicho modelo. Nota: No olvide asegurarse que cumple los supuestos y que no existen violaciones de ningún tipo.
\(\underline{\text{Solución:}}\)
El modelo de regresión final que vamos a ocupar es el que se obtuvo a través de Stepwise Selection, no sólo por que es el que tiene el mejor coeficiente de determinación, el que minimiza el cuadrado medio del error y también minimiza la información de Akaike, sino que todo modelo con una cantidad menor de variables (ajustado a la base de datos), tendrá peores resultados en cuanto a \(R^2_{adj}\), \(CME\) y \(AIC\), tal y como veremos en en siguienet ejemplo:
#Descartamos la variable ST9 que se ocupa en el modelo bwd.model
prueba1<-lm(formula = BMI18 ~ Sexo + WT18 + HT18, data = DATA)
#Estadística de Mallow
Cpprueba<-ols_mallows_cp(prueba1,modelo_completo)
#Información de Akaike
AICprueba<-ols_aic(prueba1)##
## Call:
## lm(formula = BMI18 ~ Sexo + ST9 + WT18 + HT18, data = DATA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.09151 -0.08358 0.02979 0.12202 1.51000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.597508 0.668284 62.245 <2e-16 ***
## Sexo 0.134079 0.063291 2.118 0.0360 *
## ST9 -0.003063 0.001649 -1.857 0.0655 .
## WT18 0.340215 0.002700 126.025 <2e-16 ***
## HT18 -0.242296 0.004209 -57.562 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2591 on 131 degrees of freedom
## Multiple R-squared: 0.9919, Adjusted R-squared: 0.9917
## F-statistic: 4032 on 4 and 131 DF, p-value: < 2.2e-16##
## Call:
## lm(formula = BMI18 ~ Sexo + WT18 + HT18, data = DATA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08556 -0.07552 0.02465 0.10730 1.57118
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.882836 0.656399 63.807 <2e-16 ***
## Sexo 0.122518 0.063566 1.927 0.0561 .
## WT18 0.339678 0.002709 125.395 <2e-16 ***
## HT18 -0.244859 0.004013 -61.009 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2615 on 132 degrees of freedom
## Multiple R-squared: 0.9917, Adjusted R-squared: 0.9915
## F-statistic: 5277 on 3 and 132 DF, p-value: < 2.2e-16## [1] "Estadística de Mallow Backward= 0.570700675558115"## [1] "Información de Akaike Step= 25.4785145923758"## [1] "Estadística de Mallow Prueba1= 1.90424459522274"## [1] "Información de Akaike Prueba1= 27.0140451676428"Puede que el estadístico de Fisher \(F_c\) sea mayor para el modelo prueba1 y que la estadística de Mallow se acerque más al número de variables que ocupa el modelo prueba1, pero si hablamos del cuadrado medio del error, el coeficiente de determinación ajustado y la información de Akaike, los resultados del modelo Prueba1 son peores.
En cuanto a la interpretación del modelo, (además de lo ya dicho en la sección a) del Inciso II del Ejercio 14), sólo podemos agregar que es un modelo eficiente en términos de predicciones, pues a la hora de evaluar en las observaciones la cantidad de operaciones realizadas es mucho menor al caso del modelolineal (Inciso II del Ejercici 14) ó el modelo Forward y eso, a la larga, es fundamental en el uso de un modelo, (especialmente cuando queramos ocupar bases de datos con más observaciones y analicemos los residuos).
Ahora, analicemos los residuos y el cumplimiento de principios:
Primero, si desplegamos la información básica del modelo veremos que se satisface la hipótesis de exitencia de una relación lineal entre la variable de respuesta BMI18 y las variables regresoras Sexo + ST9 + WT18 + HT18
##
## Call:
## lm(formula = BMI18 ~ Sexo + ST9 + WT18 + HT18, data = DATA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.09151 -0.08358 0.02979 0.12202 1.51000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.597508 0.668284 62.245 <2e-16 ***
## Sexo 0.134079 0.063291 2.118 0.0360 *
## ST9 -0.003063 0.001649 -1.857 0.0655 .
## WT18 0.340215 0.002700 126.025 <2e-16 ***
## HT18 -0.242296 0.004209 -57.562 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2591 on 131 degrees of freedom
## Multiple R-squared: 0.9919, Adjusted R-squared: 0.9917
## F-statistic: 4032 on 4 and 131 DF, p-value: < 2.2e-16eso con un poco más 93% y un poco menos del 94% (eso si consideramos todas las prubas parciales del modelo donde ST9 tiene un p-value del 0.0655, por lo que para niveles de significancia menores a 0.0655 podríamos aceptar la hipótesis de anular el regresor de ST9) de probabilidad al tener un F-statistic=4032 y un p-value<2.2e-16 en la prueba de hiótesis de significancia del modelo.
Ahora veremos qué podemos decir sobre las correlaciones de las variables regresoras:
#Correlaciones entre variables regresoras
D<-DATA[,c('Sexo','HT18','WT18','ST9')]
#D%>%mutate(Sexo=as.numeric(Sexo),ST9=as.numeric(ST9))
D<-lapply(D, as.numeric)
DEe<-cbind(D$Sexo,D$HT18,D$WT18,D$ST9)
colnames(DEe)<-c('Sexo','HT18','WT18','ST9')
CP1<-pcor(DEe,method = "pearson")
#Mostramos un 'HeatMap' de las correlaciones recién obtenidas para evidenciar
#de una manera práctica cúales son las variables aleatorias con mayor
#correlación.
P11<-ggcorrplot(CP1$estimate,
hc.order = TRUE,
type = "full",
lab = TRUE,
lab_size = 3,
method="square",
#colors = c("#A40000", "#04004C", "#04894C"),
title = "Correlaciones de v.a's",
ggtheme = theme_bw)
#Mostramos un 'HeatMap' de los p-values las correlaciones recién obtenidas para
#evidenciar de una manera práctica cúalespares de las variables aleatorias
#tienen el menor p-value en la pruea de correlación de Pearson,
#(mientras menor sea el p-value menor razón tendremos para suponer no
#correlación entre los pares de variables).
P22<-ggcorrplot(CP1$p.value,
hc.order = TRUE,
type = "full",
lab = TRUE,
lab_size = 3,
method="square",
#colors = c("#A40000", "#04004C", "#04894C"),
title = "Pares de p-values",
ggtheme = theme_bw)
Como puede verse en el HeatMap de correlaciones parciales, las variables con mayor correlación (en cuanto a la mangnitud) son HT18 y Sexo con un valor de -0.58. Por ende, podemos concluir que existen variables regresoras en el modelo que no son completamente independientes entre sí, (tales como HT9 y Sexo), pero al no ser una magnitud reelevante (eso es >0.7), no tenemos variables regresoras dependientes unas de otras, por lo que estamos ocupando la cantidad adecuada de regresores. En cuanto a la información de los p-value, podemos decir que en las pruebas parciales de correlación tenemos que casi todos son cero, específicamente para la prueba de correlación entre HT18 y Sexo, eso sólo confirma lo anteriormente mencionado, lo interesante en cuanto a los p-values es que para ciertas variables valores son mayores al 20%, eso conviene para poder justificar que con probabilidad positiva las correlaciones de ciertas variables son nulas, lo que vuelve a confirmar lo anteriormente dicho.
Continuamos con el primer tipo de residuo:
\[e_i=Y_i-\hat{Y_i}\]
D<-DATA[,c('Sexo','HT18','WT18','ST9')]
CME=mean(bwd.model$residuals^2)
Residuo_e=resid(bwd.model)
Residuo_t=MASS::studres(bwd.model)
D<-D%>%mutate(Residuo_E=Residuo_e,Residuo_T=Residuo_t,
Fitted=bwd.model$fitted.value)
head(D)## Sexo HT18 WT18 ST9 Residuo_E Residuo_T Fitted
## 1 0 179.0 110.2 74 -1.09151357 -5.0967223 35.49151
## 2 0 195.1 79.4 73 -0.21497785 -0.8592832 21.11498
## 3 0 183.7 76.3 64 -0.25005810 -0.9769666 22.85006
## 4 0 178.7 74.5 75 -0.11545444 -0.4484602 23.41545
## 5 0 171.5 55.7 63 0.09929878 0.3891000 18.80070
## 6 0 181.8 68.2 77 0.08514675 0.3308374 20.51485Revisaremos la ‘normalidad’ de los residuos:
Para tener un resultado mas exacto, ocuparemos el test de Shapiro-Wilk de normalidad de una muestra:
##
## Shapiro-Wilk normality test
##
## data: D$Residuo_E
## W = 0.83943, p-value = 7.079e-11Dado a que p-value = \(7.079e-11<\alpha=0.05\), rechazamos la hipótesis de normalidad de los residuos, por lo que aquí se está violando el principio de normalidad de los residuos.
Revisaremos la variabilidad constante de los residuos (homocedasticidad):
Al representar los residuos frente a los valores ajustados por el modelo (bwd.modelo$fitted.values), los primeros se tienen que distribuir de forma aleatoria en torno a cero, manteniendo aproximadamente la misma variabilidad a lo largo del eje X. Si se observa algún patrón específico, por ejemplo forma cónica o mayor dispersión en los extremos, significa que la variabilidad es dependiente del valor ajustado y por lo tanto no hay homocedasticidad.
ggplot(D, aes(Fitted,Residuo_E)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Nótese cómo la línea roja (el comportamiento de la varianza) empieza a trazar cierto patrón ‘cuadrático’ conforme aumenta el valor de los valores ajustados, ésto nos indica una variabilidad no constante de los residuos, es decir otro principio violado. Para tener una lectura más precisa ocuparemos el test studentizado de Breusch-Pagan en nuestro modelos para checar qué sucede con la homocedasticidad:
##
## studentized Breusch-Pagan test
##
## data: bwd.model
## BP = 43.867, df = 4, p-value = 6.836e-09Dado a que p-value = \(6.836e-09<\alpha=0.05\), rechazamos la hipótesis de homocedasticidad de los residuos.
Por otro lado, el supuesto de la suma de residuos igual a 0 sí se cumple:
## [1] 2.046974e-15Vamos a identificar ahora posibles valores atípicos ó influyentes. Para ello ocuparemos los residuos studentizados, ya que dentro de una mestra de distribución t de Student se considera valor atípico aque mayor ó igual a 3.
ggplot(D,aes(x=predict(bwd.model),y=abs(Residuo_T)))+
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
# se identifican en rojo observaciones con residuos estandarizados
#absolutos > 3
geom_point(aes(color = ifelse(abs(Residuo_T) >= 3, 'blue','red') ) )+
labs(title = "Distribución de los residuos studentized",
x = "predicción modelo")+
theme_bw() + theme(plot.title = element_text(hjust = 0.5),
legend.title = element_blank())
Se pueden observar 2 observaciones atípicas, siendo más precisos, hablamos de:
## Sexo HT18 WT18 ST9 Residuo_E Residuo_T Fitted
## 1 0 179.0 110.2 74 -1.091514 -5.096722 35.49151
## 134 1 162.8 97.7 44 1.509997 7.957489 35.39000Curiosamente esos dos outliers aparecieron previamente en la sección f) del inciso II del Ejercicio 14, por lo que la conclusión sobre la razón ó significado de estas observaciones atípicas es la misma que se menciona en dicho ejercicio, con la diferencia de que el Cuadrado Medio del Error en este nuevo modelo es mucho menor, es decir, se está cometiendo un error similar al primer modelo pero se está mejorando en cuanto a que la cantidad de outliers es menor y, por ende, el error provocado por la presencia de esas observaciones es menor.
Por lo tanto, sí hay violaciones en los principios fundamentales del modelo: no hay variabilidad constante, los residuos no se distribuyen normal y presencia de outliers.
\(_\blacksquare\)