## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(ham)

## ----vlogo, echo=FALSE, out.width="30%", fig.align="center"-------------------
knitr::include_graphics("logo.png") 

## ----Post0--------------------------------------------------------------------
blos1 <- Bayes(x=losmcmc)

## ----plotdx1, fig.dim = c(6, 4.5)---------------------------------------------
plot(blos1, y="dxa", parameter="muOfY")

## ----plotdx2, fig.dim = c(6, 4.5)---------------------------------------------
plot(blos1, y="dxd", parameter="muOfY")

## ----plotdx3, fig.dim = c(6, 4.5)---------------------------------------------
plot(blos1, y="dxg", parameter="muOfY")

## ----plotdx4, fig.dim = c(6, 4.5)---------------------------------------------
plot(blos1, y="dxt", parameter="muOfY")

## ----Diagnostics1-------------------------------------------------------------
blos2 <- Bayes(losmcmc, y="Dx", parameter="muOfY")
interpret(blos2$Diagnostics, digits=5)

## ----Post1--------------------------------------------------------------------
blos1 <- Bayes(x=losmcmc, y="post", parameter="muOfY", newdata=TRUE)
print(blos1$Posterior.Summary)

## ----plotPost101, fig.dim = c(6, 4.5)-----------------------------------------
plot(x=blos1, y="post", parameter="muOfY", compare=4.5, rope=c(4,5), lcol= c("blue","red"),
bcol="goldenrod", HDItext=.3, main= "Summary of average LOS (muOfY)")

## ----plotPost1, fig.dim = c(6, 4.5)-------------------------------------------
plot(x=blos1, y="post", parameter=list("sigmaOfY", "muOfY" ),math="divide",
bcol="cyan", HDItext=.3, main= "Coefficient of Variation")

## ----plotPPC1, fig.dim = c(6, 4.5)--------------------------------------------
plot(x=blos1, y="check", type="n", data=hosprog, dv="los",
parameter=c("muOfY", "sigmaOfY"), breaks=30, cex.axis=1.3, lwd=3, xlab=NULL,
pline=20, vlim=c(-2, 20), xlim=c(-2, 20), add.legend="topright",
main="Length of Stay", cex.main=1.5, xpt=5, pcol="red", lcol="orange",
cex.legend=1, bcol="cyan")

## ----anovaPlot, echo=FALSE, out.width="60%", fig.align="center"---------------
knitr::include_graphics("taov.png") 

## ----PostInt1-----------------------------------------------------------------
# These make 4 levels for the intervention program * pre/post period
#pre/post indicator
hosprog$prepost <- ifelse(hosprog$month >= 5, 1, 0)
#Control group at pre-test
hosprog$Post.Int <- 1
#Control group at post-test
hosprog[, "Post.Int"][hosprog$program == 0 & hosprog$prepost == 1 ] <- 2
#Intervention group at pre-test
hosprog[, "Post.Int"][hosprog$program == 1 & hosprog$prepost == 0 ] <- 3
#Intervention group at post-test
hosprog[, "Post.Int"][hosprog$program == 1 & hosprog$prepost == 1 ] <- 4

## ----DIDVary------------------------------------------------------------------
bvlos <- Bayes(x=losvary, newdata=TRUE)
#key parameters for post check
parls <- list(c("muOfLogY[1]", "sigmaOfLogY[1]"),  
              c("muOfLogY[2]", "sigmaOfLogY[2]"),
              c("muOfLogY[3]", "sigmaOfLogY[3]"),
              c("muOfLogY[4]", "sigmaOfLogY[4]"))

## ----plotDIDvary, fig.dim = c(6, 4.5)-----------------------------------------
plot(x=bvlos, y="post", parameter="B3DIDSig", bcol="green", compare=0,
     HDItext=.3, cex=2, main= "B3DIDSig")

## ----plotVary1, fig.dim = c(6, 4.5)-------------------------------------------
plot(x=bvlos, y="vary", type="ln", dv="los", breaks=75, pline=3, lwd=1, vlim=c(0, 15),
     data=hosprog, lcol= c("pink","red", "cyan", "blue"), xlim=c(0.5, 15), ylim=c(0, .35), parameter=parls, main="LOS and posterior pred checks", cex.main=2, cex.legend=1.5, cex.label=1.5,
     cex.axis=1.5, add.legend="topright", legend=c("Observed Data", "CTL Pre Est.", "CTL Post Est.", "Treat Pre Est.", "Treat Post Est."))

## ----plotVary2, fig.dim = c(6, 4.5)-------------------------------------------
plot(x=bvlos, y="vary", type="ln", dv="los", breaks=75, pline=1, lwd=5, vlim=c(0, 15),
     data=hosprog, lcol= c("pink","red", "cyan", "blue"), xlim=c(10, 15), ylim=c(0, .03), parameter=parls, main="LOS Zoomed In", cex.main=2, cex.legend=1.5, cex.label=1.5,
     cex.axis=1.5, add.legend="topright", legend=c("Observed Data", "CTL Pre Est.", "CTL Post Est.", "Treat Pre Est.", "Treat Post Est."))

## ----plotDIDavg1, fig.dim = c(6, 4.5)-----------------------------------------
plot(x=bvlos, y="post", parameter="B3DIDMu", bcol="yellow",
     HDItext=.3, cex=2, main= "DID for the mean")

## ----plotDIDavg2, fig.dim = c(6, 4.5)-----------------------------------------
plot(x=bvlos, y="post", parameter="B3DIDMed", bcol="tomato",
     HDItext=.3, cex=2, main= "DID for the median")

## ----plotDIDavg3, fig.dim = c(6, 4.5)-----------------------------------------
plot(x=bvlos, y="post", parameter="B3DIDMode", bcol="slategray",
     HDItext=.3, cex=2, main= "DID for the mode")

## ----Check1-------------------------------------------------------------------
bco2 <- Bayes(x=co2mcmc, y='mcmc', newdata=TRUE )

## ----plotPPC2, fig.dim = c(6, 4.5)--------------------------------------------
plot(x=bco2, y="check", type="ol", data=CO2, dv="uptake", iv="conc",
parameter=c("b0","b1"), add.data="al", cex.axis=1.3, lwd=1.5, pline=50,
vlim=c(50, 1100), xlim=c(0, 1100), ylim=c(0, 50), cex=2, cex.lab=2,
pcol="magenta", cex.main=2,cex.legend=1.2,  add.legend="topleft",
lcol="steelblue")              #vlim lets me extrapolate a little

## ----plotMulti1, fig.dim = c(6, 4.5)------------------------------------------
plot(x=co2multi, y="multi", level=2, aorder=FALSE,
subset= c("Qn2","Qn3","Qc3","Qc2","Mn3","Mn2","Mc2","Mc3"),
lcol="blue", pcol= c("red", "skyblue"), round.c=2, bcol="yellow",
xlim=c(-.1, 1), legend=NULL, add.legend="topright", lwd=3, cex.lab=1.2,
cex= 2, cex.main=1.25, cex.axis=.75, cex.legend=1.5, xlab=NULL)

## ----plotMulti2, fig.dim = c(6, 4.5)------------------------------------------
plot(x=co2multi, y="multi", level=3, aorder=FALSE, lcol="blue", pcol= c("green", "pink"),
round.c=3, bcol="lavender", xlim=c(-.1, 1), legend=NULL, add.legend="bottomleft", lwd=3, cex.lab =1.2, cex= 2, cex.main=1.25, cex.axis=.75, cex.legend=1.25, xlab=NULL)

## ----Target1------------------------------------------------------------------
btarget1 <- Bayes(x=losmcmc, y="target", type="n", parameter=c("muOfY","sigmaOfY"),
newdata=TRUE, targets=list(p=c(.35,.4,.45, .5, .55),  y=c(3,4), e= list(a=c(.35,.4,.45), b=.5))) 
print(btarget1$Target)

## ----Interpret Target1--------------------------------------------------------
interpret(btarget1$Target, digits=3)

## ----plotTarget2, fig.dim = c(7.5, 4.5)---------------------------------------
plot(x=btarget1, y="target", type="n", data=hosprog, dv="los", breaks=30,
cex.axis=1.3, lwd=1.5, pline=20, vlim=c(1, 12), xlim=c(1, 9),
parameter=c("muOfY","sigmaOfY"), add.legend="right", main="Length of Stay",
cex.main=1.5, xpt=5, pcol="black", lcol="salmon", tgtcol="blue", bcol="gray90",
cex.legend=1.25, cex.text = 1.5)

## ----plotTarget1, fig.dim = c(6, 4.5)-----------------------------------------
plot(x=btarget1, y="target", type="n", lcol="purple", tgtcol="blue", xlim=c(3.5, 5))

## ----R21----------------------------------------------------------------------
bR2 <- Bayes(x=co2mcmc, y='r2', data=CO2, iv="uptake", parameter=c("b0", "b1", "sigma"))
# R^2
print(bR2$R2.Summary$R2)
# Variance of predicted outcome
print(bR2$R2.Summary$Variance.Pred.Y)      
# Variance of residuals
print(bR2$R2.Summary$Variance.Residuals)
# A few predicted outcome values
print(head(bR2$R2.Summary$yPRED))