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]]