1 Introducción

El modelo de regresión se ocupa de caracterizar cómo el proceso generativo asociado con una variable aleatoria \(y\) varía junto con otra variable o conjunto de variables \(\boldsymbol{x} = (x_1,\ldots,x_p)\). Específicamente, un modelo de regresión especifica una forma para \(p (y \mid \boldsymbol{x})\), la distribución condicional de \(y\) dado \(\boldsymbol{x}\). La estimación de \(p (y \mid \boldsymbol{x})\) se realiza utilizando el vector de observaciones \(\boldsymbol{y} = (y_1,\ldots,y_n)\) que se recopilan bajo una variedad de condiciones \(\boldsymbol{x}_1,\ldots,\boldsymbol{x}_n\), con \(\boldsymbol{x}_i = (x_{i,1},\ldots,x_{i,p})\) para \(i=1,\ldots,n\).

Una solución a este problema es asumir que \(p (y \mid \boldsymbol{x})\) es una función suave de \(\boldsymbol{x}\), de modo que los valores de \(\boldsymbol{x}\) pueden incidir en el proceso generativo de \(y\). Un modelo de regresión lineal es un tipo particular de modelo para \(p (y \mid \boldsymbol{x})\), el cual especifica que \(\textsf{E} (y \mid \boldsymbol{x})\) tiene una forma lineal en un conjunto de parámetros \(\boldsymbol{\beta} = (\beta_1,\ldots,\beta_p)\) como sigue: \[ \textsf{E} (y \mid \boldsymbol{x}) = \int_{\mathcal{Y}} y\, p (y \mid \boldsymbol{x})\,\text{d}y = \sum_{k=1}^p \beta_k x_k = \boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}\,. \] El modelo de regresión lineal Normal especifica que la variabilidad alrededor de \(\textsf{E} (y \mid \boldsymbol{x})\) surge por medio de una distribución Normal: \[ y_i \mid \boldsymbol{x}_i,\boldsymbol{\beta},\sigma^2 \stackrel{\text {iid}}{\sim} \textsf{N}(\boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}_i,\sigma^2) \qquad\Longleftrightarrow\qquad y_i = \boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}_i + \epsilon_i\,,\quad\epsilon_i\mid\sigma^2\stackrel{\text {iid}}{\sim} \textsf{N}(0,\sigma^2) \] para \(\quad i=1,\ldots,n\). Equivalentemente, \[ \boldsymbol{y} \mid \mathbf{X},\boldsymbol{\beta},\sigma^2 \sim \textsf{N}_n(\mathbf{X}\boldsymbol{\beta},\sigma^2\mathbf{I}) \qquad\Longleftrightarrow\qquad \boldsymbol{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}\,,\quad\boldsymbol{\epsilon}\mid\sigma^2 \sim \textsf{N}_n(\boldsymbol{0},\sigma^2\mathbf{I}) \] donde \(\mathbf{X} = [\boldsymbol{x}_1,\ldots,\boldsymbol{x}_n]^{\textsf{T}}\), \(\boldsymbol{\epsilon} = (\epsilon_1,\ldots,\epsilon_n)\) y \(\mathbf{I}\) es la matriz identidad.

Esta formulación especifica completamente la distribución muestral de los datos: \[\begin{align*} p(\boldsymbol{y} \mid \mathbf{X},\boldsymbol{\beta},\sigma^2) &= \prod_{i=1}^n \textsf{N}(y_i\mid\boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}_i,\sigma^2) \\ &\propto (\sigma^2)^{-n/2}\,\exp\left\{ -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}_i)^2 \right\} \\ &= (\sigma^2)^{-n/2}\,\exp\left\{ -\frac{1}{2\sigma^2} ( \boldsymbol{y} - \mathbf{X}\boldsymbol{\beta} )^{\textsf{T}} ( \boldsymbol{y} - \mathbf{X}\boldsymbol{\beta} ) \right\} \end{align*}\]

2 Mínimos cuadrados ordinarios

¿Qué valor de \(\boldsymbol{\beta}\) es más adecuado para los datos observados \(\boldsymbol{y}\) y \(\mathbf{X}\)?

Dados \(\boldsymbol{y}\) y \(\mathbf{X}\), el exponente de la distribución muestral se maximiza cuando la suma de cuadrados residual \[ \textsf{SSR}(\boldsymbol{\beta})=\sum_{i=1}^n (y_i - \boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}_i)^2 = \boldsymbol{\beta}^{\textsf{T}}\mathbf{X}^{\textsf{T}}\mathbf{X}\boldsymbol{\beta}^{\textsf{T}} - 2\boldsymbol{\beta}^{\textsf{T}}\mathbf{X}^{\textsf{T}}\boldsymbol{y} + \boldsymbol{y}^{\textsf{T}}\boldsymbol{y}^{\textsf{T}} \] se minimiza. Esto se logra cuando \(\boldsymbol{\beta}\) asume el valor \[ \hat{\boldsymbol{\beta}}_{\text{ols}} = (\mathbf{X}^{\textsf{T}}\mathbf{X})^{-1}\mathbf{X}^{\textsf{T}}\boldsymbol{y} \] el cual se denomina estimador de mínimos cuadrados ordinarios (OLS) de \(\hat{\boldsymbol{\beta}}\).

Además, se puede demostrar que un estimador insesgado de \(\sigma^2\) está dado por \[ \hat\sigma^2_{\text{ols}} = \frac{1}{n-p}\,\textsf{SSR}(\hat{\boldsymbol{\beta}}_{\text{ols}}) \]

3 Inferencia Bayesiana

Motivación:

Intervalos de cridibilidad para \(\beta_j\), para \(j=1,\ldots,p\).
Probabilidades posteriores de la forma \(\textsf{Pr}(\beta_j > 0 \mid \mathbf{X}, \boldsymbol{y})\), para \(j=1,\ldots,p\).
El estimador ols tiende a sobreajustar los datos cuando \(p\) es grande, el estimador Bayesiano es más conservador.
Mecanismo natural para seleccionar modelos.

3.1 Modelo semi-conjugado

3.1.1 Modelo

Distribución muestral: \[ \boldsymbol{y} \mid \mathbf{X},\boldsymbol{\beta},\sigma^2 \sim \textsf{N}_n(\mathbf{X}\boldsymbol{\beta},\sigma^2\mathbf{I}) \]

Previa: \[\begin{align*} \boldsymbol{\beta} &\sim \textsf{N}_n(\boldsymbol{\beta}_0, \mathbf{\Sigma}_0) \\ \sigma^2 & \sim \textsf{GI}\left( \tfrac{\nu_0}{2}, \tfrac{\nu_0\sigma^2_0}{2} \right) \end{align*}\]

Parámetros: \((\boldsymbol{\beta},\sigma^2)\).

Hiperparámetros: \((\boldsymbol{\beta}_0, \mathbf{\Sigma}_0, \nu_0, \sigma^2_0)\).

Distribuciones condicionales completas: \[\begin{align*} \boldsymbol{\beta}\mid\text{resto} &\sim \textsf{N}_n\left( (\mathbf{\Sigma}_0^{-1} + \tfrac{1}{\sigma^2}\mathbf{X}^{\textsf{T}}\mathbf{X} )^{-1} (\mathbf{\Sigma}_0^{-1}\boldsymbol{\beta}_0 + \tfrac{1}{\sigma^2}\mathbf{X}^{\textsf{T}}\boldsymbol{y} ) , (\mathbf{\Sigma}_0^{-1} + \tfrac{1}{\sigma^2}\mathbf{X}^{\textsf{T}}\mathbf{X} )^{-1} \right) \\ \sigma^2 \mid \text{resto} &\sim \textsf{GI}\left( \frac{\nu_0 + n}{2}, \frac{\nu_0\sigma^2_0 + \textsf{SSR}(\boldsymbol{\beta})}{2} \right) \end{align*}\]

Observe que:

Si \(\mathbf{\Sigma}_0 \longrightarrow \mathbf{0}\,\), entonces \(\textsf{E}(\boldsymbol{\beta}\mid\mathbf{X}, \boldsymbol{y}) \longrightarrow \hat{\boldsymbol{\beta}}_{\text{ols}}\).
Si \(\sigma^2 \longrightarrow \infty\), entonces \(\textsf{E}(\boldsymbol{\beta}\mid\mathbf{X}, \boldsymbol{y}) \longrightarrow \boldsymbol{\beta}_0\).

3.1.2 Previa unitaria

Encontrar valores de los hiperparámetros que representen información previa sustantiva puede ser difícil; y se vuelve aun más difícil a medida que aumenta \(p\), dado que en ese caso el número de hiperparámetros aumenta drásticamente.

Si no hay información previa relevante, entonces debe ser lo menos informativa posible. En este escenario, la distribución posterior corresponde a la información posterior de una analista con poco conocimiento de la población.

Una distribución previa de información unitaria (unit information prior; Kass y Wasserman, 1995) es aquella que contiene la información asociada con una sola observación. Dado que \(\textsf{E}(\hat{\boldsymbol{\beta}}_{\text{ols}}) = \boldsymbol{\beta}\) y \((\,\textsf{Var}(\hat{\boldsymbol{\beta}}_{\text{ols}})\,)^{-1} = \tfrac{1}{\sigma^2} (\mathbf{X}^{\textsf{T}}\mathbf{X})\), se propone utilizar \[ \boldsymbol{\beta}_0 = \hat{\boldsymbol{\beta}}_{\text{ols}}\,,\qquad \mathbf{\Sigma}_0 = n\,\hat\sigma^2_{\text{ols}} (\mathbf{X}^{\textsf{T}}\mathbf{X})^{-1}\,,\qquad \nu_0 = 1\,,\qquad \sigma^2_0 = \hat\sigma^2_{\text{ols}}\,. \] Estrictamente, esta distribución no es una distribución previa, ya que utiliza los datos observados. Sin embargo, como solo se utiliza una pequeña fracción de información, \(\tfrac{1}{n}\), entonces se puede interpretar como la distribución previa de un analista con información previa débil.

Algoritmo de MCMC:

Simular \(\sigma^2 \sim p(\sigma^2\mid\boldsymbol{\beta},\mathbf{X},\boldsymbol{y})\).
Simular \(\boldsymbol{\beta} \sim p(\boldsymbol{\beta}\mid\sigma^2,\mathbf{X},\boldsymbol{y})\).

3.1.3 Previa g

Otro principio para construir una distribución previa se basa en que la estimación de los parámetros debe ser invariante a cambios en la escala de los regresores. Esta condición se satisface si \(\boldsymbol{\beta}_0 = \boldsymbol{0}\) y \(\mathbf{\Sigma}_0 = k\,(\mathbf{X}^{\textsf{T}}\mathbf{X})^{-1}\) con \(k > 0\).

La distribución previa g (g-prior; Zellner, 1986) sugiere que \(k = g\,\sigma^2\) con \(k > 0\) (si \(g = n\), entonces se tiene en cuenta la información correspondiente a una sola observación). Además, se tiene que \[ \boldsymbol{\beta}\mid\text{resto} \sim \textsf{N}_n\left( \tfrac{g}{g+1}(\mathbf{X}^{\textsf{T}}\mathbf{X})^{-1}\mathbf{X}^{\textsf{T}}\boldsymbol{y} , \tfrac{g}{g+1}\sigma^2(\mathbf{X}^{\textsf{T}}\mathbf{X} )^{-1} \right) \] Para simular fácilmente de la distribución posterior de \((\boldsymbol{\beta},\sigma^2)\) se observa que \[ p(\boldsymbol{\beta},\sigma^2\mid\mathbf{X},\boldsymbol{y}) = p(\boldsymbol{\beta}\mid\sigma^2,\mathbf{X},\boldsymbol{y})\,p(\sigma^2\mid\mathbf{X},\boldsymbol{y}) \] y \[\begin{align*} p(\sigma^2\mid\mathbf{X},\boldsymbol{y}) &\propto p(\boldsymbol{y}\mid\sigma^2,\mathbf{X})\,p(\sigma^2) \\ &=\int_{\mathbb{R}^p} p(\boldsymbol{y}\mid\boldsymbol{\beta},\sigma^2,\mathbf{X})\,p(\boldsymbol{\beta}\mid\sigma^2,\mathbf{X},\boldsymbol{y})\,\text{d}\boldsymbol{\beta} \quad p(\sigma^2) \\ &\propto \textsf{GI}\left( \sigma^2 \mid \frac{\nu_0 + n}{2}, \frac{\nu_0\sigma^2_0 + \textsf{SSR}_g}{2} \right) \end{align*}\] donde \(\textsf{SSR}_g = \boldsymbol{y}^{\textsf{T}}\left( \mathbf{I}_n - \tfrac{g}{g+1}\mathbf{X}(\mathbf{X}^{\textsf{T}}\mathbf{X})^{-1}\mathbf{X}^{\textsf{T}} \right)\boldsymbol{y}\).

Algoritmo de Monte Carlo (¡no es MCMC!):

Simular \(\sigma^2 \sim p(\sigma^2\mid\mathbf{X},\boldsymbol{y})\).
Simular \(\boldsymbol{\beta} \sim p(\boldsymbol{\beta}\mid\sigma^2,\mathbf{X},\boldsymbol{y})\).

4 Previa difusa independiente

\[ \boldsymbol{\beta}_0 = \boldsymbol{0}\,,\qquad \mathbf{\Sigma}_0 = 100\,\mathbf{I}_p\,,\qquad \nu_0 = 1\,,\qquad \sigma^2_0 = 100. \]

5 Ejemplo: Absorción de oxígeno

Estudio del efecto de dos regímenes de ejercicio sobre la absorción de oxígeno.

Seis de doce hombres fueron asignados aleatoriamente a un programa de carrera en terreno plano de 12 semanas, y los seis restantes fueron asignados a un programa de aeróbicos de 12 semanas.

Se midió el consumo máximo de oxígeno de cada sujeto (en litros por minuto) al correr en una cinta de correr inclinada, tanto antes como después del programa de 12 semanas.

El objetivo consiste en evaluar cómo el cambio en la absorción máxima de oxígenodepende del programa de entrenamiento.

\(y\) : cambio en la absorción máxima de oxígeno.
trat : programa de entrenamiento (1 si aeróbicos, 0 si correr).
edad : edad (en años).
\(n\) : tamaño de la muestra.

Modelo: \(\textsf{E}(y\mid\mathbf{X}) = \beta_1x_1 + \beta_2x_2 + \beta_3x_3 + \beta_4x_4\)

\(x_1 = 1\).
\(x_2 = 1\) si aeróbicos; \(x_2 = 0\) si correr.
\(x_3 = \text{edad}\).
\(x_4 = x_2*x_3\).

Por lo tanto:

Correr: \(\textsf{E}(y\mid\mathbf{X}) = \beta_1 + \beta_3\,\text{edad}\).
Aeróbicos: \(\textsf{E}(y\mid\mathbf{X}) = (\beta_1 + \beta_2) + (\beta_3 + \beta_4)\,\text{edad}\).

# data
trat <- c(0,0,0,0,0,0,1,1,1,1,1,1)
edad <- c(23,22,22,25,27,20,31,23,27,28,22,24)
y    <- c(-0.87,-10.74,-3.27,-1.97,7.50,-7.25,17.05,4.96,10.40,11.05,0.26,2.51)

# data
y <- as.matrix(y)
X <- cbind(1, trat, edad, trat*edad)

# dimensiones
n <- dim(X)[1]
p <- dim(X)[2]
colnames(X) <- paste0("x", 1:p)

# scatterplot
par(mfrow=c(1,1), mar=c(3,3,1,1), mgp=c(1.75,.75,0))
plot(y ~ edad, pch=16, xlab = "Edad", ylab = "Cambio en la absorción", 
     col = c("black","gray")[trat+1])
legend("topleft", legend=c("Aeróbicos","Correr"), pch=c(16,16), col=c("gray","black"))

# OLS
beta.ols <- solve(t(X)%*%X)%*%t(X)%*%y
beta.ols

##           [,1]
## x1 -51.2939459
## x2  13.1070904
## x3   2.0947027
## x4  -0.3182438

sig2.ols <- sum((y - X%*%beta.ols)^2)/(n-p)
sig2.ols

## [1] 8.542477

# fit.ols <- lm(y ~ -1 + X)
# summary(fit.ols)$coef
# summary(fit.ols)$sigm

# hiperparametros (previa g)
nu0 <- 1
s20 <- sig2.ols
g   <- n

# calculos para la posterior
Hg    <- (g/(g+1))*(X%*%solve(t(X)%*%X)%*%t(X))
SSRg  <- t(y)%*%(diag(1,n) - Hg)%*%y
Vbeta <- (g/(g+1))*solve(t(X)%*%X)
Ebeta <- Vbeta%*%t(X)%*%y

# Monte Carlo
S <- 5000
s2.post   <- matrix(NA, S, 1) 
beta.post <- matrix(NA, S, p) 
set.seed(1)
for(s in 1:S) {
      s2.post[s] <- 1/rgamma(1, (nu0 + n)/2, (nu0*s20 + SSRg)/2)
      beta.post[s,] <- c(mvtnorm::rmvnorm(1, Ebeta, s2.post[s]*Vbeta))
}

# inferencia sigma^2
quantile(s2.post, probs = c(.025,.5,.975))

##      2.5%       50%     97.5% 
##  5.373982 10.820749 26.550250

# inferencia sigma
quantile(sqrt(s2.post), probs = c(.025,.5,.975))

##     2.5%      50%    97.5% 
## 2.318185 3.289491 5.152693

# inferencia beta
apply(X = beta.post, MARGIN = 2, FUN = quantile, probs = c(.025,.5,.975))

##            [,1]      [,2]      [,3]       [,4]
## 2.5%  -74.35577 -23.88030 0.7419868 -1.7263041
## 50%   -46.89450  11.60457 1.9210013 -0.2808402
## 97.5% -19.19350  46.81701 3.0954272  1.1870935

apply(X = beta.post, MARGIN = 2, FUN = mean)

## [1] -46.9234476  11.6863701   1.9155231  -0.2774947

apply(X = beta.post, MARGIN = 2, FUN = sd)

## [1] 14.0065219 18.1531796  0.6020815  0.7480802

colMeans(beta.post > 0)

## [1] 0.0014 0.7522 0.9980 0.3494

par(mfrow=c(2,2),mar=c(2.75,2.75,.5,.5),mgp=c(1.7,.7,0))
# beta2 posterior
x<-seq(-85,130,length=200)
plot(density(beta.post[,2],adj=2), xlab=expression(beta[2]),main="",ylab="",lwd=2)
abline(v=0,col="gray")
abline(v=quantile(beta.post[,2], c(0.025,0.975)),col="blue",lty=2)
# beta4 posterior
x<-seq(-5,5,length=100)
plot(density(beta.post[,4],adj=2),xlab=expression(beta[4]),main="",ylab="",lwd=2)
abline(v=0,col="gray")
abline(v=quantile(beta.post[,4], c(0.025,0.975)),col="blue",lty=2)
# (beta2, beta4) posterior
plot(beta.post[,c(2,4)],type = "p",cex=0.1,xlab=expression(beta[2]),ylab=expression(beta[4]))
abline(h=0,col="gray") ; abline(v=0,col="gray")

# Efecto del tratamiento aerobico: Delta
#
# Delta = E(y | edad, aerobico) - E(y | edad, aerobico)
#       = ( (beta1 + beta2) + (beta3 + beta4)*edad ) - ( beta1 + beta3*edad )
#       = beta2 + beta4*edad

r.edad <- min(edad):max(edad)
n.edad <- length(r.edad)
TE <- matrix(NA, S, n.edad)
for (j in 1:n.edad)
      TE[,j] <- beta.post[,2] + beta.post[,4]*r.edad[j]

that <- colMeans(TE)
ic1  <- apply(X = TE, MARGIN = 2, FUN = function(x) quantile(x, c(0.050,0.950)))
ic2  <- apply(X = TE, MARGIN = 2, FUN = function(x) quantile(x, c(0.025,0.975)))
colo <- c("blue","black")[as.numeric((ic2[1,] < 0) & (0 < ic2[2,]))+1]

par(mfrow=c(1,1),mar=c(3,3,1.5,1),mgp=c(1.75,.75,0))
plot(NA, NA, xlab = "Edad", ylab = "Efecto tratamiento", main = "", 
     xlim = range(edad), ylim = range(TE), cex.axis = 0.75, xaxt = "n")
axis(side = 1, at = r.edad, labels = r.edad, cex.axis = 0.75)
abline(h = 0, col = "gray", lwd = 2)
abline(v = r.edad, col = "gray95", lwd = 1, lty = 2) 
for (j in 1:n.edad) {
      segments(x0 = r.edad[j], y0 = ic1[1,j], x1 = r.edad[j], y1 = ic1[2,j], lwd = 3, col = colo[j])
      segments(x0 = r.edad[j], y0 = ic2[1,j], x1 = r.edad[j], y1 = ic2[2,j], lwd = 1, col = colo[j])
      lines(x = r.edad[j], y = that[j], type = "p", pch = 16, cex = 0.8, col = colo[j])
}

6 Ejemplo: Jags

# descarcar JAGS
# http://www.sourceforge.net/projects/mcmc-jags/files

# user manuals
# https://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf
# http://www.jkarreth.net/files/bayes-cph_Tutorial-JAGS.pdf 

library(R2jags)

## Warning: package 'R2jags' was built under R version 4.1.2

## Loading required package: rjags

## Warning: package 'rjags' was built under R version 4.1.2

## Loading required package: coda

## Linked to JAGS 4.3.0

## Loaded modules: basemod,bugs

## 
## Attaching package: 'R2jags'

## The following object is masked from 'package:coda':
## 
##     traceplot

# data
trat <- c(0,0,0,0,0,0,1,1,1,1,1,1)
edad <- c(23,22,22,25,27,20,31,23,27,28,22,24)
y    <- c(-0.87,-10.74,-3.27,-1.97,7.50,-7.25,17.05,4.96,10.40,11.05,0.26,2.51)

# data
y <- c(as.matrix(y))
X <- cbind(1, trat, edad, trat*edad)

# dimensiones
n <- dim(X)[1]
p <- dim(X)[2]
colnames(X) <- paste0("x", 1:p)

# la Normal está parametrizada en términos de la precisión
# los vectores deben tener el mismo formato
# notación:
#     phi    = 1/sigma^2
#     Omega0 = X^T*X
#     a0     = nu0/2
#     b0     = nu0*sigma0^2/2

model <- function() {
      for (i in 1:n) {
            y[i] ~ dnorm(inprod(X[i,], beta), phi)
      }
      beta[1:p] ~ dmnorm(beta0[1:p], (phi/g)*Omega0[1:p,1:p])
      phi ~ dgamma(a0, b0)
}

# previa
beta0  <- c(rep(0,p))
Omega0 <- t(X)%*%X
g      <- n
nu0    <- 1
s20    <- sig2.ols
a0     <- nu0/2
b0     <- nu0*s20/2

# input
model_data <- list(y = y, X = X, n = n, p = p, g = g, beta0 = beta0, Omega0 = Omega0, a0 = a0, b0 = b0)

# parameters
model_parameters <- c("beta", "phi")

# initial values
initial_values <- list(list("beta" = beta0, "phi" = 1/s20), 
                       list("beta" = beta0, "phi" = 1/s20),
                       list("beta" = beta0, "phi" = 1/s20))

# mcmc settings
niter  <- 6000
nburn  <- 1000
nthin  <- 1
nchain <- length(initial_values)

# mcmc
set.seed(123)
fit <- jags(data = model_data, inits = initial_values, 
            parameters.to.save = model_parameters, model.file = model, 
            n.chains = nchain, n.iter = niter, n.thin = nthin, n.burnin = nburn)

## module glm loaded

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 12
##    Unobserved stochastic nodes: 2
##    Total graph size: 114
## 
## Initializing model

print(fit)

## Inference for Bugs model at "C:/Users/JUANCA~1/AppData/Local/Temp/RtmpExxIXt/model182c36456937.txt", fit using jags,
##  3 chains, each with 6000 iterations (first 1000 discarded)
##  n.sims = 15000 iterations saved
##          mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
## beta[1]  -47.327  14.044 -75.409 -56.399 -47.330 -38.473 -19.633 1.001  9900
## beta[2]   12.046  17.933 -23.380   0.634  12.030  23.391  47.656 1.001 15000
## beta[3]    1.934   0.603   0.740   1.549   1.934   2.320   3.132 1.001 11000
## beta[4]   -0.293   0.740  -1.767  -0.764  -0.287   0.180   1.163 1.001 15000
## phi        0.098   0.038   0.038   0.070   0.093   0.121   0.187 1.001  5900
## deviance  61.687   3.700  56.375  58.942  61.098  63.714  70.543 1.002  3200
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 6.8 and DIC = 68.5
## DIC is an estimate of expected predictive error (lower deviance is better).

# numero de muestras
S <- fit$BUGSoutput$n.sims

# dic
fit$BUGSoutput$DIC

## [1] 68.52862

# transformar a objecto MCMC               
fit_mcmc <- coda::as.mcmc(fit)

# diagnosticos
# mcmcplots::mcmcplot(fit_mcmc)  # despliegaa graficos en el browser
superdiag::superdiag(fit_mcmc)   # super diagnostico

## Number of chains = 3 
## Number of iterations = 5000 per chain before discarding the burn-in period
## Burn-in period = 2500 per chain
## Sample size in total = 7503 
## 
## ****************** The Geweke diagnostic: ******************
## Windows:
##            chain 1 chain 2 chain 3
## From start     0.1  0.6190  0.5577
## From stop      0.5  0.0681  0.3008
## 
## Z-scores:
##          chain 1 chain 2  chain 3
## beta[1]  -0.6470 -1.2963  2.41505
## beta[2]   0.8750  0.9512 -1.54531
## beta[3]   0.6927  1.1277 -2.35547
## beta[4]  -0.9195 -1.0987  1.63632
## deviance -0.4797 -0.3621 -0.68466
## phi       1.5208 -0.5036 -0.05521
## 
## *************** The Gelman-Rubin diagnostic: ***************
## Potential scale reduction factors:
##          Point est. Upper C.I.
## beta[1]      1.0002      1.001
## beta[2]      0.9999      1.000
## beta[3]      1.0003      1.002
## beta[4]      0.9999      1.000
## deviance     1.0006      1.001
## phi          1.0005      1.001
## 
## Multivariate psrf: 1.001
## 
## ************* The Heidelberger-Welch diagnostic ************
## Chain 1:
## epsilon=0.1, alpha=0.05                                       
##          Stationarity start     p-value
##          test         iteration        
## beta[1]  passed         1       0.8702 
## beta[2]  passed         1       0.5776 
## beta[3]  passed         1       0.8584 
## beta[4]  passed         1       0.6244 
## deviance passed         1       0.4199 
## phi      passed       752       0.1466 
##                                       
##          Halfwidth Mean      Halfwidth
##          test                         
## beta[1]  passed    -47.58947 0.551435 
## beta[2]  passed     12.45163 0.695928 
## beta[3]  passed      1.94654 0.023663 
## beta[4]  passed     -0.31152 0.028632 
## deviance passed     61.76186 0.192992 
## phi      passed      0.09557 0.002508 
## 
## Chain 2:
## epsilon=0.062, alpha=0.005                                       
##          Stationarity start     p-value
##          test         iteration        
## beta[1]  passed       1         0.8479 
## beta[2]  passed       1         0.8033 
## beta[3]  passed       1         0.8538 
## beta[4]  passed       1         0.7738 
## deviance passed       1         0.2813 
## phi      passed       1         0.1887 
##                                      
##          Halfwidth Mean     Halfwidth
##          test                        
## beta[1]  passed    -47.0904 0.545035 
## beta[2]  passed     12.2351 0.705667 
## beta[3]  passed      1.9228 0.023442 
## beta[4]  failed     -0.2998 0.029111 
## deviance passed     61.7751 0.209334 
## phi      passed      0.0984 0.002531 
## 
## Chain 3:
## epsilon=0.041, alpha=0.05                                       
##          Stationarity start     p-value
##          test         iteration        
## beta[1]  passed       252       0.09693
## beta[2]  passed         1       0.18964
## beta[3]  passed       252       0.09059
## beta[4]  passed         1       0.15292
## deviance passed         1       0.71678
## phi      passed         1       0.94419
##                                       
##          Halfwidth Mean      Halfwidth
##          test                         
## beta[1]  passed    -47.87153 0.575292 
## beta[2]  failed     12.20540 0.703709 
## beta[3]  passed      1.95735 0.024770 
## beta[4]  failed     -0.30175 0.029057 
## deviance passed     61.67226 0.188225 
## phi      passed      0.09772 0.002797 
## 
## *************** The Raftery-Lewis diagnostic ***************
## Chain 1:
## Convergence eps = 0.001
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
## 
## You need a sample size of at least 3746 with these values of q, r and s
## 
## Chain 2:
## Convergence eps = 5e-04
## Quantile (q) = 0.001
## Accuracy (r) = +/- 0.001
## Probability (s) = 0.99 
## 
## You need a sample size of at least 6629 with these values of q, r and s
## 
## Chain 3:
## Convergence eps = 0.0025
## Quantile (q) = 0.1
## Accuracy (r) = +/- 5e-04
## Probability (s) = 0.9 
## 
## You need a sample size of at least 973996 with these values of q, r and s
## 
## ************* The Hellinger distance diagnostic ************
## Between chains: 
##              Min     Max
## beta[1]  0.02536 0.03774
## beta[2]  0.03372 0.04264
## beta[3]  0.03084 0.03701
## beta[4]  0.03473 0.03922
## deviance 0.03101 0.03987
## phi      0.03015 0.05054
## 
## Within chain 1:
##              500    1000    1500    2000
## beta[1]  0.04593 0.05470 0.05812 0.05010
## beta[2]  0.05889 0.04258 0.05956 0.04329
## beta[3]  0.04703 0.06101 0.04613 0.04859
## beta[4]  0.05255 0.04471 0.05662 0.04056
## deviance 0.04367 0.04297 0.05842 0.06337
## phi      0.04802 0.05206 0.04757 0.05879
## 
## Within chain 2:
##              500    1000    1500    2000
## beta[1]  0.05859 0.08154 0.04281 0.04595
## beta[2]  0.05548 0.06650 0.05124 0.05723
## beta[3]  0.06358 0.07820 0.05078 0.04516
## beta[4]  0.06049 0.07462 0.05292 0.05854
## deviance 0.04665 0.05068 0.03630 0.06884
## phi      0.05950 0.05604 0.04386 0.07249
## 
## Within chain 3:
##              500    1000    1500    2000
## beta[1]  0.05969 0.05009 0.05010 0.05237
## beta[2]  0.05598 0.04865 0.04468 0.04368
## beta[3]  0.05238 0.04401 0.04795 0.05246
## beta[4]  0.04915 0.04757 0.05244 0.04531
## deviance 0.04905 0.06116 0.05671 0.07703
## phi      0.04095 0.07750 0.07249 0.05426

# algunos graficos
mcmcplots::denplot(fit_mcmc)

## Registered S3 method overwritten by 'mcmcplots':
##   method        from  
##   as.mcmc.rjags R2jags

mcmcplots::traplot(fit_mcmc)

mcmcplots::caterplot(fit_mcmc, parms = c("beta[2]","beta[4]"))

# mas opciones
# densidades
# fit_mcmc_gg <- ggmcmc::ggs(fit_mcmc)
# ggmcmc::ggs_density(fit_mcmc_gg)

Regresión lineal

Juan Sosa PhD

jcsosam@unal.edu.co