RJAGS 的貝氏迴歸

使用 RJAGS 的貝氏建模

Alicia Johnson

Associate Professor, Macalester College

貝氏迴歸模型

$Y_i$ = 第 $i$ 位成年人的體重(kg)
$X_i$ = 第 $i$ 位成年人的身高(cm)

$\;$

模型

$Y_i \sim N(m_i, s^2)$
$m_i = a + b X_i$

$a \sim N(0, 200^2)$
$b \sim N(1, 0.5^2)$
$s \sim \text{Unif}(0, 20)$

使用 RJAGS 的貝氏建模

先驗洞見

使用 RJAGS 的貝氏建模

來自觀察到的體重與身高的洞見

$Y_i \sim N(m_i, s^2)$
$m_i = a + b X_i$

wt_mod <- lm(wgt ~ hgt, bdims)
coef(wt_mod)
(Intercept)         hgt 
-105.011254    1.017617 
summary(wt_mod)$sigma
9.30804
使用 RJAGS 的貝氏建模

定義迴歸模型

weight_model <- "model{
    # Likelihood model for Y[i]





    # Prior models for a, b, s




}"  
使用 RJAGS 的貝氏建模

定義迴歸模型

weight_model <- "model{
    # Likelihood model for Y[i]
    for(i in 1:length(Y)) {


    }

    # Prior models for a, b, s




}"  
  • 對於 $i$ 從 1 到 507,有 $Y_i \sim N(m_i, s^2)$
使用 RJAGS 的貝氏建模

定義迴歸模型

weight_model <- "model{
    # Likelihood model for Y[i]
    for(i in 1:length(Y)) {
      Y[i] ~ dnorm(m[i], s^(-2))

    }

    # Prior models for a, b, s




}"  
  • 對於 $i$ 從 1 到 507,有 $Y_i \sim N(m_i, s^2)$
使用 RJAGS 的貝氏建模

定義迴歸模型

weight_model <- "model{
    # Likelihood model for Y[i]
    for(i in 1:length(Y)) {
      Y[i] ~ dnorm(m[i], s^(-2))
      m[i] <- a + b * X[i]
    }

    # Prior models for a, b, s




}"  
  • 對於 $i$ 從 1 到 507,有 $Y_i \sim N(m_i, s^2)$
  • $m_i = a + b X_i$
    • 「注意」:使用 <-,不要用 ~
使用 RJAGS 的貝氏建模

定義迴歸模型

weight_model <- "model{
    # Likelihood model for Y[i]
    for(i in 1:length(Y)) {
      Y[i] ~ dnorm(m[i], s^(-2))
      m[i] <- a + b * X[i]
    }

    # Prior models for a, b, s
    a ~ dnorm(0, 200^(-2))
    b ~ dnorm(1, 0.5^(-2))
    s ~ dunif(0, 20)

}"  
  • 對於 $i$ 從 1 到 507,有 $Y_i \sim N(m_i, s^2)$
  • $m_i = a + b X_i$

    • 「注意」:使用 <-,不要用 ~
  • $a \sim N(0, 200^2)$

  • $b \sim N(1, 0.5^2)$
  • $s \sim \text{Unif}(0, 20)$
使用 RJAGS 的貝氏建模

編譯迴歸模型

# 編譯模型    
weight_jags <- jags.model(textConnection(weight_model), 
    data = list(X = bdims$hgt, Y = bdims$wgt),
    inits = list(.RNG.name = "base::Wichmann-Hill", .RNG.seed = 2018))

dim(bdims)
507  25
head(bdims$hgt)

head(bdims$wgt)
174.0 175.3 193.5 186.5 187.2 181.5

65.6 71.8 80.7 72.6 78.8 74.8
使用 RJAGS 的貝氏建模

模擬迴歸模型

# 編譯模型    
weight_jags <- jags.model(textConnection(weight_model), 
    data = list(X = bdims$hgt, Y = bdims$wgt),
    inits = list(.RNG.name = "base::Wichmann-Hill", 
                 .RNG.seed = 2018))

# 模擬後驗分佈
weight_sim <- coda.samples(model = weight_jags, 
    variable.names = c("a", "b", "s"), 
    n.iter = 10000)
使用 RJAGS 的貝氏建模

使用 RJAGS 的貝氏建模

處理馬可夫鏈不穩定

  • 標準化身高自變數(減去平均數並除以標準差)

  • 拉長鏈長度

使用 RJAGS 的貝氏建模

使用 RJAGS 的貝氏建模

後驗洞見

使用 RJAGS 的貝氏建模

一起來練習吧!

使用 RJAGS 的貝氏建模

Preparing Video For Download...