[R-sig-ME] Parameter estimates and covariance structure in MCMCglmm

[Adam Lendvai]

Dear List members,

I'm doing analyses in MCMCglmm to investigate individual differences
in reaction to an experiment (i.e. differences in reaction norms). We
measured multiple variables at the same time, so we're also interested
in knowing whether there is a covariance between the reaction norms. I
have three questions regarding parameter estimates and manipulating
the variance-covariance matrix that can be found below within a
minimal reproducible example.

Your feedback would be greatly appreciated.

Best,
Adam Lendvai

# set up test data: We have 20 individuals (identified by 'band'), two
treatments, and each individual had twice the control and
# twice the other treatment during the study
test  <-  expand.grid(band = factor(1:20),  replicate = 1:2, treatment
= factor(1:2, labels = c("ctrl", "treat")))
# response variable:
set.seed(22)
test$y  <-  c(rnorm(40, 5, 5), rnorm(40, 15, 7))
test$z  <-  c(rnorm(40, 30, 5), rnorm(40, 20, 5))
# introduce more variance by defining negative slopes too
test$y[test$band == 8]  <-  c(10,13, 0, 5)
test$y[test$band == 17]  <-  c(8, 5, 1,4)
test$y[test$band == 6]  <-  c(11, 7, 10,6)


require(lattice); require(MCMCglmm)
# Step1: I'd like to estimate the individual differences in the
reaction to the treatment:
Prior1= list(R = list(V = diag(1), nu = 0.002), G = list(G1 = list(V =
diag(2), nu = 1.002)))
mcmc.test1 <- MCMCglmm(y~treatment, prior = Prior1, random =
~us(treatment):band, family = "gaussian", data = test, saveX = T,
saveZ = T, verbose = T, pr = T)
# plot the results
W1 <- cBind(mcmc.test1$X, mcmc.test1$Z)
pred1 <- W1%*%posterior.mode(mcmc.test1$Sol)
xyplot(y+pred1 at x~treatment|band, data=test, type = c("p", "a"), pch =
19, auto.key = list(text  = c("Data", "Prediction"), lines = 1, points
= F))

# the model fit is OKish, although in some cases there is a mismatch
between the data and the prediction
# in the slope parameter (e.g. individual 6,8,17), in the intercept
(e.g. 2), or both (e.g. 18)

# Question 1: is it possible to improve model fit?

# I'm interested in the between-individual variation in the slopes:
summary(mcmc.test1)$Gcovariances[4,1:3]

# Now test it against a model that does not include random slopes,
just random intercept:
Prior2  = list(R=list(V=1, nu = 0.002), G=list(G1 = list(V = 1, nu = 0.002)))
mcmc.test2 <- MCMCglmm(y~treatment, prior = Prior2, random = ~band,
family = "gaussian", data = test, saveX = T, saveZ = T, verbose = T,
pr = T)
summary(mcmc.test2)

summary(mcmc.test2)$DIC-summary(mcmc.test1)$DIC

# A substantial difference in DIC supports the random slope model:
i.e. there are individual differences
# in the reaction to the treatment

# I would be interested in quantifying the within-individual variance
for each slope, i.e. how 'consistent'
# are the individuals in their response? (e.g. individual 7 is
consistent vs. individual 5, which is less consistent).
# Howver, if I get this right, the residual error in these models is
calculated at the population level.

# Question 2: how can I estimate within-individual variance
(consistency) for each individual?


# Step2: I'd like to see if there is covariance in the slopes between
response 'y' and 'z'

# let's standardize reponse values
z = function(x) (x - mean(x, na.rm = T)) / sd(x, na.rm = T)

Prior3 = list(R = list(V = diag(2), nu = 1.002), G = list(G1 = list(V
=diag(4), nu = 3.002)))

mcmc.test3 = MCMCglmm(cbind(z(y), z(z))~(trait-1) * treatment, prior =
Prior3, random = ~us(trait:treatment):band, rcov = ~us(trait):units,
family = c("gaussian", "gaussian"), data = test)
summary(mcmc.test3)
# variance -covariance matrix:
(matrix(round(summary(mcmc.test3)$Gcov[,1], 2),4,4) -> mtr)
mtrnames = c("int.y", "int.z", "sl.y", "sl.z"); dimnames(mtr) =
list(mtrnames, mtrnames)

# I'd like to see the conditional between-individual variance in 'y'
response which is independent from variation in 'z':
(cbivy = mtr["sl.y", "sl.y"] - ((mtr["sl.z", "sl.y"])^2)/mtr["sl.z", "sl.z"])

# Question 3: How can I test this? I can use the ~idh function to
constrain all covariances to 0, but it would not
# be appropriate, because it would set all other covariances to zero.
Is it possible to set to zero just selected parts
# of the variance-covariance matrix (here: covariance of slopes)?


sessionInfo()
# R version 3.0.1 (2013-05-16)
# Platform: x86_64-apple-darwin10.8.0 (64-bit)
#
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base
#
# other attached packages:
#   [1] MuMIn_1.9.5     MCMCglmm_2.17   corpcor_1.6.6   ape_3.0-9
 coda_0.16-1     Matrix_1.0-12   tensorA_0.36
# [8] lattice_0.20-15
#
# loaded via a namespace (and not attached):
#   [1] grid_3.0.1   nlme_3.1-110 tools_3.0.1