d = read.csv("datos/kidiq.csv")
d$mom_hs = factor(d$mom_hs, labels = c("no", "si"))Para el modelo
m = lm(kid_score ~ mom_hs, data = d)
summary(m)##
## Call:
## lm(formula = kid_score ~ mom_hs, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57.55 -13.32 2.68 14.68 58.45
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 77.548 2.059 37.670 < 2e-16 ***
## mom_hssi 11.771 2.322 5.069 5.96e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.85 on 432 degrees of freedom
## Multiple R-squared: 0.05613, Adjusted R-squared: 0.05394
## F-statistic: 25.69 on 1 and 432 DF, p-value: 5.957e-07se tiene que:
Por tanto, la puntuación de los segundos es superior a la de los primeros. Pero este resultado se ha obtenido con los datos de 434 madres y sus hijos. ¿Es generalizable? ¿Si repetimos los cálculos con otras madres e hijos obtendremos las mismas conclusiones? ¿Se puede afirmar que si una madre termina secundaria influye positivamente en la capacidad cognitiva de sus hijos?
Para resolver este problema se hace lo siguiente:
\[ kid\_score_i = \beta_0 + \beta_1 mom\_hs_i + U_i, \quad U_i \sim N(0,\sigma) \]
\(U_i\) es una variable aleatoria; \(\beta_0\) y \(\beta_1\) son parámetros desconodidos.
se considera que los datos que tenemos, los 434 datos en este caso, han sido generados por el modelo anterior. Son datos PARTICULARES.
estamos interesados en los parámetros del modelo generador de los datos, ya que es el que indica conclusiones generales. Pero dichos parámetros son desconocidos y no se pueden conocer. El objetivo es inferir propiedades de los parámetros generales a partir de los datos particulares.
El modelo de la población tiene tres parámetros: \(\beta_0\), \(\beta_1\) y \(\sigma\). Lo primero que haremos es estimar un valor para dichos parámetros.
Se denomina estimador de un parámetro a cualquier expresión que asigna un valor a dicho parámetro. Y se denomina estimación al valor asignado.
Por ejemplo, podríamos utilizar los siguientes estimadores (se indica poniendo el símbolo \(\hat \square\)):
\[ \hat \beta_1 = suma(kid\_score_i), \ \hat \beta_0 = media(kid\_score_i), \ \hat \sigma = Var(kid\_score_i) \]
lo que daría lugar a las siguientes estimaciones:
# para beta_0
sum(d$kid_score)## [1] 37670# para beta_1
mean(d$kid_score)## [1] 86.79724# para sigma
var(d$kid_score)## [1] 416.5962Sin embargo, unos estimadores tienen mejores propiedades que otros. Nosotros vamos a utilizar como estimadores los valores calculados con mínimos cuadrados:
\[ \begin{bmatrix} \hat \beta_0 \\ \hat \beta_1 \end{bmatrix} = (X^T X)^{-1} (X^T y) \]
La estimación es:
coefficients(m)## (Intercept) mom_hssi
## 77.54839 11.77126El estimador de \(\sigma\) se define como (no es un resultado de mínimos cuadrados):
\[ \hat \sigma = \sqrt{ \frac{\sum e_i^2}{n-k-1} } \]
donde \(k\) es el número de regresores, 1 en este caso.
La estimación tomará el valor:
sqrt( sum(resid(m)^2)/(434-2) )## [1] 19.85253El hecho de considerar que nuestros datos han sido generados por un modelo con variables aleatorias hace que los estimadores sean variables aleatorias: si consideramos otra muestra diferente, los valores obtenidos serán diferentes, y antes de seleccionar la muestra no sabemos que valores vamos a obtener. En este apartado vamos a obtener la distribución de los estimadores.
Como hemos visto en el apartado anterior, el modelo generador de los datos es:
\[ Y_i = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_k x_{ki} + U_i, \quad U_i \sim N(0,\sigma), \quad i = 1,\cdots,n \]
En forma matricial
\[ Y = X \beta + U, \quad U \sim N(0,\sigma \mathrm{I}) \]
Como \(U \sim Normal\), debido a las propiedades de la distribución normal, \(Y \sim Normal\). Concretamente:
\[ Y \sim N(X \beta, \sigma \mathrm{I}) \]
ya que:
\[ E[Y] = E[ X\beta + U] = X \beta \]
\[ Var[Y] = Var[ X\beta + U] = \sigma^2 I \]
La ecuación de los estimadores es
\[ \hat \beta = (X^T X)^{-1} (X^T y) \]
Si consideramos variables aleatorias en lugar de datos:
\[ \hat \beta = (X^T X)^{-1} (X^T Y) \] obtenemos la distribución de los estimadores:
\[ \hat \beta \sim N(\beta, \sigma^2 (X^T X)^{-1}) \]
ya que
\[Y \sim Normal \Rightarrow \hat \beta \sim Normal\]
\[ E[\hat \beta] = E[(X^TX)^{-1}X^T Y] = (X^TX)^{-1}X^T E[Y] = (X^TX)^{-1}X^T X\beta = \beta \]
\[ Var[\hat \beta] = Var[(X^TX)^{-1}X^T Y] = (X^TX)^{-1}X^T Var[Y] ((X^TX)^{-1}X^T)^T = \sigma^2 (X^TX)^{-1} \]
Si en lugar de utilizar las matrices \((X^TX)^{-1}\) y \((X^T y)\) utilizamos matrices de varianzas y convarianzas, los estimadores son:
\[ \hat \beta_a = S_{XX}^{-1} S_{Xy} = \left( \frac{1}{n-1}X_a^T X_a \right)^{-1} \left( \frac{1}{n-1} X_a^T Y_a \right) = \left( X_a^T X_a \right)^{-1} \left( X_a^T Y_a \right) \]
donde \(\hat \beta_a = [\hat \beta_1 \ \hat \beta_2 \ \cdots \ \hat \beta_k]^T\). Por tanto hay que obtener la distribución de \(Y_a\). El modelo generador de datos es:
\[ Y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \cdots + \beta_k x_{ki} + U_i, \ i = 1,2,\cdots,n \]
Teniendo en cuenta que:
\[ \bar Y = \beta_0 + \beta_1 \bar x_{1} + \beta_2 \bar x_{2} + \cdots + \beta_k \bar x_{k} \]
y restando ambas ecuaciones se obtiene:
\[ Y_i - \bar Y = \beta_1 (x_{1i} - \bar x_{1}) + \beta_2 (x_{2i} - \bar x_{2}) + \cdots + \beta_k (x_{ki} - \bar x_{k}) + U_i, \ i = 1,2,\cdots,n \]
Estas n ecuaciones se pueden expresar en forma matricial de la misma forma que hicimos antes, obteniendo:
\[ Y_a = X_a \beta_a + U, \quad U \sim N(0,\sigma I) \]
Por tanto,
\[ Y_a \sim N \left( X_a \beta_a, \sigma^2I \right) \]
Ahora se puede demostrar que
\[ \hat \beta_a \sim N \left( \beta_a, \frac{\sigma^2}{n-1} S_{XX}^{-1} \right) \]
ya que
\[Y_a \sim Normal \Rightarrow \hat \beta_a \sim Normal\]
\[ E[\hat \beta_a] = E[(X_a^T X_a)^{-1}X_a^T Y_a] = (X_a^T X_a)^{-1}X_a^T E[Y_a] = (X_a^T X_a)^{-1}X_a^T X_a \beta_a = \beta_a \]
\[ Var[\hat \beta_a] = Var[(X_a^T X_a)^{-1}X_a^T Y_a] = (X_a^T X_a)^{-1}X_a^T Var[Y_a] ((X_a^T X_a)^{-1}X_a^T)^T = \sigma^2 (X_a^T X_a)^{-1} \]
Finalmente:
\[ S_{XX} = \frac{1}{n-1}X_a^T X_a \Rightarrow \left( X_a^T X_a \right)^{-1} = \frac{1}{n-1} S_{XX}^{-1} \]
Faltaría el estimador de \(\beta_0\) que se obtiene con la ecuación
\[ \hat \beta_0 = \bar Y - \bar x^T \hat \beta_a \]
dónde \(\bar x = [\bar x_1 \ \bar x_2 \ \cdots \ \bar x_n]^T\). Se puede demostrar que
\[ \hat \beta_0 \sim N \left( \beta_0, \sigma^2 \left( \frac{1}{n} + \frac{1}{n-1} \bar x^T S_{XX}^{-1} \bar x \right) \right) \]
ya que:
\[ \bar Y \sim N \left( \beta_0 + \bar x^T \beta_a, \frac{\sigma^2}{n} \right) \]
El modelo tiene un parámetro más, la varianza de los errores, \(\sigma^2\). Este parámetro también hay que estimarlo. Se puede demostrar que
\[ \frac{\sum e_i^2}{\sigma^2} \rightarrow \chi^2_{n-k-1} \]
donde n es el número de observaciones y k es el número de regresores. Por ello se propone el siguiente estimador
\[ \hat \sigma^2 = \frac{\sum e_i^2}{n-k-1} \]
ya que es un estimador insesgado de \(\sigma^2\). Efectivamente
\[ E[\hat \sigma^2] = E \left[ \frac{\sum e_i^2}{n-k-1} \right] = \sigma^2 \]
ya que \(E[\chi^2_n] = n\). Al término \(\sum e_i^2/(n-k-1)\) también se lo conoce como varianza residual y se representa por \(\hat s_R^2\).
\[ \hat s_R^2 = \frac{\sum e_i^2}{n-k-1} \]
A la raiz cuadrada se le conoce como residual standard error. El término (n-k-1) son los grados de libertad. La distribución en el muestreo de la varianza residual es
\[ \frac{\sum e_i^2}{\sigma^2} \rightarrow \chi^2_{n-k-1} \Rightarrow \frac{(n-k-1)\hat s_R^2}{\sigma^2} \rightarrow \chi^2_{n-k-1} \]