With repeated measures on individuals it is often of interest to see how repeatable a trait is. We can estimate the repeatability of a trait as the proportion of phenotypic variance \(V_P\) explained by individual variance \(V_{ind}\); \(R = V_{ind}/V_P = V_{ind}/(V_{ind}+V_R)\). As you already know, bayesian modelisation requires prior. Here, we create a unformative prior with one estimate for the G matrix and one estimate for the Residual matrix, in addition
# p.var <- var(gryphonRM$laydate, na.rm = TRUE)
prior3.1 <- list(G = list(G1 = list(V = 1, nu = 0.002)), R = list(
V = 1,
nu = 0.002
))
model3.1 <- MCMCglmm(laydate ~ 1,
random = ~animal, data = gryphonRM,
prior = prior3.1, verbose = FALSE
)
posterior.mode(model3.1$VCV) animal units
11.27557 21.16771 Note the use of the term animal as random allowed to partition the phenotypic variance \(V_P\) into among individual variance \(V_{ind}\) associated with animal and residual variance \(V_R\) associated with units.
Here then the repeatability of the laydate can be determined as: 22.17 (i.e., as 11.276/(11.276 + 21.168)). Just a friendly remember, we work with Monte Carlo chain with model iteration, so the point estimate can be different (but very similar) each time you run the model.
Mean lay date might change with age, so we could ask what the repeatability of lay date is after conditioning on age. This would be done by adding age into the model as a fixed effect.
model3.2 <- MCMCglmm(laydate ~ age,
random = ~animal, data = gryphonRM,
prior = prior3.1, verbose = FALSE
)
par(mar = c(1, 1, 1, 1))
plot(model3.2$Sol)
plot(model3.2$VCV)
posterior.mode(model3.2$VCV) animal units
12.06557 16.47879 The model assumption seems correct, so we can look at the different estimates. Note that the random effect structure has remained unchanged because we did not modified the prior prior3.1. The repeatability of laydate, after accounting for age effects, is now estimated as 22.17 (i.e., as 11.276/(11.276 + 21.168)). Just as we saw when estimating \(h_2\) in tutorial 1, the inclusion of fixed effects will alter the estimated effect size if we determine total phenotypic variance as the sum of the variance components. Thus, proper interpretation is vital.
posterior.mode(model3.2$Sol)(Intercept) age3 age4 age5 age6
20.228762 2.582580 4.228445 6.107462 3.037460 HPDinterval(model3.2$Sol, 0.95) lower upper
(Intercept) 19.756458 20.833641
age3 1.896463 3.228363
age4 3.573414 4.852338
age5 5.443962 6.736376
age6 2.545933 3.793747
attr(,"Probability")
[1] 0.95Here age is modeled as a 5-level factor (specified using the function as.factor() at the beginning of the analysis). We could equally have fitted it as a continuous variable, in which case, given potential for a late life decline, we would probably also include a quadratic term. In addition, using age as continuous variable can help in saving some degree of freedom in the analysis.
gryphonRM$age_c <- as.numeric(gryphonRM$age)
model3.2_2 <- MCMCglmm(laydate ~ age_c + I(age_c^2),
random = ~animal, data = gryphonRM,
prior = prior3.1, verbose = FALSE
)
par(mar = c(1, 1, 1, 1))
plot(model3.2_2$Sol)
plot(model3.2_2$VCV)
posterior.mode(model3.2_2$VCV) animal units
11.78516 16.41881 posterior.mode(model3.2_2$Sol)(Intercept) age_c I(age_c^2)
15.0732890 5.5363778 -0.7823499 HPDinterval(model3.2_2$Sol, 0.95) lower upper
(Intercept) 13.9815893 16.0702462
age_c 4.8411740 6.3412406
I(age_c^2) -0.9016996 -0.6572575
attr(,"Probability")
[1] 0.95Generally we expect that the repeatability will set the upper limit for heritability since among individual variation can be decomposed in the additive genetic variation and non additive genetic variation. In other word, the additive genetic variation is a subcomponent of the difference between individuals. Non-additive contributions to fixed among-individual differences are normally referred to as permanent environment effects. If a trait has repeated measures then it is necessary to model permanent environment effects in an animal model to prevent upward bias in \(V_A\).
To illustrate it, we first fit the animal model:
Ainv <- inverseA(gryphonped)$Ainv
model3.3 <- MCMCglmm(laydate ~ 1 + age,
random = ~animal, ginv = list(animal = Ainv),
data = gryphonRM, prior = prior3.1, verbose = FALSE
)Variance components are almost unchanged if we compare the previous model:
posterior.mode(model3.3$VCV) animal units
13.88073 16.61835 posterior.mode(model3.2$VCV) animal units
12.06557 16.47879 This suggests that most of the among-individual variance is – rightly or wrongly – being partitioned as \(V_A\) here. In fact here the partition is wrong since the simulation included both additive genetic effects and additional fixed heterogeneity that was not associated with the pedigree structure (i.e. permanent environment effects). In order to o obtain an unbiased estimate of \(V_A\), we need to fit the individual identity twice in the model: once linked to the pedigree (genetic effect) and once not linked to the pedigree (permanent environment effect).To do so, we need to duplicate the variable containing the individual identity animal and give it a new name. In addition, the prior need to be modified to integrate a seconf random effect.An more appropriate estimate of \(V_A\) is given by the model:
gryphonRM$animal_pe <- gryphonRM$animal
# p.var <- var(gryphonRM$laydate, na.rm = TRUE)
prior3.4 <- list(G = list(G1 = list(V = 1, nu = 0.002), G2 = list(
V = 1,
nu = 0.002
)), R = list(V = 1, nu = 0.002))
model3.4 <- MCMCglmm(laydate ~ 1 + age,
random = ~ animal + animal_pe,
ginv = list(animal = Ainv), data = gryphonRM, prior = prior3.4, verbose = FALSE
)
posterior.mode(model3.4$VCV) animal animal_pe units
3.500474 7.577899 16.539407 The estimate of\(V_A\) is now much lower (reduced from 13.6735 to 5.1238) due to a proper separation in the additive and permanent environment effects. We can estimate \(h^2\) and the repeatability from this model:
model3.4.VP <- model3.4$VCV[, "animal"] + model3.4$VCV[, "animal_pe"] + model3.4$VCV[, "units"]
model3.4.PE_VA <- model3.4$VCV[, "animal"] + model3.4$VCV[, "animal_pe"]
posterior.mode(model3.4.PE_VA / model3.4.VP) var1
0.4250901 posterior.mode(model3.4$VCV[, "animal"] / model3.4.VP) var1
0.1229688 Models of repeated measures can be extended to include other fixed or random effects. For example we can try including year of measurement (year) and birth year (byear) as other random effects.
# p.var <- var(gryphonRM$laydate, na.rm = TRUE)
prior3.5 <- list(G = list(G1 = list(V = 1, nu = 0.002), G2 = list(
V = 1,
nu = 0.002
), G3 = list(V = 1, nu = 0.002), G4 = list(
V = 1,
nu = 0.002
)), R = list(V = 1, nu = 0.002))
model3.5 <- MCMCglmm(laydate ~ 1 + age,
random = ~ animal + animal_pe +
year + byear, ginv = list(animal = Ainv), data = gryphonRM, prior = prior3.5,
verbose = FALSE
)
posterior.mode(model3.5$VCV) animal animal_pe year byear units
4.214827166 9.094013914 6.863855767 0.002860896 7.844281659 HPDinterval(model3.5$VCV, 0.95) lower upper
animal 1.8364038852 7.9560325
animal_pe 5.4411251224 11.3866006
year 4.7436521093 12.8796974
byear 0.0002281004 0.2656322
units 7.2322172604 8.6215165
attr(,"Probability")
[1] 0.95
This model will return additional variance components corresponding to year of measurement effects and birth year of the female effects.
\(V_{byear}\) is very low and its posterior distribution (via the function HPDinterval or plot) is very close to zero indicating its not significance. You have to remember bayesian model never estimate variable to 0 or passing zero, so you will never see a credible interval CI crossing zero for a variance. If you compared the DIC of model3.5 to a reduced model without byear, it should be very similar.
prior3.5_2 <- list(
G = list(G1 = list(V = 1, nu = 0.002), G2 = list(
V = 1,
nu = 0.002
), G3 = list(V = 1, nu = 0.002)),
R = list(V = 1, nu = 0.002)
)
model3.5_2 <- MCMCglmm(laydate ~ 1 + age,
random = ~ animal + animal_pe +
year, ginv = list(animal = Ainv), data = gryphonRM, prior = prior3.5_2,
verbose = FALSE
)
posterior.mode(model3.5_2$VCV) animal animal_pe year units
4.416027 8.590906 7.528937 7.819615 model3.5$DIC[1] 8290.961model3.5_2$DIC[1] 8290.354year effects could alternatively be included as fixed effects (try it!, you should be able to handle the new prior specification at this point). This will reduce \(V_R\) and increase the estimates of heritability and repeatability, which must now be interpreted as proportions of phenotypic variance after conditioning on both age and year of measurement effects.