Статистика Гельмана-Рубина для модели MCMCglmm в R

У меня есть многомерная модель с этой (приблизительной) формой:

library(MCMCglmm)

mod.1 <- MCMCglmm(
    cbind(OFT1, MIS1, PC1, PC2) ~ 
    trait-1 +
    trait:sex +
    trait:date,
    random = ~us(trait):squirrel_id + us(trait):year,
    rcov = ~us(trait):units, 
    family = c("gaussian", "gaussian", "gaussian", "gaussian"),
    data= final_MCMC, 
    prior = prior.invgamma, 
    verbose = FALSE,
    pr=TRUE, #this saves the BLUPs 
    nitt=103000, #number of iterations
    thin=100, #interval at which the Markov chain is stored
    burnin=3000)

В целях публикации меня попросили сообщить статистику Гельмана-Рубина, чтобы показать, что модель сошлась.

Я пытался запустить:

gelman.diag(mod.1)

Но я получаю эту ошибку:

Error in mcmc.list(x) : Arguments must be mcmc objects

Любые предложения по правильному подходу? Я предполагаю, что ошибка означает, что я не могу передать свой вывод mod.1 через gelman.diag(), но я не уверен, что я должен вместо этого поместить туда? Мои познания здесь весьма ограничены, поэтому буду признателен за любую помощь!

Обратите внимание, что я не добавил сюда данные, но я подозреваю, что ответ больше связан с синтаксисом кода, а не с данными.


r mcmc
person Blundering Ecologist    schedule 21.05.2021    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Для gelman.diag требуется mcmc.list. Если мы запускаем модели с другим набором параметров, извлеките «Sol» и поместите его в list (ниже это та же модель).

library(MCMCglmm)
model1 <- MCMCglmm(PO~1, random=~FSfamily, data=PlodiaPO, verbose=FALSE,
  nitt=1300, burnin=300, thin=1)
model2 <- MCMCglmm(PO~1, random=~FSfamily, data=PlodiaPO, verbose=FALSE,
  nitt=1300, burnin=300, thin=1 )

mclist <- mcmc.list(model1$Sol, model2$Sol)
gelman.diag(mclist)
# gelman.diag(mclist)
#Potential scale reduction factors:
#            Point est. Upper C.I.
#(Intercept)          1          1

Согласно документации, это применимо для более чем одной цепочки mcmc.

Гельман и Рубин (1992) предлагают общий подход к мониторингу сходимости выходных данных MCMC, в котором выполняется m > 1 параллельных цепочек с начальными значениями, которые чрезмерно разбросаны по отношению к апостериорному распределению.

Вход x здесь

x — объект mcmc.list с более чем одной цепочкой и начальными значениями, чрезмерно рассредоточенными по отношению к апостериорному распределению.

person akrun    schedule 21.05.2021
comment
Если я использую только 1 модель MCMCglmm(), не будет ли плохой практикой просто перейти: mclist <- mcmc.list(model1$Sol, model1$Sol)? - person Blundering Ecologist; 21.05.2021
comment
@BlunderingEcologist, согласно документации, ясно указывает, что ввод x - An mcmc.list object with more than one chain, and with starting values that are overdispersed with respect to the posterior distribution. - person akrun; 21.05.2021
comment
@BlunderingEcologist, вероятно, их запрос может состоять в том, чтобы запустить несколько цепочек для проверки конвергенции. - person akrun; 21.05.2021
comment
Я совершенно не в своей тарелке, чтобы понять документацию. Возможно, это ваш дополнительный комментарий (например, запуск второй модели с разными значениями nitt, burnin и thin?). - person Blundering Ecologist; 21.05.2021
comment
@BlunderingEcologist документация и примеры кода ограничены. Возможно, вам придется прочитать все описание, особенно их виньетки, если они доступны. - person akrun; 21.05.2021
comment
@BlunderingEcologist, может быть, лучше задать теоретический вопрос в перекрестной проверке по этому поводу. Это может помочь вам лучше - person akrun; 21.05.2021
comment
Вам нужно запустить несколько цепочек с разными значениями начальных параметров, чтобы диагностика Гельмана-Рубина имела какой-либо смысл (изменение выгорания/количества итераций/и т. д. не имеет значения) . MCMCglmm имеет аргумент start, поэтому это должно быть выполнимо (но вам нужно будет выяснить/решить, какой диапазон начальных значений подходит для вашей проблемы). - person Ben Bolker; 21.05.2021
comment
@BenBolker Спасибо, что рассказали об этом! Я должен буду прочитать больше о том, как принять решение о начале в этом случае. - person Blundering Ecologist; 21.05.2021