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)

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.

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")

Para estimar el vector de parámetros \(\hat\beta\) del modelo Definimos la matriz \(X\) a partir de los datos:

X <- cbind(1,D,C)

Calculamos \((X^tX)^{-1}\). Para multiplicar matrices en R se usa %*%; la traspuesta se calcula mediante t() y la inversa mediante solve(). Nótese que para que exista inversa, el determinante de \(X^tX\) debe ser distinto de cero:

xtx <- t(X)%*%X # Traspuesta de X por X
det(xtx)

[1] 710583801834

Por tanto, en este caso tenemos que la matriz \((X^tX)^{-1}\) es:

xtx1 <- solve(xtx) # Matriz inversa
xtx1

                           D             C
   1.200749791 -1.851980e-03 -3.199341e-03
D -0.001851980  1.026019e-05  1.323483e-07
C -0.003199341  1.323483e-07  1.179752e-05

Calculamos ahora \(X^tY\):

xty <- t(X)%*%Y
xty

         [,1]
     53078.65
D  9468343.51
C 14395833.00

Por último, calculamos el vector de coeficientes estimados \(\hat\beta=(X^tX)^{-1}X^t Y\):

beta <- xtx1%*%xty
beta

         [,1]
  141.8228390
D   0.7516564
C   1.2715591

La estimación del error estándar residual \(\hat\sigma_\varepsilon\) es entonces:

sigma2 <- (1/(n-3))*t(Y-X%*%beta)%*%(Y-X%*%beta)
sigma2

         [,1]
[1,] 727.8514

sigma2 <- as.numeric(sigma2)
sigma <- sqrt(sigma2)

Y el error estándar de los estimadores de \(\beta\):

varB <- sigma2*xtx1
sb <- sqrt(diag(varB))
sb

                      D           C 
29.56293959  0.08641696  0.09266521

La variabilidad residual es \(VR=(Y-XB)^t \cdot (Y-XB)\):

VR <- as.numeric(t(Y-X%*%beta)%*%(Y-X%*%beta))
VR

[1] 60411.66

La variabilidad explicada es \(VE=\left(XB-\overline{y}\right)^t\left(XB-\overline{y}\right)\):

VE <- as.numeric(t(X%*%beta-mean(Y))%*%(X%*%beta-mean(Y)))
VE

[1] 190054.5

y la variabilidad total es la suma de ambas:

VT <- VE+VR
VT

[1] 250466.1

El coeficiente de determinación es el cociente entre la variabilidad explicada y la varibilidad total:

R2 <- VE/VT
R2

[1] 0.7588031

Por último, el valor de \(F_{obs}\):

p <- 2
Fobs <- (VE/p)/(VR/(n-p-1))
Fobs

[1] 130.5586

y su p-valor asociado:

pvalue <- 1-pf(Fobs,p,n-p-1)
pvalue

[1] 0

Utilizando todos estos cálculos matriciales hemos obtenido los mismos resultados que si hubiésemos aplicado lm() (pues esta función contiene una implementación eficiente de todos estos cálculos):

modelo <- lm(Y~D+C,data=fito)
summary(modelo)


Call:
lm(formula = Y ~ D + C, data = fito)

Residuals:
   Min     1Q Median     3Q    Max 
-64.47 -14.37  -0.18  17.66  60.33 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 141.82284   29.56294   4.797 7.00e-06 ***
D             0.75166    0.08642   8.698 2.63e-13 ***
C             1.27156    0.09267  13.722  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 26.98 on 83 degrees of freedom
Multiple R-squared:  0.7588,    Adjusted R-squared:  0.753 
F-statistic: 130.6 on 2 and 83 DF,  p-value: < 2.2e-16

Si repetimos esta simulación muchas veces, generando los datos siempre a partir del mismo modelo (\(Y_i = 130 + 0.75\cdot D_i + 1.3\cdot C_i + \varepsilon\)), obtenemos muchas estimaciones distintas de los parámetros. El siguiente código realiza 10 simulaciones y muestra los coeficientes estimados:

simula <- function(){ 
  # Esta función simula una muestra y estima los parámetros beta y sigmaE
  D <- runif(n,120,230)
  C <- runif(n,210,330)
  sigmaE=25
  eps <- rnorm(n,0,sigmaE)
  Y <- 130+0.75*D+1.3*C+eps
  fito <- data.frame(Y,D,C)
  modelo <- lm(Y~D+C,data=fito)
  estimacion <- c(b=coef(modelo), s=summary(modelo)$sigma)
  estimacion
}
# Se repite la simulación 10 veces:
data.frame(t(replicate(10,simula())))

   b..Intercept.       b.D      b.C        s
1       104.5992 0.7064713 1.423895 24.00497
2       187.8914 0.6141341 1.155648 22.01456
3       145.0013 0.6992625 1.289776 28.13526
4       104.9130 0.8697813 1.309946 23.85009
5       133.8804 0.7522024 1.285287 22.54462
6       102.8798 0.7888521 1.381973 23.77345
7       161.7706 0.6530456 1.237563 27.75090
8       149.0331 0.6315697 1.314691 23.58824
9       116.5137 0.7860573 1.318801 26.88984
10      115.8226 0.8778812 1.274648 25.92411

Cada fila es el resultado de estimar los parámetros a partir de una muestra simulada de 86 observaciones. Estas simulaciones permiten hacernos una idea de como cambian los valores estimados de los coeficientes de una muestra a otra. Como vemos los valores estimados oscilan alrededor de los verdaderos valores de los parámetros. En el modelo que genera los datos los parámetros son \(\beta_0=130\), \(\beta_D=0.75\), \(\beta_C=1.3\) y \(\sigma_\varepsilon=25\). Como vemos en la tabla anterior el coeficiente que se estima con mayor variabilidad es el \(\beta_0\), pues su valor en el modelo es 130 y hemos obtenido valores estimados que van desde 102.87 hasta 187.89; el valor de \(\beta_D\) del modelo es 0.75 y las estimaciones obtenidas oscilan entre 0.614 y 0.878; el parámetro \(\beta_C\) vale 1.3 y las estimaciones oscilan desde 1.139 hasta 1.387. Por último el valor de \(\sigma_\varepsilon\) es 25 y sus estimaciones oscilan desde 22.01 hasta 28.14.

En realidad 10 simulaciones son pocas. Si usamos el mismo código para generar, por ejemplo, 10.000 simulaciones:

sims <- data.frame(t(replicate(10000,simula())))

podemos calcular, para cada parámetro, la media y los percentiles 2.5 y 97.5 de sus estimaciones:

library(flextable)


Attaching package: 'flextable'

The following object is masked from 'package:purrr':

    compose

sims %>% 
  pivot_longer(everything(),names_to = "parametro",values_to = "Valor estimado") %>% 
  group_by(parametro) %>% 
  summarize(Media=mean(`Valor estimado`),
            p2.5=quantile(`Valor estimado`,probs=0.025),
            p97.5=quantile(`Valor estimado`,probs=0.975)) %>% 
  flextable() %>% 
  colformat_double(digits=2) %>% 
  autofit()

parametro	Media	p2.5	p97.5
b..Intercept.	129.49	77.95	181.23
b.C	1.30	1.15	1.46
b.D	0.75	0.58	0.92
s	24.94	21.24	28.89

Como vemos, las medias de las 10.000 estimaciones prácticamente coinciden con los valores de los parámetros. Esto significa que las estimaciones son centradas: a veces dan un valor mayor, a veces un valor menor que el parámetro, pero no tienen desviaciones sistemáticas. Asimismo, vemos que el 95% de los valores estimados de \(\beta_0\) se encuentran entre 77.95 y 181.23; el 95% de los valores estimados de \(\beta_C\) están entre 1.15 y 1.46; el 95% de los valores estimados de \(\beta_D\) están entre 0.58 y 0.92 y los de \(sigma_\varepsilon\) entre 21.24 y 28.89. Podemos expresar este resultado diciendo también que el 95% de los valores estimados de \(\beta_0\) se encuentran a poco más de 51 unidades de distancia de su media; el 95% de los valores estimados de \(\beta_C\) se encuentran a 0.15 unidades de distancia de su media; y el 95% de los valores de \(\beta_D\) se encuentran a 0.17 unidades de distancia de su media. Teniendo en cuenta que las medias de las estimaciones prácticamente coinciden con los parámetros, lo que esta simulación nos enseña es que cuando tomamos una muestra (de tamaño 86) de datos generados de acuerdo al patrón señalado, los valores estimados en general no coinciden con los valores de los parámetros; pero en un 95% de las ocasiones, el valor de \(\hat\beta_0\) estimado no se aparta más de 51 unidades del verdadero valor de \(\beta_0\); el valor de \(\hat\beta_1\) estimado no se aparta más de 0.15 unidades del verdadero valor de \(\beta_1\); y el valor estimado \(\hat\beta_2\) no se aparta más de 0.17 unidades del verdadero valor de \(\beta_2\)

Ahora bien, en este caso hemos podido calcular estos márgenes de error porque jugamos con la ventaja de que conocemos exactamente los valores de los parámetros del modelo que genera los datos, y con ellos hemos podido simular 10.000 muestras de 86 observaciones cada una que nos han permitido calcular estas cantidades. Pero cuando trabajamos con datos reales para construir un modelo, no conocemos los valores de los parámetros, por lo que no tenemos ninguna referencia para simular. Pero realmente no nos hace falta. Más arriba ya hemos calculado los errores estándar \(\hat{\sigma}_{\beta_{i}}\) de los estimadores a partir de nuestra única muestra original (recuérdese, los errores estándar de los estimadores son las raíces cuadradas de la diagonal de la matriz \(\textrm{Var}\left(\hat{\beta}\right)=\hat{\sigma}^{2}_\varepsilon\left(X^{t}X\right)^{-1}\)):

sb

                      D           C 
29.56293959  0.08641696  0.09266521

Bajo las hipótesis del modelo (distribución normal de los residuos), la distribución de las estimaciones de los \(\beta_i\) es normal y el intervalo de confianza correspondiente es de la forma:

\[\left[\hat{\beta}_{i}\pm t_{n-p,1-\alpha/2}\hat{\sigma}_{\beta_{i}}\right]\] En este caso \(t_{86-2,0.975}=1.989\). Si multiplicamos este valor por los errores estándar calculados obtenemos el ancho de cada intervalo (\(t_{n-p,1-\alpha/2}\hat{\sigma}_{\beta_{i}}\)):

sb*1.989

                    D          C 
58.8006868  0.1718833  0.1843111

que como podemos ver coincide bastante aproximadamente con los márgenes de error que nos daba la simulación (usando 10.000 muestras). Esto ocurre precisamente porque los residuos siguen una distribución normal de varianza constante \(\sigma_\varepsilon\), y por tanto los \(\hat\beta_i\) sigue una distribución \(t\) de Student. Por tanto cuando se den estas condiciones, podemos estimar el margen de error de nuestras estimaciones utilizando una única muestra.

La siguiente función permite obtener fácilmente intervalos de confianza al 95% para los valores de los coeficientes del modelo:

confint(modelo)

                 2.5 %      97.5 %
(Intercept) 83.0233412 200.6223369
D            0.5797765   0.9235362
C            1.0872517   1.4558665

y mediante la función predict podemos predecir el valor de la variable respuesta (Y) para los valores D=150, C=250:

predict(modelo,newdata=data.frame(D=150,C=250), interval="prediction")

       fit      lwr      upr
1 572.4611 518.1717 626.7504

mientras que si deseamos obtener el intervalo de confianza para el valor medio de todos los Y que podrían observarse para esos valores de C y D, utilizaríamos:

predict(modelo,newdata=data.frame(D=150,C=250), interval="confidence")

       fit      lwr      upr
1 572.4611 564.2158 580.7063

El problema de la multicolinealidad

Observemos con atención la forma que hemos utilizado para generar nuestros datos originalmente:

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.

Los valores de \(D\) (las concentraciones de dinofitas medidas en cada una de las 86 observaciones que componen la muestra simulada) se han generado aleatoriamente entre 120 y 230 mediante runif(n,120,230). Asimismo los valores de \(C\) (las concentraciones de criptofitas) se han generado aleatoriamente entre 210 y 330 mediante runif(n,210,330). Ambos conjuntos de valores son independientes; no hemos incluido en el programa ninguna instrucción para que los valores de las criptofitas dependan de los valores de la Dinofitas o al revés. De hecho si representamos ambas variables conjuntamente podemos observar que no muestran ningún tipo de asociación entre ellas:

fito %>% 
  ggplot(aes(x=C,y=D)) +
  geom_point()

La estimación del modelo global para predecir la clorofila (\(Y\)) en función de \(D\) y \(C\) fue:

coef(modelo)

(Intercept)           D           C 
141.8228390   0.7516564   1.2715591

Es decir, el estimador \(\hat\beta_D\) vale 0.752 y el estimador de \(\hat\beta_C\) vale 1.272. La interpretación de estos valores es simple: por cada unidad que se incrementa la concentración de Dinofitas, la concentración de clorofila se incrementa en 0.752 unidades; por cada unidad que se incrementa la concentración de Criptofitas, la clorofila se incrementa en 1.272 unidades. Por cierto el estimador \(\hat\beta_0=141.823\) no puede interpretarse en este contexto; su valor correspondería al valor basal de la concentración de algas cuando no hay criptofitas ni dinofitas (esto es, cuando C=D=0); ahora bien, en los datos esta condición no se ha explorado (no hay en la muestra ningún valor 0 en C ni en D), por lo que no puede decirse nada de lo que pasaría en tal caso. El valor de \(\hat\beta_0\) debe considerarse en este caso un simple coeficiente de ajuste sin mayor interpretación.

Si hubiésemos ajustado un modelo lineal para predecir el valor de clorofila en función sólo del valor de Dinofitas obtenemos:

lm(Y~D,data=fito)


Call:
lm(formula = Y ~ D, data = fito)

Coefficients:
(Intercept)            D  
   486.6539       0.7374

El coeficiente de D (el \(\hat\beta_D\)) que en el modelo anterior era 0.752 ahora es 0.737 (hay solo 15 milésimas de diferencia), es decir, no ha habido un gran cambio.

Si ajustamos la clorofila solo en función de las criptofitas obtenemos:

lm(Y~C,data=fito)


Call:
lm(formula = Y ~ C, data = fito)

Coefficients:
(Intercept)            C  
    277.498        1.262

El coeficiente de C (el \(\hat\beta_C\)) que en el modelo anterior era 1.271 es ahora 1.262 (9 milésimas de diferencia). Tampoco en este coeficiente se ha registrado un gran cambio.

Teniendo en cuenta que los valores de \(\hat\beta_D\) y \(\hat\beta_C\) miden el efecto que tiene el cambio de una unidad de Dinofitas y Criptofitas, respectivamente, sobre la concentración de clorofila, el hecho de que estos parámetros apenas cambien si se estiman en modelos por separado o se estiman conjuntamente significa que la presencia de criptofitas no altera el efecto de las dinofitas sobre la clorofila, ni la presencia de dinofitas altera tampoco el efecto de las criptofitas. En realidad era de esperar que ocurriese así porque en la forma que hemos generado los datos no hemos introducido ninguna interacción o dependencia entre ambas variables.

Ahora bien, imaginemos que en la muestra se recoge la concentración de otro componente del plancton (las haptofitas), y se sabe que la concentración de haptofitas se relaciona con la concentración de criptofitas mediante la ecuación:

\[H = 1.5\cdot C + \delta\]

siendo \(\delta \approx N\left(0,15\right)\). En otras palabras la concentración de haptofitas es aproximadamente igual a vez y media la concentración de criptofitas, con una pequeña variabilidad. Podemos añadir esta variable a nuestra simulación (dejando el resto de las cosas igual) del siguiente modo:

set.seed(seed=20087)
n <- 86
D <- runif(n,120,230)
C <- runif(n,210,330)
sigmaE <- 25
eps <- rnorm(n,0,sigmaE)
Y <- 130+0.75*D+1.3*C+eps
H <- 1.5*C+rnorm(n,0,15)
fito <- data.frame(Y,D,C,H)

Si representamos gráficamente las haptofitas frente a las criptofitas y las dinofitas obtenemos:

library(cowplot)
figDH <- fito %>% 
  ggplot(aes(x=D,y=H)) +
  geom_point()

figCH <- fito %>% 
  ggplot(aes(x=C,y=H)) +
  geom_point()

plot_grid(figDH, figCH)

Vemos, pues, que nuestra simulación ha generado los datos conel patrón esperado, en el que las haptofitas no tienen relación con las Dinofitas, pero sí con las criptofitas.

Si queremos predecir la clorofila a partir de estas variables tenemos varias posibilidades:

Predecir la clorofila a partir de solo uno de los componentes del fitoplancton (D, C o H)
Predecir la clorofila a partir de D y C.
Predecir la clorofila a partir de D y H.
Predecir la clorofila a partir de C y H.
Predecir la clorofila a partir de D, C y H.

A continuación ajustaremos todos estos modelos:

Modelos con solo una variable explicativa:

modeloD <- lm(Y~D,data=fito)
summary(modeloD)


Call:
lm(formula = Y ~ D, data = fito)

Residuals:
     Min       1Q   Median       3Q      Max 
-105.325  -35.214   -0.125   37.706  130.811 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 486.6539    27.9839  17.390  < 2e-16 ***
D             0.7374     0.1553   4.748 8.35e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 48.48 on 84 degrees of freedom
Multiple R-squared:  0.2116,    Adjusted R-squared:  0.2022 
F-statistic: 22.55 on 1 and 84 DF,  p-value: 8.346e-06

modeloC <- lm(Y~C,data=fito)
summary(modeloC)


Call:
lm(formula = Y ~ C, data = fito)

Residuals:
   Min     1Q Median     3Q    Max 
-84.22 -24.24   0.43  25.77  74.33 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 277.4980    34.5132   8.040 5.04e-12 ***
C             1.2619     0.1273   9.909 8.83e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 37.08 on 84 degrees of freedom
Multiple R-squared:  0.5389,    Adjusted R-squared:  0.5335 
F-statistic: 98.19 on 1 and 84 DF,  p-value: 8.83e-16

modeloH <- lm(Y~H,data=fito)
summary(modeloH)


Call:
lm(formula = Y ~ H, data = fito)

Residuals:
    Min      1Q  Median      3Q     Max 
-84.114 -26.739  -1.852  27.020  83.852 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 295.36015   31.99475   9.232 2.03e-14 ***
H             0.79764    0.07869  10.136 3.10e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 36.62 on 84 degrees of freedom
Multiple R-squared:  0.5502,    Adjusted R-squared:  0.5448 
F-statistic: 102.7 on 1 and 84 DF,  p-value: 3.099e-16

Sabemos (porque lo hemos hecho así) que la variable \(H\) no intervino en la generación de datos; sin embargo el modelo estimado con esta variable es \(Y=295.36+0.798\cdot H\), siendo el coeficiente de \(H\) significativo, lo que implica que \(H\) tiene cierta capacidad para predecir \(Y\) (el \(R^2\) es 0.55, lo que significa que \(H\) explica el 55% de la variabilidad de la concentración de clorofila \(Y\)). Es lógico que esto suceda pues \(H\) está muy asociada con \(C\) que es la que hemos usado para generar \(Y\).

Modelos con dos variables explicativas

Ya hemos estimado antes el modelo usando las variables \(C\) y \(D\):

modeloDC <- lm(Y~D+C,data=fito)
summary(modeloDC)


Call:
lm(formula = Y ~ D + C, data = fito)

Residuals:
   Min     1Q Median     3Q    Max 
-64.47 -14.37  -0.18  17.66  60.33 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 141.82284   29.56294   4.797 7.00e-06 ***
D             0.75166    0.08642   8.698 2.63e-13 ***
C             1.27156    0.09267  13.722  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 26.98 on 83 degrees of freedom
Multiple R-squared:  0.7588,    Adjusted R-squared:  0.753 
F-statistic: 130.6 on 2 and 83 DF,  p-value: < 2.2e-16

Vemos que con estos datos los coeficientes estimados (\(\hat\beta_D=0.752\) y \(\hat\beta_C=1.272\)) son muy parecidos a los que utilizamos para generar el modelo (\(\beta_D=0.75\) y \(\beta_C=1.3\)), lo que indica que estos datos permiten hacer una muy buena estimación del modelo.

Si estimamos ahora el modelo utilizando como variables explicativas la \(D\) y la \(H\) obtenemos:

modeloDH <- lm(Y~D+H,data=fito)
summary(modeloDH)


Call:
lm(formula = Y ~ D + H, data = fito)

Residuals:
    Min      1Q  Median      3Q     Max 
-61.860 -18.299  -0.807  14.633  68.468 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 162.31295   27.75356   5.848 9.50e-08 ***
D             0.74464    0.08512   8.749 2.08e-13 ***
H             0.80067    0.05710  14.022  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 26.57 on 83 degrees of freedom
Multiple R-squared:  0.766, Adjusted R-squared:  0.7603 
F-statistic: 135.8 on 2 and 83 DF,  p-value: < 2.2e-16

Obsérvese que ahora el efecto estimado de \(D\) ha cambiado un poquito (ahora es \(\hat\beta_D=0.745\)), a la vez que el efecto de \(H\) (su coeficiente estimado es \(\hat\beta_H=0.801\)) resulta significativo. Los coeficientes se diferencian también un poquito más respecto a los obtenidos en las regresiones simples. Nótese además que aunque para generar los datos utilizamos las variables \(D\) y \(C\), el coeficiente de determinación ajustado de la regresión cuando se usan estas dos variables resulta ser \(\bar{R}^2=0.753\), mientras que cuando se usan \(D\) y \(H\) dicho coeficiente es algo mayor(\(\bar{R}^2=0.76\)). Así pues, sorprendentemente, se obtiene un peor ajuste con las variables que en realidad se usaron para generar los datos.

Por último, si estimamos el modelo usando como variables explicativas \(C\) y \(H\) obtenemos:

modeloCH <- lm(Y~C+H,data=fito)
summary(modeloCH)


Call:
lm(formula = Y ~ C + H, data = fito)

Residuals:
    Min      1Q  Median      3Q     Max 
-80.423 -25.227  -3.372  24.650  79.974 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 282.8245    34.2072   8.268  1.9e-12 ***
C             0.4748     0.4597   1.033   0.3046    
H             0.5119     0.2876   1.780   0.0787 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 36.61 on 83 degrees of freedom
Multiple R-squared:  0.5559,    Adjusted R-squared:  0.5452 
F-statistic: 51.95 on 2 and 83 DF,  p-value: 2.346e-15

Ninguna tiene un p-valor menor que el 5%, por lo que ninguna de las dos variables resulta significativa para predecir el valor de \(Y\), aún cuando ambas lo eran consideradas por separado. Es más, los coeficientes estimados (0.475 para el efecto de C y 0.512 pare el efecto de H) son muy distintos de los estimados individualmente para cada variable. Lo que está ocurriendo aquí es que, como sabemos, el valor de \(Y\) realmente depende de \(C\) (la hemos generado así), pero \(H\) mantiene una estrecha relación lineal con \(C\); de alguna manera \(C\) y \(H\) nos están dando la misma información (conociendo \(C\) se puede aproximar muy bien \(H\) y viceversa). Al introducir ambas variables en el mismo modelo para predecir \(Y\) es como si el efecto de una variable se repartiera entre dos y por tanto dicho efecto no es significativo en ninguna de ellas.

Nótese además que cuando se estima el modelo sólo con la \(C\) el error estándar del estimador es 0.127; cuando se estima el modelo solo con la \(H\) el error estándar del estimador es 0.079; pero cuando se estima el modelo simultáneamente con \(C\) y \(H\) los errores estándar son, respectivamente 0.46 y 0.29; en otras palabras, al poner juntas las dos variables sus errores estándar respectivos casi se cuadruplican, lo que significa que los efectos de ambas variables (sus coeficientes de regresión) se estiman peor. La causa de que el error estándar se haya cuadruplicado se puede entender fácilmente si se observa como se calculan los errores estándar de los estimadores. Recordemos que los errores estándar se obtenían como la raíz cuadrada de la diagonal de la matriz:

\[\textrm{Var}\left(\hat{\beta}\right)=\hat{\sigma}^{2}_\varepsilon\left(X^{t}X\right)^{-1}\] En el caso que nos ocupa, la matriz \(X\) tiene tres columnas: una primera columna de unos, una columna con los valores de \(C\) y otra columna con los valores de \(H\). Ahora bien, los valores de \(H\) son casi proporcionales a los de \(C\) (pues \(H=1.5\cdot C + \varepsilon\)), lo que significa que la columna de las \(H\) es casi proporcional a la columna de las \(C\). Ello se traduce en que al calcular \(X^{t}X\) se obtiene una matriz con una columna casi proporcional a otra; su determinante será entonces un valor próximo a cero (recuérdese que si en una matriz dos columnas son linealmente dependientes, su determinante es cero; si la dependencia lineal no es exacta como ocurre en este caso, el determinante puede que no sea cero, pero será un valor pequeño) y al calcular su inversa \(\left(X^{t}X\right)^{-1}\) los valores de esta matriz (y en consecuencia los errores estándar de los estimadores de los \(\beta\)) serán grandes. Este fenómeno se denomina multicolinealidad.

Modelos con tres variables explicativas

Si estimamos el modelo con las tres variables obtenemos:

modeloDCH <- lm(Y~D+C+H,data=fito)
summary(modeloDCH)


Call:
lm(formula = Y ~ D + C + H, data = fito)

Residuals:
    Min      1Q  Median      3Q     Max 
-63.544 -16.071   0.962  15.138  63.984 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 147.46860   28.93318   5.097 2.18e-06 ***
D             0.74776    0.08429   8.872 1.30e-13 ***
C             0.54112    0.33043   1.638   0.1053    
H             0.47506    0.20671   2.298   0.0241 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 26.31 on 82 degrees of freedom
Multiple R-squared:  0.7734,    Adjusted R-squared:  0.7651 
F-statistic: 93.29 on 3 and 82 DF,  p-value: < 2.2e-16

Como vemos, resultan significativas la \(D\) y la \(H\), aún cuando en origen el modelo que con el que generamos los datos dependía de \(D\) y de \(C\). En realidad, la estimación en este caso indica que una vez que se ha tenido en cuenta el efecto de \(D\) y \(H\), la \(C\) no aporta más información para predecir \(Y\). De hecho, si miramos el valor del \(R^2\) ajustado con las tres variables, resulta ser 0.7651, mientras que el ajustado solo con la \(D\) y la \(H\) era 0.7603. Así pues, El modelo con \(D\), \(C\) y \(H\) tiene la misma capacidad predictiva que el que tiene \(D\) y \(H\), que es algo más sencillo, por lo que este segundo sería el modelo preferido.

En la práctica se suele usar el criterio de información de Akaike (AIC) para comparar modelos distintos que pretenden predecir (o explicar) la misma variable respuesta \(Y\). El criterio de Akaike combina la verosimilitud del modelo (la verosimilitud mide cuál es la probabilidad de los datos dado el modelo; el método de máxima verosimilitud encuentra el modelo que maximiza la probabilidad de observar los datos disponibles; en el caso de que los residuos sean normales, este método es equivalente al de mínimos cuadrados) con el número de parámetros que se incluyen en el mismo (número de coeficientes). Teniendo en cuenta que un modelo tiende a ajustar mejor cuanto más parámetros tiene, el criterio de Akaike penaliza el hecho de usar muchos parámetros, de tal manera que un modelo no resulte mejor que otro solo por estar sobreparametrizado. R permite calcular el valor del criterio de Akaike mediante la función AIC. Debe tenerse en cuenta que el AIC depende también del número de datos (a más datos mayor AIC). Por ello, si se va a utilizar este criterio para comparar varios modelos, éstos tienen que haberse ajustado al mismo conjunto de datos; en tal caso resultará preferible aquel modelo con menor valor de AIC. En general se acepta que modelos entre los cuales la diferencia de AIC sea inferior a 3 unidades pueden considerarse equivalentes. Entre los modelos que hemos ajustado a nuestros datos, los valores de AIC son los siguientes:

AIC(modeloC,modeloD,modeloH,modeloDC,modeloDH,modeloCH,modeloDCH)

          df      AIC
modeloC    3 869.4711
modeloD    3 915.6078
modeloH    3 867.3473
modeloDC   4 815.7522
modeloDH   4 813.1515
modeloCH   4 868.2489
modeloDCH  5 812.3839

Vemos que el modelo con menor AIC es el que tiene las tres variables (812.38). No obstante, el modelo solo con D y H es prácticamente equivalente, pues la diferencia en AIC es menor de una unidad. El modelo con las variables D y C (que sabemos fue el que generó los datos) quedaría descartado pues su AIC es más de tres unidades mayor que el del modelo con las tres variables. Vemos, por tanto, que la presencia de colinealidad entre las variables puede conducir a no identificar el modelo correcto para los datos disponibles.

Validación del modelo lineal

Los constrastes de hipótesis y los intervalos de confianza que se construyen en el modelo lineal dependen de las hipótesis sobre el que éste se sustenta, a saber:

La relación entre la variable respuesta \(Y\) y las variables explicativas \(X_i\) es lineal.
Los residuos siguen una distribución \(N\left(0,\sigma_\varepsilon\right)\).
Los residuos son independientes.
Los residuos son homoscedásticos (la varianza \(\sigma_\varepsilon\) es constante).