[R] Specifying priors in a multi-response MCMCglmm

Michelle Kline michelle@@nn@kline @ending from gm@il@com
Tue May 1 23:15:06 CEST 2018


Hi all,

I previously emailed about a multinomial model, and after seeking some
additional help, realized that since my response/outcome variables are not
mutually exclusive, I need to use a multi-response model that is *not*
multinomial. I'm now trying to figure out how to specify the priors on the
multi-response model. Any help would be much appreciated.

 My data look like this:

     X other focal village present           r teaching Opp_teacher
Dir_teacher Enh_teacher SocTol_teacher Eval_teacher Total_teacher
f_Age f_Ed Age Ed1   61 10202 10213       0      15 0.250000000
2           0           0           0              0            2
       2     1    0  48  82   63 10203 10213       0      19
0.500000000        6           0           0           4
0            6            10     1    0  27 103   64 10204 10213
0       1 0.250000000        0           0           0           0
         0            0             0     1    0  25  94   69 10206
10213       0       6 0.250000000        2           0           0
      1              0            1             2     1    0  20 115
72 10207 10213       0       4 0.250000000        0           0
   0           0              0            0             0     1    0
18  86   80 10210 10213       0       4 0.250000000        0
0           0           0              0            0             0
 1    0  30 127   83 10211 10213       0       8 0.062500000        0
         0           0           0              0            0
    0     1    0  73  38   85 10212 10213       0      11 0.125000000
      8           0           1           1              0
8            10     1    0   9  19  132 10403 10213       0       1
0.000976563        0           0           0           0
0            0             0     1    0  10  010 241 11703 10213
0       3 0.015625000        1           0           0           0
         0            1             1     1    0  49  8

Columns Opp_teacher through Eval_Teacher are count data different possible
teaching behaviors that I have observed, with each row being a dyad. The
teaching types are not mutually exclusive. They can co-occur. This is why I
am using a multi-response model but not a multi-nomial model. Focals as
well as others can appear in more than one dyad, so I have included those
as random effects. The fixed effects in the model are r (relatedness) and
present (# observations together).

I've specified my model as follows:

m3.random.present.r <- MCMCglmm(cbind(Opp_teacher , Dir_teacher,
Enh_teacher, SocTol_teacher, Eval_teacher) ~ +present + r + trait -1,
               random = ~ other + focal,
               prior = prior.m3,
               burnin = burn,
               nitt = iter,
               family =c("poisson","poisson","poisson","poisson","poisson"),
               data = data,
               pr=TRUE,
               pl=TRUE,
               DIC = TRUE,
               verbose = FALSE)

The prior, prior.m3 is as follows:

prior.m3 <- list(R = list(V = diag(2), nu = 2),
                     G = list(G1 = list(V = diag(2), nu = 5),
                              G2 = list(V = diag(2), nu = 5),
                              G3 = list(V = diag(2), nu = 5),
                              G4 = list(V = diag(2), nu = 5),
                              G5 = list(V = diag(2), nu = 5)))

This is based on Hadfield's Course Notes, as well as some advice found in this
post
<https://stackoverflow.com/questions/40617099/mcmcglmm-binomial-model-prior>.
It's consistent with how I've specified priors for simpler models (with
single outcome variables), but I am obviously missing something that must
change with respect to the G-structures when using multiple responses,
because running the model results in the following error:

Error in MCMCglmm(cbind(Opp_teacher, Dir_teacher, Enh_teacher,
SocTol_teacher, : prior$G has the wrong number of structures

I am not sure what this error message refers to. My understanding is that
there should be 5 G-structures listed because I have 5 dependent variables.
(Trial & error suggests this isn't the meaning of the error message - a
different number of G-structures does not change the result). This suggests
the problem has to do with the rest of the G-structure code: I've set `V =
diag(2)` because there are two random effects.

I can't come up with any other rationale, despite having scoured the
internet for additional help.
Thanks,

Michelle


-- 
Michelle A. Kline, PhD

Assistant Professor
Department of Psychology
Simon Fraser University

	[[alternative HTML version deleted]]




More information about the R-help mailing list