# [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

________________________________
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
________________________________
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]]

```