MCMCglmm mô hình đa thức R

Question

Tôi đang cố gắng để tạo ra một mô hình sử dụng gói MCMCglmm trong R.

Các dữ liệu được cấu trúc như sau, trong đó hai vật giống nhau, tiêu cự, khác là tất cả các hiệu ứng ngẫu nhiên, predict1-2 là yếu tố dự báo biến, và phản ứng 1-5 là các biến kết quả mà chụp # hành vi quan sát của phân nhóm khác nhau:

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 
… 

vì vậy, một mô hình chỉ với một kết quả (giảng dạy) như sau:

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) 

Ghi chú khóa học của Hadfield (Ch 5) đưa ra một ví dụ về mô hình đa thức chỉ sử dụng một biến kết quả duy nhất với 3 cấp độ (sừng cừu của 3 loại). Điều trị tương tự có thể tìm thấy ở đây: http://hlplab.wordpress.com/2009/05/07/multinomial-random-effects-models-in-r/ Điều này không hoàn toàn phù hợp với những gì tôi đang làm, nhưng chứa thông tin cơ bản hữu ích.

Một tham chiếu khác (Hadfield 2010) đưa ra ví dụ về MCMCglmm đa phản hồi theo cùng định dạng nhưng sử dụng cbind() để dự đoán vectơ phản hồi, chứ không phải là một kết quả. Cùng một mô hình với nhiều phản ứng sẽ trông như thế này:

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) 

Tôi có hai câu hỏi lập trình ở đây:

1. Làm thế nào để xác định một trước cho mô hình này? Tôi đã xem xét các tài liệu được đề cập trong bài viết này nhưng không thể tìm ra nó.

2. Tôi đã chạy phiên bản tương tự chỉ với hai biến trả lời, nhưng tôi chỉ có một độ dốc - nơi tôi nghĩ tôi sẽ có độ dốc khác nhau cho mỗi biến số resp. Tôi đang đi sai ở đâu, hoặc tôi đã hiểu nhầm mô hình?

Answer 1

trả lời cho câu hỏi đầu tiên của tôi, dựa trên bài HLP và một số sự giúp đỡ từ một colleage/số liệu thống kê tư vấn:

# 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 

Và đối với mô hình của tôi, sử dụng multinomial5 gia đình và 5 kết quả biến, trước là:

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

Đối với câu hỏi thứ hai của tôi, tôi cần phải thêm một thuật ngữ tương tác với các tác động cố định trong mô hình này:

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

Kết quả cho cả hai hiệu ứng chính cho biến Response và ước tính sau cho tương tác Response/Predictor (tác dụng của biến dự báo trên mỗi biến phản hồi).