Как смоделировать интересующие вас количества, используя пакеты arm или rstanarm в R?

Question

Я хотел бы знать, как имитировать представляющие интерес количества из регрессионной модели, оцененной с использованием либо arm или rstanarm пакеты в R. Я новичок в байесовских методах и R и использую Zelig пакет на некоторое время. Я задавал подобный вопрос раньше, но я хотел бы знать, возможно ли имитировать эти величины, используя апостериорное распределение, оцененное этими пакетами.

В Zelig Вы можете установить значения, которые вы хотите для независимых значений, и он вычисляет результаты для переменной результата (ожидаемое значение, вероятность и т. д.). Пример:

# Creating a dataset:
set.seed(10)
x <- rnorm(100,20,10)
z <- rnorm(100,10,5)
e <- rnorm(100,0,1)
y <- 2*x+3*z+e
df <- data.frame(x,z,e,y)

# Loading Zelig
require(Zelig)

# Model
m1.zelig <- zelig(y ~ x + z, model="ls", data=df)
summary(m1.zelig)

# Simulating z = 10
s1 <- setx(m1.zelig, z = 10)
simulation <- sim(m1.zelig, x = s1)
summary(simulation)

Так Zelig сохраняет х в среднем (20,56) и моделирует интересующее количество с z = 10. В этом случае у приблизительно 71.

Та же модель с использованием arm:

# Model
require(arm)
m1.arm <- bayesglm(y ~ x + z, data=df)
summary(m1.arm)

И используя rstanarm:

# Model
require(rstanarm)
m1.stan <- stanlm(y ~ x + z, data=df)
print(m1.stan)

Есть ли способ смоделировать z = 10 и x равно его среднему значению с апостериорным распределением, оцененным этими двумя пакетами, и получить ожидаемое значение y? Большое спасибо!

Answer 1

В случае с байесглом вы можете сделать

sims <- arm::sim(m1.arm, n = 1000)
y_sim <- rnorm(n = 1000, mean = sims@coef %*% t(as.matrix(s1)), sd = sims@sigma)
mean(y_sim)

Для (неизданного) rstanarm это было бы похоже

sims <- as.matrix(m1.stan)
y_sim <- rnorm(n = nrow(sims), mean = sims[,1:(ncol(sims)-1)] %*% t(as.matrix(s1)), 
               sd = sims[,ncol(sims)])
mean(y_sim)

В общем, для Стэна вы можете передать s1 как row_vector и использовать его в блоке сгенерированных количеств файла.stan, например

generated quantities {
  real y_sim;
  y_sim <- normal_rng(s1 * beta, sigma);
}

в этом случае заднее распределение y_sim появится, когда вы print задняя сводка.