set.seed(seed=20087)
n <- 86
D <- runif(n,120,230) # Se generan 86 valores aleatorios entre 120 y 230
C <- runif(n,210,330) # Se generan 86 valores aleatorios entre 210 y 330
sigmaE <- 25
eps <- rnorm(n,0,sigmaE) # Se genera el error residual para cada observación
Y <- 130+0.75*D+1.3*C+eps # Se calcula el valor de Y a partir de D, C y el error residual eps
fito <- data.frame(Y,D,C) # Se guardan todos los datos en un data.frame.Modelos Lineales
Definición
Un vector de variables aleatorias: \[\mathbf{Y}=\left(\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array}\right)\]
sigue un modelo lineal si:
\[\mathbf{Y} \cong \mathbf{N_n} \left(\mathbf{X} \beta, \mathbf{I_n}\sigma_\varepsilon \right)\]
donde:
- \(\mathbf{X}=\left(\begin{array}{ccccc} 1 & x_{11} & x_{21} & \ldots & x_{p1}\\ 1 & x_{12} & x_{22} & \vdots & x_{p2}\\ 1 & \vdots & \vdots & \ddots\\ 1 & x_{1n} & x_{2n} & \ldots & x_{pn} \end{array}\right)\)
es una matriz de dimensión \(n\times p\), cuya columna \(j\)-ésima (\(j\ge2\)) contiene los valores observados en la variable \(X_j\),
- \(\beta=\left(\begin{array}{c} \beta_{0}\\ \beta_{1}\\ \vdots\\ \beta_{p} \end{array}\right)\) es un vector de parámetros \(p\)-dimensional,
- \(\mathbf{I_{n}}\) es la matriz identidad de orden \(n\),
- \(\sigma_\varepsilon\) es una constante.
De manera alternativa y equivalente, el modelo lineal puede expresarse como: \[\mathbf{Y=}\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\dots+\beta_{p}X_{p}+\varepsilon\] siendo las \(X_j\) las columnas de la matrix \(\mathbf{X}\) anterior, y \[\mathbf{\varepsilon}=\left(\begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{array}\right)\]
un vector de variables aleatorias independientes tales que \(\varepsilon_i\approx N \left(0,\sigma_\varepsilon\right)\;\; \forall i\).
la variable \(\mathbf{Y}\) se conoce habitualmente como variable respuesta o dependiente y las variables \(\mathbf{X}\) como variables explicativas o independientes.
En la práctica la matriz \(\mathbf{X}\) (exceptuando la primera columna de unos) representa los datos observados en \(n\) sujetos sobre los que se han medido las \(p\) variables explicativas \(X_1, X_2, \dots, X_p\). De esta forma el valor \(x_{ij}\) representa el valor de la variable \(X_j\) medido en el sujeto \(i\). A su vez, el vector \(\left(y_1, y_2, \dots, y_n\right)\) representa el conjunto de valores observados en la variable respuesta \(\mathbf{Y}\), siendo \(y_i\) el valor observado en el sujeto \(i\). El objetivo del modelo es predecir o explicar el valor de la \(Y\) a partir de los valores de las \(X_j\). Los términos residuales \(\varepsilon_i\) representan la parte de \(\mathbf{Y}\) que no es posible predecir a partir de las \(X_j\). Por ello suelen denominarse residuos o errores de predicción.
Bajo esta definición general se engloban los modelos clásicos de:
Regresión: La variable \(\mathbf{Y}\) y todas las \(X_j\) son continuas
Análisis de la varianza: La variable \(\mathbf{Y}\) es continua y todas las \(X_j\) son dicotómicas (toman solo los valores 0 ó 1).
Análisis de la covarianza: La variable \(\mathbf{Y}\) es continua y entre las \(X_j\) hay variables continuas y variables dicotómicas.
Cuando en este modelo la distribución normal se cambia por otra distribución y la relación entre la \(\mathbf{Y}\) y las \(\mathbf{X}\) es lineal a través de una función de enlace \(f\), nos encontramos ante los modelos lineales generalizados:
Regresión logística binaria: La variable \(\mathbf{Y}\) es dicotómica y las \(X_j\) pueden ser dicotómicas y/o continuas. Si \(\mathbf{X}_i\) es la \(i\)-ésima fila de la matriz \(\mathbf{X}\), se asume que la distribución de probabilidad de \(Y_i\) es de Bernoulli, siendo: \[P\left(Y_i=1\right)=f(\mathbf{X_i} \beta)=\frac{e^{\mathbf{X_i}\beta}}{1+e^{\mathbf{X_i}\beta}}\]
Regresión de Poisson: La variable \(\mathbf{Y}\) es una variable de recuento y las \(X_j\) pueden ser dicotómicas y/o continuas. En este modelo se asume que la distribución de \(Y_i\) es de Poisson, \[Y_i \cong \mathbf{P} \left(f(\mathbf{X}_i \beta)\right)\]
Cuando la relación entre la \(\mathbf{Y}\) y las \(\mathbf{X}\) es de la forma: \[Y=f_{1} \left(X_{1}\right)+f_{2}\left(X_{2}\right)+\dots f_{p}\left(X_{p}\right)+\varepsilon\] donde las \(f_i\) son funciones de variación suave (smoothers) nos encontramos ante los modelos aditivos generalizados (GAM).
Estimación del modelo lineal
El modelo lineal para predecir (o explicar) la variable respuesta \(\mathbf{Y}\) en función de la matriz de variables explicativas \(\mathbf{X}\) puede estimarse por dos métodos:
Mínimos cuadrados: Este método puede aplicarse siempre y consiste en encontrar los valores de los parámetros \(\beta_i\) que minimizan la distancia entre los valores observados y los predichos.
Máxima verosimilitud: Este método sólo es de aplicación cuando las observaciones son independientes y las distancias de los puntos al modelo (residuos) siguen una distribución normal. Cuando estas condiciones se cumplen, el método de máxima verosimilitud y el de mínimos cuadrados producen los mismos estimadores.
Sólo en el caso de que los residuos del modelo sean normales e independientes puede llevarse a cabo adecuadamente la inferencia clásica sobre los parámetros del modelo.
Los estimadores de los parámetros, utilizando cualquiera de los dos métodos anteriores, se obtienen mediante:
\[ \hat{\beta} = \left(X^{t}X\right)^{-1}X^{t}Y\]
\[\hat{\sigma}^{2}_\varepsilon=\frac{1}{n-p-1}\sum_{i=1}^{n}\left(y_{i}-\hat{y_{i}}\right)^{2}=\frac{1}{n-p-1}\sum_{i=1}^{n}\left(y_{i}-X_{i}\hat{\beta}\right)^{2}=\frac{1}{n-p-1}\left(\mathbf{Y-X}\hat{\beta}\right)^{t}\left(\mathbf{Y-X}\hat{\beta}\right)\]
La matriz de covarianzas de los estimadores de \(\beta_i\) es:
\[\textrm{Var}\left(\hat{\beta}\right)=\hat{\sigma}^{2}_\varepsilon\left(X^{t}X\right)^{-1}\]
La raiz cuadrada de la diagonal de esta matriz nos proporciona los llamados errores estándar de los \(\beta_i\), que permiten construir contrastes e intervalos de confianza para estos parámetros. Concretamente, si llamamos \(\hat{\sigma}_{\beta_{i}}\) al error estándar del parámetro \(\beta_i\), puede obtenerse un intervalo de confianza a nivel \(1-\alpha\) para dicho parámetro mediante: \[\left[\hat{\beta}_{i}\pm t_{n-p,1-\alpha/2}\hat{\sigma}_{\beta_{i}}\right]\]
Descomposición de la variabilidad.
La variabilidad total presente en la variable respuesta puede descomponerse en la parte explicada por las variables explicativas, más la variabilidad residual:
VARIABILIDAD TOTAL = VARIABILIDAD EXPLICADA + VARIABILIDAD RESIDUAL
donde:
VARIABILIDAD TOTAL: VT = \(\sum_{i=1}^{n}\left(y_{i}-\overline{y}\right)^{2}\)
VARIABILIDAD EXPLICADA: VE = \(\sum_{i=1}^{n}\left(\hat{y}_{i}-\overline{y}\right)^{2}\)
VARIABILIDAD RESIDUAL: VR = \(\sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2}\)
Es decir:
\[\sum_{i=1}^{n}\left(y_{i}-\overline{y}\right)^{2}=\sum_{i=1}^{n}\left(\hat{y}_{i}-\overline{y}\right)^{2}+\sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2}\]
A partir de esta descomposición es posible definir coeficientes para medir la calidad de ajuste del modelo (coeficientes de determinación y de determinación corregido), y construir un contraste (conocido como análisis de la varianza del modelo) que permite decidir si las variables explicativas son realmente explicativas; dicho de otra forma, permite decidir si la variabilidad explicada por dichas variables representa una proporción significativa de la variabilidad total presente en el modelo.
Coeficiente de determinación \(R^2\)
Este coeficiente mide simplemente la proporción de la variabilidad en la variable respuesta que es explicada (o está determinada) por las variables explicativas. Cuanto más se aproxime su valor a 1, mejor es el modelo; cuanto más se aproxime a cero tanto peor:
\[R^2=\frac{VE}{VT}\]
Coeficiente de determinación corregido (o ajustado) \(\tilde{R}^2\)
El coeficiente de determinación que se acaba de definir presenta el problema práctico de que su valor se va aproximando a 1 a medida que se añaden variables explicativas \(X_i\) al modelo, aunque la aportación de cada una sea irrelevante. Interesa por ello redefinir el coeficiente de determinación de forma que se vea “penalizado” si su valor aumenta por introducir muchas variables que no expliquen nada en lugar de pocas variables que expliquen mucho. Para ello observemos en primer lugar que el coeficiente de determinación \(R^2\) puede expresarse como: \[R^{2}=\frac{VE}{VT}=1-\frac{VR}{VT}\]
Para penalizar el exceso de variables explicativas, el coeficiente de determinación corregido se define de manera similar mediante:
\[\tilde{R}^{2}=1-\frac{VR/\left(n-p\right)}{VT/\left(n-1\right)}=1-\frac{n-1}{n-p}\frac{VR}{VT}\]
De esta forma, el término \(\frac{n-1}{n-p}\) es mayor que 1 (tanto mayor cuanto mayor sea el número de variables explicativas \(p\)). Dicho término multiplica a la proporción de varianza no explicada \(\frac{VR}{VT}\), por lo cual su efecto es “inflar” la variabilidad no explicada y por tanto reducir la variabilidad explicada artificialmente por un exceso de variables. Si después de estos ajustes el coeficiente de determinación así obtenido es próximo a 1, ello significaría que aún penalizando la introducción de variables, la varibilidad explicada es alta, y por tanto las variables consideradas puede considerarse que tienen buena capacidad explicativa o predictiva de la variable respuesta.
Análisis de la varianza del modelo lineal
Se denomina así al contraste:
\[H_0:\beta_{1}=\beta_{2}=\dots=\beta_{p}=0\]
que de forma matricial puede expresarse como:
\[H_0:\left(\begin{array}{cccc} 1 & 0 & \ldots & 0\\ 0 & 1 & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & 1 \end{array}\right)\left(\begin{array}{c} \beta_{1}\\ \beta_{2}\\ \vdots\\ \beta_{p} \end{array}\right)=\left(\begin{array}{c} 0\\ 0\\ \vdots\\ 0 \end{array}\right)\]
Este contraste tiene como objetivo decidir si, consideradas conjuntamente, las variables \(X_i\) tienen capacidad explicativa/predictiva sobre la variable respuesta \(\mathbf{Y}\). En caso de aceptar la hipótesis nula, ello significaría que el posible efecto de las \(X_i\) sobre la \(\mathbf{Y}\) puede atribuirse simplemente al azar; en caso de rechazarla, ello implicaría que las \(X_i\) (todas o algunas de estas variables) tienen realmente efecto sobre la \(\mathbf{Y}\), más allá de lo que cabría esperar por azar.
Este es un contraste lineal de la forma:
\[H_0:\mathbf{A\beta}=0\]
Cuando \(H_0\) es cierta el test estadístico:
\[F=\frac{1}{p\hat{\sigma}^{2}}\hat{\beta}^{t}\mathbf{A}^{t}\left(\mathbf{A}\left(\mathbf{X}^{t}\mathbf{X}\right)^{-1}\mathbf{A}^{t}\right)^{-1}\mathbf{A\hat{\beta}}\]
sigue una distribución \(F\) de Fisher con \(p\) y \(n-p-1\) grados de libertad. El p-valor del contraste es
\[p-valor=P\left(F_{p,n-p-1}>F_{obs}\right)=1-P\left(F_{p,n-p-1} \le F_{obs}\right)\]
siendo \(F_{obs}\) el valor observado del test estadístico \(F\). En este caso puede probarse además que:
\[F_{obs}=\frac{VE/p}{VR/(n-p-1)}\]
lo que nos indica que este contraste está, en definitiva, comparando la variabilidad explicada por el modelo con la variabilidad residual (tras ajustar por el número de casos y por el número de variables); sólo si la variabilidad explicada es significativamente mayor que la residual es cuando el valor de \(F_obs\) resultaría lo suficientemente grande para rechazar \(H_0\).
Nótese que este test simplemente decide si hay o no efecto de las \(X_i\) sobre la \(\mathbf{Y}\). No nos dice nada sobre la intensidad de dicho efecto, que deberá cuantificarse en términos del coeficiente de determinación.
Predicción
Una vez ajustado un modelo lineal, puede utilizarse para realizar predicciones sobre los valores de la variable respuesta \(Y\). Hay dos tipos de predicciones que se pueden realizar:
Predicción del valor medio de \(Y\) dado el valor de \(X\).
Predicción del valor individual de \(Y\) dado el valor de \(X\).
En ambos casos la predicción es análoga; si \(X=x_0\) el valor predicho es:
\[\hat{Y} = X\hat{\beta} = x_0\hat{\beta}\]
Aunque el valor predicho sea el mismo, el error estándar difiere entre ambas predicciones. De manera intuitiva, es obvio que el valor medio tiene menor dispersión que los valores individuales y, en efecto, así ocurre. La varianza de la predicción del valor medio dado que \(X=x_0\) es:
\[Var\left(\left.\hat{y}\right|X=x_{0}\right)=\hat{\sigma}^{2}x_{0}\left(X^{t}X\right)^{-1}x_{0}^{t} \]
y la varianza de la predicción del valor individual dado que \(X=x_0\) es:
\[Var\left(\left.\widetilde{y}\right|X=x_{0}\right)=\hat{\sigma}^{2}\left(1+x_{0}\left(X^{t}X\right)^{-1}x_{0}^{t}\right)\]
Los errores estándar son, respectivamente, las raíces cuadradas de estas varianzas. Llamando:
\[ \hat{y}_0 = x_0 \hat{\beta}\]
podemos construir intervalos de confianza para las predicciones mediante:
- Intervalo de confianza a nivel \(1-\alpha\) para la predicción del valor medio de \(Y\) cuando \(X=x_0\):
\[\left[\hat{y}_{0}\pm t_{n-p-1,1-\alpha/2}\hat{\sigma}\sqrt{x_{0}\left(X^{t}X\right)^{-1}x_{0}^{t}}\right]\]
- Intervalo de confianza a nivel \(1-\alpha\) para la predicción del valor individual de \(Y\) cuando \(X=x_0\):
\[\left[\hat{y}_{0}\pm t_{n-p-1,1-\alpha/2}\hat{\sigma}\sqrt{1+x_{0}\left(X^{t}X\right)^{-1}x_{0}^{t}}\right]\]
Ejemplo de estimación con R
Se sabe que la concentración total de clorofila en el fitoplancton marino depende de la composición del fitoplancton. Simularemos la siguiente situación: se recogen muestras de fitoplancton en \(n=86\) localizaciones diferentes a lo largo de una costa. Las muestras de fitoplancton contienen diversos grupos fitoplanctónicos, entre los que se estudian particularmente las Dinofitas y las Cryptofitas. En la muestra \(i\) llamaremos:
\(Y_i\) a la concentración total de clorofila.
\(D_i\) a la concentración de Dinofitas. Supondremos que \(D_i \approx U(120,230)\)
\(C_i\) a la concentración de Cryptofitas; supondremos que \(C_i \approx U(210,330)\)
Supondremos además que las variables \(D_i\) y \(C_i\) son independientes y que sus valores se relacionan con \(Y_i\) de acuerdo a la ecuación:
\[Y_i = \beta_0 + \beta_D\cdot D_i + \beta_C\cdot C_i + \varepsilon_i\] con \(\beta_0=130\), \(\beta_D=0.75\) y \(\beta_C=1.3\). Es decir, la concentración de clorofila se obtiene a partir de las concentraciones de Dinofitas y Criptofitas mediante:
\[Y_i = 130 + 0.75\cdot D_i + 1.3\cdot C_i + \varepsilon_i\] siendo \(\varepsilon_i \approx N(0,95)\).
Simulamos las 86 observaciones y las guardamos en un dataframe. El uso de set.seed(2468) al comienzo de este código garantiza obtener siempre la misma simulación; modificando el número (o suprimiendo el comando set.seed) se generan simulaciones distintas)
Gráficamente:
library(tidyverse)
fito %>%
pivot_longer(-Y,names_to="Variable",values_to="Valor") %>%
ggplot(aes(x=Valor,y=Y)) +
geom_point() +
facet_wrap(~Variable,scales="free_x") +
geom_smooth(method="lm")