MCMCglmm биномиальная модель до

Question

Я хочу оценить биномиальную модель с пакетом RMCMCglmm. Модель должна включать перехват и наклон - как фиксированные, так и случайные. Как мне указать принятый prior? (Обратите внимание, here is a similar question, но в гораздо более сложной обстановке.)

Предположим, что данные имеют следующий вид:

y   x cluster 
1 0 -0.56047565  1 
2 1 -0.23017749  1 
3 0 1.55870831  1 
4 1 0.07050839  1 
5 0 0.12928774  1 
6 1 1.71506499  1 

---

В самом деле, данные были получены с помощью

set.seed(123) 
nj <- 15 # number of individuals per cluster 
J <- 30 # number of clusters 
n <- nj * J 
x <- rnorm(n) 
y <- rbinom(n, 1, prob = 0.6) 
cluster <- factor(rep(1:nj, each = J)) 

dat <- data.frame(y = y, x = x, cluster = cluster) 

Answer 1

Информация в вопросе о модели, предлагаемая для указания fixed = y ~ 1 + x и random = ~ us(1 + x):cluster. С us() вы позволить случайные эффекты коррелируют (раздел 3.4 см и в таблице 2, в Hadfield's 2010 jstatsoft-article)

Прежде всего, как вы есть только один зависимую переменную (y), в G участие в предшествующем уровне (ср уравнение 4 и раздел 3.6 в Hadfield's 2010 jstatsoft-article) для дисперсии (ов) случайных эффектов должен иметь только один элемент списка, называемый G1. Этот элемент списка не является фактическим предыдущим распределением - это было указано Hadfield как inverse-Wishart distribution. Но с G1 вы задаете параметры этого распределения обратного Whishart, которые являются масштаб матрица ( в Википедии и обозначениях V в MCMCglmm нотации) и степени свободы ( в Википедии нотации и nu в MCMCglmm обозначениях). Поскольку у вас есть два случайных эффекта (перехват и наклон), V должен быть матрицей 2 x 2. Частым выбором является двумерная тождественная матрица diag(2). Хэдфилда часто использует nu = 0.002 для степеней свободы (см his course notes)

Теперь, вы также должны указать R части в предшествующем для остаточной дисперсии. Здесь снова Haidfield задает обратное-Whishart-распределение, оставляя пользователя указывать его параметры. Поскольку мы имеем только одну остаточную дисперсию, V должен быть скаляром (скажем, V = 0.5). Дополнительным элементом для R является fix. С помощью этого элемента вы указываете, должна ли остаточная дисперсия фиксироваться на определенное значение (чем вы должны писать fix = TRUE или fix = 1) или нет (тогда fix = FALSE или fix = 0). Обратите внимание, что вы не фиксируете остаточную дисперсию как 0.5 по fix = 0.5! Поэтому, когда вы найдете в примечаниях к курсу Хадфилда fix = 1, прочитайте его как fix = TRUE и посмотрите, какое значение V было исправлено.

Все togehter мы создали предшествующее следующим образом:

prior0 <- list(G = list(G1 = list(V = diag(2), nu = 0.002)), 
       R = list(V = 0.5, nu = 0.002, fix = FALSE)) 

С этого до мы можем запустить MCMCglmm:

library("MCMCglmm") # for MCMCglmm()  
set.seed(123) 

mod0 <- MCMCglmm(fixed = y ~ 1 + x, 
      random = ~ us(1 + x):cluster, 
      data = dat, 
      family = "categorical", 
      prior = prior0) 

Розыгрыши из Гиббса-пробоотборник для фиксированных эффектов найдено в mod0$Sol, ничьи для параметров дисперсии в mod0$VCV.

Обычно бином модель требует остаточной дисперсии, чтобы быть фиксированной, поэтому мы устанавливаем остаточная дисперсия будет зафиксирована на уровне 0.5

set.seed(123) 

prior1 <- list(G = list(G1 = list(V = diag(2), nu = 0.002)), 
      R = list(V = 0.5, nu = 0.002, fix = TRUE)) 

mod1 <- MCMCglmm(fixed = y ~ 1 + x, 
      random = ~ us(1 + x):cluster, 
      data = dat, 
      family = "categorical", 
      prior = prior1) 

Различие можно видеть при сравнении с mod0$VCV[, 5]mod1$VCV[, 5]. В последнем случае все записи: 0.5, как указано.