[R-sig-ME] Question about the predict function in MCMCglmm

[Jarrod Hadfield]

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

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20161103/037234d2/attachment.pl>