[R-sig-ME] how to specify priors in MCMCglmm?

[Jarrod Hadfield]

Hi Darren,

The prior you specified for the first model has the correct form:

prior = list(R = list(V = 1, n = 0, fix = 1), G = list(G1 = list(V =  
1,n = 0)))

However, the prior is improper because n=0. This does not matter for  
the residual variance because it is fixed at 1 (V=1), but for the  
group variance it does matter. The posterior is tending to zero where  
two things may happen. It can get trapped at low values and give odd  
looking traces, or it can reach very low values at which point it will  
be picked up during matrix inversion and you'll get the error message  
" ill-conditioned G/R structure". If you make the prior proper the  
posterior is guaranteed to be proper also. In this case n=1 is enough  
to ensure propriety, and the chain should not get trapped at zero.   
You ca use higher n but this means the prior will be more informative.  
For more general models where V is a matrix (rather than a scalar)   
n>=dim(V) to ensure propriety.

For the second model you need to add an additional prior element in  
the G-structure because of the additional term (+atime). Something like

prior2 = list(R = list(V = 1, n = 0, fix = 1), G = list(G1 = list(V =  
1,n = 1), G2 = list(V = 1,n = 1)))

but a more sensible value for V (or a check on the sensitivity to  
different choices of V) would be recommended.

Hope this helps,

Jarrod


On 7 May 2009, at 01:33, Darren Norris wrote:

> Dear all,
>
> I would like to know how to specify priors for a MCMCglmm model with  
> a binary response. I have read the help (?MCMCglmm) and searched the  
> best I can through nabble R and this list.
> Unfortunately if the answer is there I am too stupid to understand it.
>
>> From the MCMCglmm tutorial I see how to code a univariate version -  
>> but
> then I get lost about how to extend the prior specification to more  
> complex models.
> I also do no understand the prior terminology or how / why to  
> specify what I need -  so if at all possible I need an idiots guide  
> how to specify the priors.
>
> I work for a conservation NGO and unfortunately we're not affiliated  
> with an academic institution so don't have access to books or  
> journal articles or any academic statistical support.
> Any guidance would be much appreciated.
> Many thanks,
> Darren
>
> R version 2.8.1 (2008-12-22) i386-pc-mingw32
> locale:
> LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United  
> Kingdom.1252;LC_MONETARY=English_United Kingdom. 
> 1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
> other attached packages:
> [1] MCMCglmm_1.09      gtools_2.5.0       combinat_0.0-6      
> orthopolynom_1.0-1 polynom_1.3-4      pscl_1.03         [7]  
> mvtnorm_0.9-4      ape_2.2-3          coda_0.13-4         
> Matrix_0.999375-18 lattice_0.17-17    tensorA_0.31      [13]  
> corpcor_1.5.2      MASS_7.2-45
> loaded via a namespace (and not attached):
> [1] gee_4.13-13 grid_2.8.1  nlme_3.1-89
>
>
> #make the data frame- using "sample" so data is not the same but the  
> important difference is between probability of events in "cD" and  
> the other two ("aD" and "bD")
> # we coducted the work in 3 different places, each place has 3  
> habitat types.
> #Each habitat type has 2 groups of samples (3 samples per group)  
> each group in a habitat has a different type of  bait.
> #I am trying to model how the occurence of events was influenced by  
> habitat and bait. We checked samples repeatedly (atime).
> #Specifying atime and agroup as random effects to deal with  
> potential pseudo replication.
>
> aD<-sample(0:1,108,replace="True",prob=c(0.6,0.4));
> bD<-sample(0:1,108,replace="True",prob=c(0.6,0.4));
> cD<-rep(0,108);
> aevent<-c(aD,bD,cD);
> aplace<-rep(c(8,8,8,8,8,8,9,9,9,9,9,9,10,10,10,10,10,10),9);
> ahabitat<-c(rep(1,108),rep(2,108),rep(3,108));
> atime<-rep(c(0,1,2,4,8,24),54);
> abait<-rep(c(1,1,1,1,1,1,2,2,2,2,2,2),27);
> agroup<-paste(aplace,ahabitat,abait);
> MyDf<-data.frame(aplace,agroup,aevent,ahabitat,atime,abait)
>
> #-using prior code examples from the MCMCglmm tutorial (section 1.2:  
> produced April 17 2009), but priors are not right.
> library(MCMCglmm);
> prior = list(R = list(V = 1, n = 0, fix = 1), G = list(G1 = list(V =  
> 1,n = 0)));
> MMCMm1 <- MCMCglmm(aevent ~ as.factor(ahabitat), random = ~agroup,  
> family = "categorical",prior=prior,data = MyDf, verbose = FALSE)
>
>
> # I would like to specify, but don't know how to specify priors:
> MMCMm2 <- MCMCglmm(aevent ~ as.factor(ahabitat)*as.factor(bait),  
> random = ~agroup + atime, family = "categorical",prior=prior,data =  
> MyDf, verbose = FALSE)
>
> _______________________________________________
> 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.