[R-sig-ME] set up of priors in MCMCglmm for different error families

[albrechj at staff.uni-marburg.de]

Dear list,

I want to us MCMCglmm to fit three phylogenetic mixed models. One for  
each of three response variables. Basically, the model formula is as  
follows:

response ~ trait1 * trait2 + (1|plant_species) + (1|site) + (1|year);

The data is based on observations from 15 plant species in 13 study  
sites across two years, but the dataset is not completely full  
factorial with respect to the random grouping factors. The dataset is  
relatively small (n = 121).

The error families of the three response variables are gaussian,  
multinomial2 (i.e., binomial), and ztpoisson (zero-truncated poisson).

My question is now which priors should I use for these three models. I  
have searched the course notes and on the list for suggestions, and  
ended up with flat *improper* priors for the residual variance and  
parameter expanded priors for the (co)variance in the random terms.

However, at this point I'am stuck and I don't know whether the priors  
are specified correctly.

Could anyone give advice on the correct specification of the priors  
for these models? I have attached the r script below (I hope it is not  
too confusing, because I have written this script to fit the models on  
multiple cores parallel).


Thank you.

Jörg


The code for model fitting:
################################################################################
### modelfitting
################################################################################
if(! file.exists(paste(compiled_data_victoria_dir,
    "phylogenetic_mcmcglmm_results.rda",sep="/"))) {
### set up the phylogenetic covariance matrix
# based on a tree including 15 plant species
Ainv<-inverseA(phylo1, nodes="TIPS")$Ainv

### set up the model formulas
responses <- c("plant_d", "cbind(removal, total_rem)", "richness")
### paste response and predictors into one character string
fixed1 <- paste(responses, "s.fruit_density_con * s.phen_d", sep=" ~ ")
### formula for random effects
random1 <- "~ plant_species + site + year"
### set up error families
families1 <- c("gaussian", "multinomial2", "ztpoisson")

### set up priors
priors1 <- list(list(R=list(V=1, nu=0),
                      G=list(G1=list(V=1, nu=0, alpha.mu=0, alpha.V=1e3),
                             G2=list(V=1, nu=0, alpha.mu=0, alpha.V=1e3),
                             G3=list(V=1, nu=0, alpha.mu=0, alpha.V=1e3))),
                 list(R=list(V=1, nu=0),
                      G=list(G1=list(V=1, nu=1e3, alpha.mu=0, alpha.V=1),
                             G2=list(V=1, nu=1e3, alpha.mu=0, alpha.V=1),
                             G3=list(V=1, nu=1e3, alpha.mu=0, alpha.V=1))),
                 list(R=list(V=1, nu=0),
                      G=list(G1=list(V=1, nu=0, alpha.mu=0, alpha.V=1e3),
                             G2=list(V=1, nu=0, alpha.mu=0, alpha.V=1e3),
                             G3=list(V=1, nu=0, alpha.mu=0, alpha.V=1e3))))

### make a list containing the formulas for fixed and random terms,
# error family, and prior information for for each modelfit
mod_setup <- list()
for (i in 1:length(mod_formulas1)) {
mod_setup[[i]] <- list(fixed = fixed1[i], random = random1,
                        family = families1[i], prior = priors1[[i]])
}

### function to fit the models
run1 <- function(mod_setup)  {
   MCMCglmm(fixed = as.formula(mod_setup$fixed),
            random = as.formula(mod_setup$random),
            ginverse = list(plant_species = Ainv), prior = mod_setup$prior,
            data = data_complete, nodes = "TIPS", family = mod_setup$family,
            nitt = 1.5e6, thin = 1e3, burnin = 5e5)
}

### parallel fitting of models on multiple cores
system.time(result1 <- mclapply(mod_setup, run1, mc.cores=length(mod_setup)))
### save the result
save(list=c("result1"), file=paste(compiled_data_victoria_dir,
      "phylogenetic_mcmcglmm_results.rda",sep="/"))
} else {
### load the result if allready exists
   base::load(file=paste(compiled_data_victoria_dir,
              "phylogenetic_mcmcglmm_results.rda",sep="/"))
}
################################################################################
### End of script
################################################################################