Подберите байесовскую линейную регрессию и прогнозируйте ненаблюдаемые значения

Question

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

Мои данные выглядят следующим образом:

ngroups <- 2
group <- 1:ngroups
nobs <- 100
dta <- data.frame(group=rep(group,each=nobs),y=rnorm(nobs*ngroups),x=runif(nobs*ngroups))
head(dta)

Answer 1

У JAGS есть мощные способы сделать вывод об отсутствующих данных, и как только вы освоите их, это легко! Я настоятельно рекомендую вам ознакомиться с превосходной книгой Марка Кери, в которой содержится прекрасное введение в программирование на языке BUGS (JAGS достаточно близок к ошибкам, которые передает почти все).

Самый простой способ сделать это, как вы говорите, настроить модель. Ниже я приведу полный проработанный пример того, как это работает. Но вы, похоже, спрашиваете, как получить интервал прогнозирования без повторного запуска модели (ваша модель очень большая и требует больших вычислительных ресурсов?). Это также может быть сделано.
Как прогнозировать - трудный путь (без повторного запуска модели) Для каждой итерации MCMC смоделируйте ответ для желаемого значения x на основе апостериорных прорисовок этой итерации для ковариатных значений. Итак, представьте, что вы хотите предсказать значение для X=10. Затем, если итерация 1 (после записи) имеет наклон =2, пересечение = 1 и стандартное отклонение =0,5, нарисуйте значение Y из

Y=rnorm(1, 1+2*10, 0.5)  

И повторите для итерации 2, 3, 4, 5... Это будут ваши последующие результаты для ответа при X=10. Примечание: если вы не отслеживали стандартное отклонение в своей модели JAGS, вам не повезло, и вам нужно снова подгонять модель.

Как предсказать - простой способ - с обработанным примером. Основная идея состоит в том, чтобы вставить (в ваши данные) значения x, ответы которых вы хотите предсказать, с соответствующими значениями y NA. Например, если вам нужен интервал прогнозирования для X=10, вам просто нужно включить точку (10, NA) в ваши данные и установить монитор трассировки для значения y.

Я использую JAGS от R с пакетом rjags. Ниже приведен полный обработанный пример, который начинается с моделирования данных, затем добавляет к ним некоторые дополнительные значения x, определяет и запускает линейную модель в JAGS с помощью rjags и суммирует результаты. Y[101:105] содержит отрывки из задних интервалов прогнозирования для X[101:105]. Обратите внимание, что Y[1:100] просто содержит значения y для X[1:100]. Это наблюдаемые данные, которые мы передали в модель, и они никогда не меняются при обновлении модели.

library(rjags)
# Simulate data (100 observations)
my.data <- as.data.frame(matrix(data=NA, nrow=100, ncol=2))
names(my.data) <- c("X", "Y")
# the linear model will predict Y based on the covariate X

my.data$X <- runif(100) # values for the covariate
int <- 2     # specify the true intercept
slope <- 1   # specify the true slope
sigma <- .5   # specify the true residual standard deviation
my.data$Y <- rnorm(100, slope*my.data$X+int, sigma)  # Simulate the data

#### Extra data for prediction of unknown Y-values from known X-values
y.predict <- as.data.frame(matrix(data=NA, nrow=5, ncol=2))
names(y.predict) <- c("X", "Y")
y.predict$X <- c(-1, 0, 1.3, 2, 7)

mydata <- rbind(my.data, y.predict)


set.seed(333)
setwd(INSERT YOUR WORKING DIRECTORY HERE)
sink("mymodel.txt")
cat("model{

    # Priors

    int ~ dnorm(0, .001)
    slope ~ dnorm(0, .001)
    tau <- 1/(sigma * sigma)
    sigma ~ dunif(0,10) 

    # Model structure

    for(i in 1:R){
    Y[i] ~ dnorm(m[i],tau)
    m[i] <- int + slope * X[i]
    }
    }", fill=TRUE)
sink()
jags.data <- list(R=dim(mydata)[1], X=mydata$X, Y=mydata$Y)

inits <- function(){list(int=rnorm(1, 0, 5), slope=rnorm(1,0,5),
                         sigma=runif(1,0,10))}

params <- c("Y", "int", "slope", "sigma")

nc <- 3
n.adapt <-1000
n.burn <- 1000
n.iter <- 10000
thin <- 10
my.model <- jags.model('mymodel.txt', data = jags.data, inits=inits, n.chains=nc, n.adapt=n.adapt)
update(my.model, n.burn)
my.model_samples <- coda.samples(my.model,params,n.iter=n.iter, thin=thin)
summary(my.model_samples)