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
------------------------------------------------------------------------
*From:* Jarrod Hadfield [j.hadfield at ed.ac.uk]
*Sent:* Wednesday, November 02, 2016 11:27 AM
*To:* Joelle Mbatchou; r-sig-mixed-models at r-project.org
*Subject:* Re: [R-sig-ME] Question about the predict function in MCMCglmm

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