[R] Specifying priors in a multi-response MCMCglmm

[Bert Gunter]

1. (Mainly) Statistical issues are generally off topic on this list.
You might want to try the r-sig-mixed-models list instead.

2. However, I think a better answer is to seek local statistical
expertise in order to have an extended discussion about your research
intent in order to avoid producing yet more irreproducible
psychological research.

Cheers,
Bert
Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Tue, May 1, 2018 at 2:15 PM, Michelle Kline
<michelle.ann.kline using gmail.com> wrote:
> 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]]
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.