In n8thangreen/STIecoPredict: Multilevel Regression With Post-Stratification Using Natsal Data

Following explicit-mrp-Natsal-MAIN-'...' scripts ouput.

Some validation plots

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")
}

n8thangreen/STIecoPredict documentation built on June 7, 2020, 12:50 p.m.