Modelos Estadísticos. Grado Biotecnología

Introducción

Recuperamos aquí la información presentada en la primer unidad del curso referida a los modelos estadísticos. Cuando el investigador se plantea un diseño experimental y comienza con la recogida de datos es porque persigue el estudio de o verificación de un objetivo planteado sobre la población bajo estudio. Estos objetivos se suelen establecer en base a teorías o hipótesis que se desean verificar sobre le funcionamiento de la población bajo ciertas condiciones experimentales. Por ejemplo: * Teorías que establezcan la posible relación entre dos características de la población. * Teorías que plateen la idea de comportamientos distintas para una característica de la población en función de una variable que clasifica a los sujetos bajo estudio en diferentes grupos.

Un primer paso en la modelización estadística es el planteamiento de una expresión matemática que represente el comportamiento general de la población bajo estudio teniendo en cuenta el diseño experimental establecido y el objetivo u objetivos que se desean verificar. Esto es lo que se conoce como componente sistemática del modelo y se centra únicamente en la parte controlado de nuestro diseño experimental. Por ejemplo, si nos planteamos como objetivo conocer la suma de dos números todo el mundo conoce cual es la función matemática que permite obtener la suma de forma única. Esto se denomina función determinista pues siempre proporciona el mismo resultado si los valores de entrada son iguales. Para poder establecer la parte sistemática del modelo es necesario fijar la variable que se asocia al objetivo o hipótesis planteada (variable respuesta) y la variables o variables (variables predictoras) que se supone están relacionadas con ella a través de la función matemática planteada.

Por otro lado, dado que cada uno de los individuos de la población puede mostrar un comportamiento propio, los datos muestrales de la variable respuesta pueden mostrar ciertas fluctuaciones para un mismo valor de la variable predictora. Estas fluctuaciones pueden ser más o menos pronunciadas dependiendo de la variables que están siendo medidas. Es necesario introducir en el modelo este comportamiento fluctuante que denominamos componente aleatoria del modelo. Dichas fluctuaciones se denominan errores aleatorios y el establecimiento del modelo pasa por asumir la distribución de probabilidad que puede generarlos. Los modelos estadísticos no tratan de explicar este comportamiento aleatorio sino que intentan reducir su impacto a la hora de obtener como se relaciona la respuesta con las predictoras. Cuanto mejor sea el modelo menor será el impacto de esta componente aleatoria y resultará más fácil contestar a los objetivos planteados en el diseño experimental.

Ejemplo Supongamos un diseño experimental donde tenemos una variable respuesta $Y$ (las variables siempre se escriben en mayúsculas) y un conjunto de variables predictoras, $X_1,X_2,...,X_p$. El modelo estadístico general viene dado por la expresión: \[Y = f(X_1,X_2,...,X_p) + \epsilon\] donde $f$ es la función matemática que asumimos en la relación entre las variables predictoras y la respuesta, y $\epsilon$ es la variable aleatoria para los errores cuya distribución de probabilidad se asume que es $F$. La parte sistemática de dicho modelo se corresponde con $f(X_1,X_2,...,X_p)$ y representa el valor esperado de la respuesta ($\widehat{Y}$) para un conjunto de valores de las variables predictoras $X_1 = x_1,…, X_p =x_p $, es decir, $\widehat{Y} = E(Y | X_1 = x_1,X_2 = x_2,...,X_p = x_p)$.

Notación

Antes de continuar introducimos algo de notación. Dada una muestra de $n$ sujetos de los que conocemos sus valores en la variable respuesta $Y_i, i=1,...,n$ y los valores de las variables predictoras:

\[X = \left(\begin{array}{cccc} x_{11} & x_{12} & ... & x_{1p} \\ x_{21} & x_{22} & ... & x_{2p} \\ ... & ... & ... & ... \\ x_{n1} & x_{n2} & ... & x_{np} \\ \end{array}\right)\]

donde cada fila representa las observaciones del $i$-esimo sujeto en cada variable. Si evaluamos la función $f$ para cada fila de $X$, obtendremos el valor esperado de la respuesta para cada sujeto $\widehat{y}_i$, y el residuo para dicho sujeto como la diferencia entre el valor real y el valor ajustado, es decir,
\[r_i = y_i - \widehat{y}_i \text{, } i=1,...,n\] La matriz $X$ se conoce como matriz de diseño, $\widehat{Y} = (\widehat{y}_1,...,\widehat{y}_1)$ es el vector de valores ajustados, y $R = (r_1,...,r_n)$ es el vector de residuos del modelo.

El modelo más sencillo

El modelo más simple que podemos considerar es aquel que resulta al asumir que$f$ es una función lineal sobre todas las variables predictoras, es decir: \[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_p X_p + \epsilon\] donde $\beta_0, \beta_1,...,\beta_p$ son los coeficientes desconocidos del modelo. $\beta_0$ es la interceptación del modelo y el resto de coeficientes son las pendientes asociadas a cada variable predictora.

Cuando tenemos una única variable predictora el modelo resultante viene dado por: \[Y = \beta_0 + \beta_1 X + \epsilon\] que se puede representar gráficamente de forma muy sencilla a través de un diagrama de dispersión de la respuesta con respecto a la predictora, lo que nos permite de una forma muy rápida comprobar la posible relación entre ambas. Veamos un par de ejemplos:

# Cargamos las librerías
library(tidyverse)
library(forcats)
library(broom)

Ejemplo 1.

Treinta aleaciones del tipo 90/10 Cu-Ni, cada una con un contenido específico de hierro son estudiadas bajo un proceso de corrosión. Tras un período de 60 días se obtiene la pérdida de peso (en miligramos al cuadrado por decímetro y día) de cada una de las aleaciones debido al proceso de corrosión. Los resultados obtenidos aparecen a continuación:

# carga de datos
hierro <- c(0.01, 0.48, 0.71, 0.95, 1.19, 0.01, 0.48, 1.44, 0.71, 1.96, 0.01, 1.44, 1.96)
peso <- c(127.6, 124, 110.8, 103.9, 101.5, 130.1, 122, 92.3, 113.1, 83.7, 128, 91.4, 86.2)
ejemplo01 <- data.frame(hierro,peso)

Para este diseño experimental consideramos como variable respuesta la pérdida de peso recogida al final del experimento, y como variable predictora el contenido específico de hierro al comienzo de la experiencia: \[ peso = \beta_0 + \beta_1 hierro + \epsilon\]

A continuación podemos ver el gráfico de dispersión resultante:

# gráfico de dispersión
ggplot(ejemplo01,aes(hierro,peso)) + 
  geom_point() + 
  labs(x = "Contenido en hierro", y = "Pérdida de peso") + 
  theme_bw()

En el gráfico podemos ver que conforme aumenta el contenido de hierro disminuye la pérdida de peso, es decir, la variable respuesta y la predictora están relacionadas de forma inversa. A la vista de la configuración de los puntos podríamos considerar un modelo de tipo lineal que relacione los valores de ambas variables. Para representar dicho modelo utilizamos la opción geom_smooth(method = "lm") en la función gráfica:

# Comportamiento lineal
ggplot(ejemplo01,aes(hierro,peso)) + 
  geom_point() + 
  labs(x = "Contenido en hierro", y = "Pérdida de peso") + 
  theme_bw() + 
  geom_smooth(method = "lm", se = FALSE) # Por el momento añadimos la opción se = FALSE

El modelo obtenido permite establecer la relación entre ambas variables, y permite aproximar el valor de la respuesta (línea azul) para cada valor de la variable predictora. Por ejemplo podemos ver que para un contenido de hierro de 0.25 tendríamos una pérdida de peso próxima a 125.

También podríamos considerar un tipo de modelo no lineal: \[ peso = \beta_0 + \beta_1 f(hierro) + \epsilon\]

# Comportamiento no lineal
ggplot(ejemplo01,aes(hierro,peso)) + 
  geom_point() + 
  labs(x = "Contenido en hierro", y = "Pérdida de peso") + 
  theme_bw() + 
  geom_smooth(method = "loess", se = FALSE)

Se puede ver que la tendencia global (descendente) es similar a la del modelo lineal pero dicha curva se acerca más a los puntos observados. Sin embargo, este modelo es más complejo (desde el punto de vista estadístico) ya que involucra tanto los coeficientes del modelo como el establecimiento de la función $f$. Más adelante estableceremos los criterios para elegir entre uno y otro modelo.

ggplot(ejemplo01,aes(hierro,peso)) + 
  geom_point() + 
  labs(x = "Contenido en hierro", y = "Pérdida de peso") + 
  theme_bw() + 
  geom_smooth(method = "lm", se = FALSE) +
  geom_smooth(method = "loess", se = FALSE ,colour = "red")

Ejemplo 2.

Se propone a una empresa que fabrica vasos de cristal un nuevo proceso de control de calidad. Hasta ahora la empresa seleccionaba una caja de vasos al final de la fabricación y observaba si había alguno roto. Esto provocaba un gran gasto ya que en caso de encontrar algún defecto la caja se desembala y los vasos vuelven a la cadena de embalaje. Ahora se propone seleccionar vasos antes de embalar y determinar así el porcentaje de defectos. Se desea saber si ambos porcentajes están relacionados. Los datos aparecen a continuación:

# carga de datos
cajas <- c(3.0, 3.1, 3.0, 3.6, 3.8, 2.7, 3.1, 2.7, 2.7, 3.3, 3.2, 2.1, 3.0, 2.6)
vasos <- c(3.1, 3.9, 3.4, 4.0, 3.6, 3.6, 3.1, 3.6, 2.9, 3.6, 4.1, 2.6, 3.1, 2.8)
ejemplo02 <- data.frame(cajas,vasos)

En este caso tratamos de utilizar un procedimiento menos agresivo (detección de vasos defectuosos) que el que se venia utilizando (cajas devueltas con vasos rotos). Pero para poder hacer esto debemos tener certeza de que los resultados con este nuevo procedimiento son similares a los que obtendríamos con el método antiguo. Utilizamos como variable predictora el porcentaje de vasos defectuosos. A continuación podemos ver el gráfico de dispersión resultante:

# gráfico de dispersión
ggplot(ejemplo02,aes(vasos,cajas)) + 
  geom_point() + 
  labs(x = "% vasos defectuosos", y = "% de cajas con vasos defectuosos") + 
  theme_bw()

En este caso podemos ver como el diagrama de puntos está menos concentrado que en el ejemplo anterior, es decir, para un mismo valor de la predictora hay mucha variabilidad en la respuesta (valor de 3.6). Esto implica que aunque podamos usar un modelo de tipo lineal su capacidad para explicar el comportamiento de la respuesta será más reducido. Como antes representamos un modelo lineal y otro no lineal con ecuaciones:

\[ cajas = \beta_0 + \beta_1 vasos + \epsilon\] \[ cajas = \beta_0 + \beta_1 f(vasos) + \epsilon\]

# Comportamiento lineal
ggplot(ejemplo02,aes(vasos,cajas)) + 
  geom_point() + 
  labs(x = "% vasos defectuosos", y = "% de cajas con vasos defectuosos") + 
  theme_bw() + 
  geom_smooth(method = "lm", se = FALSE) + 
  geom_smooth(method = "loess", se = FALSE, colour = "red")

¿Cuál de las dos tendencias ajustadas se comporta mejor? ¿en que basas tus conclusiones?

Tipos de modelos

Como hemos podido ver en los ejemplos anteriores la suposición de un modelo que relacione la variable respuesta con las predictoras no es única. Los modelos a considerar dependen de:

Tipo de la variable respuesta: Se pueden modelizar tanto variables respuestas de tipo numérico como categórico. Cuando tenemos variables de tipo categórico deberemos transformar nuestro modelo para obtener relaciones sencillas. En los ejemplos anteriores las dos respuestas eran de tipo numérico.
Tipo/s de la/s varaible/s preditora/s: Se pueden dar todas las combinaciones posibles entre variables de estos tipos (numéricas, categóricas, y numéricas y categóricas). Si intervienen variables de tipo categórico hablaremos de los posibles efectos de interacción, que son aquellos que aparecen cuando el modelo presenta una pendiente diferente para una predictora en función de los diferentes grupos de sujetos de conforma la variable categórica.
Función de relación: Se pueden considerar una gran cantidad de funciones $f$, aunque nosotros nos centraremos en las situaciones más sencillas (generalmente de tipo lineal). En algunos casos partiremos de modelos teóricos establecidos y será necesario obtener funciones $f$ que nos permitan una fácil interpretación de resultados.
Hipótesis sobre la componente aleatoria del modelo. Las hipótesis que establezcamos sobre los errores aleatorios (distribución $F$) es la parte más importante del modelo, ya que los procedimiento para la obtención de los coeficientes del modelo dependen de ellas. Deberemos verificar que los residuos obtenidos podrían ser generados por la función $F$, lo que en ocasiones nos llevará a desechar el modelo considerado y buscar uno nuevo.

En los temas siguientes iremos presentando las diferentes opciones de modelización en función de los puntos anteriores. Nuestro objetivo no es mostrar el desarrollo teórico de todos ellos, sino más bien su utilización en la práctica. En concreto veremos los modelos lineales, los modelos lineales generalizados, los modelos de supervivencia, los modelos de suavizado y los modelos con efectos aleatorios.

Fases en la construcción de un modelo

La construcción de un modelo estadístico es un proceso continuo de búsqueda de la mejor relación que podamos obtener entre la variable respuesta y las predictoras. Por tanto, nunca estaremos seguros de haber encontrado el mejor modelo posible sino el menos malo entre todos los que hemos probado.

Los pasos para la obtención de un modelo estadístico son:

Selección de variable respuesta y predictoras: En base al diseño experimental se deben establecer la variable respuesta y las posibles variables predictoras, lo que permitirá decidir a priori el tipo de modelo más adecuado. Para la elección del modelo será necesario realizar una representación gráfica que nos permita visualizar de forma rápida el comportamiento de todas las variables involucradas.
Hipótesis del modelo: Las hipótesis del modelo hacen referencia al comportamiento de la parte aleatoria del modelo. De forma general, se pide que los errores obtenidos sean independientes (hipótesis de independencia), que su media sea cero (hipótesis de linealidad), que su variabilidad sea constante (hipótesis de homocedasticidad), y deberemos asumir una forma para la distribución $F$ de los errores. La hipótesis de independencia depende del diseño experimental establecido, mientras que las otra tres dependen del modelo asumido para explicar la respuesta.
Estimación del modelo: Esta fase nos permite obtener las estimaciones de los coeficientes del modelo propuesto y los residuos correspondientes. También podremos realizar el correspondiente análisis inferencial de dichos parámetros (intervalos de confianza y contrastes de hipótesis) lo que nos permitirá construir unas bandas de confianza para la ecuación del modelo obtenida.
Bondad de ajuste del modelo: Junto con la estimación del modelo podemos obtener las denominadas medidas de bondad del ajuste que nos permiten concluir sobre si el modelo obtenido tiene capacidad explicativa. La capacidad explicativa representa lo bien que las variables predictoras pueden explicar el comportamiento de la respuesta a través del modelo asumido.
Comparación de modelos: Cuando el modelo contiene diferentes variables predictoras se establecen procedimientos para seleccionar el conjunto de ellas que mejor explican el comportamiento de la respuesta, es decir, no todas tienen que estar representadas en el modelo final. Un principio básico en la modelización estadística es la obtención de modelos parsimoniosos, es decir, construir el modelo lo más sencillo posible que explique el comportamiento de la respuesta en función de las variables explicativas despreciando todas aquellas predictoras que influyan poco o nada en el comportamiento de la respuesta.
Diagnóstico del modelo: El diagnóstico del modelo incluye todos los procedimientos que tratan de verificar si los residuos del modelo obtenido cumple con las hipótesis de partida. El incumplimiento de alguna de estas hipótesis obliga a modificar el modelo propuesto o desecharlo y proponer uno nuevo.
Predicción: Una vez validado el modelo propuesto (cumplimiento de hipótesis) solo nos resta utilizarlo para predecir el comportamiento de la respuesta. “Predecir”" significa utilizar un conjunto de valores de las variables predictoras que sustituiremos en el modelo obtenido para predecir el valor de la respuesta. Los valores a considerar son todos aquellos que queden dentro del rango de valores posibles de cada predictora (aunque no los hallamos observado de forma directa en el diseño experimental). Cuando los valores de las predictoras son los obtenidos en la muestra se denomina ajuste en lugar de predicción.

En los temas siguientes

En los temas siguientes iremos tratando cada uno de los posibles modelos desarrollando todo el proceso de construcción, validación y predicción en cada uno de ellos. Obviaremos en la medida de lo posible los desarrollamos matemáticos involucrados, y nos centraremos en la resolución práctica de cada tipo de modelo. Se introducirán las herramientas de R necesarias para resolver cada una de las fases de construcción del modelo.

El material de los temas siguientes se puede completar con Mayoral and J. (2001), Krzanowski (1998), Ugarte, Militino, and Arnholt (2008), Maindonald and Braun (2007), Dobson (2002), Wood (2006), y Faraway (2006)

Bibliografía

Dobson, A. 2002. An Introduction to Generalized Linear Models. CRC Press.

Faraway, JJ. 2006. Extending the Linear Model with R. Chapman & Hall.

Krzanowski, Wojtek J. 1998. An Introduction to Statistical Modelling. Arnold.

Maindonald, J., and J. Braun. 2007. Data Analysis and Graphics Using R. an Example-Based Approach. Second. Cambridge University Press.

Mayoral, A.M., and Morales J. 2001. Modelos Lineales Generalizados. Servicio de publicaciones de la UMH.

Ugarte, MD, AF Militino, and A Arnholt. 2008. Probability and Statistics with R. CRC Press.

Wood, SN. 2006. Generalized Additive Models. Chapman & Hall.

Introducción a los Modelos Estadísticos