MCMCglmm modello multinomiale in R

Question

Sto cercando di creare un modello utilizzando il pacchetto MCMCglmm in R.

I dati sono strutturati come segue, dove diade, focali, altri sono tutti effetti casuali, predict1-2 sono predittore variabili, e la risposta 1-5 sono variabili di risultato che catturano ° comportamenti osservati su diversi sottotipi:

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 
… 

Quindi, un modello con un solo risultato (insegnamento) è la seguente:

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) 

Le note del corso di Hadfield (Cap 5) forniscono un esempio di un modello multinomiale che utilizza solo una singola variabile di risultato con 3 livelli (corna di pecora di 3 tipi). Un trattamento simile può essere trovato qui: http://hlplab.wordpress.com/2009/05/07/multinomial-random-effects-models-in-r/ Questo non è giusto per quello che sto facendo, ma contiene informazioni di base utili.

Un altro riferimento (Hadfield 2010) fornisce un esempio di MCMCglmm a più risposte che segue lo stesso formato ma utilizza cbind() per prevedere un vettore di risposte, piuttosto che un singolo risultato. Lo stesso modello con più risposte sarebbe simile a questa:

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) 

Ho due domande di programmazione qui:

1. Come faccio a specificare una preventiva per questo modello? Ho esaminato i materiali citati in questo post ma non riesco a capirlo.

2. Ho eseguito una versione simile con solo due variabili di risposta, ma ottengo solo una pendenza - dove ho pensato di ottenere una pendenza diversa per ciascuna variabile resp. Dove sto andando male o ho frainteso la modella?

Answer 1

risposta alla mia prima domanda, sulla base del palo HLP e qualche aiuto da un colleage/statistiche consulente:

# 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 

E per il mio modello, utilizzando i multinomial5 famiglia e 5 di risultato variabili, le prima è:

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

per la mia seconda domanda, ho bisogno di aggiungere un termine di interazione per gli effetti fissi in questo modello:

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

Il risultato fornisce entrambi gli effetti principali per le variabili di risposta e le stime posteriori per l'interazione risposta/predittore (l'effetto della variabile predittore su ciascuna variabile di risposta).