Mcmcglmm мультиномиальная модель в R

Я пытаюсь создать модель, используя пакет MCMCglmm в R.

Данные структурированы следующим образом, где dyad, focal, other-все случайные эффекты, predict1 - 2-переменные предиктора, а response 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)
В примечаниях к курсу Хэдфилда (Ch 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. Где я ошибаюсь, или я неправильно понял модель?

1 7
glm bayesian

1 ответ:

Ответ на мой первый вопрос, основанный на посте 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,
 ...

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

7