[R-sig-ME] Plotting MCMCglmm regression with credible intervals

[Jarrod Hadfield]

Hi Jon,

You can use predict on newdata:

newdata<-df
newdata$x<-seq(min(df$x), max(df$x), length=100)

p1<-predict(model1, newdata=newdata, interval="confidence")

plot(df$x,df$y)
lines(p1[,1]~newdata$x,col="red")
lines(p1[,2]~newdata$x,col="red",lty=2)
lines(p1[,3]~newdata$x,col="red", lty=2)

The predict function with default arguments is quite slow with 
non-Gaussian data because it uses numerical integration to integrate 
over the random effects (if they are to be marginalised). However, the 
integrals can be approximated using Taylor expansion (approx="taylor") 
or if the response involves a logit link (as here) it can be approximate 
using approx="diggle" or approx="mculloch". These are much faster, and 
the latter two especially can be very accurate. However, if you try to 
use these approximations on your model MCMCglmm will tell you they are 
not implemented for this distribution.  This isn't true, and I will put 
a new version on CRAN ASAP that allows you to use these approximations. 
I can also send you a version in the mean time if you let me know your OS.

Cheers,

Jarrod



On 09/02/2018 20:26, jonnations wrote:
> Hi all,
>
> I am working on a logistic regression model using MCMCglmm. I would like to
> plot the mean with 95% credible intervals on my regression plot, however
> this is not so straightforward as the MCMCglmm output is in a very
> different format than most software. I think I can manually bang it out
> with a page of code, but I was curious if anyone knows of a simple, elegant
> way to do this.
>
> Below is a sample script. Imagine this plot with 2 additional lines (95%
> CI's) on either side of the mean, shaded in between.
>
> If r-sig-mixed-models is not the correct venue for this kind of question, I
> apologize.
>
> Best,
> Jon
>
> library(MCMCglmm)
> #generate binary data
> intercept1 = -6.0
> b      = 2.75
> x     = rnorm(100, 1, 3)
> linpred   = intercept1 + x * b
> prob1      = exp(linpred) / (1 + exp(linpred))
> run1     = runif(100, 0, 1)
> y     = ifelse(run1 < prob1, 1, 0)
> df <- data.frame(y = y, x = x)
>
> model1<-MCMCglmm(y ~ x, data = df, family = "categorical",
>                   verbose = F, nitt  =2000, burnin = 500, thin = 1)
>
> x <- df$x
> y <- df$y
> mod.x <- mean(model1$Sol[,1])
> mod.y <- mean(model1$Sol[,2])
> plot(x,y)
> curve(plogis(mod.x+mod.y*x),col="red",add=TRUE)
>


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.