Извлечение коэффициентов из объекта mcmc

Извлечение данных из объектов всегда было для меня одним из самых запутанных аспектов R. Я установил модель байесовской линейной регрессии с использованием rjags и получил следующий объект mcmc:

summary(m_csim)
Iterations = 1:150000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 150000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean        SD  Naive SE Time-series SE
BR2     0.995805 0.0007474 1.930e-06      3.527e-06
BR2adj  0.995680 0.0007697 1.987e-06      3.633e-06
b[1]   -5.890842 0.1654755 4.273e-04      1.289e-02
b[2]    1.941420 0.0390239 1.008e-04      1.991e-03
b[3]    1.056599 0.0555885 1.435e-04      5.599e-03
sig2    0.004678 0.0008333 2.152e-06      3.933e-06

2. Quantiles for each variable:

            2.5%       25%       50%       75%    97.5%
BR2     0.994108  0.995365  0.995888  0.996339  0.99702
BR2adj  0.993932  0.995227  0.995765  0.996229  0.99693
b[1]   -6.210425 -6.000299 -5.894810 -5.784082 -5.55138
b[2]    1.867453  1.914485  1.940372  1.967466  2.02041
b[3]    0.942107  1.020846  1.057720  1.094442  1.16385
sig2    0.003321  0.004082  0.004585  0.005168  0.00657

Чтобы извлечь средние значения коэффициентов, я сделал b = colMeans(mod_csim)[3:5]. Я хочу рассчитать достоверные интервалы, поэтому мне также нужно извлечь квантили 0,025 и 0,975. Как мне это сделать программно?


r mcmc jags
person FranGoitia    schedule 05.11.2018    source источник
comment
Обычно для этого есть методы. Не уверен насчет rjags, но coef и fixef обычно используются либо в объекте модели, либо в сводке объекта модели. Есть и хорошие пакеты, например tidybayes.   -  person Axeman    schedule 05.11.2018
comment
Вам будет легче помочь, если вы включите простой воспроизводимый пример с образцом ввода и желаемым выводом, которые можно использовать для тестирования и проверки возможных решений.   -  person MrFlick    schedule 05.11.2018
comment
Я бы начал с str(m_csim), чтобы понять, что это такое. Затем помните, что сводные функции много раз выполняют свои собственные вычисления перед печатью, а также запускают str(summary(m_csim)).   -  person Rui Barradas    schedule 05.11.2018

Ответы (3)


arrow_upward
1
arrow_downward

Вы, вероятно, ищете

model_summary_object <- summary(m_csim)
model_summary_object$quantiles[,c('2.5%','97.5%')]
person Pepacz    schedule 10.09.2019

arrow_upward
0
arrow_downward

Надеюсь, я не перехожу границы своих знаний, но хочу ответить с точки зрения "вообще", а не конкретно для рягов. m_csim — это объект, и для него можно использовать множество методов. Вы использовали метод сводки, чтобы увидеть что-то. Как прокомментировали люди, возможно, есть метод coef. Однако, как прокомментировал кто-то еще (пока я отвечал!), использование str() для просмотра того, что содержит объект, - лучший способ увидеть, какая информация содержится в объекте и как к ней обращаться. Я был бы очень удивлен, если бы использование str() не показало, как найти не только коэффициенты, но и достаточно информации о доверительных интервалах, чтобы позволить вам найти желаемый CI.

person meh    schedule 05.11.2018

arrow_upward
0
arrow_downward

Вероятно, вы можете извлечь квантили напрямую. Как указывали другие, вы можете вызвать str(m_csim), и вы можете ограничить вывод вызова str с помощью str(m_csim, max.level=1) и продолжать добавлять единицу к аргументу max.level=, пока не увидите что-то похожее на квантили.

Что мне нравится делать, так это преобразовывать вывод MCMC в data.frame, чтобы с ним было легче работать. Я использую jagsUI, а не rjags, но часто делаю что-то вроде

mcmc_df <- as.data.frame(as.matrix(MY_MCMC_OBJECT$samples))

Примечание: с rjags может быть немного по-другому, но я уверен, что вы сможете найти это, немного покопавшись.

Преимущество: затем я могу получить доступ к одному вектору для каждого mcmc_df$PARAMETER и создать матрицу квантилей с

mcmc_quants <- apply(mcmc_df, 2, quantile, p=c(0.025, 0.25, 0.5, 0.75, 0.975))

или любые квантили, которые вы хотите.

person Matt Tyers    schedule 07.11.2018