Estimar un modelo de regresión lineal univariada con muestreo de Gibbs (estadística bayesiana)

 1.1.           Simulando datos

Primero simulamos los datos, asumiendo un beta1 real de 0.5. Puntualmente se generan dos variables aleatorias x e y, con una relación de 0.5, más un término de error. El código se desarrolló en el marco de un curso sobre “Modelos para el análisis de políticas macroeconómicas”. En esta sección se puede indicar que BETA y SIG2 son los parámetros a estimar en el procedimiento. Como créditos se aclara que el código original en Matlab fue elaborado por el profesor Daniel Parra.

rm(list=ls())
set.seed(1)
 

# Generar datos
# Modelo
# y_{t} = BETA * x_{t} + v_{t}
# v_{t} ~ N(0,SIG2)
 
BETA <- 0.5
SIG2 <- 0.006 # Varianza del error
k <- 1  # Número de variables explicativas
t <- 500  # Longitud de la muestra
e1 <- rnorm(t, mean = 0, sd = sqrt(SIG2))  # Generar errores distribuidos N(0, SIG2)
x <- rnorm(t, mean = 0, sd = 1)  # Variable explicativa: extracciones i.i.d. de N(0,1)
 
y <- x * BETA + e1 # Calcular la variable dependiente

1.2. Enfoque frecuentista

Utiizando el enfoque tradicional del modelo de regresión podemos recuperar, a partir de las variables simuladas, los coeficientes a estimar –tenga pendiente que estos procesos están automatizados en R mediante el comando lm tradicional (lm(y~x))-:

 # Estimación de OLS (para comparación)
bols <- solve(t(x) %*% x) %*% t(x) %*% y
sigols <- t(y - x * bols) %*% (y - x * bols) / (t - k)

Observando los resultados se puede apreciar como el enfoque logra aproximar bien los datos.

        coef  real         mco
1      beta 0.500 0.496877122
2 sig.error 0.006 0.006136154

1.3. Enfoque bayesiano

En primer lugar, es importante recordar que el enfoque bayesiano está interesado en estimar una posterior que resulta del producto de la prior y una función de verosimilitud. La prior se estima o supone basado en el conocimiento previo del analista, mientras que la mle debe estimarse mediante algún algoritmo. En el siguiente ejemplo se utiliza el algoritmo de Gibbs. En esta línea se asume una distribución normal del beta y una gamma del error del modelo y se asumen parámetros de estas definiciones, como la media del set de parámetros que se estimaran mediante la estimación bayesiana.

b0 <- 0.2  # Media de la distribución a priori para beta
P0 <- 10   # Varianza de la distribución a priori para beta
 
t0 <- 0  # Parámetro de forma a priori para sigma
R0 <- 0  # Parámetro de escala a priori para sigma: cero implica un prior "plano"
# es decir, sin información a priori

Posteriormente se establecen los parámetros para el rastreador de Gibbs. La idea de estos es tomar x cantidad de simulaciones, pero omitir las primeras dado que estas pueden estar incididas por las condiciones de inicio establecidas en el algoritmo. Puntualmente, se generan 100,000 simulaciones (n), pero de estas, posteriormente se desechan (Período de quema: las extracciones se descartan) las primeras 90,000. Por tanto, es sobre las ultimas 10 extracciones sobre las que se basara la inferencia. Note el condicional if que se encuentra en el bucle (for) de abajo, este guarda solo las últimas 10,000 iteracciones.

nburn <- 90000  # Período de quema: las extracciones se descartan
nkeep <- 10000  # Número de extracciones en las que se basará la inferencia
n <- nkeep + nburn  # Número total de iteraciones del muestreador de Gibbs

Luego se generan los objetos donde se guardará la distribución de betas y errores (parámetros) que serán obtenidas a partir de Gibbs.

sig2 <- 1  # Valor inicial para el muestreador de Gibbs
 
Beta <- matrix(0, nkeep, k)
Sig2 <- numeric(nkeep)

Finalmente se utiliza el algoritmo de Gibbs para obtener la distribución de parámetros, y en función de la media de esta distribución el valor estimado. Puntualmente, aquí se define un algoritmo iterativo de las simulaciones, desde 1 hasta n. En el primer paso se utiliza el prior para estimar el beta, más una corrección conocida como factor de incertidumbre. Luego con este beta, estimado en la primera iteración se utiliza para estimar la varianza. En tal sentido, sig2 guarda un posible candidato aleatorio para la varianza. Como la varianza no puede ser negativa se usa una distribución gamma. posterior usa la información del prior, con la información que recibo de los datos.

for (iter in 1:n) {
 
  # Generar una extracción para beta condicional a sig2
  P1 <- solve(solve(P0) + (1/sig2) * t(x) %*% x)
  b1 <- P1 %*% (solve(P0) %*% b0 + (1/sig2) * t(x) %*% y)
  beta <- b1 + chol(P1) %*% rnorm(k, mean = 0, sd = 1)
 
  # Generar una extracción para sig2 condicional a beta
  t1 <- t0 + t
  R1 <- R0 + t(y - x * beta) %*% (y - x * beta)
  shpe <- t1
  scle <- 1 / R1
  sig2inv <- rgamma(1, shape = shpe, rate = scle)
  sig2 <- 1 / sig2inv
 
  if (iter > nburn) {
    Beta[iter - nburn,] <- beta
    Sig2[iter - nburn] <- sig2
  }
}

Estos pasos están relacionados con la actualización de parámetros en el muestreo de Gibbs.

 -          P0 es la matriz de varianza-covarianza a priori para el parámetro beta.

-          sig2 es la varianza del error en el modelo.

-        P1 es la inversa de la suma de la inversa de la matriz prior P0 y el producto (1/sig2) * (x'*x). se utiliza para actualizar la matriz de covarianza-condicional de la distribución posterior de los parámetros beta dados los valores actuales de otros parámetros y los datos observados. La actualización se realiza en la fase de muestreo de Gibbs, donde se obtienen nuevas muestras de los parámetros del modelo. La matriz resultante P1 es la varianza-covarianza condicional de la distribución posterior de beta dada la varianza del error sig2 y los datos observados.

 - b0 es el valor medio a priori para el parámetro beta. Esta línea actualiza el vector de medias condicionales de la distribución posterior de beta dada la varianza del error sig2 y los datos observados.

-y es el vector de observaciones de la variable dependiente.

-          (1/sig2) * (x'*y) es la contribución de los datos observados al término de actualización.

En resumen, estas líneas de código están realizando la actualización de parámetros en el contexto del muestreo de Gibbs para obtener una nueva muestra de la distribución conjunta de beta condicional en los valores actuales de los demás parámetros y los datos observados.

 1.4.           Resultados

En primer lugar, se observa que pese a tener un prior algo alejado del verdadero valor del beta, el algoritmo es capaz de empujar sus resultados hacia el valor verdadero.

hist(Beta, main = expression(beta), col = "lightblue", xlim = c(0.32, 0.6))
abline(v = BETA, col = "red", lwd = 2)
abline(v = b0, col = "black", lwd = 2)
legend("topright", legend = c("Estimado", "Verdadero", "Prior"), col = c("blue", "red", "black"), lty = 1, cex = 0.8)

Posteriormente, podemos comparar la distribución obtenida, vs. El valor verdadero del parámetro, y las estimaciones mco.

 hist(Beta, main = expression(beta), col = "lightblue")
abline(v = BETA, col = "red", lwd = 2)
abline(v = mean(Beta), col = "green", lwd = 3)
abline(v = bols, col = "black", lwd = 2)
legend("topright", legend = c("Estimado", "Verdadero", "Media posterior", "MCO"), col = c("blue", "red", "green"), lty = 1, cex = 0.8)

 Y las tablas de resultados:

        coef  real         mco       gibbs
1      beta 0.500 0.499524587 0.499512621
2 sig.error 0.006 0.006116056 0.006568336

 

Curso de Análisis estadístico en R

La siguiente entrada contiene un curso de Métodos Estadístico de Investigación en R, usando los fundamentos del tidyverso para la manipulación y comunicación de datos. 

Parte I. Análisis exploratorio en R

Clase 1.0. pdf. video. Introducción a R

Clase 1.1. pdf. video. Estructura de datos (tomado del curso de Derek Corcoran)

Clase 2.1. pdf. video. Análisis exploratorio: posición central y dispersión

Clase 2.2. pdf. video. Análisis exploratorio: posición (quintiles), distribuciones y asociación

Clase 2.3. pdf. video. Análisis exploratorio: datos perdidos y valores atípicos

Clase 2.4. pdf. video. Tablas en R

Parte II. Manipulación y comunicación de datos en R

Clase 3.1. pdf. video. Manipulando datos usando dplyr I (select, filter, mutate)

Clase 3.2. pdf. video. Manipulando datos usando dplyr I (summarize, by_group, accros, if, at)

Clase 3.3. pdf. video. Ejercicios usando dplyr

Clase 4.1. pdf. video. Comunicación de datos con ggplot2 y rmakdowm

Clase 4.2. pdf. video. Ejercicios con ggplot2

Clase 5.1. pdf. video. Gestión de datos relacionales en R (join, reshape, long, wide)

Parte III. Inferencia estadística y nociones de ciencia de datos en R

Clase 6.1. pdf. video. Repaso de inferencia estadística

Clase 6.2. pdf. video. Análisis de independencia estadística en R

Clase 6.3. pdf. video. Ejercicios aplicados a independencia y datos relacionales 

Clase 7.1. pdf. video. Métodos de clasificación en R (Statistical Learning)

Clase 8.1. pdf. video. Modelos de regresión y aprendizaje supervisado

Parte IV. Análisis de series temporales

Clase 9.1.   pdf. video. Exploración y agregación

Clase 10.1. pdf. video. Estacionariedad y modelos ARIMA

Clase 10.2. pdf. video. Ejericicios

Clase 11.1. pdf. video. Pronóstico

Clase 12.1. pdf. video. Modelos de volatilidad (GARCH)

Clase 12.2. pdf. video. Markov GARCH

Test de independencia estadística en R

 1. Análisis preliminar

En la siguiente entrada se muestran algunos ejemplos de cómo realizar test de independencia estadística en R, para responder a uno de los ejercicios de la asignatura de Seminario de Investigación (clase 6 de la parte de R). Los ejercicios se realizan con la base utown del libro Principles of Econometrics 4e, by Carter Hill, William Griffiths, and Guay Lim. Los mismos se encuentran disponibles en el paquete PoEdata.

library(PoEdata)

library(tidyverse)
             
data(utown)        
head(utown,5)        

##     price  sqft age utown pool fplace
## 1 205.452 23.46   6     0    0      1
## 2 185.328 20.03   5     0    0      1
## 3 248.422 27.77   6     0    0      0
## 4 154.690 20.17   1     0    0      0
## 5 221.801 26.45   0     0    0      1

1. Análisis de independencia

1.1. Usando ggplo2, presente un histograma de precios y muestre un summary (min, max,median, q1, q5) condicionado a si la casa tiene o no piscina (pool).

 theme_set(theme_minimal())

#histograma
utown%>%
ggplot(aes(price)) +
geom_histogram(col = "black",fill = "darkblue", alpha = 0.6)

# Summary de precios por grupo
tapply(utown$price, utown$pool, summary)

## $`0`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   134.3   214.8   244.3   246.5   276.8   345.2
##
## $`1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   160.3   219.5   249.1   252.0   285.4   340.7

1.2. Obtenga un boxplot donde relacione el tener piscina o no, con el precio de la casa, y establezca la independencia entre las variables.

utown%>%
  ggplot(aes(y=price, x=factor(pool), fill=factor(pool))) +
  geom_boxplot(col = "black", alpha = 0.6)

1.3. A partir de la comparación de promedios, que conclusión podemos obtener y que limitaciones puede enfrentar esta conclusión sobre la posible obtención de un efecto causal.

Al considerar los intervalos del gráfico anterior se puede observar que los promedios son relativamente parecidos, sin embargo, esta comparación no se puede interpretar como un efecto causal, dado que esta simple comparación omite otros factores relevantes que inciden sobre el precio de las casas.

1.4. Presente una tabla de frecuencia, donde los porcentajes se presenten en relación a las tablas, interprete sus resultados y posteriormente realice el test chi2 para determinar la posible independencia entre las variables.

Necesitamos una tabla de frecuencia de precio con pool, pero price es una variable continua, por lo que, debemos discretizarla, antes de realizar la tabla de frecuencia. Una alternativa es crear rangos de precios (cut(price, breaks = 5)), utilizar los quintiles (ejemplo debajo), o separarla en grupos alrededor de alguna medida como la media (price>mean(price)).

# discretizar la variable precio
utown <- utown %>%
  mutate(quintiprecio = ntile(price, 5))

utown %>%
select(quintiprecio, pool) %>%
table() %>%
prop.table(.,2) 

##             pool
## quintiprecio         0         1
##            1 0.2085427 0.1666667
##            2 0.2010050 0.1960784
##            3 0.1959799 0.2156863
##            4 0.2022613 0.1911765
##            5 0.1922111 0.2303922

Los valores son bastante diferentes, eso es evidencia de independencia, pero se debe realizar el test chi2 de independencia entre variables discretas (p.valor=0.5508, no se rechaza ho de independencia):

# discretizar la variable precio
utown %>%
select(quintiprecio, pool) %>%
table() %>%
chisq.test()

##
##  Pearson's Chi-squared test
##
## data:  .
## X-squared = 3.0422, df = 4, p-value = 0.5508

2. Análisis de independencia: Correlación

2.1. Sabiendo que tenemos una medida del size del terreno de la casa (sqft), crear los quintiles de esta variable y muestre una tabla resumen donde indique el precio max, min, mediano, q1 y q3 del precio de las casas para cada quintil.

utown %>%
  mutate(quintiprecio = ntile(price, 5)) %>%
  select(price,quintiprecio) %>%
  group_by(quintiprecio) %>%
  summarise(min = min(price),
    q1 = quantile(price,0.2),
    mean = mean(price),
    q3 = quantile(price,0.6),
    max = max(price))

## # A tibble: 5 x 6
##   quintiprecio   min    q1  mean    q3   max
##          <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1            1  134.  178.  190.  197.  208.
## 2            2  209.  214.  222.  225.  235.
## 3            3  235.  239.  246.  249.  258.
## 4            4  258.  263.  272.  274.  287.
## 5            5  287.  293.  309.  311.  345.

2.2. Obtenga un scatter donde relacione el size de la casa con el precio, interprete los resultados, coloreando los puntos según la casa tenga o no piscina (pool).

utown %>%
ggplot(aes(sqft, price, color =factor(pool)))+
geom_point()+
geom_smooth()

2.3. Realice el test de correlación de pearson y determine si rechaza o no ho.

Se rechaza ho.

cor.test(utown$price,utown$sqft, method = "pearson")

##
##  Pearson's product-moment correlation
##
## data:  utown$price and utown$sqft
## t = 23.367, df = 998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5530750 0.6333235
## sample estimates:
##       cor
## 0.5946785

2.4. Utilice el qqplot para realizar un test de normalidad sobre el price y sqft.

car::qqPlot(utown$price)

2.5. En base al test de normalidad, que podemos concluir sobre las medidas de resumen y el análisis de correlación realizado.

Como para las colas se rechaza que la variable se comporte como una normal, se entiende que es preferible usar alguna estimación no parámetrica para obtener el coeficiente de correlación, tal como el método de Spearman.

3. Test de media

3.1. Realice un test de media para testear si las casas con chimeneas son estadísticamente más grandes que aquellas casas que no tienen.

t.test(sqft~fplace,
       data=utown,
       var.equal=F,
       alternative = "greater",
       conf.level=0.95)

##
##  Welch Two Sample t-test
##
## data:  sqft by fplace
## t = -3.1236, df = 989.19, p-value = 0.9991
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -0.8778771        Inf
## sample estimates:
## mean in group 0 mean in group 1
##        24.91187        25.48674

3.2. Realice boxplot donde compare el precio de la casa según el size esté por encima o por debajo de la mediana.

utown%>%
  mutate(sqftmedian = sqft>median(sqft)) %>%
  ggplot(aes(y=price, x=factor(sqftmedian), fill=factor(sqftmedian))) +
  geom_boxplot(col = "black", alpha = 0.6)

3.3. Realice un test de media para testear si las casas con edad por encima de la mediana, tienen a tener un precio menor al resto de casa.

utown <- utown %>%
  mutate(medianEdad = as.numeric(age>median(age)))

t.test(price~medianEdad,
       data=utown,
       var.equal=F,
       alternative = "less",
       conf.level=0.95)

##
##  Welch Two Sample t-test
##
## data:  price by medianEdad
## t = 1.2334, df = 994.12, p-value = 0.8911
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##      -Inf 7.685081
## sample estimates:
## mean in group 0 mean in group 1
##        249.2916        246.0001

Una demostración del teorema del límite central en R

"El Teorema Central del Límite indica que, en condiciones muy generales, la distribución de la suma de variables aleatorias tiende a una distribución normal cuando la cantidad de variables es muy grande. Es decir, garantiza una distribución normal cuando n es suficientemente grande" (Gómez y Benlloch, nd).

En el siguiente ejemplo se muestra como la distribución de medias muéstrales tiende hacia una distribución normal, aunque las muestras no procedan de una distribución normal. Esta aproximación será más cercana en la medida en que aumente la cantidad de muestras usadas. 

Usando la data Wage1 del texto de Econometría Wooldridge, se toma una muestra de tamaño 70 sobre la variable wage, para obtener el estadístico de la media a partir de esa muestra.

library(tidyverse)

library(wooldridge)

 attach(wage1)

# semilla

set.seed(1)

# toma una muestra n=70, con reemplazo.

sample(wage, size = 70, replace = TRUE) %>%

  mean()

[1] 5.252571

Obtenemos la distribución en el muestreo del estadístico al repetir el proceso 10 veces usando la función replicate, es decir, repetimos el procedimiento anterior 10 veces, obteniendo 10 medias (una para cada muestra), que al asignarle probabilidades se obtiene una distribución en el muestreo del estadístico, sobre las que dibujamos un histograma.

 dist_mues_media1 <- replicate(10, sample(wage, size = 70, replace = TRUE) %>% mean())

 dist_mues_media1 %>%  data.frame() %>%

 ggplot(aes(.)) +

  geom_histogram(bins=10)

Ahora repetimos el ejercicio anterior, pero obteniendo 100 muestras de igual tamaño, para verificar como la distribución en el muestreo converge a una distribución normal.

 dist_mues_media2 <- replicate(100, sample(wage, size = 70, replace = TRUE) %>% mean())

 dist_mues_media2 %>%  data.frame() %>%

  ggplot(aes(.)) +

  geom_histogram(bins=10)

Finalmente, se verifica como el promedio (centro de la distribución anterior) se va acercando al promedio de la población global (o "parámetro poblacional") en la medida en que crece la cantidad de muestras usadas para crear la distribución en el muestreo.

> mean(wage)

[1] 5.896103

> mean(dist_mues_media1)

[1] 6.015486

> mean(dist_mues_media2)

[1] 5.853953

Referencias

- Benlloch, M.; Gómez, M. (nd). Utilización Práctica del Teorema Central del Límite. Universidad Politécnica de Valencia.

- Data Camp (2020). Introduction to statistics in R. 

Curso de inferencia estadística

La siguiente entrada presenta las clases de un curso de 20 horas de fundamentos de inferencia estadística, impartido a estudiantes de economía de la Universidad Autónoma de Santo Domingo (UASD) en noviembre de 2019. El programa del curso pretende ser un primer reforzamiento de las estadísticas universitarias orientado a continuar estudios más avanzados así como sentar bases sobre los aspectos estadísticos básicos usados en la materia de econometría.

Clase 1. Teoría de probabilidad. [Semana 1] se explica el papel de la inferencia estadística y la teoría de la probabilidad en la obtención de estimadores, además de nociones de la teoría de medida (espacios medibles), álgebra de conjuntos y análisis de independencia basados en correlaciones, productos de probabilidad y análisis condicional (de media y probabilidad).

Clase 2. Variables y vectores aleatorios. [Semana 2 y 3]  En esta segunda clase se explica que es una variable aleatoria y sus funciones, tanto de distribución como de probabilidad. Se muestra como obtener probabilidades a partir de estas funciones y calcular los momentos de una variable aleatoria. Además, se estudian probabilidades condicionales a partir de funciones de densidad conjunta, se estudian marginales, momentos condicionales, matrices de covarianzas y correlaciones. [ejemplos de clases: 2.2 | 2.3]

Clase 3. Estimadores puntuales. [Semana 4 y 5] La tercera semana trata sobre muestreo y distribución en el muestreo de los estimadores con métodos computacionales (bootstrap, ejemplo en R), además se tratan las propiedades de los estimadores, conjuntamente con métodos de estimación puntual como máxima verosimilitud y método de los momentos.

Clase 4. Test de hipótesis.

En construcción