Plotting MCMCglmm regression with credible intervals
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.