Skip to content

MCMCglmm: Fixing the priors in multivariate response models without random effects

3 messages · Iker Vaquero Alba, Jarrod Hadfield

Original
Iker Vaquero Alba
# 
?? Hello all,
?? I have managed to pretty much understand the MCMCglmm function, at least to the point of being able to write a model with just a single response categorical variable and run it. It can be done without any need to specify any priors. However, when I try to run a more complicated model with a bivariate response, the problems start. 
?? This is the model I am trying to run and the error message I get:
testmodel1<-MCMCglmm(cbind(natapshort,nataplong)~gender+age+religion+sexor+selfattr+partnerattr+gender:age+gender:religion+gender:sexor+gender:selfattr+gender:partnerattr+age:religion+age:sexor+age:selfattr+age:partnerattr+religion:sexor+religion:selfattr+religion:partnerattr+sexor:selfattr+sexor:partnerattr+selfattr:partnerattr,random=NULL,rcov=~us(trait):units,family=c("categorical","categorical"),data=extphen,nitt=10000,singular.ok=TRUE)

????????????????????? MCMC iteration = 0

? Acceptance ratio for latent scores = 0.000154

????????????????????? MCMC iteration = 1000

? Acceptance ratio for latent scores = 0.211168
Error in MCMCglmm(cbind(natapshort, nataplong) ~ gender + age + religion +? : 
? Mixed model equations singular: use a (stronger) prior


So I need to define my own priors. As I am still in the testing stage, I tried some simple ones found in the CourseNotes. Also from the CourseNotes and other sources, I understand that the term G refers to random effects, so I should not include it. Then I should include B and R, am I right? I haven't been able to find very clear information about what G, B and R refer to.
?? I have tried this:? 
?? prior<- list(B= list(B1 = list(V = diag(2), n = 1.002)),R = list(V = diag(2), n = 1.002))
??? And this is what I get:
?? testmodel1<-MCMCglmm(cbind(natapshort,nataplong)~gender+age+religion+sexor+selfattr+partnerattr+gender:age????? +gender:religion+gender:sexor+gender:selfattr+gender:partnerattr+age:religion+age:sexor+age:selfattr+age:partnerattr????? +religion:sexor+religion:selfattr+religion:partnerattr+sexor:selfattr+sexor:partnerattr+selfattr:partnerattr,random=NULL,????? rcov=~us(trait):units,family=c("categorical","categorical"),data=extphen,nitt=10000,prior=prior,singular.ok=TRUE)
Error in priorformat(if (NOpriorG) { : 
? V is the wrong dimension for some prior$G/prior$R elements

After getting this error, I have tried adding another five B terms (as I have 6 explanatory variables), but the result is tha same. 


If I try just with G and B:
prior<- list(G = list(G1 = list(V = diag(2), n = 1.002)),B = list(V = diag(2), n = 1.002))
Error in MCMCglmm(cbind(natapshort, nataplong) ~ gender + age + religion +? : 
? either both or neither R and G structures need a prior

?? I am getting crazy. Could anybody shed some light on the priors' issue and help me find some simple ones that can make my model work? I don't even know where to look any more, I have read plenty of sources and documents, but I haven't got any clear conclusion yet.
?? Thank you very much.?? Best wishes,?? Iker
 ?__________________________________________________________________

?? Iker Vaquero-Alba
?? Visiting Postdoctoral Research Associate
?? Laboratory of Evolutionary Ecology of Adaptations 
?? Joseph Banks Laboratories
?? School of Life Sciences
?? University of Lincoln?? Brayford Campus, Lincoln
?? LN6 7DL
?? United Kingdom

?? https://eric.exeter.ac.uk/repository/handle/10036/3381
2 days later
Jarrod Hadfield
# 
Hi,

The issue is that both outcomes are binary and you are trying to  
estimate an unstructured residual covariance matrix. The diagonal  
elements (the variances) are not identifiable and so need to be  
constrained. The simplest method is to constrain the matrix to a  
correlation matrix using corg(trait):units.

Its hard to say without knowing what the data are, but I would think  
you need to fit trait in the fixed effect part of the model together  
with an interaction between trait and other predictors.

Also, I would recommend using family="threshold" rather than  
family="categorical" for bivariate problems. I've given the reasons  
for this in older posts. For example, the probit section of:

https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q1/021875.html

Regarding the correct dimension of the prior for the fixed effects, B  
should be equal to the number of fixed effects fitted. I can't see how  
many you have, but definitely more than 2: it looks closer to 20.

Cheers,

Jarrod





Quoting Iker Vaquero Alba <karraspito at yahoo.es> on Fri, 18 Sep 2015  
22:26:20 +0000 (UTC):

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Iker Vaquero Alba
    

    

      
      #
    

  


  

      
?

     De: Jarrod Hadfield <j.hadfield at ed.ac.uk>
 Para: Iker Vaquero Alba <karraspito at yahoo.es> 
CC: "r-sig-mixed-models at r-project.org" <r-sig-mixed-models at r-project.org> 
 Enviado: Lunes 21 de septiembre de 2015 11:48
 Asunto: Re: [R-sig-ME] MCMCglmm: Fixing the priors in multivariate response models without random effects
   
Hi,

The issue is that both outcomes are binary and you are trying to? 
estimate an unstructured residual covariance matrix. The diagonal? 
elements (the variances) are not identifiable and so need to be? 
constrained. The simplest method is to constrain the matrix to a? 
correlation matrix using corg(trait):units.
Dear Jarrod, 

?? Thank you so much for your invaluable help. I'm not sure I understand why "both outcomes are binary". My two dependent variables, as well as some of the explanatory ones, are multinomial, at least as far as I can understand: they are valuations from 1 to 5 given to different attributes by a group of volunteers. 


Its hard to say without knowing what the data are, but I would think? 
you need to fit trait in the fixed effect part of the model together? 
with an interaction between trait and other predictors.
The data are as follows: 
Response variables: 
natapshort and nataplong: both are categorical and multinomial variables, consisting on values from 1 to 5, depending on the importance several volunteers give to natural appearance on their potential partners, either for a short term relationship or for a longer commitment.Explanatory variables:
gender: factor with 3 levels (male, female, other). Gender of the participants.age: factor with 6 levels. Age of the participants, classified by age ranges.religion: factor with 2 levels (yes or no).sexor: factor with 3 levels. Sexual orientation of the participants.selfattr: factor with 5 levels. Value (from 1 to 5) that each participants gives to their own attractiveness.partnerattr: factor with 5 levels. Value (from 1 to 5) that each participant considers as the minimum attractiveness for a person to be considered as a potential partner.

When you say I need to "fit trait in the fixed effect part of the model together with an interaction between trait and other predictors.", you mean physically include the term "trait" in the model, as I have seen in many places, or you mean just including 2-way interactions between explanatory variables (that has been done already)??
Also, I would recommend using family="threshold" rather than? 
family="categorical" for bivariate problems. I've given the reasons? 
for this in older posts. For example, the probit section of:

https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q1/021875.html
Ok, thank you very much for that, I'll have a look to understand the reasons.

 Regarding the correct dimension of the prior for the fixed effects, B? should be equal to the number of fixed effects fitted. I can't see how? 
many you have, but definitely more than 2: it looks closer to 20.
Cheers,

Jarrod
Quoting Iker Vaquero Alba <karraspito at yahoo.es> on Fri, 18 Sep 2015? 
22:26:20 +0000 (UTC):

Just another question: should it be equal to the number of single fixed effects fitted (six in this case), or each two-way interaction (for example, gender:selfattr) qualifies as a new fixed effect (in which case I would have 21)?
Also, I have just read in a very useful document I've found from Tufts University, that when writing the priors, you need to fit an R structure for each fixed effect and a G structure for each random effect. According to that, I would have to fit 21 (6?) R structures and no G structures at all. But you say I have to fit 21 B structures as well. So, if R is the structure of fixed effects and G the structure of random effects, what is B?
Thank you very much again for your patience and huge help, and sorry for the endless questions.
Best wishes,Iker

__________________________________________________________________

?? Iker Vaquero-Alba
?? Visiting Postdoctoral Research Associate
?? Laboratory of Evolutionary Ecology of Adaptations 
?? Joseph Banks Laboratories
?? School of Life Sciences
?? University of Lincoln?? Brayford Campus, Lincoln
?? LN6 7DL
?? United Kingdom

?? Animal sexual signals: Do they maximise or optimise information content?