Skip to content
Prev 5258 / 20628 Next

MCMCglmm: ixn between binary family-level factor and us() random effect

Hi Jarrod,

Thanks again, for the syntax and for the modelling advice.

I think (hope) I can answer some of your concerns about my model, particularly about the residual and family-level variances being confounded. In the most general model I'm specifying 8 separate variances, for mothers, fathers, daughters and sons at the within-family (R) and between-family (G) levels:

    random=~us(fampos.f):famid,
    rcov=~idh(fampos.f):units,

There's only one mother per family so individual level = family level and two variances can't be estimated. The same applies to fathers. For both mothers and fathers the variance at one level has to be fixed at (close to) zero:

    R=list(nu=nu,V=diag(c(V,V,0.0001,0.0001)),fix=3)
# the variances on the diagonal are for daughters, sons, mothers, fathers

I didn't think it mattered which of the R or G variances was fixed (unless the R and G priors differ?). I chose to fix the R variance because I don't know how to fix the variances in the G structure without also fixing the covariances that fall within the partition (in MLWIN it makes no difference which is fixed). I calculate the overall 4x4 covariance matrix by adding the daughter and son residual variances to their family variances and converting to a correlation matrix (for each row of the MCMC output).

I'm not totally confident about this method but I've used it for a previous analysis in MLWIN - this gave sensible estimates for family correlations which I replicated by regressing offspring z-scores on parent z-scores. I haven't yet had time yet to check the MCMCglmm results from the current analysis in a similar way - you might save me the trouble by telling me it's invalid.

Cheers,
Paul

PS Here's the call to MCMCglmm:

  MCMCglmm(
    fixed=fixed.form,
    rcov=~idh(fampos.f):units,
    random=~us(fampos.f):famid,
    data=mfs.fam,nitt=nitt,thin=thin,burnin=burnin,
    prior=
      list(
        R=list(nu=nu,V=diag(c(V,V,0.0001,0.0001)),fix=3),
        G=list(G1=list(nu=nu,V=diag(V,4)))))
[1] "Daughter" "Son"      "Mother"   "Father"