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

[jonnations]

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)

-- 
Jonathan A. Nations
PhD Candidate
Esselstyn Lab <http://www.museum.lsu.edu/esselstyn>
Museum of Natural Sciences <http://sites01.lsu.edu/wp/mns>
Louisiana State University

	[[alternative HTML version deleted]]