Mínimos cuadrados y la enfermedad COVID-19 (I)

El objetivo de esta entrada (y la que sigue) es hacer algunas consideraciones sobre la diferencia entre estadística descriptiva e inferencia estadística, e ilustrarlas con algún ejemplo aplicado a la epidemia de Covid-19.

1. Modelo empírico de la regresión simple

Supongamos que se recogen los datos experimentales \((x_i,y_i)\), \(i=1,\ldots,n\), donde \(x_i\) e \(y_i\) representan respectivamente valores numéricos de dos variables de interés, \((X,Y)\). En muchas ocasiones un simple análisis gráfico induce a pensar que los puntos \((x_i,y_i)\), \(i=1,\ldots,n\), (nube de puntos) se sitúan de manera aproximada en torno a una curva. Supongamos en esta sección que dicha curva es una recta.

Describiremos a continuación el modelo empírico de la regresión simple. Una vez recogidos los datos este modelo asigna una distribución de probabilidad uniforme sobre los puntos observados:

$$(1)\hspace{4cm}\mathbb{P}(X^*=x_i,\, Y^*=y_i) =\frac{1}{n} ,\;\;i=1,\ldots,n.$$

Son importantes las siguientes consideraciones acerca del modelo (1). \((X^*, Y^*)\) es un vector aleatorio cuya distribución no debe confundirse con la distribución del vector aleatorio \((X,Y)\) que da lugar a los datos observados \((x_i, y_i)\). Técnicamente hablando la distribución de probabilidad de \((X^*, Y^*)\) en (1) debe entenderse como una probabilidad condicionada a los datos observados, esta distribución se denomina distribución empírica de los datos. Las distribuciones empíricas son ampliamente usadas en Estadística, como ejemplo más sobresaliente podríamos citar las técnicas Bootstrap que supusieron una auténtica revolución en los métodos estadísticos desarrollados en las últimas dos décadas del Siglo XX.

Por simplicidad supondremos que los \(x_i’\)s son distintos entre sí y que también lo son entre sí los \(y_i’\)s. Como simple ejercicio podemos obtener a partir del modelo (1):

Las funciones de probabilidad marginales de \(X^*\) e \(Y^*\):

$$\mathbb{P}(X^*=x_i)=\frac{1}{n},\;\;\; \mathbb{P}(Y^*=y_i) = \frac{1}{n},\;\; i=1,\ldots, n.$$

Media empírica, varianza empírica, … etc:

$$\mathbb{E}X^*=\frac1n\sum_{i=1}^nx_i :=\overline{x};\;\;\;
\mathbb{E} Y^* =\frac{1}{n}\sum_{i=1}^n y_i:=\overline{y}$$
$$\text{Var}X^* =\frac{1}{n}\sum_{i=1}^n (x_i-\overline{x})^2:=s_x^2;\;\;\;
\text{Var}Y^* =\frac{1}{n}\sum_{i=1}^n (y_i-\overline{y})^2:=s_y^2;\;\;\;$$
$$\text{Cov}(X^*,Y^*)=\frac{1}{n}\sum_{i=1}^n(x_i-\overline{x})(y_i \overline{y}):=s_{xy}$$

El siguiente paso consiste en modelar la idea anteriormente expresada de que la nube de puntos está en el entorno de una recta, desconocida, \(r(x):=r(x;a,b) = bx+a\). El modelo empírico de la regresión simple es:

$$Y^* = bX^*+a+\epsilon$$

que intenta expresar que los puntos observados en el modelo \((X^*, Y^*)\) se encuentran aproximadamente en la recta \(r(x)=bx+a\) y la razón por la que no están exactamente en la recta es por la presencia de una variable aleatoria \(\epsilon\), denominada error aleatorio.

Los valores \(a\) y \(b\) son desconocidos y se denominan parámetros. En cualquier modelo estadístico los parámetros deben ser inferidos (aproximados o estimados) a partir de los datos observados.

La idea consiste en elegir los valores de \(a\) y \(b\) tales que el error, \(\epsilon = Y^*-bX^*-a\), sea «pequeño». Las diferentes interpretaciones o criterios sobre lo que quiere decir «pequeño» da lugar a diferentes procedimientos de estimación de los parámetros.

Según el criterio de mínimos cuadrados ordinarios los valores de \(a\) y \(b\) se determinan de tal forma que se haga mínimo el momento de segundo orden del error, es decir, debemos resolver el problema de optimización

$$\min_{a,b} \mathbb{E}(Y^*-bX^*-a)^2= \min_{a,b}\,\frac{1}{n}\sum_{i=1}^n (y_i-bx_i-a)^2.$$

La función

$$D(a,b) = \frac{1}{n}\sum_{i=1}^n (y_i-bx_i-a)^2$$

se denomina desviación cuadrática media y es un polinomio en dos variables de grado 2. La minimización se consigue resolviendo el sistema de derivadas parciales igualadas a cero, conocido como sistema de ecuaciones normales,

$$\frac{\partial D}{\partial a}=0,\;\;\frac{\partial D}{\partial b}=0$$

cuya resolución nos proporciona los estimadores de mínimos cuadrados de los parámetros \(a\) y \(b\):

$$\widehat{a} = \overline{y} – \widehat{b}\,\overline{x};\;\;\;\;\widehat{b} = \displaystyle\frac{s_{xy}}{s^2_x}.$$

Es fácil además comprobar que \(\widehat{a}\) y \(\widehat{b}\) proporcionan un mínimo global único pues la función \(D(a,b)\) es estrictamente convexa si, como es usual en la práctica, \(s_x^2>0\).

Las cantidades \(e_i:=y_i -\widehat{b}x_i-\widehat{a}\), \(i=1,\ldots,n\) se denominan residuos y no deben confundirse con el error aleatorio definido anteriormente. El error aleatorio no es observable mientras que los residuos pueden calcularse perfectamente dados los datos y las estimaciones \(\widehat{a}\) y \(\widehat{b}\). El modelo empírico induce sobre los residuos una distribución uniforme. El valor mínimo de la función desviación es por tanto:

$$D(\widehat{a}, \widehat{b}) = \frac{1}{n}\sum_{i=1}^n e_i^2.$$

Llegados a este punto, tenemos la siguiente descomposición:

$$Y^* = (\widehat{b}X^* + \widehat{a}) + e = \text{parte explicada por el modelo} + \text{parte residual}.$$

Un modelo es tanto mejor cuanto mayor sea la parte explicada por el mismo. Las mediciones de las partes explicadas y residuales se hacen a través de sus respectivas varianzas. Puede probarse que el residuo es incorrelado con la parte explicada, es decir (bajo el modelo empírico)

$$\text{Cov}(\widehat{b}X^* + \widehat{a}, e)=0,$$

y entonces

$$\text{Var}Y^* =\text{Var} (\widehat{b}X^* + \widehat{a}) +\text{Var}( e)$$

que en la «jerga estadística» se explica diciendo que la variabilidad total es igual a la variabilidad explicada más la variabilidad residual.

Resulta interesante la proporción de variabilidad que explica un modelo frente a la variabilidad total, que tras algunos cálculos resulta ser

$$\frac{\text{Var} (\widehat{b}X^* + \widehat{a})}{\text{Var}Y^* } = \frac{s^2_{xy}}{s_x^2 s_y^2}:=r^2$$

donde \(r:=s_{xy}/(s_x^2 s_y^2)^{1/2}\) se denomina coeficiente de correlación lineal de Pearson y es frecuentemente utilizado en la práctica como una medida de la asociación lineal entre dos variables. Cuanto más próximo es el valor absoluto de \(r\) a \(1\) mayor es la relación lineal entre ambas variables. La proximidad de \(r\) a cero indica ausencia de asociación de tipo lineal, aunque pueda haber otro tipo de asociación de índole no lineal.

2. Qué es y qué no es el modelo empírico. Estadística Descriptiva e Inferencia Estadística

Es preciso hacer algunas consideraciones sobre el modelo empírico dado en (1). La más importante es que el modelo empírico asigna una probabilidad uniforme a los datos observados como lógicamente se desprende de la definición del modelo y no hace ninguna suposición sobre el vector aleatorio que generó los datos, el vector \((X,Y)\), de ahí el cambio de notación por \((X^*, Y^*)\). La asignación de una probabilidad uniforme a los datos observados se realiza a posteriori, una vez son observados los datos, por eso se comentó anteriormente que el modelo empírico aquí descrito, en rigor, debería estudiarse como un modelo condicionado a los datos observados.

El estudio del modelo empírico cae dentro del marco de la Estadística Descriptiva, cuyo objetivo es simplemente intentar describir mediante resúmenes numéricos (medias, varianzas, cuantiles, etc) o representaciones gráficas (histogramas, boxplots, scatter-plots, etc) el comportamiento de los datos observados en la experimentación. Bajo la perspectiva de la Estadística Descriptiva, los datos \((x_i,y_i)\) no son el resultado o la realización de un vector aleatorio \((X,Y)\), por lo tanto no tiene sentido bajo este marco hablar de distribuciones asociadas al muestreo de las características numéricas calculadas como medias, varianzas, etc. Una consecuencia inmediata de esto es que la inferencia no es posible en este marco.

Por el contrario, la Inferencia Estadística tiene por objeto el estudio de las distribuciones poblacionales, en nuestro caso correspondería al estudio de propiedades del vector aleatorio \((X,Y)\) que puedan inferirse a partir de las observaciones. Bajo esta otra perspectiva, los datos observados \((x_i,y_i)\), \(i=1,\ldots,n\) son considerados como realizaciones de \(n\) vectores aleatorios \((X_i,Y_i)\). El objetivo de la Inferencia Estadística consiste en proponer modelos razonables para la distribución del vector \((X,Y)\), esos modelos pueden depender de parámetros los cuales deben ser estimados a partir de los datos observados, estos estimadores poseen una distribución inducida por el modelo considerado, (a diferencia de lo que ocurría en los estudios descriptivos). La presencia de estas distribuciones de probabilidad inducidas por el modelo permiten la construcción de intervalos de confianza (las famosas horquillas que citan periodistas y políticos), contrastar hipótesis y realizar predicciones, siempre que el modelo sea correcto.

Un error frecuente consiste en confundir los vectores aleatorios \((X^*,Y^*)\) y \((X,Y)\), de hecho la diferenciación que hemos hecho aquí en la notación no es habitual en los libros de texto y ni siquiera en las aplicaciones. Generalmente el contexto determina si estamos tratando con el modelo empírico o con el modelo poblacional.

Como comentario final, vale la pena mencionar que los modelos empíricos no están alejados, al menos teóricamente, de los modelos poblacionales. Cuando el tamaño muestral, \(n\), tiende a infinito las distribuciones empíricas aproximan razonablemente bien a las distribuciones poblacionales de ahí lo que se oye en múltiples ocasiones «la Estadística funciona para tamaños muestrales grandes», pero estas son cuestiones que se alejan del propósito de esta nota.

3. Ejemplo: ¿Ha evolucionado de igual forma el número de pacientes detectados en Andalucía que en el resto de España?

Vamos a intentar dar una respuesta a esta pregunta con el modelo empírico de la regresión simple. Podríamos razonar de la siguiente manera, si la evolución en Andalucía fuera similar a la del resto de España, los datos deberían guardar una cierta proporcionalidad debido a la diferencia de tamaño de las poblaciones, además, posiblemente haya algún desplazamiento debido a posibles diferencias en el comienzo de la epidemia en cada sitio. En términos un poco más precisos, la idea intuitiva anterior podría traducirse en que salvo diferencias en traslación y escala los datos correspondientes deberían ser similares, es decir deberían situarse alrededor de una recta.

Procedemos entonces a ajustar el modelo de regresión simple empírico. Los datos que hemos usado fueron descargados el 09 de Abril de 2020 de la página:

https://www.juntadeandalucia.es/institutodeestadisticaycartografia/salud/index.htm

La nube de puntos considerada es:

$$x_i =\text{número acumulado de pacientes detectados en Andalucía hasta el día \(i\)}$$

$$y_i = \text{diferencia entre el acumulado en España y el correspondiente dato en Andalucía}$$

donde \(i=0,\ldots,30\), y los datos corresponden a las fechas comprendidas entre 10/03/2020 y el 09/04/2020) .

Hemos utilizado el módulo OLS de la librería statsmodels.api en Python 3.6. y hemos obtenido los siguientes valores

$$\widehat{a}=852.74560141,\quad \widehat{b}=15.00957533,\quad r^2=0.9987$$

En la Figura 1 está representada la nube de puntos y la recta estimada.

Podemos ver cómo los puntos se alinean en torno a la recta, lo que, según la idea intuitiva expresada anteriormente, indicaría que no existen diferencias en la evolución del número de pacientes detectados entre Andalucía y el resto de España. Sin embargo, debemos ser cautelosos a la hora de realizar este tipo de afirmaciones; no hemos realizado un proceso de inferencia, simplemente estamos describiendo el comportamiento de los datos observados.

Por otro lado, observamos un valor inusualmente alto del coeficiente de correlación de Pearson, \(r^2\). No es habitual en la práctica estadística observar valores cercanos a 1. Una posible explicación que se me ocurre está relacionada con las Leyes de los Grandes Números. Estas leyes establecen que las sumas de las variables (valores acumulados), convenientemente normalizadas, se estabilizan a medida que aumenta el tamaño de la muestra. Cuando los valores de dos variables se estabilizan, como en nuestro caso, tienden a agruparse en torno a una línea recta.

4. Conclusiones

En un delicioso artículo (ver https://rssdss.design.blog/2020/03/31/all-models-are-wrong-but-some-are-completely-wrong/) Martin Goodson, Director de la sección de datos científicos de la Royal Statistical Society recomienda seis simples reglas a seguir. En nuestra traducción libre son:

Regla 1. Los científicos y periodistas deberían expresar el nivel de incertidumbre asociado

Regla 2. Los periodistas deben obtener opiniones de otros expertos antes de publicar.

Regla 3. Los científicos deberían describir claramente las variables críticas y las hipótesis de sus modelos.

Regla 4. Ser tan trasparentes como sea posible.

Regla 5. Los decisores deberían usar múltiples modelos para informar de sus políticas.

Regla 6. Indicar cuando un modelo ha sido propuesto por alguien sin formación en enfermedades infecciosas.

Seamos críticos con nosotros mismos, ¿hemos cumplido con estas reglas?

Regla 1. Tal como hemos comentado en la sección 2 los modelos presentados tiene un carácter descriptivo y no inferencial, por lo que no sería estrictamente correcto hacer predicciones con los mismos. Las medidas del error que pueden usarse en este contexto son las desviaciones cuadráticas medias y los coeficientes de determinación.

Regla 2. No somos periodistas.

Regla 3. Hemos formulado nuestros modelos, sus variables y descrito sus limitaciones.

Regla 4. Hemos citado las fuentes donde hemos recopilado los datos y el software utilizado, cualquier otra persona podría fácilmente reproducir nuestros resultados.

Regla 5. No somos políticos que debamos tomar decisiones.

Regla 6. Rotundamente, los autores de esta nota no sabemos absolutamente nada de enfermedades contagiosas.

Nota:

Originalmente estas notas estaban destinadas a nuestros alumnos de Segundo Curso de la Facultad de Matemáticas de la Universidad de Sevilla para la asignatura Teoría de la Probabilidad. Desde el IMUS surgió la idea de hacerlas más visibles, y por eso las hemos incluido aquí como entradas.