[R-sig-ME] Plotting MCMCglmm regression with credible intervals
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sat Feb 10 10:14:37 CET 2018
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.
More information about the R-sig-mixed-models
mailing list