Introducción al pensamiento bayesiano

Introducción

Introduciremos el pensamiento bayesiano con el análisis de una proporciòn poblacional. Antes de tomar una muestra, uno puede tener una buena idea de cuàl es el valor de la proporción poblacional. Esta información se vuelca en una distribuciòn a priori, propia del teorema de Bayes:

$$g[p|data] \propto g(p)*L(data|p)$$

Donde vemos que la distribución a posteriori según los datos observados $g[p|data]$ es proporcional a la verosimilitud de la muestra para dicho valor p y nuestras creencias $g(p)$.

Ilustraremos las diferentes maneras de llevar estos conocimientos a una distribución a priori.

Aprendiendo acerca del sueño de los estudiantes

Estamos interesados en estudiar los tiempos de reposo de los estudiantes universitarios ^1. Siendo qu ese recomiendan ocho horas de sueño, ¿qué proporción de estudiantes duerme al menos la cantidad de horas recomendadas?

Sea p dicha proporción. Estamos interesados en conocer el comportamiento de p. Desde un punto de vista bayesiano, p es desconocido, pero tiene una distribución de probabilidad que refleja nuestras conjeturas a priori acerca de su valor. Una muestra aleatoria puede ayudarnos a conocer algo de p, pero tenemos mucha información disponible antes de eso.

1. Disponemos de una nota de una revista de medicina donde dice que los estudiantes duermen 6 horas por día.

2. Otra nota refuerza esta opinión con un estudio que dice que durante la semana el 70% de los estudiantes duermen de 5 a 6 horas, un 28% de 7 a 8 hs, y un 2% las 9 horas saludables.

3. También conocemos que es inviavle dormir menos de cuatro horas de forma sostenida.

Basada esn esta información, podríamos esbozar una distribución de p antes de observar una muestra. p debiera ser menor que 0.5. Probablemente el mejor valor sobre el cual arriesgar sería 0.3, pero con la mayoría de las posibilidades entre 0 y 0.5.

Luego de este ejercicio previo, vemos que disponemos de una muestra de 27 estudiantes con la cual validaremos y perfeccionaremos nuestro conocimiento a priori de p. Observamos 11 estudiantes que durmieron al menos 8 horas. Nos proponemos dos desafíos:

1. estimar la proporción p con toda la información disponible.

2. predecir el número de estudiantes que duerma al menos 8 horas en una nueva muestra de 20 estudiantes.

Supongamos que nuestra función de probabilidad a priori de p es g(p). Si consideramos un ”éxito” el dormir al menos 8 horas y tomamos una muestr aaleatoria con s ”éxitos” y f “fracasos”, entonces la función de verosimilitud de p en dicha muestra es:

$$L(p) \propto p^s(1-p)^f \qquad 0 < p <1$$

La distribución a posteriori de p, es obtenida proporcionalmente con el teorema de Bayes, tal como se mostró al inicio de la nota:

$$g[p|data] \propto g(p)*L(p)$$

Veremos tres métodos comunes para obtener las distribuciones a priori g(p) y a posteriori g(p|data).

Utilización de una función a priori discreta

Creemos que nuestra distribución a priori de p es así:

p	.05	.15	.25	.35	.45	.55	.65	.75	.85	.95
peso	1	5.2	8	7.2	4.6	2.1	0.7	0.1	0	0

Convertiremos los pesos de nuestras crrencias en probabilidades dividiendolos por su suma. Definiremos en R a p como el vector de proporciones y prior lso correspondientes pesos normalizados, o probabilidades asociadas. El comando plot nos permitirá graficar nuestra distribución a priori.

p = seq(0.05, 0.95, by = 0.1)
prior = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
prior = prior/sum(prior)
plot(p, prior, type = "h", ylab="Probabilidad a priori")

^1: Muestra obtenida de un colegio americano de estudiantes

Cuidado con las sumarizaciones

Me ha gustado este ejemplo del que surge una advertencia para aquel analista que realice conclusiones sobre totales en una tabla de doble entrada.
El siguiente cuadro presenta los resultados de dos tipos de tratamiento para cálculos renales. Las tasas de éxito totales (78% para cirugía abierta y 83% para tratamiento por ultrasonido). No obstante, cuando se comparan los métodos según el tamaño del cálculo, en ambos casos la tasa de éxito favorece a la cirugía abierta.

	Pequeña			Grande			Total
	cir abierta	ultrasonido		cir abierta	ultrasonido		cir abierta	ultrasonido
éxito	81	234	éxito	192	55	éxito	273	289
fracaso	6	36	fracaso	71	25	fracaso	77	61
Tasa de éxito	93%	87%	Tasa de éxito	73%	69%	Tasa de éxito	78%	83%

Los peligros de malinterpretar resultados están a la vista. Veamos cómo hacerlo en R:

piedras <- array(c(81,6,234,36,192,71,55,25), dim=c(2,2,2),


dimnames=list(Resultado=c("Exito","Fracaso"),


Método=c("cir abierta","ultrasonido"),


Tamaño=c("<2cm", ">=2cm")))


mosaicplot(piedras, sort=3:1)







  




El diagrama de mosaico muestra que la tasa general de éxito, para la cirugía abierta está sesgada hacia los resultados por cirugía abierta en piedras grandes, miestras que el ultrasonido está sesgado hacia los resultados con càlculos pequeños.

Introducción a los métodos de reservamiento

En esta nota explicaremos brevemente el proceso de un siniestro y como surge de ello la necesidad de reservar.

Proceso del siniestro

Asumamos que estamos ante una situación de "riesgo", asociado con una cobertura de seguro. Sus características esenciales son la existencia de una persona o corporación que se ve afectada por acontecimientos fortuitos dentro de un período definido. Un ejemplo sería un individuo con seguro de daños en su coche privado sufriendo un choque que le genera un costo de reparación.
Si existe cobertura de seguro, la situación será una pérdida asegurada, que se convertirá en un siniestro para la aseguradora.

Normalmente habrá un retraso entre el hecho causante del siniestro y su liquidación monetaria, que puede tratarse de días o años, según se trate el tipo de daño cubierto. Otras fechas importantes participan en el proceso. Después de la ocurrencia del siniestro, viene la fecha en la cual el asegurador se entera del evento. Este período dependerá del tipo de riesgo asociado. Por ejemplo, cuando un barco está dañado en el puerto, pero el daño se hace evidente sólo cuando está el dique seco en una fecha posterior.

Reserva de casos

La demora entre el evento y la liquidación significa que el asegurador debe establecer "reservas" respecto de los siniestros pendientes de liquidación (SPL). Dichas reservas son recursos necesarios para equiparar y cubrir los costos que van surgiendo de las reclamaciones no firmemente establecidas en el momento del evento dañoso.
Las reservas de casos cubren así eventos ya ocurridos, que es distinto a cubrir siniestros futuros derivados de los riesgos cubiertos por la aseguradora hasta la extinción de la cobertura de la póliza. Estos pasivos se cubren con la reserva de prima no devengada o reserva de riesgos en curso, que no trataremos aquí.
Nos preocupan la constitución de reservas para siniestros ocurridos hasta una fecha particular (la fecha de "valuación").
Aquí se hace una clasificación entre los siniestros ocurridos y reportados a la aseguradora, y los siniestros ocurridos pero aún no reportados a la aseguradora (IBNR).
Muchas veces se incluye en el segundo item, la adición de reservas por siniestros no correctamente evaluados y necesita una adición futura de reserva (se conoce como IBNER).
Podemos ver a las reservas como un intento de atribuír valor financiero a estos pagos aún no realizados. Por supuesto que la presición no existe acá y las reservas deben ser estimadas. Distintas suposiciones sobre las futuras influencias en los resultados de estos siniestros llevarán a una menor o mayor estimación de las reservas, lo que implica un esfuerzo financiero diferente que se verá en la hoja de balance de la compañía de seguros.
Es por ello que la estimación de las reservas debe realizarse cuidadosamente y tener en claro cuáles son sus supuestos de base.
Enumeramos a continuación razones por las cuales las reservas son necesarias:
- evaluar la situación financiera de la aseguradora. Los movimientos de dichas reservas son esenciales para evaluar su progreso.
- calcular el costo final de un negocio en el sentido de que es necesario conocer los costos futuros aún no liquidados de siniestros ocurridos en el pasado.
- evaluar la solvencia del asegurador, para ver si sus activos disponibles pueden hacer frente a sus compromisos futuros por siniestros que le han contraído obligaciones.
- determinar el valor del patrimonio neto de los accionistas, sobretodo con la intención de calcular ganancias y repartición de dividendos sin afectar la solvencia de la aseguradora.
- cerrar o ceder un rama de negocio con un reasegurador, poniendo un valor económico a la cartera de siniestros actual.

Hasta principio de los años '70, el enfoque de reservas era normalmente sobre siniestros conocidos. La práctica consistía en evaluar individualmente cada siniestro periódicamente y estimar su valor vigente. De allí el nombre "reserva de casos" para las reservas de siniestros conocidos. Todas estas estimaciones individuales se agregaban para formar una gran reserva total por créditos no pagados. Con el paso del tiempo y el aumento de la capacidad informática de someter los datos a un escrutinio estadístico más profundo, se comenzó a ver que estos métodos no eran los más adecuados.
También a veces el volumen de siniestros y evoluciones hace que sea imposible mantener una evaluación individual y periódica como se plantea. También el IBNR no puede tener una estimación individual, debido a que no se conoce aún el siniestro  (se necesita otro tipo de técnica de reservamiento).
El estimador de casos - liquidador o gestor de siniestros - requiere de acceso a los detalles de la póliza cubierta, a los siniestros y su historial, antes de poner realizar una estimación de reserva. De manera similar, los métodos de estimación estadística también precisan de acceso a bases de información relacionados con el grupo de siniestros a reservar.
Los datos más relevantes para realizar una reserva refieren a los costos incurridos, el número de siniestros reportados, y de ellos discriminar el grupo de siniestros terminados y los que siguen vigentes. También se tiene en cuenta otro tipo de información más "blanda" quedan característica a subgrupos de siniestros.
También es necesario establecer un aserie de supuestos para poder proyectar los costos futuros. Aquí es esencial el uso de la información blanda. Otro punto importante a tener en cuenta es la naturaleza y calidad de los datos sobre los que el método se basará.
Por todas estas razones, es crucial que la persona encargada de confeccionar las reservas esté familiarizado con el negocio subyacente.

Gráficos de autocorrelación

Los gráficos de auto-correlación comenzaron a utilizarse por Box y Jenkins en análisis de series de tiempo.
Son comúnmente utilizados para chequear la aleatoriedad de los datos, que se comprueba analizando la relación entre pares de datos de la muestra, a diferentes rezagos (lags). Si existe una correlación significativa en algún nivel de rezago, se puede decir que los datos no son completamente aleatorios.
También se utilizan estos tipos de gráficos para comenzar a construir modelos auto-rregresivos y de medias móviles.
Hay tres tipos básicos de gráficos de auto-correlación:

• Gráfico de auto-correlación puro
• Gráfico de auto-correlación parcial
• Gráfico de auto-correlación cruzada

Gráfico de auto-correlación puro

Mide el nivel de correlación entre una variable aleatoria $Yt$ y $Y(t+h)$, h períodos después. El rezago h se va variando en el gráfico para analizar un amplio espectro.

El siguiente gráfico de auto-correlación fue ejecutado con el comando acf (o plot.acf, es lo mismo) del lenguaje R. :

a<-acf(presidents, na.action = na.pass,lwd=3,col="red",main="")
title(main="Ranking Trimestral de Aprobación",font.main= 2, line=2)
title(main="Presidentes de Estados Unidos 1945-1974",font.main= 1, line=1)

Esta autocorrelación en la muestra de satisfacción de presidentes de EEUU  muestra que la serie de tiempo no es aleatoria debido a que los valores observados superan los umbrales de aleatoriedad.

Explicamos en más detalle los elementos del gráfico:

En el eje vertical se tiene la medida del coeficiente de auto-correlación.
$$R_h=C_h/C_0$$
donde $C_h$ es la función de auto-covarianzas.
$$C_h=\frac{1}{N}\sum\limits_{t=1}^{N-h}{(Y_t-\bar{Y})(Y_{t+h}-\bar{Y})}$$
donde $C_0$ es la función de varianza.
$$C_0=\frac{1}{N}\sum\limits_{t=1}^{N}{(Y_t-\bar{Y})^2}$$
Al ser la varianza por definición mayor al valor de la covarianza, $R_h$ se encuentar entre -1 y 1.

Hay que prestar atención que a veces la función de auto-covarianzas viene definida de forma insesgada:
$$C_h=\frac{1}{N-h}\sum\limits_{t=1}^{N-h}{(Y_t-\bar{Y})(Y_{t+h}-\bar{Y})}$$
Aunque esta variante cuenta con menos propiedades algebraicas que al dividir por N.

En el eje horizontal se detallan los valores del rezago o lag h (h=0,1,2,3,..).

La función acf maneja este límite con el argumento lag.max.

Las líneas punteadas en azul destacan un intervalo de confianza para $R_h$. En el gráfico se tomaron al 95% (valor predeterminado de la función R acf, que por lo que estuve averiguando no se puede modificar en el gráfico)

Hay dos formas distintas de generar las bandas de confianza:

1. Si el gráfico de auto-correlación se está utilizando para detectar no-aleatoriedad (del tipo dependencia del tiempo), se recomienda la fórmula
   $$\pm\frac{z_{1-\alpha/2}}{\sqrt{N}}$$
   donde N es el tamaño de la muestra y z el percentil elegido de una distribución normal para un nivel de confianza $\alpha$.
2. Los gráficos de auto-correlación, como ya hemos mencionado, sirven para estimar modelos de ajuste ARIMA.  En dicho caso, se asume un promedio móvil en la serie y los intervalos se calculan con la siguiente formula:
   $$\pm{z_{1-\alpha/2}}\sqrt{\frac{1}{N}({1+2\sum\limits_{i=1}^{h}R_i^2)}}$$
   donde h es el lag, N es el tamaño de la muestra y z el percentil elegido de una distribución normal para un nivel de confianza $\alpha$. En este caso, el intervalo de confianza crece cuando el lag h crece.
   Se trata de la curva verde en el gráfico. Puede lograrse en R escribiendo los siguientes comandos:
   

   z <- qnorm(0.975)
   y <- c(z*sqrt(1/a$n.used),z*sqrt(1/a$n.used*(1+2*cumsum(a$acf^2))))
   lines(a$lag,y[-length(y)],col="green")

El gráfico de auto-correlación puede brindar respuesta a las siguientes incógnitas:

1. ¿son aleatorios los datos en el sentido del tiempo?
2. ¿son los datos un ruido blanco?
3. ¿son los datos una serie sinusoidal?
4. ¿se observan auto-rregresiones significativas?
5. ¿qué modelo parece seguir la serie de tiempo observada?
6. ¿es el modelo Y = constante + error válido o insuficiente?
7. ¿es válida la fórmula $s_\bar{Y}=s/\sqrt{N}$?

Gráfico de correlaciones parciales

Cuando se mide la auto-correlación de una variable aleatoria Y entre dos momentos t y t+h, pasa por momentos intermedios que condicionan el resultado final y estamos midiendo la resultante de varios lags $Y(t+1)$,$Y(t+2)$,....,$Y(t+h)$.
La auto-correlación parcial estima la relación entre $Y(t)$ y $Y(t+h)$, sin los efectos añadidos por los lags de 1 a  h-1.
Los gráficos de auto-correlación parcial son útiles para identificar comportamientos auto-regresivos en los procesos. La auto-correlación parcial de un proceso AR(p) es cero a partir de un rezago igual o superior a p+1. Si el análisis de la muestra indica que un modelo auto-regresivo puede ser utilizado, entonces el se examina el gráfico de auto-correlación parcial para conocer el orden.
El gráfico se lee de la misma manera que el gráfico de auto-correlación: se busca el punto a partir del cual las auto-correlaciones parciales comienzan a no ser significativas respecto a la hipótesis nula (auto-correlación cero). Un intervalo de confianza cercano al 95% es tomar bandas de $2\sqrt{N}$.

Realicemos el gráfico de auto-correlación parcial del ejemplo anterior. El código es el mismo sólo que el comando de R cambia de acf a pacf.

pacf(presidents, na.action = na.pass,lwd=3,col="red",main="")
title(main="Ranking Trimestral de Aprobación",font.main= 2, line=2)
title(main="Presidentes de Estados Unidos 1945-1974",font.main= 1, line=1)

Como puede verse, el gráfico de auto-correlación puede resultar engañoso. La gran cantidad de rezagos afectados lo es debido únicamente al primer lag y la alta intensidad que tiene hizo que queden afectados varios niveles. Se trata entonces de un modelo de proceso auto-regresivo de orden 1 AR(1).

Dejo el siguiente link a una aplicación online muy interesante donde se muestra justamente le efecto contagio del primer coeficiente de un proceso AR(1) en un gráfico de auto-correlación.

Resumimos las preguntas que se pueden formular con este gráfico:

1. ¿es un modelo auto-regresivo AR conveniente para el modelo?
2. ¿de qué orden p es el modelo AR?

Gráfico de auto-correlación cruzada

Mide el nivel de correlación entre una variable aleatoria $Yt$ y otra variable $X(t+h)$, h períodos después. Desde el punto de vista gráfico, no hay diferencias respecto al gráfico de auto-correlación puro.

G. E. P., and Jenkins, G. (1976), Time Series Analysis: Forecasting and Control, Holden-Day

Cuando los supuestos subyacentes no se cumplen

En nuestra nota anterior hemos descripto los cuatro supuestos básicos que subyacen en una mayoría de modelos estadísticos. Ahora investigaremos qué sucede cuando alguno(s) de ello(s) no se cumple(n).

Los cuatro supuestos básicos eran:

1. aleatoriedad
2. posición fija
3. dispersión fija
4. distribución fija 

Consecuencias de no aleatoriedad

Si el supuesto aleatoriedad no se cumple, entonces:

1. Todos los tests estadísticos quedan invalidados.
2. Los rangos de variación pierden significado.
3. El cálculo de muestra mínima o suficiente carece de sentido. En un escenario totalmente no aleatorio, tomar muestras carece de sentido.
4. El modelo respuesta = constante + error pierde sentido.
5. La estimación de sus parámetros a través de muestras se torna sospechoso.

Un caso específico de no aleatoriedad es la existencia de autocorrelación. La autocorrelación es la relación que existe entre la variable aleatoria $Y_t$ y la misma variable k períodos antes $Y_{t-k}$. La autocorrelación sería una no-aleatoriedad en base al tiempo de observación. La autocorrelación se detecta con el gráfico de lags que hemos visto o con un correlograma.

Si los datos están afectados por autocorrelación, entonces:

1. Los datos adyacentes están relacionados.
2. No se tienen n muestras independientes.
3. Los outliers son más difíciles de detectar.

Consecuencias de no existir un parámetro de posición fijo

La estimación usual de posición es la media muestral:
$${\bar{Y} = \frac{1}{N}}{\sum\limits_{j=1}^{N} Y_i}$$
Si el supuesto de posición fija no se cumple,

1. La posición estimada puede tener tendencia.
   
2. La fórmula usual de incertidumbre de la media:  $${s(\bar{Y}) = \frac{1}{\sqrt{N(N-1)}}}\sqrt{\sum\limits_{j=1}^{N} {(Y_i - \bar{Y})}^2}$$ puede ser inválida y su valor numérico sería menor.
   
3. La estimación de otra medida de posición simple sería pobre y sesgada. 

Consecuencias de no existir un parámetro de dispersión fijo

La estimación usual de dispersión es el desvío estándar:
$${s(Y) = \frac{1}{\sqrt{N-1}}}\sqrt{\sum\limits_{j=1}^{N} {(Y_i - \bar{Y})}^2}$$
Si el supuesto de dispersión fija no se cumple,

1. La varianza puede tener tendencia.
2. La estimación simple de la varianza pierde significado y puede estar sesgada.

Consecuencias de no conocer la distribución subyacente

Rutinariamente se utiliza la media (promedio) para estimar el "medio" de una distribución de probabilidad. Así también se utiliza la fórmula de desvío estándar para estimarle un rango de variación:
$${s(\bar{Y}) = \frac{1}{\sqrt{N(N-1)}}}\sqrt{\sum\limits_{j=1}^{N} {(Y_i - \bar{Y})}^2}$$
No es muy conocido que dicha fórmula tiene supuestos de una distribución normal subyacente (aparte de un buen tamaño de muestra), y si no se cumple, estos márgenes quedan invalidados.
Para algunos tipos de distribuciones, sobre todo las distribuciones con una distribución asimétrica importante, la media como medida de posición resumen es una elección pobre. Cada distribución tiene una medida de posición óptima para elegir; se recomienda elegir aquella que menos variabilidad y ruido tenga. Esta elección óptima podría ser, por ejemplo, la mediana, el rango medio, la moda, etc... Esto implica tener que estimar la distribución subyacente, y una vez conocida la distribución de probabilidad, elegir el estimador de posición.
Trataremos estas técnicas de estimación en una nota futura. lo importante es estar conscientes que si el histograma y gráfico p-p Plot no muestran la distribución esperada (habiendo ya afirmado la componente determinística del modelo), es mejor dedicar tiempo a estimar la distribución subyacente.

Otras consecuencias de subestimar la distribución subyacente son:

 En la distribución

1. La  distribución subyacente puede ir cambiando.
2. Si sucede lo del punto 1, una estimación por única vez de la distribución subyacente es insuficiente y quedará obsoleta.
3. La distribución subyacente puede ser marcadamente no normal
4. La verdadera distribución de probabilidad para el error puede mantenerse desconocida y deberemos conformarnos con una aproximación. 

 En el modelo

1. El modelo puede ir cambiando.
2. Si sucede lo del punto 1, una estimación por única vez de la distribución subyacente es insuficiente y quedará obsoleta.
3. El modelo respuesta = constante + error puede quedar inválido.
4. Si no se tiene identificada la distribución correcta, una mejora necesaria al modelo puede quedar oculta.

 En el proceso

1. El proceso puede estar fuera de control y no detectarse.
2. El proceso puede ser impredecible y no ser consciente de ello. 
3. El proceso puede ser no modelable.

Los 4 supuestos subyacentes básicos

Existe una base de supuestos estadísticos generales sobre los que mas o menos se posan gran parte de los fenómenos aleatorios que queremos analizar.
En esta nota trataremos de conceptualizarlos y propondremos un set de herramientas para una primer mirada exploratoria.

Supuestos subyacentes

Existen cuatro supuestos muy típicos en los análisis estadísticos, es decir, que los los datos "se comportan como...". Y son los siguientes:

1. Aleatoriedad, o ruido blanco para series de tiempo
2. Proveniencia de una distribución dada
3. La distribución tiene una medida de posición fija dada
4. La distribución tiene una medida de dispersión fija dada

La posición fija mencionada en el punto 3 puede diferir según se trate del tipo de proceso analizado, pero para problemas univariados, utilizar este supuesto implica convertir un problema típico

respuesta = componente determinístico + componente aleatorio

a lo siguiente

respuesta = constante + error

la constante sería nuestra medida de posición fija a determinar. Así pues, podemos imaginar el proceso a la mano de estar operando bajo condiciones constantes que producen una sola columna de datos con las propiedades que:

• los datos no están correlacionados entre sí;
• la componente aleatoria tiene una distribución fija;
• el componente deterministico consiste en sólo una constante, y
• el componente aleatorio tiene una variación fija.

El poder universal y la importancia de este modelo univariado es que puede ser fácilmente extendido a problemas más generales donde el componente determinístico no es una constante, sino en efecto una función de varias variables.

El punto clave es que, independientemente de cuántos factores haya, e independientemente de lo complicado que sea la función, y si el analista tiene éxito en la elección del modelo, entonces las diferencias (residuos) entre los datos de respuesta y los valores pronosticados por el modelo deben comportarse como un proceso univariado, el cual tendrá:

• Aleatoriedad, o ruido blanco para series de tiempo
• una distribución dada
• una medida de posición fija dada (cero en este caso)
• una medida de dispersión fija dada

Así, los residuos se comportan de acuerdo al modelo ideal que se asumió. Entonces las pruebas de los supuestos subyacentes sobre los residuos se convierten en una herramienta para la validación y la calidad de ajuste del modelo elegido.
Por otro lado, si los residuos del modelo escogido violan uno o más de los supuestos univariantes anteriores, entonces el modelo escogido es insuficiente y existe una oportunidad de hallar un modelo mejor.

Importancia

La predictibilidad es por lejos el objetivo más importante de un análisis estadístico. Pero para que se cumpla dicho objetivo debemos tener el proceso "bajo control estadístico", condición que se cumple cuando las cuatro condiciones mencionadas se comprueban y tenemos suficientes seguridades de que se seguirán cumpliendo en el futuro.

Técnicas gráficas para comprobar los cuatro supuestos básicos

Existe un innumerable conjunto de herramientas numéricas y gráficas para testear los cuatro supuestos. Nosotros nos centraremos en un conjunto reducido de técnicas gráficas, que aprenderemos a desarrollar con lenguaje R. Considero que las mismas cumplirán de manera simple y eficiente el análisis EDA de cualquier serie univariada. Utilizaremos el siguiente set de gráficos que pasaremos a llamar "4-plot":

• gráfico de secuencias
• gráfico de rezagos o diagrama de fase
• histograma
• gráfico de probabilidad normal o P-P Plot

Ya he tratado el histograma y otras no mencionadas (diagrama de hojas, densigrama, cajas) en una de mis primeras notas de EDA. El conjunto elegido es simple de entender y se presta muy bien a comprobar los 4 supuestos mencionados.

Construcción de los gráficos con R

A continuación detallo el código necesario para construir los gráficos. Dejé el código comentado y muy simple para que pueda ser reutilizado facilmente:

## 4-Plot
## Generando los gráficos de 4-plot para los datos.

## La variable con datos es y. Es un array univariado. No funciona con data.frame  

## El data.frame se puede convertor a array con el comando y <- as.array(y)  

##variable de conteo
n <- length(y)
t <- 1:n
## parámetros de gráfico
## mfrow: nro de gráficos c(nr,nc) 2x2=4, al ser mfrow y no mfcol, se enumeran por fila
## oma: márgenes externos entre gráfico y texto c(bottom, left, top, right)
## mar: márgenes internos entre área de trazado y gráfico c(bottom, left, top, right)
par(mfrow = c(2, 2),                 #2 filas por 2 columnas#
      oma = c(0, 0, 2, 0),           #márgenes externos# 
      mar = c(5.1, 4.1, 2.1, 2.1))   #márgenes internos#
 
##gráfico de secuencias
plot(t,y,ylab="Y",xlab="orden",type="l")
abline(h=mean(y), add=TRUE, col="blue")
abline(h=-2*sd(y), add=TRUE, col="blue", lty="dashed")
abline(h=2*sd(y), add=TRUE, col="blue", lty="dashed")
##gráfico de fases o lags
plot(y[-1],y[-n],xlab="Y[i-1]",ylab="Y[i]")
##gráfico de histograma
hist(y,main="",xlab="Y",ylab="f(Y)",freq=FALSE)
curve(dnorm(x,mean=mean(y),sd=sd(y)), add=TRUE, col="blue")
##gráfico QQPlot normal
qqnorm(y,main="",xlab="teórico",ylab="muestra")
qqline(y,col="blue")
mtext("4-Plot", line = 0.5, outer = TRUE) ##gráfico QQPlot normal

Veamos cómo funciona a través de ejemplos:

 

Ejemplo 1: números aleatorios distribuidos normalmente

Sería este el caso ideal donde se cumplen los 4 supuestos básicos. Nos servirá para presentar los gráfico en un estado "ideal".

Generamos primero una sucesión de 200 números distribuidos normalmente

y <- rnorm(200,0,1)

Luego ejecutamos el código R provisto en la sección anterior. El resultado debe ser similar a este:

Veamos  qué podemos aprender. Primero vayamos investigando cada uno de los 4 supuestos básicos:

• una medida de posición fija: si la hipótesis de cumple, el gráfico de secuencias (esquina superior izquierda) debe ser uniforme y sin tendencia. Puede apreciarse ello en el gráfico.
• una medida de dispersión fija: si la hipótesis de cumple, el gráfico de secuencias debe tener la misma amplitud sobre todo el eje x.
• aleatoriedad o ruido blanco: en caso de no haber un comportamiento correlacionado o determinista entre los datos. El gráfico de lags no debe mostrara patrones identificables, tal como se muestra en segundo gráfico arriba a la derecha.
• una distribución dada: en dicho caso, el histograma (esq inf izq) debe ajustar a la curva supuesta y el P-P Plot (esq inf derecha) mostrarse en línea recta. En este caso (y en el código R) la distribución subyacente es la Normal y se aprecia en los gráficos un buen ajuste a los datos.

Ejemplo 2: reflexiones de un láser

El siguiente ejemplo toma una muestra del año 1969 de 200 disparos de láser, computando su desvío en cada uno de ellos

  [1] -213 -564  -35  -15  141  115 -420 -360  203 -338 -431  194 -220 -513  154
 [16] -125 -559   92  -21 -579  -52   99 -543 -175  162 -457 -346  204 -300 -474
 [31]  164 -107 -572   -8   83 -541 -224  180 -420 -374  201 -236 -531   83   27
 [46] -564 -112  131 -507 -254  199 -311 -495  143  -46 -579  -90  136 -472 -338
 [61]  202 -287 -477  169 -124 -568   17   48 -568 -135  162 -430 -422  172  -74
 [76] -577  -13   92 -534 -243  194 -355 -465  156  -81 -578  -64  139 -449 -384
 [91]  193 -198 -538  110  -44 -577   -6   66 -552 -164  161 -460 -344  205 -281
[106] -504  134  -28 -576 -118  156 -437 -381  200 -220 -540   83   11 -568 -160
[121]  172 -414 -408  188 -125 -572  -32  139 -492 -321  205 -262 -504  142  -83
[136] -574    0   48 -571 -106  137 -501 -266  190 -391 -406  194 -186 -553   83
[151]  -13 -577  -49  103 -515 -280  201  300 -506  131  -45 -578  -80  138 -462
[166] -361  201 -211 -554   32   74 -533 -235  187 -372 -442  182 -147 -566   25
[181]   68 -535 -244  194 -351 -463  174 -125 -570   15   72 -550 -190  172 -424
[196] -385  198 -218 -536   96

  En este caso, los gráficos resultaron ser:

Verifiquemos las hipótesis:

• una medida de posición fija: el gráfico de secuencias no presenta tendencias
• una medida de dispersión fija: el gráfico de secuencias presenta aproximadamente la misma variación a través de todo su recorrido.
• aleatoriedad o ruido blanco: El gráfico de lags muestra un patrón muy específico, confirmando que la serie no es aleatoria.
• una distribución dada: siendo que se comprobó que no existe aleatoriedad, los gráficos de histograma y P-P plot no tienen un significado útil.

 Entonces concluimos que el modelo respuesta = constante + error no es aplicable para este problema y debemos buscar un mejor modelo.

Daré a continuación un modelo al cual ajustan mejor los datos, sin explicar en detalle cómo se llegó a esta nueva propuesta. Sólo diré que El gráfico de lags sugiere una función sinuidal u oscilante

siendo C una constante de nivel, alfa la amplitud, omega la frecuencia y tita la fase de la función seno.

Analizaremos los errores de medición y - y modelado para ver si el muevo modelo pudo capturar la parte determinista con éxito

t <- 1:length(y)
 
C <- -178.786
AMP <- -361.766 
FREC <- 0.302596 
PHASE <- 1.46536 

e <- C + AMP*sin(2*pi*FREC*t+PHASE) - y

En la última línea hemos restado la parte determinista a los datos reales. Los gráficos de residuales e son:

Verifiquemos nuevamente las hipótesis sobre los residuos o parte aleatoria:

• una medida de posición fija: el gráfico de secuencias no presenta tendencias
• una medida de dispersión fija: el gráfico de secuencias presenta saltos durante su recorrido, indicando que aún es mejorable el modelo en este aspecto.
• aleatoriedad o ruido blanco: El gráfico de lags no muestra patrones.
• una distribución dada: los gráficos de histograma y p-p Plot muestran una distribución aproximadamente normal centrada en cero (insesgado). Quizás haya una leve asimetría hacia la derecha por la caída izquierda del p-p Plot.

El nuevo modelo cumple bastante bien con los cuatro supuestos subyacentes. Aunque aún puede mejorarse. por ejemplo,

• los gráficos muestran la existencia de 4 o 5 valores extraños, 
• hay una asimetría muy leve que se puede tratar,
• también hay oscilaciones en la amplitud vistas en el diagrama de secuencias que se podrían capturar volviendo a cambiar la constante C por un modelo más complejo.

Así y todo, los resultados son lo suficientemente satisfactorios para realizar inferencias. Así, adoptamos el nuevo modelo.

Conclusiones

La herramienta 4-plot desarrollada es muy efectiva para validar los cuatro supuestos subyacentes. Incluso nos aporta más información para resolver los hipótesis no cumplidas y pensar en un nuevo modelo.

 Este es el listado de preguntas que se pueden responder con este set de gráficos:

1. Está el proceso en control; es estable y predecible?
2. El proceso muestra tendencias en su posición?
3. El proceso muestra tendencias en su variación?
4. Los datos son aleatorios?
5. Las observaciones están relacionadas con las observaciones adyacentes?
6. Si los datos son una serie de tiempo, se trata de un ruido blanco?
7. Si los datos son una serie de tiempo y no se trata de un ruido blanco, es un modelo autoregresivo o sinusoidal?
8. Los datos o residuos siguen una distribución normal?
9. Si no es normal, qué forma tiene la distribución?
10. Es el modelo  Y = C + e  válido y suficiente?
11. Si el modelo básico es insuficiente, hay posibilidad de mejorarlo?
12. Es la varianza muestra un buen estimador de variación?
13. Es la media muestral un buen estimador de posición?
    
14. Si no lo fuera, cuál sería mejor?
15. Existen valores extraños (outliers)? 

El Análisis Exploratorio de Datos (EDA)

EDA es un enfoque de análisis de datos. ¿Qué otros enfoques existen? Podemos destacar tres enfoques populares:

1. Enfoque Clásico
2. Enfoque Exploratorio (EDA)
3. Enfoque Bayesiano

Estos tres paradigmas tienen el objetivo común de abordar un problema de manera científica. Sus diferencias radican en el orden de los pasos a seguir:

1. Enfoque Clásico:  secuencia  Problema>Datos>Modelo>Análisis>Conclusiones.
   La revisión de los datos impone un modelo (normalidad, linearidad, etc...), y el análisis, estimación y testeo subsiguiente está concentrado en los parámetros de dicho modelo seleccionado.
2. Enfoque Exploratorio (EDA):  secuencia  Problema>Modelo>Análisis>Datos>Conclusiones.
   Los datos no imponen un modelo; en su lugar se prueban modelos en base a análisis y se selecciona el más apropiado.
3. Enfoque Bayesiano: :  secuencia  Problema>Datos>hipótesis a Priori>Modelo>Análisis>Conclusiones.
   El analista intenta incorporar su experiencia con funciones independientes de los datos recolectados que las asume como ciertas, y luego se plantean modelos con los datos para arribar a conclusiones.

Enfocándonos en los primeros dos enfoques, apreciamos las siguientes diferencias en los siguientes temas:

El modelo

El enfoque clásico impone modelos (tanto determinísticos como probabilísticos) sobre los datos observados. El enfoque se encuentra en el modelo.
El enfoque Exploratorio (EDA) no impone un modelo a los datos, sino que son estos últimos los que sugieren los modelos. El enfoque se encuentra en los datos.

Herramientas utilizadas

Sin ser listas exhaustivas y determinantes de ambos enfoques, podemos decir que los enfoques utilizan con mayor frecuencia las siguientes herramientas:

• El enfoque clásico: Herramientas cuantitativas para validar los modelos a probar: ANOVA, Pruebas t, F y Chi cuadrado, 
• El enfoque Exploratorio (EDA) utiliza más herramientas gráficas: dispersión, histogramas, cajas, residuales, etc.

Rigor científico

Sin duda el enfoque clásico responde al paradigma clásico del método científico, por lo que es objetivo y formal y riguroso.

El análisis EDA no tiene los mismos niveles de rigurosidad y formalismo. Aborda los problemas inferencialmente y sus conclusiones son sugestivas e indicativas, muy en base a la experiencia del analista.

Tratamiento de los datos

Las estimaciones del enfoque clásico se realizan generalmente comparando laa calidad del modelo evaluando algunas propiedades de los datos (media, mediana, desvío estándar, nivel de asimetría, percentiles de cola, etc). Por lo tanto hay cierta pérdida de información.

El enfoque EDA, por otro lado, hace uso de toda la información disponible y no hay pérdida de información.

Hipótesis planteadas

El enfoque toma como punto de partida un modelo hipotético a contrastar con los datos de la muestra. Pero estos contrastes se basan a su vez en distribuciones teóricas que asumen hipótesis sobre los datos que pueden no llegar a cumplirse. Entonces, la validación del método científico termina valiendose de tests que tienen hipótesis subyacentes no necesariamente probadas, y las conclusiones pueden terminar siendo dudosas.

En EDA, varios de los tests de comprobación son no paramétricos, y otros meramente visuales.

Finalmente decidirá el analista qupe modelo serpa el mejor, independientemente del enfoque utilizado.

Como resumen, podemos decir que le enfoque clásico es una reducción numérica de un set de dato. Es así un enfoque pasivo. Su propósito es arribar con sus modelos a unas pocas medidas clave (por ej, media y desvío estándar). En contraste, EDA tiene el objetivo de explorar los procesos científicos que se hallan detrás de los datos. Utiliza los datos como una "ventana" que mira los procesos subyacentes que los han generado.