##read  bugs data set into R (using stata version of data)
library(foreign)
pp.dat<-read.dta(file.choose())


##copy and paste 'stats' for p from bugs into R
pp.out<-read.table('clipboard')
pp.dat$p<-pp.out[,2]

###in sample predicted prob
plot(pp.dat$p~pp.dat$age,xlim=c(min(pp.dat$age), max(pp.dat$age)))


######### simulated pred probs for given values of covariates
pp.chain<-boa.importBUGS("pp.turn", "C:/bugs.out")
new.age<-seq(18,82)
new.education<-rep(median(pp.dat$educ), 65)
new.closing<-rep(median(pp.dat$closing), 65)
new.govelec<-rep(median(pp.dat$govelec), 65)
new.south<-rep(median(pp.dat$south), 65)
constant<-rep(1,65)


pp.df<-cbind(constant,new.education,new.age, new.south, new.govelec, new.closing)
pp.chain<-as.matrix(pp.chain)

Xb<-t(pp.df%*% t(pp.chain))


mean.Xb<-apply(Xb, 2, mean)




pp.age<-exp(mean.Xb)/(1+exp(mean.Xb))

plot(new.age, pp.age, pch=19)

sd.Xb<-apply(Xb,2, sd)
lb<-mean.Xb-(1.96*sd.Xb)
ub<-mean.Xb+(1.96*sd.Xb)
pp.lb<-exp(lb)/(1+exp(lb))
pp.ub<-exp(ub)/(1+exp(ub))

segments(new.age,pp.lb, new.age,pp.ub, lty=2)


########specific values of South dummy var--different groups
pp.df.south<-cbind(constant,new.education,new.age,rep(1,65), new.govelec, new.closing)
pp.df.north<-cbind(constant,new.education,new.age, rep(0,65), new.govelec, new.closing)

Xb.s<-t(pp.df.south%*% t(pp.chain))
Xb.n<-t(pp.df.north%*% t(pp.chain))


mean.Xb.s<-apply(Xb.s, 2, mean)
mean.Xb.n<-apply(Xb.n, 2, mean)



pp.age.s<-exp(mean.Xb.s)/(1+exp(mean.Xb.s))
pp.age.n<-exp(mean.Xb.n)/(1+exp(mean.Xb.n))

plot(new.age, pp.age.s, pch=19)
points(new.age, pp.age.n, pch=19, col='red')