[R-sig-ME] Running multinomial models with random effects

[MORAN LOPEZ, TERESA]

Dear all,
I am a phD student working on animal behavior under different  
predation risk. We have conducted an experiment in which acorn  
selection by jays was evaluated in savanna vs forest-type  landscapes.  
  We placed 8 feeders in two different forests and 8 feeders in two  
different savannas.
My response variable is factorial, acorn choice (big, small, both).   
In my design I have spatial pseudorreplication within Areas.  I have  
not found multinom function which allows fitting random effects.  
However, MCMCglmm package allows incorporating them using Bayesian  
Methods. I am new using Bayesian Methods so it has taken me a lot of  
effort to partially understand these models implementation. However,  
once I have run my model I am not sure about output interpretation and  
most importantly how to convert fix effect parameters into  
probabilities.

I have tried three different things: to convert my data in binomial  
(ignoring  infrequent acorn choice “both”) and run mixed model lmer.   
Pros: I can include random effects , Cons: I lose an infrequent but  
informative level of the response.

  Keep categorical responses and run multinom function without random  
effects.  Pros: I can fit my response variable as categorical and keep  
all levels. Cons: I am ignoring spatial autocorrelation (however, when  
running binomial models acorn choice variation between areas was much  
more lower than between risk levels).
  +
Run MCMCglmm model with random effects. Pros:  I can keep both spatial  
autocorrelation and all levels of the response factor. Cons: Not sure  
about  model implementation and interpretation. Low sample size.


Which option is the most accurate?

Thanks a lot! I am really stucked with my data.

Here I post my MCMCglmm model which is the one I have more doubts

#head(choice)-dataset with choice as factor  (b=big, s=small, bt=both)
    RISK       Area             CHOICE
LOW        Sopie                  b
  LOW       Sopie                  b
  LOW      Sopie                   b
  LOW  Anchurones                  s
  LOW  Anchurones                  bt
  LOW  Anchurones                  b

  I have read package tutorial and Hadfield notes (though I have not  
read all the chapter of Hadfield course notes). Besides I have used  
information posted  Jaeger lab blog:  
http://hlplab.wordpress.com/2009/05/07/multinomial-random-effects-models-in-r/

I am not very sure about random effects priors but I believe the rest  
of them are ok.
k <- length(levels(choice$CHOICE))
I<- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
#Constraints to variance-covariance matrix

prior = list(R = list(fix=1, V =(1/k) * (I + J), n=k), G = list(G1 =  
list(V =diag(2), n=2)))

#R priors, Hadfield manual suggests that structure. Are they right?   
NOT SURE ABOUT RANDOM PRIORS.

M<- MCMCglmm(CHOICE~ -1+ trait + RISK,
               random = ~ us(trait):Area,
               rcov = ~ us(trait):units,
               prior = prior,
               family = "categorical",
               data = choice)

#-1+trait following Hadfield suggestions in order to estimate  
interception if every outcome.
us(trait):-----  since we are not sure of independence assumptions. Is  
model function allright?

My model runs but I find it very difficult to interpret my summary

summary(M)

Iterations = 3001:12991
  Thinning interval  = 10
  Sample size  = 1000

  DIC: 44.21802

  G-structure:  ~us(trait):Area

                          post.mean l-95% CI u-95% CI eff.samp
CHOICE.bt:CHOICE.bt.Area    14.622   0.1881    64.46    7.278
CHOICE.s:CHOICE.bt.Area      3.274 -11.1431    19.29  116.593
CHOICE.bt:CHOICE.s.Area      3.274 -11.1431    19.29  116.593
CHOICE.s:CHOICE.s.Area       5.071   0.1416    15.82  127.395

  R-structure:  ~us(trait):units

                           post.mean l-95% CI u-95% CI eff.samp
CHOICE.bt:CHOICE.bt.units    0.6667   0.6667   0.6667        0
CHOICE.s:CHOICE.bt.units     0.3333   0.3333   0.3333        0
CHOICE.bt:CHOICE.s.units     0.3333   0.3333   0.3333        0
CHOICE.s:CHOICE.s.units      0.6667   0.6667   0.6667        0

  Location effects: CHOICE ~ -1 + trait + RISK

                post.mean l-95% CI u-95% CI eff.samp pMCMC
traitCHOICE.bt   -4.5042 -11.8926   0.3824    15.18 0.052 .
traitCHOICE.s    -2.0187  -6.3370   1.3569    84.55 0.200
RISKRAÑA          1.9095  -2.3541   8.2421    56.13 0.336
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

posterior.mode(M$Sol)

traitCHOICE.bt   traitCHOICE.s       RISKRAÑA
      -2.594132      -1.836527             1.554387

#Since there is not any intercept I can´t apply plogis() in order to  
get probabilities of choice in different areas.


Sorry that I have that many  questions and thanks a lot!!!