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

 

Usar la función summarise+apply en dplyr

En la siguiente entrada, se simula una data, que contienen diferentes tipos de series, con el objetivo de mostrar, como al combinar las funciones summarise (del paquete dplyr) y la función apply (base de R), podemos obtener resúmenes estadísticos de un conjunto de variables, incluyendo test formales como el Dickey Fuller, o estadísticos como el coeficiente de autocorrelación. Puntualmente, tentemos una base de datos con 30 variables y deseamos obtener una tabla donde a cada variable le calculemos el promedio y la desviación estándar.

library(dplyr)
library(moments)
library(tseries)

Generamos de forma aleatoria una data llamada my_data, esta es un ejemplo de datos aleatorios para ilustrar el proceso.

set.seed(123)
my_data <- data.frame(
  serie1 = rnorm(100),
  serie2 = rnorm(100),
  fecha = as.Date("2024-01-01") + 0:99,
  serie3 = rnorm(100),
  categoria = sample(c("A", "B", "C"), 100, replace = TRUE)
)
 
          serie1      serie2      fecha      serie3 categoria
1   -0.560475647 -0.71040656 2024-01-01  2.19881035         B
2   -0.230177489  0.25688371 2024-01-02  1.31241298         A
3    1.558708314 -0.24669188 2024-01-03 -0.26514506         C
4    0.070508391 -0.34754260 2024-01-04  0.54319406         C
5    0.129287735 -0.95161857 2024-01-05 -0.41433995         C
…

Ahora creamos un filtro para seleccionar únicamente las series numéricas dentro de mi base de datos, dato que puede contener fechas y caracteres. Es decir, las series 1,2 y 3.

numeric_columns <- my_data %>%
  select_if(is.numeric) %>%
  colnames()
 
umeric_columns
[1] "serie1" "serie2" "serie3"

Ahora, usamos select para recuperar las variables numéricas y luego a cada columna de aplicamos la opción deseada combinando summarise y apply como se muestra en los siguientes ejemplos.

my_data %>%
  dplyr::select(all_of(numeric_columns)) %>%
  summarise(
    varianza = apply(., 2, var),
    media = apply(., 2, mean)
  )
 
varianza       media
1 0.8332328  0.09040591
2 0.9350631 -0.10754680
3 0.9022700  0.12046511

Finalmente, calculamos para cada columna (numérica, porque solo va a trabajar con las condiciones impuestas por alguna operación lógica), se estimarán los estadísticos indicados, y los organizara en un cuadro. Usamos la función select del paquete dplyr, esta se indica el nombre del paquete, dado que puede coincidir con funciones de otros paquetes y generar errores. Además, usamos mutate, para colocar el nombre en cada fila, permitiendo la lectura de los datos.

Error in select(., all_of(numeric_columns)) :
  unused argument (all_of(numeric_columns))

Posteriormente, se solicita un resumen de la data, pero realizada a cada una de las columnas usando apply. En el caso de funciones menos directas, como el p valor del test DF, se inserta una fusión que permita extraerlo.  

stats <- my_data %>%
  dplyr::select(all_of(numeric_columns)) %>%
  summarise(
    varianza = apply(., 2, var),
    asimetria = apply(., 2, skewness),
    curtosis = apply(., 2, kurtosis),
    autocorr_lag1 = apply(., 2, function(x) {
      acf(x, plot = FALSE)$acf[2]
    }),
    p_valor_dickey_fuller = apply(., 2, function(x) {
      adf.test(x)$p.value
    })
  )
 
# Imprimir las estadísticas
print(stats)

Creando variables por grupos en dplyr (group_by + mutate)

 Simulemos una base de hogares, donde se identifica el hogar, el sexo (1 mujer) y provincia y edad para cada miembro.
 
# Definir la lista de relaciones familiares por hogar
hogar <- c(1,1,1,1,2,2,2,3,4,4,4)
 prov <- c("a","a","a","a","c","c","c","a","b","b","b")
   id <- 1:length(hogar)
mujer <- c(1,0,1,1,0,0,2,0,1,0,0)
 edad <- c(28,29,10,8,33,27,13,24,32,11,6)
 
hogar_data <- data.frame(id, hogar, prov, mujer, edad)
hogar_data
  id hogar prov mujer edad
1   1     1    a     1   28
2   2     1    a     0   29
3   3     1    a     1   10
4   4     1    a     1    8
5   5     2    c     0   33
6   6     2    c     0   27
7   7     2    c     2   13
8   8     3    a     0   24
9   9     4    b     1   32
10 10     4    b     0   11
11 11     4    b     0    6
 
Si combinanos group_by con summarise, obtendríamos un collapse de la base de datos, con la edad promedio del grupo indicado en group_by. Es decir, pasamos a tener tantas observaciones como hogares en nuestra base de datos. Noten que en este caso solo quedan aquellas variables a las que le hemos indicado la forma de agregación.
 
hogar_data |>
  group_by(hogar) |>
  summarise(mean(edad))
 
# A tibble: 4 x 2
  hogar `mean(edad)`
  <dbl>        <dbl>
1     1         18.8
2     2         24.3
3     3         24 
4     4         16.3
 
Ahora, en lugar de obtener el resumen anterior, quisiéramos agregar una variable que indique la edad promedio por hogar. Podemos combinar group_by + mutate. Ahora agregamos una variable con la edad promedio de cada hogar.
 
hogar_data |>
  group_by(hogar) |>
  mutate(edad_mhogar = mean(edad))
 
# A tibble: 11 x 6
# Groups:   hogar [4]
      id hogar prov  mujer  edad edad_mhogar
   <int> <dbl> <chr> <dbl> <dbl>       <dbl>
 1     1     1 a         1    28        18.8
 2     2     1 a         0    29        18.8
 3     3     1 a         1    10        18.8
 4     4     1 a         1     8        18.8
 5     5     2 c         0    33        24.3
 6     6     2 c         0    27        24.3
 7     7     2 c         2    13        24.3
 8     8     3 a         0    24        24 
 9     9     4 b         1    32        16.3
10    10     4 b         0    11        16.3
11    11     4 b         0     6        16.3
 
Contar número de menores de edad en cada hogar. Contar la cantidad de menores de edad en cada hogar.
 
hogar_data |>
  group_by(hogar) |>
  mutate(num_menores = sum(edad < 18))
 
# A tibble: 11 x 6
# Groups:   hogar [4]
      id hogar prov  mujer  edad num_menores
   <int> <dbl> <chr> <dbl> <dbl>       <int>
 1     1     1 a         1    28           2
 2     2     1 a         0    29           2
 3     3     1 a         1    10           2
 4     4     1 a         1     8           2
 5     5     2 c         0    33           1
 6     6     2 c         0    27           1
 7     7     2 c         2    13           1
 8     8     3 a         0    24           0
 9     9     4 b         1    32           2
10    10     4 b         0    11           2
11    11     4 b         0     6           2
 
Contar el numero de menores de edad que son mujeres.
 
hogar_data |>
  group_by(hogar) |>
  mutate(num_menores = sum(edad < 18 & mujer ==1))
 
# A tibble: 11 x 6
# Groups:   hogar [4]
      id hogar prov  mujer  edad num_menores
   <int> <dbl> <chr> <dbl> <dbl>       <int>
 1     1     1 a         1    28           2
 2     2     1 a         0    29           2
 3     3     1 a         1    10           2
 4     4     1 a         1     8           2
 5     5     2 c         0    33           0
 6     6     2 c         0    27           0
 7     7     2 c         2    13           0
 8     8     3 a         0    24           0
 9     9     4 b         1    32           0
10    10     4 b         0    11           0
11    11     4 b         0     6           0

20 ejemplos de Indexación de un vector en R

Entenderemos indexación acceder a una parte de un objeto. Recordemos que, en R, tenemos objetos atómicos (de un solo tipo de clases, por ejemplo, solo números, o solo texto) y recursivos. Estos combinaban diferentes objetos con formas distintas, por ejemplo, texto y números. Según chatGPT “en R, "indexar" se refiere a acceder a elementos específicos de un vector, matriz, lista o data frame utilizando su posición numérica o su nombre. Indexar te permite extraer o modificar valores individuales o conjuntos de valores de una estructura de datos”.

 INDEXACIÓN NUMÉRICA

En el caso de la indexación numérica, accedemos a parte de ese vector indicando a cuál posición deseamos acceder. Debajo colocamos diferentes ejemplos usando el vector edad.

Ej. 1: La primera posición, colocamos corchetes [] posterior al vector que contienen los elementos al cual queremos acceder (edad en este caso).

edad <- c(17,22,31,14,12)
edad[1]
[1] 17

Ej. 2: Si el número se coloca en negativo, indica que no queremos acceder a esa parte del vector. Regresa todo, menos la posición indicada con menos.

edad[-1]
[1] 22 31 14 12

Ej. 3: Igual podemos acceder a múltiples posiciones de un vector, colocando un vector de números entre los corchetes. Si repites el 1 varias veces te regresara la posición 1 cuantas veces como repitas el 1.

> edad[c(1,4)]
[1] 17 14

Ej. 4: Acceder a la última posición de un vector usando la función length.

edad[length(edad)]
[1] 12

Ej. 5: Igual podemos usar el operador dos puntos para generar una secuencia. Sí le indicamos esta secuencia a R, entonces regresa las posiciones del vector que están dentro de la misma.

edad[1:3]
[1] 17 22 31

Ej. 6: imaginemos ahora que queremos acceder a las dos ultimas posiciones del vector, podemos indicar la longitud de la misma:

edad[c(4,5)]
[1] 14 12

Ej. 7: No obstante, el ejemplo anterior deja de funcionar si agregamos datos al vector anterior:

edad <- c(edad, 14,24)
[1] 17 22 31 14 12 14 24

edad[(length(edad)-1):length(edad)]
[1] 14 24

Ej. 8: Entonces, si queremos eliminar las dos ultimas posiciones de un vector, colocamos el signo de menos delante de la secuencia.

edad[-c((length(edad)-1):length(edad))]
[1] 17 22 31 14 12

INDEXACIÓN LÓGICA 

En este caso colocamos un vector lógico entre corchetes, en lugar de números. Este permite acceder a elementos que cumplen con alguna condición, sin importar donde estén colocados los mismos.

Ej. 9: Llamar a los menores de edad.

 edad[edad<18]
[1] 17 14 12 14

Ej. 10: Llamar a los que tienen edades entre 5-15 años.

edad[edad>=5 & edad<=15]
[1] 14 12 14

Ej. 11: Podemos indexar a partir de la información de un segundo vector. Supongamos tenemos el vector mujer que indica el sexo de las personas. Entonces podemos recuperar las edades de las mujeres.

mujer<-c(1,0,1,0,1,0,0)
edad[mujer==1]
[1] 17 31 12

Tener pendiente que, si usas esta variable directamente para indexar el vector, este te regresara la primera posición tantas veces como se repita en el vector numérico mujer. Recuerde que este es un vector numérico.

edad[c(1,1,1)]
[1] 17 17 17

edad[mujer]
[1] 17 17 17

Por tanto, si se quiere usar directamente, se debe definir como un vector lógico.

mujer1<-c(T,F,T,F,T,F,F)
edad[mujer1]

OTROS EJEMPLOS

Ej. 12: Indexar utilizando nombres asignados a los datos.

nombres <- c("Ana", "Pedro", "Juana", "Johan", "Lia", NA, "Julio")
names(edad) <- nombres
edad
  Ana Pedro Juana Johan   Lia  <NA> Julio
   17    22    31    14    12    14    24
 
edad["Juana"]
Juana
   31

Ej. 13:  Combinar vectores. Mujeres menores de edad. Si colocamos el 0, se elegirían mujeres o menores de edad.

edad[mujer==1 & edad<=17]
Ana Lia
 17  12

Como se puede ver hasta aquí, solo es cuestión de usar adecuadamente los operadores lógicos, que vimos en la clase 1.

Ej. 14: acceder a una secuencia de posiciones.

Esta forma de generar secuencias igual nos permite acceder por ejemplo a todas las posiciones pares, 2,4,6,…. Pensar en alguna aplicación, imagine tiene una data mensual, y desea acceder a los meses 1,4,7,10… Para esto usamos la función seq.

1:5
[1] 1 2 3 4 5
 
seq(1,5, by = 1)
[1] 1 2 3 4 5

En el siguiente ejemplo accederá a las posiciones entre 1 y 7, iniciando en 1 y avanzando de 2 en e2. 1,3,4…

edad[seq(1, length(edad), by = 2)]
  Ana Juana   Lia Julio
   17    31    12    24

Ej. 15: Ordenar de menor a mayor. Como orden regresa la posición que ocupa cada elemento en un vector determinado, caminando de menor a mayor. Como podemos ver en el ejemplo siguiente, se verifica la salida de order, 5, 4, 6, … indica que el menor valor se encuentra en la posición 5 del vector, el segundo menor, en la posición 4, y así de forma sucesiva.

edad
  Ana Pedro Juana Johan   Lia  <NA> Julio
   17    22    31    14    12    14    24
 
order(edad)
[1] 5 4 6 1 2 7 3

Ahora, podemos obtener un vector ordenado de menor a mayor edad.

edad[order(edad)]
  Lia Johan  <NA>   Ana Pedro Julio Juana
   12    14    14    17    22    24    31

 Igual, podemos obtener los nombres, ordenados del menor al mayor.

nombres[order(edad)]
[1] "Lia"   "Johan" NA      "Ana"   "Pedro" "Julio" "Juana" 

Ej. 16: acceder al mínimo y máximo usando order. Solo debemos acceder a las posiciones extremas.

Persona con menor edad.

edad[order(edad)[1]]
Lia
 12

Persona con mayor edad.

edad[order(edad)[length(edad)]]
Juana
   31

Ej. 17: Indexación con which. Este comando regresa la posición donde la variable cumple con una determinada condición. Por ejemplo, identificar la posición donde se encuentra la persona con la edad mínima.

which(edad==min(edad))
Lia
  5

Entonces, podemos usar esta coordenada para poder recuperar la edad mínima disponible en la data. Antes vimos en qué posición estaba, ahora vemos cual es su edad.

edad[which(edad==min(edad))]
Lia
 12

Ej. 18: Podemos seleccionar múltiples nombres usando el operador contenido en. El mismo permite indicar un conjunto de nombres, para posteriormente, buscar estos nombres en el vector de edad.

personas_seleccionadas <- c("Ana", "Lia", "Julio")
edad[names(edad) %in% personas_seleccionadas]
Ana   Lia Julio
   17    12    24

Ej. 19: Acceder a las edades correspondientes a los nombres que no son NA.

edad[!is.na(names(edad))]
  Ana Pedro Juana Johan   Lia Julio
   17    22    31    14    12    24

Ej. 20: Igual podemos usar la indexación para modificar valores de un vector. Podemos verificar que un nombre es NA, entonces podemos identificar ese NA y cambiarlo por Walker.

nombres[is.na(nombres)] <- "Walker"
names(edad) <- nombres
edad
   Ana  Pedro  Juana  Johan    Lia Walker  Julio
    17     22     31     14     12     14     24

Ahora vamos a sustituir la edad de Walker, por 16 años. Debido a que el se cambio la edad cuando jugaba pelota. Le quitaron dos años.

edad[nombres == "Walker"] <- 16
edad
   Ana  Pedro  Juana  Johan    Lia Walker  Julio
    17     22     31     14     12     16     24 

Espero les haya sido de utilidad… Gracias por visitarnos.