Skip to content

Specifying priors in a multi-response MCMCglmm

5 messages · Michelle Kline, Bert Gunter, Rolf Turner

#
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
#
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 at gmail.com> wrote:
#
Hi Bert,

That was distinctly unhelpful, and your outward hostility to a field you
obviously don't understand reveals a regrettable level of ignorance.

By the way, my research is Anthropology despite my job title.

Michelle
On Tue, May 1, 2018 at 2:48 PM, Bert Gunter <bgunter.4567 at gmail.com> wrote:

            

  
    
#
On 02/05/18 09:53, Michelle Kline wrote:

            
Not if you actually follow Bert's advice.
I didn't see any hostility to any field.  Bert, like many of us, objects 
to people blithely and arrogantly applying possibly deep statistical 
techniques of which they apparently have little understanding, with an 
attitude of "Oh, it's just statistics.  How hard can it be?"

I don't think anyone who knows anything about anything would classify 
Bert as ignorant.
OK.  Your job title is misleading.  Let's say "in order to avoid 
producing yet more irreproducible anthropological research" then.

cheers,

Rolf Turner
#
Okay, you've all successfully convinced me to leave this list. Bye!
On Tue, May 1, 2018 at 3:13 PM, Rolf Turner <r.turner at auckland.ac.nz> wrote: