Question about the predict function in MCMCglmm

Hi Jarrod,

Thanks for the quick reply! The V matrix I provide has 1 on the diagonal; I implemented what you suggested:

pred1 <- predict(mMCMC, marginal = ~idU, type = "response")
pred2 <- predict(mMCMC, marginal = NULL, posterior = "mean", type = "response")
pred3 <- predict(mMCMC, marginal = NULL, posterior = "all", type = "response")

beta_hat <- as.matrix(mMCMC$Sol[,1:2])
beta_hat_mean <- colMeans(beta_hat)  # Posterior mean for beta
U <- as.matrix(mMCMC$Sol[,-(1:2)])
U_mean <- colMeans(U)  # Posterior mean for random effects u
X <- as.matrix(mMCMC$X)
k <- ((16*sqrt(3))/(15*pi))^2
sig_sq <- mMCMC$VCV[,"units"] + tcrossprod(mMCMC$VCV[,"idU"], diag(V))  #Variance of latent variables for each MCMC iteration

pred1_hand <- colMeans(plogis(tcrossprod(beta_hat, X) / sqrt(1 + k^2 * sig_sq)))
pred2_hand <- plogis( (tcrossprod(beta_hat_mean, X) + U_mean)/sqrt(1 + k^2 * mean(mMCMC$VCV[,"units"])) )
pred3_hand <- colMeans(plogis( (tcrossprod(beta_hat, X) + U)/sqrt(1 + k^2 * mMCMC$VCV[,"units"]) ))

head(cbind(pred1, pred1_hand, pred2, pred2_hand, pred3, pred3_hand))

pred1 pred1_hand      pred2 pred2_hand     pred3 pred3_hand
1 0.08467802  0.1772529 0.21552374 0.11490073 0.1622893  0.1869202
2 0.23766175  0.3115590 0.17228113 0.24496249 0.2938466  0.3114496
3 0.19338132  0.2764658 0.07845114 0.43098835 0.4278642  0.4364018
4 0.15685811  0.2458417 0.10811526 0.06985046 0.1201754  0.1443062
5 0.14048409  0.2314295 0.35230597 0.17401086 0.2186460  0.2418515
6 0.17927898  0.2648570 0.09011039 0.11425514 0.1730352  0.1961313

I do get that pred3 and pred3_hand are similar though not equal. I understand it's because in the pred3 the expectation of logit^-1(Xb+u+e) is taken over the distribution of the residual e whereas in pred3_hand, I am getting an approximation of that.

I have looked at the code for predict to get a better understanding and I saw was that when 'posterior=mean' is specified in pred2, the posterior mean for sigma_u^2 is computed (code line 93-94):
object$VCV <- matrix(colMeans(object$VCV), 1, ncol(object$VCV))
it <- 1

but only the first MCMC iteration is kept for b and u (code line 109-110):
object$Sol <- object$Sol[it, , drop = FALSE]

Is there a reason why the posterior mean is not obtained for b and u as well before computing the latent variable for each sample (like I used for pred2_hand)? I understand that the predictions obtained won't depend on sigma_u^2 since we fix u and only consider residual e (that has variance of 1) as random when taking expectation of logit^-1(Xb+u+e).

Cheers,
Joelle

Hi,

Also

tcrossprod(beta_hat, X) + U

should be

(tcrossprod(beta_hat, X) + U)/sqrt(1 + k^2 * rowSums(mMCMC$VCV[,"units"]))

for pred2_hand

Cheers,

Jarrod


On 02/11/2016 16:07, Jarrod Hadfield wrote:
tcrossprod(beta_hat, X) + U