MCMCglmm wielomianu model R

Question

Próbuję stworzyć model przy użyciu pakietu MCMCglmm w R.

Dane są pogrupowane w następujący sposób, gdzie dyad, ogniskowe, inne są wszystkie losowe efekty, predict1-2 są predyktorem zmienne i odpowiedź 1-5 są zmiennymi rezultatu że capture # obserwowanych zachowań różnych podtypów:

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 
… 

więc model z jednym tylko wynik (nauczanie) przedstawia się następująco:

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) 

Notatki z kursu Hadfielda (rozdział 5) podają przykład modelu wielomianowego, który wykorzystuje tylko jedną zmienną wyniku z 3 poziomami (rogi owcze 3 typów). Podobne leczenie można znaleźć tutaj: http://hlplab.wordpress.com/2009/05/07/multinomial-random-effects-models-in-r/ To nie jest w porządku dla tego, co robię, ale zawiera przydatne informacje podstawowe.

Kolejne odniesienie (Hadfield 2010) podaje przykład MCMCglm o wielu odpowiedziach, który ma ten sam format, ale wykorzystuje metodę cbind() do przewidywania wektora odpowiedzi zamiast pojedynczego wyniku. Ten sam model z wieloma odpowiedziami będzie wyglądać następująco:

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) 

Mam dwa pytania jest programowanie tutaj:

1. Jak mogę określić przed dla tego modelu? Przyjrzałem się materiałom wymienionym w tym poście, ale po prostu nie mogę tego rozgryźć.

2. Uruchomiłem podobną wersję tylko z dwiema zmiennymi odpowiedzi, ale otrzymałem tylko jedno nachylenie - gdzie sądziłem, że powinienem uzyskać inne nachylenie dla każdej zmiennej resp. Gdzie się mylę lub źle zrozumiałem model?

Answer 1

odpowiedź na moje pierwsze pytanie, w oparciu o stanowisko HLP i trochę pomocy z Colleage/status konsultanta:

# 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 

a dla mojego modelu, z wykorzystaniem multinomial5 rodzinnych i 5 wyników leczenia zmiennych, przed brzmi:

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

na moje drugie pytanie, muszę dodać określenie interakcji z trwałych efektów w tym modelu:

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

Wynik daje zarówno główne efekty dla zmiennych odpowiedzi, jak i oszacowania a posteriori dla interakcji między reakcjami/predyktorami (wpływ zmiennej predykcyjnej na każdą zmienną odpowiedzi).