Following explicit-mrp-Natsal-MAIN-'...' scripts ouput.
This is taken from Gelman (xxxx). We predict using the logistic fit without one of the covariates and then plot a separate line for each.
gor.names <- c("North East", #1 "North West", #2 "notused", #3 "Yorkshire and The Humber", #4 "East Midlands", #5 "West Midlands", #6 "South West", #7 "East", #8 "London", #9 "South East", #10 "Wales", #11 "Scotland") #12 pred.logit <- fixef(fit)["sexWomen"]*(as.character(sim_prop_la.fit$sex)=="Women") + fixef(fit)["studentTRUE"]*(sim_prop_la.fit$student) + ranef(fit)$ethnic2[as.character(sim_prop_la.fit$ethnic2),1] + ranef(fit)$age[as.character(sim_prop_la.fit$age),1] + fixef(fit)["I(`Average Score` > 28)TRUE"]*(sim_prop_la.fit$`Average Score`>28) pred.logit.Natsal <- fixef(fit)["sexWomen"]*(as.character(Natsal2$sex)=="Women") + fixef(fit)["studentTRUE"]*(Natsal2$student) + ranef(fit)$ethnic2[as.character(Natsal2$ethnic2),1] + ranef(fit)$age[as.character(Natsal2$age),1] + fixef(fit)["I(`Average Score` > 28)TRUE"]*(Natsal2$`Average Score`>28) pred.logit.upperCI <- (fixef(fit)["sexWomen"]+2*se.fixef(fit)["sexWomen"])*(as.character(sim_prop_la.fit$sex)=="Women") + (fixef(fit)["studentTRUE"]+2*se.fixef(fit)["studentTRUE"])*(sim_prop_la.fit$student) + (fixef(fit)["I(`Average Score` > 28)TRUE"]+2*se.fixef(fit)["I(`Average Score` > 28)TRUE"])*(sim_prop_la.fit$`Average Score`>28) + ranef(fit)$ethnic2[as.character(sim_prop_la.fit$ethnic2),1]+2*se.ranef(fit)$ethnic2[as.character(sim_prop_la.fit$ethnic2),1] + ranef(fit)$age[as.character(sim_prop_la.fit$age),1]+2*se.ranef(fit)$age[as.character(sim_prop_la.fit$age),1] pred.logit.lowerCI <- (fixef(fit)["sexWomen"]-2*se.fixef(fit)["sexWomen"])*(as.character(sim_prop_la.fit$sex)=="Women") + (fixef(fit)["studentTRUE"]-2*se.fixef(fit)["studentTRUE"])*(sim_prop_la.fit$student) + (fixef(fit)["I(`Average Score` > 28)TRUE"]-2*se.fixef(fit)["I(`Average Score` > 28)TRUE"])*(sim_prop_la.fit$`Average Score`>28) + ranef(fit)$ethnic2[as.character(sim_prop_la.fit$ethnic2),1]-2*se.ranef(fit)$ethnic2[as.character(sim_prop_la.fit$ethnic2),1] + ranef(fit)$age[as.character(sim_prop_la.fit$age),1]-2*se.ranef(fit)$age[as.character(sim_prop_la.fit$age),1] ### old version # pred.logit <- (fixef(fit)["smokenowTRUE"]*sim_prop_la$smokenow) + # fixef(fit)["studentTRUE"]*sim_prop_la$student + # ranef(fit)$ethnic2[as.character(sim_prop_la$ethnic2),1] + # ranef(fit)$'sex:age'[paste(sim_prop_la$sex, sim_prop_la$age,sep=":"), 1] # # pred.logit.Natsal <- (fixef(fit)["smokenowTRUE"]*Natsal$smokenow) + # fixef(fit)["studentTRUE"]*Natsal$student + # ranef(fit)$ethnic2[as.character(Natsal$ethnic2),1] + # ranef(fit)$'sex:age'[paste(Natsal$sex, Natsal$age,sep=":"), 1] # # pred.logit.upperCI <- ((fixef(fit)["smokenowTRUE"]+2*se.fixef(fit)["smokenowTRUE"])*sim_prop_la$smokenow) + # (fixef(fit)["studentTRUE"]+2*se.fixef(fit)["studentTRUE"])*sim_prop_la$student + # (ranef(fit)$ethnic2[as.character(sim_prop_la$ethnic2),1]+2*se.ranef(fit)$ethnic2[as.character(sim_prop_la$ethnic2),1]) + # (ranef(fit)$'sex:age'[paste(sim_prop_la$sex, sim_prop_la$age,sep=":"), 1]+2*se.ranef(fit)$'sex:age'[paste(sim_prop_la$sex, sim_prop_la$age,sep=":"), 1]) # # pred.logit.lowerCI <- ((fixef(fit)["smokenowTRUE"]-2*se.fixef(fit)["smokenowTRUE"])*sim_prop_la$smokenow) + # (fixef(fit)["studentTRUE"]-2*se.fixef(fit)["studentTRUE"])*sim_prop_la$student + # (ranef(fit)$ethnic2[as.character(sim_prop_la$ethnic2),1]-2*se.ranef(fit)$ethnic2[as.character(sim_prop_la$ethnic2),1]) + # (ranef(fit)$'sex:age'[paste(sim_prop_la$sex, sim_prop_la$age,sep=":"), 1]-2*se.ranef(fit)$'sex:age'[paste(sim_prop_la$sex, sim_prop_la$age,sep=":"), 1]) # par(mfrow=c(2,1)) # par(mfrow=c(4,5)) for (i in c(1,2,4,5,6,7,8,9,10)){ ## logistic curve plot(NA, ylab="", xlab="", ylim=c(0,0.7), main=gor.names[i]) lines(na.omit(pred.logit[sim_prop_la.fit$gor==i][order(pred[sim_prop_la.fit$gor==i])]), sort(pred[sim_prop_la.fit$gor==i])) lines(na.omit(pred.logit.lowerCI[sim_prop_la.fit$gor==i][order(pred[sim_prop_la.fit$gor==i])]), sort(pred[sim_prop_la.fit$gor==i]), col="darkgrey") lines(na.omit(pred.logit.upperCI[sim_prop_la.fit$gor==i][order(pred[sim_prop_la.fit$gor==i])]), sort(pred[sim_prop_la.fit$gor==i]), col="darkgrey") ymin <- par("usr")[3]+0.05; ymax <- par("usr")[4]-0.05 points(pred.logit.Natsal[Natsal$gor==i], (Natsal$cttestly[Natsal$gor==i] * (ymax-ymin)+runif(sum(Natsal$gor==i))/100) + ymin, col=rgb(0.1,0.1,0.1,0.5)) # hist(pred.logit.Natsal[Natsal$gor==i & Natsal$cttestly==0], breaks=20, col=rgb(0.1,0.1,0.1,0.5), main="", prob=T, xlim=c(-1,3),ylim=c(0,0.5), xlab="") # hist(pred.logit.Natsal[Natsal$gor==i & Natsal$cttestly==1], breaks=10, col=rgb(0.8,0.8,0.8,0.5), add=T, prob=T) ## densities # plot(density(pred.logit.Natsal[Natsal$gor==i & Natsal$cttestly==0], prob=T, breaks=30), col="black", type="l", main=gor.names[i], xlim=c(-1,3),ylim=c(0,0.7), xlab="") # lines(density(pred.logit.Natsal[Natsal$gor==i & Natsal$cttestly==1], breaks=10), col="darkgrey") }
pred.logit <- (fixef(fit)["smokenowTRUE"]*sim_prop_la$smokenow) + fixef(fit)["studentTRUE"]*sim_prop_la$student + ranef(fit)$gor[as.character(sim_prop_la$gor),1] + ranef(fit)$'sex:age'[paste(sim_prop_la$sex, sim_prop_la$age,sep=":"), 1] pred.logit.Natsal <- (fixef(fit)["smokenowTRUE"]*Natsal$smokenow) + fixef(fit)["studentTRUE"]*Natsal$student + ranef(fit)$gor[as.character(Natsal$gor),1] + ranef(fit)$'sex:age'[paste(Natsal$sex, Natsal$age,sep=":"), 1] # par(mfrow=c(2,1)) # par(mfrow=c(3,4)) for (i in levels(sim_prop_la$ethnic2)[1:6]){ ## logistic curve plot(pred.logit[sim_prop_la$ethnic2==i][order(pred[sim_prop_la$ethnic2==i])], sort(pred[sim_prop_la$ethnic2==i]), type="l", ylab="", xlab="", ylim=c(0,0.65), main=i) points(pred.logit.Natsal[Natsal$ethnic2==i], Natsal$cttestly[Natsal$ethnic2==i]*0.6+runif(sum(Natsal$ethnic2==i))/100, col=rgb(0.1,0.1,0.1,0.5)) # hist(pred.logit.Natsal[Natsal$ethnic2==i & Natsal$cttestly==0], breaks=20, col=rgb(0.1,0.1,0.1,0.5), main="", prob=T, xlim=c(-1,3),ylim=c(0,0.5), xlab="") # hist(pred.logit.Natsal[Natsal$ethnic2==i & Natsal$cttestly==1], breaks=10, col=rgb(0.8,0.8,0.8,0.5), add=T, prob=T) ## densities plot(density(pred.logit.Natsal[Natsal$ethnic2==i & Natsal$cttestly==0]), col="black", type="l", main=i, xlim=c(0,2.5),ylim=c(0,0.7), xlab="") lines(density(pred.logit.Natsal[Natsal$ethnic2==i & Natsal$cttestly==1]), col="darkgrey") }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.