Полиномиальная модель MCMCglmm в R

Я пытаюсь создать модель с помощью пакета MCMCglmm в R.

Данные структурированы следующим образом, где диада, фокус, прочее - все случайные эффекты, предикторы 1-2 - переменные-предикторы, а ответы 1-5 - переменные результата, которые фиксируют количество наблюдаемых поведений различных подтипов:

 dyad focal other r    present  village  resp1 resp2 resp3 resp4 resp5 
 1    10101 14302 0.5  3        1        0     0     4     0     5
 2    10405 11301 0.0  5        0        0     0     1     0     1
 …

Итак, модель, имеющая только один результат (обучение), выглядит следующим образом:

 prior_overdisp_i <- list(R=list(V=diag(2),nu=0.08,fix=2), 
 G=list(G1=list(V=1,nu=0.08), G2=list(V=1,nu=0.08), G3=list(V=1,nu=0.08), G4=list(V=1,nu=0.08)))

 m1 <- MCMCglmm(teaching ~ trait-1 + at.level(trait,1):r + at.level(trait,1):present, 
 random= ~idh(at.level(trait,1)):focal + idh(at.level(trait,1)):other + 
 idh(at.level(trait,1)):X + idh(at.level(trait,1)):village, 
 rcov=~idh(trait):units, family = "zipoisson", prior=prior_overdisp_i, 
 data = data, nitt = nitt.1, thin = 50, burnin = 15000, pr = TRUE, pl = TRUE, verbose = TRUE, DIC  = TRUE)

В примечаниях к курсу Хэдфилда (глава 5) приводится пример полиномиальной модели, в которой используется только одна переменная результата с 3 уровнями (овечьи рога 3 типов). Аналогичное лечение можно найти здесь: http://hlplab.wordpress.com/2009/05/07/multinomial-random-effects-models-in-r/ Это не совсем подходит для того, что я делаю, но содержит полезную справочную информацию.

В другой ссылке (Hadfield 2010) приводится пример MCMCglmm с несколькими ответами, который следует тому же формату, но использует cbind () для прогнозирования вектора ответов, а не одного результата. Та же модель с несколькими ответами будет выглядеть так:

 m1 <- MCMCglmm(cbind(resp1, resp2, resp3, resp4, resp5) ~ trait-1 + 
 at.level(trait,1):r + at.level(trait,1):present, 
 random= ~idh(at.level(trait,1)):focal + idh(at.level(trait,1)):other + 
 idh(at.level(trait,1)):X + idh(at.level(trait,1)):village, 
 rcov=~idh(trait):units, 
 family = cbind("zipoisson","zipoisson","zipoisson","zipoisson","zipoisson"), 
 prior=prior_overdisp_i, 
 data = data, nitt = nitt.1, thin = 50, burnin = 15000, pr = TRUE, pl = TRUE, verbose = TRUE, DIC  = TRUE)

У меня есть два вопроса по программированию:

  1. Как указать априор для этой модели? Я просмотрел материалы, упомянутые в этом посте, но никак не могу понять.

  2. Я запустил аналогичную версию только с двумя переменными ответа, но я получил только один наклон - где я думал, что должен получить другой наклон для каждой переменной resp. В чем я ошибаюсь, или что я неправильно понял модель?


r bayesian mcmc glm
person M.A.Kline    schedule 07.10.2014    source источник
comment
Вы проверяли, действительно ли имеет смысл fix = 2 в R=list(V=diag(2),nu=0.08,fix=2)? В моем понимании предыдущей спецификации MCMCglmm fix следует читать как логическое значение: fix = 0 - значение по умолчанию, не фиксирующее отклонение до V, а fix = 1 означает исправление отклонения на значении V. Так что fix = 2 (или подобное) imo не должно иметь никакого значения. (Но на странице 103 своего курса Хэдфилд использует эту спецификацию: cran.r-project.org/pub/R/web/packages/MCMCglmm/vignettes/)   -  person Qaswed    schedule 28.11.2016
comment
@Qaswed Я возвращаюсь к этим данным через пару лет и снова смотрю на эти модели. Насколько я понимаю, компонент исправления связан с тем, для какой части модели предназначена априорная ... поскольку есть категориальный компонент (прогнозирование нулей) и непрерывный компонент (прогнозирование ненулевых значений). Это характерно для моделей Ципуассона, которые технически сами по себе полиномиальны. Предупреждение: я могу запутаться!   -  person M.A.Kline    schedule 08.03.2018

Ответы (1)


arrow_upward
7
arrow_downward

Ответ на мой первый вопрос, основанный на сообщении HLP и некоторой помощи от консультанта по работе с коллегами / статистикой:

# values for prior
k <- 5 # originally: length(levels(dative$SemanticClass)), so k = # of outcomes for SemanticClass     aka categorical outcomes 
I <- diag(k-1) #should make matrix of 0's with diagonal of 1's, dimensions k-1 rows and k-1 columns
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1)) # should make k-1 x k-1 matrix of 1's

А для моей модели, использующей семейство multinomial5 и 5 переменных результата, априор таков:

prior = list(
             R = list(fix=1, V=0.5 * (I + J), n = 4),
             G = list(
               G1 = list(V = diag(4), n = 4))

Что касается моего второго вопроса, мне нужно добавить термин взаимодействия к фиксированным эффектам в этой модели:

 m <- MCMCglmm(cbind(Resp1, Resp2...) ~ -1 + trait*predictorvariable,
 ...

Результат дает как основные эффекты для переменных ответа, так и апостериорные оценки для взаимодействия ответа / предиктора (влияние переменной-предиктора на каждую переменную ответа).

person M.A.Kline    schedule 19.11.2014