Modelos Estadísticos. Grado Biotecnología

Modelos de Regresión Múltiple y Polinómicos

En este tema se presentan los modelos de regresión múltiple (RLM) y de regresión polinómica (RP).

Modelos de RLM

Los modelos de regresión lineal múltiple surgen cuando tratamos de explicar el comportamiento de una variable predictora de tipo continuo a través de un conjunto de variables predictoras de tipo continuo mediante una función lineal. Si \(X_1, X_2,...,X_p\) son las variables predictoras el modelo viene dado por:

\[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_p X_p + \epsilon\] Las hipótesis de este modelo es que los errores se distribuyen de forma independiente mediante una distribución Normal de media cero y varianza constante \(\sigma^2\).

Los parámetros desconocidos de este modelo son \((\beta_0, \beta_1,... ,\beta_p, \sigma^2)\) donde:

  • \(\beta_0\) se conoce como interceptación y representa el valor de la respuesta cuando la variable predictora toma el valor cero, interpretándose como un efecto común en la relación entre la predictora y la respuesta.
  • Los \(\beta_i\) son las pendientes de la recta asociadas con cada predictora y representa el aumento o disminución del valor de la respuesta cuando aumentamos en una unidad el valor de la predictora. En este tipo de modelos dicho parámetro se conoce también como el efecto de la predictora sobre la respuesta.
  • \(\sigma^2\) es la varianza residual del modelo.

Dada un muestra de \(n\) sujetos de la variable respuesta \((y_1,...,y_n)\) y de las variables predictoras \((x_{11},...,x_{n1}), (x_{12},...,x_{n2}),...,(x_{1p},...,x_{np})\), el modelo de regresión lineal múltiple se puede escribir como:

\[Y = \left(\begin{array}{c} y_1 \\ y_2 \\ ...\\ y_n\\ \end{array} \right) = \left(\begin{array}{cccc} 1 & x_{11} & ... & x_{1p}\\ 1 & x_{21} & ...& x_{2p}\\ ...& ... & ...& ...\\ 1 & x_{n1}& ...& x_{np}\\ \end{array} \right) \left(\begin{array}{c} \beta_0 \\ \beta_1 \\ ....\\ \beta_p\\ \end{array} \right) + \left(\begin{array}{c} e_1 \\ e_2 \\ ...\\ e_n \end{array} \right) = X \beta + \epsilon\]

donde \(X\) se denomina matriz del diseño, representando el efecto común (columna de 1’s) y el efecto de las predictoras (cada columna con los valores de la variable), y los \(e_i\) representan los errores aleatorias para cada uno de los sujetos de la muestra.

Modelos de RP

Los modelos de regresión lineal múltiple surgen cuando tratamos de explicar el comportamiento de una variable predictora de tipo continuo a través de una variable predictora de tipo continuo mediante una función polinómica lineal. Si \(X\) es la variable predictora y queremos un polinomio de grado \(k\) el modelo viene dado por:

\[Y = \beta_0 + \beta_1 X + \beta_2 X^2 + ... + \beta_k X^k + \epsilon\] Las hipótesis de este modelo es que los errores se distribuyen de forma independiente mediante una distribución Normal de media cero y varianza constante \(\sigma^2\).

Los parámetros desconocidos de este modelo son \((\beta_0, \beta_1,... ,\beta_k, \sigma^2)\) donde:

  • \(\beta_0\) se conoce como interceptación y representa el valor de la respuesta cuando la variable predictora toma el valor cero, interpretándose como un efecto común en la relación entre la predictora y la respuesta.
  • Los \(\beta_i\) son las pendientes de la recta asociadas con cada predictora y representa el aumento o disminución del valor de la respuesta cuando aumentamos en una unidad el valor de la predictora. En este tipo de modelos dicho parámetro se conoce también como el efecto de la potencia de la predictora sobre la respuesta.
  • \(\sigma^2\) es la varianza residual del modelo.

Dada un muestra de \(n\) sujetos de la variable respuesta \((y_1,...,y_n)\) y de la variable predictora \((x_{11},...,x_{n1})\), el modelo de regresión polinómico se puede escribir como:

\[Y = \left(\begin{array}{c} y_1 \\ y_2 \\ ...\\ y_n\\ \end{array} \right) = \left(\begin{array}{cccc} 1 & x_{11} & ... & x^k_{11}\\ 1 & x_{21} & ...& x^k_{21}\\ ...& ... & ...& ...\\ 1 & x_{n1}& ...& x^k_{n1}\\ \end{array} \right) \left(\begin{array}{c} \beta_0 \\ \beta_1 \\ ....\\ \beta_k\\ \end{array} \right) + \left(\begin{array}{c} e_1 \\ e_2 \\ ...\\ e_n \end{array} \right) = X \beta + \epsilon\]

donde \(X\) se denomina matriz del diseño, representando el efecto común (columna de 1’s) y el efecto del grado del polinomio (cada columna con los valores de la variable), y los \(e_i\) representan los errores aleatorias para cada uno de los sujetos de la muestra.

Ejemplos

# Cargamos las librerías
library(tidyverse)
library(forcats)
library(broom)
library(reshape2)
library(lmtest)
library(mgcv)
library(MASS)

Ejemplo 3. Queremos estudiar la relación existente entre la concentración de madera contenida en la pulpa a partir de la que se elabora papel, y la resistencia (en términos de tensión que soporta) del papel resultante. El objetivo del análisis es describir la tendencia observada.

madera <- c(1.0,1.5,2.0,3.0,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0)
resistencia <- c(6.3,11.1,20.0,24.0,26.1,30.0,33.8,34.0,38.1,39.9,42.0,46.1,53.1,52.0,52.5,48,42.8,27.8,21.9)
ejemplo03 <- data.frame(madera,resistencia)

Estos datos ya lo hemos analizado en el tema anterior. Hemos detectado que no se pueden modelizar con un modelo de regresión lineal simple, y que necesitamos un modelo polinómico de grado 2. El modelo (en expresión reducida) vendrá dado por: \[ resistencia \sim madera + madera^2 \]

Ejemplo 6. Para estimar la producción en madera de un bosque se suele realizar un muestreo previo en el que se realizan una serie de medidas no destructivas. Disponemos de mediciones para 20 árboles, así como el volumen de madera que producen una vez cortados. Las variables consideradas son: HT o altura en pies, DBH el diámetro del tronco a 4 píes de altura (en pulgadas), D16 el diámetro del tronco a 16 pies de altura (en pulgadas), y VOL el volumen de madera conseguida (en píes cúbicos). El objetivo del análisis es determinar cuál es la relación entre dichas medidas y el volumen de madera, con el fin de poder predecir este último en función de las primeras.

dbh <- c(10.2,13.72,15.43,14.37,15,15.02,15.12,15.24,15.24,15.28,13.78,15.67,15.67,15.98,16.50,16.87,17.26,17.28,17.87,19.13)
d16 <- c(9.30,12.10,13.30,13.40,14.20,12.80,14.00,13.50,14.00,13.80,13.60,14.00,13.70,13.90,14.90,14.90,14.30,14.30,16.90,17.30)
ht <- c(89.00,90.07,95.08,98.03,99.00,91.05,105.60,100.80,94.00,93.09,89.00,102.00,99.00,89.02,95.09,95.02,91.02,98.06,96.01,101.00)
vol <- c(25.93,45.87,56.20,58.60,63.36,46.35,68.99,62.91,58.13,59.79,6.20,66.16,62.18,57.01,65.62,65.03,66.74,73.38,82.87,95.71)
ejemplo06 <- data.frame(dbh,d16,ht,vol)

El objetivo es tratar de explicar el volumen de madera en función del conjunto propuesto de variables predictoras. El modelo con todas las variables predictoras se puede escribir: \[ vol \sim dbh + d16 + ht\]

Representamos la asociación de cada predictora con la respuesta:

# Represntación gráfica de todas las relaciones lineales con respecto a la respuesta
library(reshape2)
datacomp = melt(ejemplo06, id.vars='vol')
ggplot(datacomp) +
  geom_jitter(aes(value,vol, colour=variable),) + 
  geom_smooth(aes(value,vol, colour=variable), method=lm, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Volumen") +
  theme_bw()

En el gráfico se aprecia una tendencia con mucha pendiente para las variables predictoras dbh y d16. Existe una pendiente inferior con respecto a la variable ht. Esto parece indicar que las dos primeras están más relacionadas con la respuesta que la última.

Ejemplo 7. Se ha llevado a cabo un experimento para estudiar la concentración presente de un fármaco en el hígado después de sufrir un tratamiento. Se piensa que las variables que pueden influir en la concentración son el peso del cuerpo, el peso del hígado y la dosis de fármaco administrada.

p.cuerpo <- c(176,176,190,176,200,167,188,195,176,165,158,148,149,163,170,186,146,181,149)
p.higado <- c(6.5,9.5,9.0,8.9,7.2,8.9,8.0,10.0,8.0,7.9,6.9,7.3,5.2,8.4,7.2,6.8,7.3, 9.0,6.4) 
dosis <- c(.88,.88, 1.0,.88,1.0,.83,.94,.98,.88,.84,.80,.74,.75,.81,.85,.94,.73,.90,.75)
concen <- c(.42,.25,.56,.23,.23,.32,.37,.41,.33,.38,.27,.36,.21,.28,.34,.28,.30,.37,.46)
ejemplo07 <- data.frame(p.cuerpo,p.higado,dosis,concen)

El objetivo es tratar de explicar la concentración del fármaco en función del conjunto propuesto de variables predictoras. El modelo con todas las variables predictoras se puede escribir: \[ concen \sim p.cuerpo + p.higado + dosis\]

Representamos la asociación de cada predictora con la respuesta:

# Represntación gráfica de todas las relaciones lineales con respecto a la respuesta
library(reshape2)
datacomp = melt(ejemplo07, id.vars='concen')
ggplot(datacomp) +
  geom_jitter(aes(value,concen, colour=variable),) + 
  geom_smooth(aes(value,concen, colour=variable), method=lm, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Concentración") +
  theme_bw()

Las tres variables predictoras muestran un grado de relación similar con la respuesta (pendientes de las rectas similares). En todos ellos se aprecia una observación un poco más alejada del resto (concentración > 0.6) que podría ser influyente en la obtención del modelo correspondiente.

Un aspecto ha tener muy en cuenta es que los gráficos realizados sólo muestran la relacíón entre cada predictora y la respuesta, cuando nuestro modelo contiene a todas las predictoras. En la práctica resulta imposible hacer un gráfico que represente todas las posibles asociaciones.

Modelo saturado y anidado

En modelos donde hay más de un efecto sobre la predictora, es decir, tenemos diferentes predictoras o un modelo polinómico debemos introducir dos conceptos que resultan muy relevantes, y que utilizaremos de forma muy habitual en la selección del mejor modelo.

El modelo saturado es aquel que contiene todos los efectos asociados con las diferentes predictoras consideradas. Para los tres ejemplos considerados tendríamos:

\[\left\{ \begin{array}{lc} \text{Ejemplo 4} & resistencia \sim madera + madera^2\\ \text{Ejemplo 6} & vol \sim dbh + d16 + ht\\ \text{Ejemplo 7} & concen \sim p.cuerpo + p.higado + dosis\\ \end{array} \right. \]

Los modelos anidados son todos los modelos que podemos considerar y que no contienen todos los efectos asociados con las predictoras. Si tenemos un modelo con dos predictoras \(X_1\), y \(X_2\) lo modelos anidados serán:

\[\left\{ \begin{array}{lc} \text{con } X_1 & Y \sim X_1\\ \text{con } X_2 & Y \sim X_2\\ \text{Sin ninguna} & Y \sim 1\\ \end{array} \right. \]

Todos ellos están “anidados” dentro del modelo saturado y reflejan diferente información. El primero refleja que la respuesta sólo está relacionada con \(X_1\), el segundo que la respuesta está relacionada con \(X_2\), y el último refleja que no hay ninguna predictora relacionada con la respuesta.

Debemos tener en cuenta que al incluir más de una predictora debemos decidir si todas ellas son relevantes para explicar el comportamiento de la respuesta, o bien si podemos prescindir de algunas de ellas.

La consideración de los diferentes modelos anidados varía en función del modelo con el que trabajemos. En el caso de los de RLM el orden de los modelos anidados no es relevante, pero sin embargo si lo es los modelos polinómicos. No tiene sentido considerar un modelo en el que sólo se incluya el efecto del polinomio de grado 2 pero que no se incluya el de grado 1. Pos su propia construcción cuando consideramos un modelo polinómico de grado \(k\) se deben considerar obligatoriamente todos los grados desde \(1\) hasta \(k-1\). Si consideramos un modelo polinómico de grado 4, el orden de los modelos anidados viene dado por:

\[\left\{ \begin{array}{ll} \text{saturado } & Y \sim X + X^2 + X^3 + X^4\\ \text{grado 3 } & Y \sim X + X^2 + X^3 \\ \text{grado 2 } & Y \sim X + X^2 \\ \text{grado 1 } & Y \sim X \\ \text{sin efectos } & Y \sim 1 \\ \end{array} \right. \]

Siempre trataremos de seleccionar le modelo más pequeño, es decir, con menos efectos pero con igual predictivo que el modelo saturado.

Estimación del modelo

En esta fase procedemos con la estimación del modelo, es decir, obtener estimaciones puntuales,intervalos de confianza, y contrastes de hipótesis sobre cada uno de los parámetros del modelo. Todo el desarrollo teórico se puede consultar en los libros recomendados en la bibliografía. Aquí nos ocupamos únicamente de como obtener y analizar dicho proceso de estimación.

Ejemplo 4

Para ajustar modelos polinómicos utilizaremos la función I() para identificar el cálculo del polinomio.

# Ajuste del modelo
ajuste <- lm(resistencia ~ madera + I(madera^2),data = ejemplo03)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(ajuste)
fitcoef
# Intervalo de confianza
confint(ajuste)
##                   2.5 %     97.5 %
## (Intercept) -13.8812496  0.5328664
## madera        9.6382023 13.8898090
## I(madera^2)  -0.7655346 -0.5035638

El modelo ajustado viene dado por: \[\widehat{resistencia} = -6.674 + 11.76 * madera - 0.63*madera^2\]

En este caso no se pueden interpretar las significatividades de los coeficientes por debajo del grado 2, ya que dichos coeficientes cambiarían si finalmente nos quedáramos con un modelo anidado de grado inferior. En el punto de comparación de modelos veremos como comparar y elegir entre el modelo de grado 2 y el modelo de grado 1. En el gráfico siguiente se comparan los modelos ajustados para el grado 1 (azul) que vimos en el tema anterior, y el grado 2 (rojo) que acabamos de ajustar. Se puede ver que el modelo de grado 2 recoge mejor el comportamiento de la respuesta con respecto a la predictora.

Ejemplo 6

El ajuste para el modelo propuesto para estos datos viene dado por:

# Ajuste del modelo
ajuste <- lm(vol ~ dbh + d16 + ht, data = ejemplo06)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(ajuste)
fitcoef
# Intervalo de confianza
confint(ajuste)
##                    2.5 %      97.5 %
## (Intercept) -282.1716025 -107.287034
## dbh            3.4599983   16.147066
## d16          -11.3069840    3.560317
## ht             0.6212243    2.637807

El modelo obtenido es: \[\widehat{vol} = -194.73 + 9.80 * dbf - 3.87 * d16 + 1.62 *ht\] Ente tipo de modelos la magnitud del coeficiente no implica directamente que la predictora sea más relevante,salvo que todas las predictoras estén medidas en la misma escala. Sin embargo, si todas las predictoras son positivas, si podemos interpretar los signos de cada coeficiente para identificar que predictoras influyen en el crecimiento o decrecimiento de la respuesta. Se puede ver que el incremento de dbf o ht aumenta el volumen, mientras que el aumento de d16 disminuye el volumen final. Algo similar ocurre con los pvalores. Estos pvalores no sirven para estudiar si la predictora influye en la respuesta de forma individual, es decir, si el coeficiente asociado es distinto de cero, pero como veremos más adelante nuestra propuesta de modelo no puede depender de estudiar esos pvalores individuales. A la vista de dichos pvalores podríamos creer que nuestro mejor modelo vendría dado por: \[vol \sim dbf +d16,\]

Esta forma de proceder no es adecuada porque no tiene en cuenta la posible relación que existe entre las predictoras. Esta asociación entre las predictoras, y no de las predictoras con la respuesta, provoca en ocasiones que el pvalor individual de una de ellas resulte no significativo cunado a lo mejor es la variable con mayor capacidad explicativa sobre la respuesta. Se trataría entonces de construir todos los posibles modelos que surgen al considerar todas las combinaciones de predictoras y elegir aquél con mayor capacidad predictiva. Ese proceso puede resultar muy complejo si el número de predictoras es muy amplio. Más adelante veremos los procedimientos se selección automática de variables que nos permiten seleccionar el conjunto de predictores que mejor explican el comportamiento de la respuesta.

Ejemplo 7

El ajuste para el modelo propuesto para estos datos viene dado por:

# Ajuste del modelo
ajuste <- lm(concen ~ p.cuerpo + p.higado + dosis, data = ejemplo07)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(ajuste)
fitcoef
# Intervalo de confianza
confint(ajuste)
##                   2.5 %       97.5 %
## (Intercept) -0.14882688  0.680670424
## p.cuerpo    -0.03824316 -0.004249511
## p.higado    -0.02239941  0.050995526
## dosis        0.93271260  7.423510209

El modelo obtenido es: \[\widehat{concen} = 0.27 - 0.02 * p.cuerpo + 0.01 * p.higado + 4.18 *dosis\]

¿Cómo interpretamos los coeficientes y los pvalores obtenidos para este modelo?

Bondad de ajuste del modelo

La bondad de ajuste para este tipo de modelos utiliza los mismos procedimientos que en el caso de los modelos de regresión lineal simple: el \(R^2\) ajustado y el test \(F\) de regresión.

Analizamos los resultados obtenidos para cada uno de los ejemplos.

Ejemplo 4

En este caso comparamos los resultados con los obtenidos en el tema anterior cuando ajustábamos un modelo de regresión lineal simple.

Se puede ver como el \(R^2\) pasa de 0.26 a 0.89 indicando una clara mejoría en la bondad del ajuste. El pvalor también resulta significativo. Dado que el modelo de grado 1 no cumplía con las hipótesis del modelo, la propuesta de modelo de grado 2 parece mejorar hasta el momento los resultados del modelo inicial.

Ejemplo 6

# Ajuste del modelo
ajuste <- lm(vol ~ dbh + d16 + ht, data = ejemplo06)
gof <- glance(ajuste)
gof

El \(R^2\) ajustado es de 0.76 mostrando una gran capacidad explicativa. Además, dado que el pvalor del test \(F\) resulta significativo hay evidencias de que al menos existe un efecto asociado con una predictora que ayuda a explicar el comportamiento de la respuesta. El problema de este test es que resulta imposible identificar cual de esos efectos es significativo. Atendiendo a los pvalores individuales obtenidos en el punto anterior podemos ajustar un nuevo modelo (sin la predictora d16) y comparar la capacidad explicativa.

# Ajuste del modelo
ajuste <- lm(vol ~ dbh + ht, data = ejemplo06)
gof <- glance(ajuste)
gof

El \(R^2\) ajustado es prácticamente el mismo y el pvalor resulta significativo. Tenemos un modelo con un efecto menos, y por tanto más simple, pero con la misma capacidad explicativa. En el punto siguiente estudiaremos con detalle la selección del mejor modelo, pero esto ya nos da una indicación del modelo que deberíamos considerar para explicar el comportamiento del volumen.

Ejemplo 7

El ajuste para el modelo propuesto para estos datos viene dado por:

# Ajuste del modelo
ajuste <- lm(concen ~ p.cuerpo + p.higado + dosis, data = ejemplo07)
gof <- glance(ajuste)
gof

¿qué conclusiones podemos extraer de este análisis? ¿Y si ajustamos el modelo sin el peso del hígado?

# Ajuste del modelo
ajuste <- lm(concen ~ p.cuerpo + dosis, data = ejemplo07)
gof <- glance(ajuste)
gof

Comparación de modelos

En el caso de tener más de una variable predictora se hace necesario introducir nuevas herramientas para la selección del mejor modelo. Algunas de ellas están basadas en los estadísticos que presentamos en el tema de regresión lineal simple (\(AIC\) y \(BIC\)), pero otras son necesarias para valorar los efectos presentes en un modelo.

Tabla Anova

La tabla ANOVA es un herramienta que se basa en la variabilidad explicada por cada predictora. Recordemos que la variabilidad total se ha descompuesto entre variabilidad del modelo y variabilidad residual. Sin embargo, la variabilidad del modelo se puede descomponer en términos de la variabilidad explicada por cada predictora, y comparándola con la residual para saber que variable posee una mayor capacidad explicativa. Se compara mediante un test \(F\) la variabilidad aportada por cada predictora con respecto a la residual, obtenido un test que nos permite identificar los efectos significativos, es decir, aquellos que son capaces de explicar el comportamiento de la respuesta. Para realizar este análisis utilizamos la función anova(ajuste).

Selección de variables o efectos

Los métodos de selección de variables o de efectos se utilizan para seleccionar el mejor modelo de entre todos los que se pueden construir con las combinaciones de todas las variables predictoras. Dichos procedimientos se basan en comparar el valor del \(AIC\) de todos los modelos y quedarse con el que tiene un valor más pequeño. Aunque en la literatura existen diferentes formas de proceder, aquí se presenta únicamente el método secuencial hacia delante y atrás que toma como modelo inicial el saturado, es decir, aquel que contiene todos los efectos o variables predictoras presentes. Para realizar este análisis utilizamos la función step(ajuste, direction = "both"). LA forma de proceder es:

  1. Establecemos el modelo saturado y calculamos su AIC
  2. (Selección inicial) Calculamos todos los modelos quitando un efecto y valoremos su AIC. Si el AIC de alguno de estos modelos es superior al del paso 1 deberíamos eliminar dicho efecto, mientras que sin todos son mayores que el del paso 1, el modelo resultante es el modelo saturado.
  3. (hacia atrás) Si hemos eliminado un efecto recalculamos todos los modelos posibles y valoramos de nuevo el AIC de todos ellos. Procedemos como en el paso 2 comparando con el modelo saturado obtenido en el paso anterior.
  4. (hacia delante) Si hemos eliminado un nuevo efecto comprobamos si podemos volver a incluir el efecto eliminado en el paso 1 (es decir comparamos su valor de AIC). Si al incluir dicho efecto el AIC disminuye volveríamos a incluirlo.
  5. Repetimos iterativamente los pasos 3 y 4 hasta alcanzar un único modelo.

Análisis de los ejemplos

Pasamos al estudio de cada uno de los ejemplos

Ejemplo 4

La tabla ANOVA para este modelo se obtiene mediante

Al tratarse de un modelo polinómico, hemos de tener en cuenta que nos se puede estudiar el efecto asociado al grado 1 hasta descartar que el grado 2 no es necesario. En este caso el pvalor asociado con el efecto de grado 2 indica que es significativo y que debe estar presente en el modelo. El modelo final viene dado por \[ resistencia \sim madera + madera^2\]

El procedimiento de selección de variables viene dado por:

## Start:  AIC=59.21
## resistencia ~ madera + I(madera^2)
## 
##               Df Sum of Sq     RSS     AIC
## <none>                      312.64  59.212
## - I(madera^2)  1    2060.8 2373.46  95.726
## - madera       1    2689.2 3001.82 100.188
## 
## Call:
## lm(formula = resistencia ~ madera + I(madera^2), data = ejemplo03)
## 
## Coefficients:
## (Intercept)       madera  I(madera^2)  
##     -6.6742      11.7640      -0.6345

La primera tabla nos proporciona el valor del AIC par el modelo saturado (fila ), la segunda el valor del AIC si elimináramos el grado 2, y la tercera el valor del AIC si elimináramos el grado 1. Esta última no tiene sentido mirarla porque estamos en un modelo polinómico y sólo se puede interpretar si eliminamos el efecto asociado con el grado 2. En cuanto al grado 2 podemos ver que si lo eliminásemos del modelo el AIC pasaría de 59.212 (modelo saturado) a 95.726, lo que implica que dicho efecto es imprescindible para mejorar la capacidad predictiva del modelo. En este caso sólo se realiza una iteración porque no podemos eliminar ninguna variable. Al final aparecen los coeficientes del modelo ajustado: \[ resistencia = -6.6742 + 11.7640 * madera - -0.6345 * madera^2\]

Ejemplo 6

La tabla ANOVA para este modelo se obtiene mediante

# Ajuste del modelo
ajuste <- lm(vol ~ dbh + d16 + ht, data = ejemplo06)
anova(ajuste)

La tabla ANOVA resultante muestra que el efecto asociado con d16 no resulta significativo, es decir, dicha variable podría ser eliminada del modelo. El modelo resultante contendría únicamente a las predictoras dbh y ht.

Veamos que resultados proporciona el procedimiento de selección de variables:

## Start:  AIC=91.55
## vol ~ dbh + d16 + ht
## 
##        Df Sum of Sq    RSS     AIC
## - d16   1     99.45 1403.6  91.021
## <none>              1304.1  91.551
## - dbh   1    874.86 2179.0  99.818
## - ht    1    956.71 2260.8 100.555
## 
## Step:  AIC=91.02
## vol ~ dbh + ht
## 
##        Df Sum of Sq    RSS     AIC
## <none>              1403.6  91.021
## + d16   1     99.45 1304.1  91.551
## - ht    1    864.90 2268.5  98.623
## - dbh   1   2545.28 3948.9 109.709
## 
## Call:
## lm(formula = vol ~ dbh + ht, data = ejemplo06)
## 
## Coefficients:
## (Intercept)          dbh           ht  
##    -189.857        6.782        1.506

El proceso de selección realiza dos iteracciones. En la primera se elimina el efecto asociado con d16 (AIC = 91.021 más pequeño que el del modelo sin eliminar ninguna variable 91.551). En el segundo paso valoramos si podemos incluir el efecto eliminado en el paso anterior. Como era de esperar el AIC al incluir dicho efecto es superior al de no incluirlo, y por tanto no concluimos el efecto asociado con dicha variable. El modelo resultante viene dado por \[vol = -189.857 + 6.782 * dbh + 1.506 * ht\] ¿Cómo interpretamos los coeficientes del modelo?

Ejemplo 7

La tabla ANOVA para este modelo se obtiene mediante

# Ajuste del modelo
ajuste <- lm(concen ~ p.cuerpo + p.higado + dosis, data = ejemplo07)
anova(ajuste)

Podemos ver que el único efecto significativo es el asociado con la variable dosis. Los otros dos efectos muestran resultados no significativos, de forma que el modelo final contendría únicamente la variable dosis. Veamos el procedimiento de selección de variables.

## Start:  AIC=-93.78
## concen ~ p.cuerpo + p.higado + dosis
## 
##            Df Sum of Sq      RSS     AIC
## - p.higado  1  0.004120 0.093729 -94.924
## <none>                  0.089609 -93.778
## - p.cuerpo  1  0.042408 0.132017 -88.416
## - dosis     1  0.044982 0.134591 -88.049
## 
## Step:  AIC=-94.92
## concen ~ p.cuerpo + dosis
## 
##            Df Sum of Sq      RSS     AIC
## <none>                  0.093729 -94.924
## + p.higado  1  0.004120 0.089609 -93.778
## - p.cuerpo  1  0.039851 0.133580 -90.192
## - dosis     1  0.043929 0.137658 -89.621
## 
## Call:
## lm(formula = concen ~ p.cuerpo + dosis, data = ejemplo07)
## 
## Coefficients:
## (Intercept)     p.cuerpo        dosis  
##     0.28552     -0.02044      4.12533

El proceso muestra que el modelo resultante debe contener las variables p.cuerpo y dosis. En el primer paso podemos ver como eliminar la variable p.higado disminuye el valor de AIC, mientras que al intentar incluirlo en el modelo sin él, dicho valor aumenta. El modelo obtenido viene dado por: \[concen = 0.28552 -0.02044*p.cuerpo + 4.12533 *dosis\] ¿Cómo interpretamos los coeficientes del modelo?

Diagnóstico del modelo

El diagnóstico en este tipo de modelos utiliza los mismo procedimientos que los modelos de regresión lineal simple. Hay que tener en cuenta que los procedimientos de diagnóstico se centran en el estudio de los residuos del modelo que únicamente dependen del valor de la respuesta y los valores ajustados del modelo. Pasamos al análisis de los ejemplos directamente. Para el cálculo de los residuos utilizaremos los modelos seleccionados en el punto anterior.

Ejemplo 3

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(resistencia ~ madera + I(madera^2),data = ejemplo03)
fitinf <- fortify(ajuste)

Comenzamos con el análisis gráfico de los residuos

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","madera", "I(madera^2)",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method= mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

No se aprecian problemas respecto a la linealidad y varianza constante, aunque si se ve cierto agrupamiento de las observaciones.

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) + 
  stat_qq() + 
  geom_abline() +
  theme_bw()

El gráfico de normalidad parece indicar que no hay problemas con la normalidad de los residuos. Pasamos ahora a la verificación de las hipótesis a través del los diferentes tests estadísticos asociados con las hipótesis del modelo.

# Varianza constante
bptest(ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 5.6722, df = 2, p-value = 0.05865
# Normalidad
shapiro.test(fitinf$.stdresid)
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid
## W = 0.92299, p-value = 0.1286
# Independencia
dwtest(ajuste)
## 
##  Durbin-Watson test
## 
## data:  ajuste
## DW = 0.69747, p-value = 2.382e-05
## alternative hypothesis: true autocorrelation is greater than 0

Tanto el test de homogeneidad como el de normalidad proporcionan resultados no significativos indicando el cumplimiento de las hipótesis. Por otro lado, el test de correlación de los residuos resulta significativo. Sin embargo, en los modelos polinómicos es habitual encontrarnos con este problema debido a la propia estructura del modelo (curvas en lugar de rectas). Hay dependencia entre las observaciones situadas en los extremos del modelo ajustado, debido a que valores muy distintos de la predictora proporcionan la misma respuesta.

Por último realizamos un análisis de influencia para conocer si hay alguna observación que se pueda considerar anómala:

# Calculamos los valores absolutos
abs(fitinf$.cooksd)
##  [1] 0.045065495 0.018793975 0.159479502 0.002960342 0.030045238
##  [6] 0.019860882 0.010837085 0.044040594 0.018027389 0.021084854
## [11] 0.015798264 0.001346582 0.070761814 0.047760688 0.098745245
## [16] 0.063130767 0.057866905 0.285372258 0.640308010
# Valoramos si hay alguna observación con distancia mayor que 1
abs(fitinf$.cooksd) > 1
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

No se registran observaciones influyentes.

Ejemplo 6

Consideramos el modelo tras la selección del mejor modelo.

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(vol ~ dbh + ht, data = ejemplo06)
fitinf <- fortify(ajuste)

Comenzamos con el análisis gráfico de los residuos

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","dbh", "ht",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method= mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

No se observan problemas con las hipótesis del modelo pero si se aprecian dos observaciones (una por arriba y otra por abajo) con residuos muy altos, lo que podría indicar observaciones anómalas. En concreto la de abajo tiene un residuo próximo a cuatro. La forma de proceder sería eliminar dichas observaciones y volver ajustar el modelo.

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) + 
  stat_qq() + 
  geom_abline() +
  theme_bw()

El análisis de normalidad presenta el mismo problema con las observaciones alejadas.

Seleccionamos las observaciones alejadas, las eliminamos y ajustamos de nuevo el modelo.

# Localización de las observaciones anómalas
sel <- abs(fitinf$.stdresid) < 1.9
# Seleccionamos el nuevo banco de datos
ejemplo06new <- ejemplo06[sel,]
# Obtenemos el ajuste de nuevo con todas las predictoras
ajuste <- lm(vol ~ dbh + d16 + ht, data = ejemplo06new)
# Calidad del ajuste
glance(ajuste)
# Selección del mejor modelo
step(ajuste)
## Start:  AIC=37.01
## vol ~ dbh + d16 + ht
## 
##        Df Sum of Sq     RSS    AIC
## <none>               90.191 37.008
## - dbh   1    70.506 160.697 45.405
## - d16   1    91.504 181.695 47.615
## - ht    1   219.553 309.744 57.217
## 
## Call:
## lm(formula = vol ~ dbh + d16 + ht, data = ejemplo06new)
## 
## Coefficients:
## (Intercept)          dbh          d16           ht  
##   -138.3318       3.5098       4.4057       0.8737
# Cargamos resultados
fitinf <- fortify(ajuste)

El modelo resultante no prescinde de ninguna variable predictora y tiene buena capacidad explicativa con test \(F\) de la regresión significativo. El modelo obtenido es: \[vol = -138.3318 + 3.5098 * dbh + 4.4057 * d16 + 0.8737 * ht\]

Realizamos el análisis de los residuos.

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","dbh", "d16", "ht",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method= mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

En este caso los gráficos no presentan ningún problema con las hipótesis. Veamos el gráfico de normalidad:

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) + 
  stat_qq() + 
  geom_abline() +
  theme_bw()

No se detectan problemas con la normalidad.

En cuanto a los tests:

# Varianza constante
bptest(ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 3.55, df = 3, p-value = 0.3143
# Normalidad
shapiro.test(fitinf$.stdresid)
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid
## W = 0.96895, p-value = 0.7778
# Independencia
dwtest(ajuste)
## 
##  Durbin-Watson test
## 
## data:  ajuste
## DW = 1.8533, p-value = 0.2844
## alternative hypothesis: true autocorrelation is greater than 0

Se verifican todas las hipótesis del modelo. EL modelo con las tres predictoras eliminando las dos observaciones alejadas es el más adecuado y el que deberemos utilizar para la predicción.

Ejemplo 7

Consideramos el modelo tras la selección del mejor modelo.

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(concen ~ p.cuerpo + dosis, data = ejemplo07)
fitinf <- fortify(ajuste)

Realizamos el análisis de los residuos.

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","p.cuerpo", "dosis", ".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method= mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

No se aprecian problemas con las hipótesis pero si una valor ajustado muy grande (valor de 0.55 aproximadamente) con respecto al resto. Habrá que tenerlo en cuenta por si se trata de un valor influyente. NO se trata de una valor anómalo porque no tiene un residuo grande.

Veamos el estudio de normalidad:

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) + 
  stat_qq() + 
  geom_abline() +
  theme_bw()

El gráfico de normalidad no parece detectar problemas con esta hipótesis. Realizamos los tests de hipótesis y el estudio de influencia.

# Varianza constante
bptest(ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 0.40345, df = 2, p-value = 0.8173
# Normalidad
shapiro.test(fitinf$.stdresid)
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid
## W = 0.9549, p-value = 0.4766
# Independencia
dwtest(ajuste)
## 
##  Durbin-Watson test
## 
## data:  ajuste
## DW = 1.7624, p-value = 0.2748
## alternative hypothesis: true autocorrelation is greater than 0

Los tests resultan no significativos lo que nos permite concluir que se verifican las hipótesis del modelo.

# Calculamos los valores absolutos
abs(fitinf$.cooksd)
##  [1] 4.591374e-02 2.003294e-02 1.820157e+00 3.363485e-02 2.116883e-01
##  [6] 3.592818e-03 2.387308e-02 6.482268e-02 6.709915e-04 5.105119e-05
## [11] 5.941172e-02 3.062382e-02 2.002660e-01 1.382759e-03 2.131548e-03
## [16] 4.558724e-02 2.356329e-03 4.443919e-02 2.112934e-01
# Valoramos si hay alguna observación con distancia mayor que 1
abs(fitinf$.cooksd) > 1
##  [1] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

El estudio de influencia detecta la observación situada en la posición 3 como influyente. Veamos sus resultados

fitinf

Dicha observación coincide con la detectada en el gráficos de residuos vs ajustados, como se puede comprobar en la columna .fitted, con una valor ajustado de 0.53. Eliminamos dicha observación y realizamos el análisis completo de nuevo. Procedemos de forma similar al ejemplo anterior.

Seleccionamos la observación influyente y la eliminamos de los datos.

ejemplo07new <- ejemplo07[-3,]
# Obtenemos el ajuste de nuevo con todas las predictoras
ajuste <- lm(concen ~ p.cuerpo + p.higado + dosis, data = ejemplo07new)
# Calidad del ajuste
glance(ajuste)
# Selección del mejor modelo
step(ajuste)
## Start:  AIC=-88.25
## concen ~ p.cuerpo + p.higado + dosis
## 
##            Df  Sum of Sq      RSS     AIC
## - dosis     1 0.00097916 0.086696 -90.043
## - p.cuerpo  1 0.00105871 0.086776 -90.026
## - p.higado  1 0.00142114 0.087138 -89.951
## <none>                   0.085717 -88.247
## 
## Step:  AIC=-90.04
## concen ~ p.cuerpo + p.higado
## 
##            Df  Sum of Sq      RSS     AIC
## - p.cuerpo  1 0.00035562 0.087052 -91.969
## - p.higado  1 0.00082681 0.087523 -91.872
## <none>                   0.086696 -90.043
## 
## Step:  AIC=-91.97
## concen ~ p.higado
## 
##            Df  Sum of Sq      RSS     AIC
## - p.higado  1 0.00050917 0.087561 -93.864
## <none>                   0.087052 -91.969
## 
## Step:  AIC=-93.86
## concen ~ 1
## 
## Call:
## lm(formula = concen ~ 1, data = ejemplo07new)
## 
## Coefficients:
## (Intercept)  
##      0.3228
# Cargamos resultados
fitinf <- fortify(ajuste)

Al eliminar la observación influyente el modelo obtenido no contiene ninguna variable predictora, es decir, no es posible explicar el comportamiento de la concentración en función de las variables consideradas. En esta situación habría que revisar los datos experimentales para verificar la observación influyente, con le objetivo de determinar si es incorrecta o si de verdad hay que considerarla en el modelo. Sin variables es imposible realizar la predicción del modelo ya que está sería siempre constante.

Predicción

El proceso de predicción en este tipo de modelos es similar al de los modelos de regresión lineal simple. Necesitamos un grid de valores para todas las variables predictoras del modelo y calcular el valor de predicción correspondiente. Lo que resulta más difícil es representar las bandas de confianza al disponer únicamente de dos dimensiones. En los modelos polinómicos si podemos construir dichas bandas ya que en realidad tenemos una única variable predictora, pero no es así en los modelos de regresión lineal múltiple con más de una predictora.

Ejemplo 3

Queremos predecir la resistencia en función de la concentración de madera contenida en la pulpa. Para ello construimos una secuencia de valores dentro del rango de valores del contenido de hierro.

Predecimos la media:

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(resistencia ~ madera + I(madera^2),data = ejemplo03)
# Secuencia de valores de predicción
newdata <- data.frame(madera = seq(min(ejemplo03$madera),max(ejemplo03$madera), length = 20))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence")) # Opción interval = "confidence" 
newdata
# Gráfico resultante (predicción y bandas de predicción)
ggplot(newdata, aes(x = madera, y = fit)) +
  geom_line() + 
  geom_ribbon(aes(ymax=upr, ymin=lwr), alpha=1/5) +
  geom_point(data = ejemplo03, aes(x = madera, y = resistencia)) +
  labs(x = "Concentración de madera", y = "Resistencia") +
  theme_bw()

Predecimos una observación:

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(resistencia ~ madera + I(madera^2),data = ejemplo03)
# Secuencia de valores de predicción
newdata <- data.frame(madera = seq(min(ejemplo03$madera),max(ejemplo03$madera), length = 20))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="prediction")) # Opción interval = "confidence" 
newdata
# Gráfico resultante (predicción y bandas de predicción)
ggplot(newdata, aes(x = madera, y = fit)) +
  geom_line() + 
  geom_ribbon(aes(ymax=upr, ymin=lwr), alpha=1/5) +
  geom_point(data = ejemplo03, aes(x = madera, y = resistencia)) +
  labs(x = "Concentración de madera", y = "Resistencia") +
  theme_bw()

En ambos casos se puede ver el comportamiento de las bandas de predicción. Como es habitual, en la predicción de la media tenemos una mayor precisión lo que motiva que se queden más puntos fuera de las bandas de confianza.

Ejemplo 6

Queremos predecir el volumen de madera conseguida (en píes cúbicos) en función de las tres variables predictoras eliminando las observaciones anómalas. En este caso no nos planteamos realizar el gráfico ya que necesitaríamos al menos cuatro dimensiones. Utilizamos cinco valores para cada predictora y obtenemos los intervalos de predicción asociados a cada combinación.

Predecimos la media:

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(vol ~ dbh + d16 + ht ,data = ejemplo06new)
# Valores de predicción para cada predictora
# Como no podemos representar los datos tomamos los valores de los precentiles 25, 50 y 75 de cada predictora
dbhsec <- quantile(ejemplo06new$dbh, probs = c(0.25,0.5,0.75))
d16sec <- quantile(ejemplo06new$d16, probs = c(0.25,0.5,0.75))
htsec <- quantile(ejemplo06new$ht, probs = c(0.25,0.5,0.75))
# obtenemos todas las combinaciones de las secuencias obtenidas
newdata <- data.frame(cbind(dbhsec,d16sec,htsec))
colnames(newdata) <- c("dbh","d16","ht")
# Predicción
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence")) # Opción interval = "confidence" 
newdata

Para representar los efectos de las predictoras sobre la variable odemos utilizar la libreria visreg, y más concretamente al función visreg2d que permite representar la predicción dela respuesta en función de parejas de preditoras.

library(visreg)
visreg2d(ajuste, xvar='dbh', yvar='d16', scale='response')

visreg2d(ajuste, xvar='dbh', yvar='ht', scale='response')

visreg2d(ajuste, xvar='d16', yvar='ht', scale='response')

Bibliografía

Copyright © 2018 Javier Morales. Universidad Miguel Hernández de Elche.