[R-sig-ME] Question about the predict function in MCMCglmm
Joelle Mbatchou
jmbatchou at uchicago.edu
Thu Nov 3 00:19:48 CET 2016
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
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list