[R-sig-ME] Question about the predict function in MCMCglmm
Joelle Mbatchou
jmbatchou at uchicago.edu
Thu Nov 3 22:09:26 CET 2016
Hi Jarrod,
Thank you for the clarification! I tried both of the approximations specified in the course notes and I do get closer predictions using the approximation from McCulloch and Searle (2001) but only when marginalization is done with respect to the residuals.
When marginalizing over the random effects, I get that the approximation gives quite different results from those of the predict function:
######## Using predict.MCMCglmm() -- I modified function for pred2
pred1 <- as.vector(predict(mMCMC, marginal = ~idU, type = "response"))
pred2 <- as.vector(predict_MCMC(mMCMC, marginal = NULL, posterior = "mean", type = "response"))
pred3 <- as.vector(predict(mMCMC, marginal = NULL, posterior = "all", type = "response"))
sig_sq <- mMCMC$VCV[,"units"] + mMCMC$VCV[,"idU"]
######## Using approximation from Diggle et al. (2004)
pred1_hand_Da <- colMeans(plogis(Xb / sqrt(1 + ck^2 * sig_sq)))
pred2_hand_Da <- as.vector(plogis( (Xb_mean + U_mean)/ sqrt(1 + k^2 * 1) ))
pred3_hand_Da <- colMeans(plogis( (Xb + U)/sqrt(1 + k^2 * 1) ))
######## Using approximation from McCulloch and Searle (2001)
pred1_hand_MS <- colMeans(plogis( Xb - 0.5 * sig_sq * tanh(Xb * (1 + 2 * exp(-0.5 * sig_sq))/6) ))
pred2_hand_MS <- as.vector(plogis( (Xb_mean + U_mean - 0.5 * 1 * tanh( (Xb_mean + U_mean) * (1 + 2 * exp(-0.5 * 1))/6)) ))
pred3_hand_MS <- colMeans(plogis( (Xb + U) - 0.5 * 1 * tanh((Xb + U) * (1 + 2 * exp(-0.5 * 1))/6) ))
#################
###### Comparison
## Marginalizing over u and e
pred1 pred1_hand_Da pred1_hand_MS
1 0.08467802 0.03564130 0.3434171
2 0.23766175 0.15399697 0.4817053
3 0.19338132 0.11215077 0.4522571
4 0.15685811 0.08232281 0.4229498
5 0.14048409 0.07032301 0.4079111
6 0.17927898 0.10013013 0.4415578
## Taking posterior mean of b and u then marginalizing over e
pred2 pred2_hand_Da pred2_hand_MS
[1,] 0.07628308 0.11490073 0.07631337
[2,] 0.20802371 0.24496249 0.20949077
[3,] 0.41934658 0.43098835 0.42044645
[4,] 0.03839855 0.06985046 0.03826145
[5,] 0.13323146 0.17401086 0.13387440
[6,] 0.07569975 0.11425514 0.07572540
## Marginalizing over e then taking the posterior mean
pred3 pred3_hand_Da pred3_hand_MS
1 0.1622893 0.1869202 0.1626918
2 0.2938466 0.3114496 0.2942767
3 0.4278642 0.4364018 0.4282031
4 0.1201754 0.1443062 0.1205276
5 0.2186460 0.2418515 0.2191097
6 0.1730352 0.1961313 0.1734039
As you can see, I can approximate pred2 and pred3 pretty well in 3rd column. I am trying to reason why the pred1 approximation is so off... It seems to me that having a variance different from 1 in the normal pdf used in the integration of logit^-1(xb+u+e) should not make the approximation that much worse. Or is that not the case? I compared the results of 'normal.multilogistic(x,v)' and 'plogis( x - 0.5 * v * tanh(x * (1 + 2 * exp(-0.5 * v))/6) )' for a few (x,v) values from data:
> x <- Xb[,1]; v <- sig_sq[,1] # Looking at first individual across all MCMC iterations
> difference <- sapply(1:length(x), function(i) abs(normal.multilogistic(x[i],v[i])-plogis( x[i] - 0.5 * v[i] * tanh(x[i] * (1 + 2 * exp(-0.5 * v[i]))/6) )))
> summary(difference)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000004 0.0065390 0.0614500 0.2588000 0.4812000 0.9545000
> summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-11.760 -6.042 -4.865 -5.346 -4.084 -2.678
> summary(x[difference>.1]) ## Not much difference
Min. 1st Qu. Median Mean 3rd Qu. Max.
-11.760 -8.156 -6.268 -6.755 -5.422 -3.826
> summary(v)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.146 5.414 9.154 13.670 16.490 70.350
> summary(v[difference>.1])
Min. 1st Qu. Median Mean 3rd Qu. Max.
9.699 13.030 17.630 23.380 33.950 70.350
So the difference observed seems due to very high values of Var(u+e). I guess I answered my own question (still including code above if others have similar Q)... So in that case the other approximation from Diggle et al. (2004) would produce overall more accurate results.
One last Q: Is there an intuitive explanation as to why one might prefer using pred2 vs pred3, that is getting posterior means of b and u then marginalizing logit^-1(Xb+u+e) (i.e. posterior = "mean") versus first marginalizing logit^-1(Xb+u+e) and then taking the mean over all iterations (i.e. posterior = "all")?
Thank you for the very quick replies and going through all my messages!
Cheers,
Joelle
________________________________
From: Jarrod Hadfield [j.hadfield at ed.ac.uk]
Sent: Thursday, November 03, 2016 1:15 AM
To: Joelle Mbatchou; r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Question about the predict function in MCMCglmm
Oh dear - it is a bug with the predict function. The default posterior="all" which gets the posterior predicted mean is fine, but as you say making the prediction on the posterior mean/mode of the parameters is not. Having
object$Sol<-matrix(colMeans(object$Sol), 1, ncol(object$Sol))
below (L66)
object$VCV<-matrix(colMeans(object$VCV), 1, ncol(object$VCV))
and
object$Sol<-matrix(posterior.mode(object$Sol, ...), 1, ncol(object$Sol))
below (L70)
object$VCV<-matrix(posterior.mode(object$VCV, ...), 1, ncol(object$VCV))
fixes it. I will update MCMCglmm today.
Cheers,
Jarrod
On 02/11/2016 23:19, Joelle Mbatchou wrote:
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<mailto: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<mailto: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