[R-sig-ME] Interpreting a MCMCglmm model with a bivariate response variable

Jarrod Hadfield j.hadfield at ed.ac.uk
Wed Sep 23 16:38:47 CEST 2015


Hi,

a) Section 8 of the course notes details the use of 'triat'.

b) DIC is not implemented for multivariate threshold models. DIC, as  
its currently `focussed' is not useful anyway, I think.

Cheers,

Jarrod




The muklQuoting Iker Vaquero Alba <karraspito at yahoo.es> on Wed, 23 Sep  
2015 14:32:00 +0000 (UTC):

>
>    Hi, Jarrod and everyone else.
>    Thank you very much for your advise. Actually, I just realized  
> after sending the email, I was going to send the correct one. The  
> model:
> 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=~corg(trait):units,family=c("threshold","threshold"),data=extphen,nitt=100000,singular.ok=TRUE)
>  And the summary outcome:
>  summary(testmodel1)
>
>  Iterations = 3001:99991
>  Thinning interval  = 10
>  Sample size  = 9700
>
>  DIC:
>
>  R-structure:  ~corg(trait):units
>
>                             post.mean l-95% CI u-95% CI eff.samp
> natapshort:natapshort.units    1.0000   1.0000   1.0000        0
> nataplong:natapshort.units     0.3469   0.2243   0.4744     7885
> natapshort:nataplong.units     0.3469   0.2243   0.4744     7885
> nataplong:nataplong.units      1.0000   1.0000   1.0000        0
>
>  Location effects: 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
>
>                        post.mean   l-95% CI   u-95% CI eff.samp   pMCMC  
> (Intercept)            3.170e+00 -1.940e+00  8.301e+00     8794 0.21897  
> genderM               -1.048e-01 -2.006e+00  1.915e+00     9408 0.90948  
> genderO               -5.316e+02 -1.827e+05  1.854e+05     9864 0.99052  
> age                   -1.006e+00 -3.414e+00  1.266e+00     8847 0.39959  
> religionY             -2.823e+00 -5.121e+00 -5.111e-01     8322 0.01691 *
> sexorHOM               1.338e+03 -1.663e+05  1.738e+05     9700 0.99608  
> sexorOT                8.310e+02 -1.580e+05  1.586e+05     9700 0.99732  
> selfattr              -8.544e-01 -2.236e+00  5.512e-01     9092 0.22454  
> partnerattr            4.296e-01 -9.066e-01  1.784e+00     8818 0.52598  
> genderM:age            3.254e-02 -4.556e-01  5.027e-01     9700 0.88557  
> genderO:age            2.113e+02 -7.761e+04  7.336e+04     9700 0.99918  
> genderM:religionY     -4.660e-01 -1.141e+00  2.360e-01     9503 0.17258  
> genderO:religionY      3.408e+02 -1.551e+05  1.567e+05     9700 0.99485  
> genderM:sexorHOM      -7.522e-01 -1.687e+00  1.528e-01     9700 0.11155  
> genderO:sexorHOM       8.643e+01 -1.468e+05  1.512e+05     9700 0.99258  
> genderM:sexorOT       -5.950e-01 -1.661e+00  4.325e-01     9326 0.26742  
> genderO:sexorOT        6.955e+02 -1.589e+05  1.664e+05     9700 0.99340  
> genderM:selfattr       3.038e-01 -1.054e-01  6.976e-01     9334 0.13526  
> genderO:selfattr      -1.968e+02 -7.999e+04  7.118e+04     9341 0.99340  
> genderM:partnerattr   -2.686e-01 -6.974e-01  1.595e-01     9345 0.20866  
> genderO:partnerattr    1.064e+00 -1.365e+00  3.685e+00     8894 0.41155  
> age:religionY          3.473e-01 -3.822e-01  1.100e+00     9215 0.34763  
> age:sexorHOM          -6.704e+02 -8.690e+04  8.316e+04     9700 0.99608  
> age:sexorOT           -4.160e+02 -7.930e+04  7.899e+04     9700 0.99732  
> age:selfattr           3.904e-01 -2.322e-01  9.908e-01     9257 0.20351  
> age:partnerattr       -5.805e-02 -6.213e-01  5.200e-01     8974 0.85155  
> religionY:sexorHOM     3.131e-01 -8.014e-01  1.381e+00     9700 0.57505  
> religionY:sexorOT     -3.484e-03 -1.612e+00  1.651e+00     7850 0.99278  
> religionY:selfattr     1.074e-01 -3.688e-01  5.790e-01     8506 0.65649  
> religionY:partnerattr  5.330e-01  1.455e-01  8.974e-01     8817 0.00495 **
> sexorHOM:selfattr      3.787e-01 -3.850e-01  1.134e+00     9101 0.31588  
> sexorOT:selfattr       3.805e-01 -2.086e-01  9.754e-01     9163 0.20866  
> sexorHOM:partnerattr   4.977e-01 -2.861e-01  1.323e+00     8472 0.22639  
> sexorOT:partnerattr   -7.482e-02 -6.965e-01  5.284e-01     9023 0.80000  
> selfattr:partnerattr   5.051e-02 -1.438e-01  2.540e-01     8700 0.62041  
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
>
>  Cutpoints:
>                            post.mean l-95% CI u-95% CI eff.samp
> cutpoint.traitnatapshort.1     1.024   0.7313    1.335    872.3
> cutpoint.traitnatapshort.2     2.362   2.0725    2.698    756.2
> cutpoint.traitnatapshort.3     4.002   3.6299    4.396    884.4
> cutpoint.traitnataplong.1      1.187   0.9010    1.489    828.4
> cutpoint.traitnataplong.2      2.496   2.1718    2.812    727.6
> cutpoint.traitnataplong.3      3.959   3.5890    4.350    872.9
> Some additional questions, if that's fine:
> 1. Do I have to actually include the word 'trait' in the model,  
> something like 'cbind(natapshort, nataplong) ~ trait + gender +  
> age..."? What is the reason for that? Also, do I have to include  
> interactions between 'trait' and other variables?
> 2. Why the summary of my model does not give me the DIC value? It  
> appears empty.
>
> Thank you very much.
> 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
>
>
>       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: Miércoles 23 de septiembre de 2015 15:25
>  Asunto: Re: [R-sig-ME] Interpreting a MCMCglmm model with a  
> bivariate response variable
>
> Hi Iker,
>
> You need to follow the advice given to your previous post. With 
> unconstrained residual variance the model is largely generating 
> nonsense. Use `corg' instead of `us'. Also, depending on what the 
> outcomes are you almost certainly need to have `trait' in the fixed 
> effect specification.
>
> Jarrod
>
>
>
>
>   Quoting Iker Vaquero Alba <karraspito at yahoo.es> on Wed, 23 Sep 2015 
> 13:28:09 +0000 (UTC):
>
>
>
>>
>>    Hello all,
>>    I am implementing a model with a multiple (bivariate) response 
>> variable using MCMCglmm. Both response variables and all the 
>> explanatory variables are categorical variables, with between 2 and 
>> 6 levels. The model is as follows:
>>   
>> 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("threshold","threshold"),data=extphen,nitt=100000,singular.ok=TRUE)
>>    And this is the summary of the model after all the iterations:
>>   
>> summary(testmodel1)  Iterations =3001:99991 Thinninginterval  = 
>> 10 Sample size  = 9700   DIC:   R-structure:  
>> ~us(trait):units                            post.mean l-95% CI u-95% 
>> CI eff.sampnatapshort:natapshort.units    114108   992.1   213000    
>> 7.105nataplong:natapshort.units      33245   310.8    66300    
>> 4.656natapshort:nataplong.units      33245   310.8    66300    
>> 4.656nataplong:nataplong.units       82964   671.5   155869    
>> 2.218  Location effects: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 
>>                        post.mean   l-95% CI  u-95% CI eff.samp   
>> pMCMC   (Intercept)           8.934e+02 -6.995e+02 2.596e+03   
>> 275.73 0.22660   genderM              -8.936e+01 -7.437e+0
>> 2 5.066e+02  4236.67 0.75794   genderO              -1.292e+03 
>> -1.934e+05 1.765e+05  9700.00 0.99052   
>> age                  -3.493e+02 -1.170e+03 3.615e+02   505.92 
>> 0.31918   religionY            -7.361e+02 -1.598e+03 1.481e+01    
>> 33.71 0.03402 * sexorHOM             -1.235e+03 
>> -1.808e+05 1.679e+05  9700.00 0.99814   sexorOT              
>>   2.193e+03 -1.589e+05  1.687e+05 10583.09 0.97814   
>> selfattr             -2.367e+02 -7.424e+02 1.706e+02   314.82 
>> 0.24948   partnerattr           1.391e+02 -2.667e+02 6.147e+02   
>> 966.40 0.49546   genderM:age           2.696e+01 -1.313e+02  
>> 1.748e+02  3786.10 0.69670   genderO:age          -1.055e+04 
>> -8.325e+04 7.163e+04  1722.63 0.78474   
>> genderM:religionY    -1.295e+02 -3.725e+02 8.194e+01   200.08 
>> 0.20495   genderO:religionY    -1.016e+04 -1.589e+05 1.505e+05  
>> 8731.86 0.89052   genderM:sexorHOM     -2.245e+02 
>> -5.713e+02 4.443e+01   105.67 0.10495   
>> genderO:sexorHOM     -8.104e+03 -1.620e+05 1.385e+05  5318.22 
>> 0.90474   genderM:sexorOT      -1.52
>> 0e+02 -5.124e+02 1.856e+02   423.76 0.33402   
>> genderO:sexorOT       2.628e+03 -1.654e+05 1.658e+05  9700.00 
>> 0.97670   genderM:selfattr      9.029e+01 -3.152e+01 2.334e+02   
>> 119.78 0.12907   genderO:selfattr      6.281e+03 
>> -6.511e+04 8.524e+04  3504.67 0.88412   
>> genderM:partnerattr  -7.284e+01 -2.160e+02 6.729e+01   263.29 
>> 0.25052   genderO:partnerattr   2.536e+02 -5.113e+02 1.121e+03  
>> 1291.76 0.49113   age:religionY         8.732e+01 
>> -1.283e+02 3.457e+02   727.46 0.42289   
>> age:sexorHOM          2.809e+02 -8.592e+04 8.847e+04  9700.00 
>> 0.99711   age:sexorOT          -1.246e+03 -8.447e+04 7.941e+04  
>> 9370.57 0.97526   age:selfattr          1.195e+02 
>> -6.636e+01 3.452e+02   212.35 0.19567   
>> age:partnerattr      -8.598e+00 -1.963e+02 1.714e+02  9700.00 
>> 0.92227   religionY:sexorHOM    8.506e+01 -2.392e+02 4.612e+02  
>> 2059.92 0.59959   religionY:sexorOT     1.420e+01 
>> -5.170e+02 5.464e+02  9700.00 0.96268   
>> religionY:selfattr    2.782e+01 -1.198e+02 1.833e+02  3520.80 
>> 0.68701   religionY:part
>> nerattr 1.407e+02  1.423e+01  2.886e+02   22.99 0.00928 
>> **sexorHOM:selfattr     1.160e+02 -1.141e+02 3.707e+02   394.74 
>> 0.28495   sexorOT:selfattr      1.006e+02 -8.528e+01 3.050e+02   
>> 305.50 0.24577   sexorHOM:partnerattr  1.231e+02 
>> -1.246e+02 3.990e+02   415.43 0.31072   
>> sexorOT:partnerattr  -1.401e+00 -2.007e+02  1.956e+02 9700.00 
>> 0.99237   selfattr:partnerattr  5.483e+00 -6.017e+01 7.207e+01  
>> 3007.45 0.85464   ---Signif. codes:  0?***? 0.001 ?**? 0.01 ?*? 0.05 
>> ?.? 0.1 ? ? 1  Cutpoints:                           post.mean l-95% 
>> CI u-95% CI eff.sampcutpoint.traitnatapshort.1     235.2   62.34    
>> 376.2    8.822cutpoint.traitnatapshort.2     633.4  202.16    
>> 944.0    3.578cutpoint.traitnatapshort.3   1139.5   364.35   
>> 1683.7   5.832cutpoint.traitnataplong.1      293.4   54.85    
>> 433.7    5.203cutpoint.traitnataplong.2      651.8  223.70    
>> 961.6    2.604cutpoint.traitnataplong.3    1023.1   344.82   
>> 1483.0   2.353
>>   
>> So, my question is: in that summary, where are the effect sizes, are 
>> they the "post. mean" column? And have they been transformed in some 
>> way? Because obviously, for response variables that can only take 
>> values 1,2,3,4 or 5, I would expect to see those as the effect size.
>> Also, is there any way of knowing to what extent are those results 
>> due to each specific response variable, and the degree of covariance 
>> between both? Is it possible to get all that information from that 
>> summary output I have copied above?
>>    Thank you very much.   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
>>
>>
>>     [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
>
>
> --
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
>
>
>
>
>



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.



More information about the R-sig-mixed-models mailing list