[R-sig-ME] Multivariate Response Model in MCMCglmm - Estimating Posterior Correlation Between Traits
Jarrod Hadfield
j.hadfield at ed.ac.uk
Mon Jun 3 18:05:37 CEST 2013
Hi,
1) If you only have one observation per individual you should drop the
random effects (random=~us(trait):Bird) and have
rcov=~us(trait):units. The posterior correlation can be found as
posterior.cor(mc2$VCV). You could estimate both if the Bird effects
are structured by a pedigree, but I think this is not the case.
2) The effective sample sizes are good enough that the MCMC
approximation to the posterior should be fairly accurate for many
summary statistics (e.g. posterior mean). Following the suggestion in
1) you should get (many) more effective samples per iteration.
3) You need to add one to each of the variances prior to calculating
the correlation in the model as you have specified it. Better to
follow suggestion 1)
4) This is probably because you are fitting a non-identifiable set of
parameters - try running with suggestion 1).
Sorry I took so long to answer this - I am in the field.
Cheers,
Jarrod
Quoting Xavier Harrison <xav.harrison at gmail.com> on Fri, 17 May 2013
17:04:40 +0100:
> 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]]
>
> _______________________________________________
> 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.
More information about the R-sig-mixed-models
mailing list