Hi All

I am trying to fit a bivariate response model to estimate the covariance between two traits (Mass and Number of Offspring), after controlling for several variables that I know affect these traits. I know that Mass is affected by body size (skull length) and day of season (cycle day); and suspect that number of offspring (JuvFuture) is affected by NAO (njs). I have strong reason to expect a positive covariance between Mass and Offspring Number, as the animals are capital breeders and fuel reproductive effort from fat stores. The dataset is small - 213 individuals, each measured only once, which I suspect may cause problems with this kind of model. Nevertheless, the modelled covariates come out as significant and in the direction of effect I would expect. My questions largely pertain to the estimation of the posterior correlation between traits once controlling for the covariates. Questions are below, and prior, model and calculation code are below that. 

1 ) As I only have a single observation per individual, am I right in thinking that I cannot estimate both the G and R matrix simultaneously, and therefore must fix the residual matrix to a small positive value like 1? 

2) Should I be worried that the effective sample size for the 'JuvFuture' variance is approx. half of the total expected effective sample size? Is this an issue of my own (potentially poor) model specification?

3) The 95% credible intervals for the mass/offspring covariance cross zero, which I take to mean that one cannot accurately state that there is positive covariance between the two traits. However, when I scale the covariance to a correlation (code below) I get a positive mean value with credible intervals that approach, but do not cross zero. Is this sufficient to conclude significant posterior correlation between traits once controlling for the covariates in the model?

4) Assuming all so far is ok / not heretical, I have been testing out sensitivity of the model to priors. If I try the parameter expanded priors Jarrod Hadfield has previously suggested for such models (list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*a)) for the G structure, I get some wild results, with a mean correlation of 0.99 and then a CI of -0.8 to 0.999, so very asymmetric and suspiciously high. For a previous, unrelated model I ran, these priors worked well. Are the current priors ('prior 3' below), overwhelming the data somehow to 'force' a positive correlation? The symmetry of the CIs around the mean for prior3 are more reassuring than the hugely asymmetric CIs on the parameter expanded priors, but that could just reflect a lack of understanding on my part. 

 Any help greatly appreciated. 

Xav



#Prior
prior3 <- list( R=list(V=diag(2),fix=1), G=list(G1=list(V=diag(2),nu=0.5)) )

#Model
mc2<-MCMCglmm(cbind(Mass,JuvFuture)~trait-1 +at.level(trait,1):poly(cycleday,2)+at.level(trait,1):Skull + at.level(trait,2):njs ,random=~us(trait):Bird  ,rcov=~idh(trait):units,prior=prior3,data=allstage,family=c("gaussian","poisson"),verbose=F,pr=T,nitt=250000,burnin=50000,thin=50)

#Estimate Mode and 95% CI of Posterior Correlation
posterior.mode(mc2$VCV[,2] / (sqrt(mc2$VCV[,1]) * sqrt(mc2$VCV[,4])))
HPDinterval(mc2$VCV[,2] / (sqrt(mc2$VCV[,1]) * sqrt(mc2$VCV[,4])))
	#Gives mean = 0.47; 95%CI = 0.06 - 0.85

#Output
 Iterations = 50001:249951
 Thinning interval  = 50
 Sample size  = 4000 

 DIC: 1334.423 

 G-structure:  ~us(trait):Bird

                         post.mean   l-95% CI  u-95% CI eff.samp
Mass:Mass.Bird           9009.4427 7341.21278 1.080e+04     4000
JuvFuture:Mass.Bird        21.6487   -1.59616 4.626e+01     2988
Mass:JuvFuture.Bird        21.6487   -1.59616 4.626e+01     2988
JuvFuture:JuvFuture.Bird    0.2813    0.05117 5.705e-01     2206

 R-structure:  ~idh(trait):units

                post.mean l-95% CI u-95% CI eff.samp
Mass.units              1        1        1        0
JuvFuture.units         1        1        1        0

 Location effects: cbind(Mass, JuvFuture) ~ trait - 1 + at.level(trait, 1):poly(cycleday, 2) + at.level(trait, 1):Skull + at.level(trait, 2):njs 

                                      post.mean  l-95% CI  u-95% CI eff.samp  pMCMC    
traitMass                             -134.0018 -532.3044  311.4352     4000 0.5155    
traitJuvFuture                          -0.7929   -1.0669   -0.5184     3119 <3e-04 ***
at.level(trait, 1):poly(cycleday, 2)1 2515.7841 2256.6187 2781.0072     4231 <3e-04 ***
at.level(trait, 1):poly(cycleday, 2)2  407.4622  136.2618  660.5247     4000 0.0015 ** 
at.level(trait, 1):Skull                21.9526   16.9723   26.6290     4000 <3e-04 ***
at.level(trait, 2):njs                  -0.9975   -1.3345   -0.7135     3841 <3e-04 ***
---






	[[alternative HTML version deleted]]

