## ----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") ## ----group-------------------------------------------------------------------- gr1 <- group(x="program", y="los", z="month", data=hosprog, dist="t", cluster=TRUE) print(gr1$Group.CI) ## ----group2------------------------------------------------------------------- gr2 <- group(x="age", y="los", z="month", data=hosprog, dist="t", increment=3, rolling=6, quarts=TRUE) print(gr2$Roll.CI$Rolling) ## ----plotGroup1, fig.dim = c(6, 4.5)------------------------------------------ plot(x=gr1, y="group", order="numeric", lwd=4, gcol= "blue", pcol="red", overall=TRUE, oband=TRUE, ocol="gray", tcol="green", tgt=4.5, cex=2, cex.axis=1, cex.lab=1.1, cex.text=2, cex.main=1.25, adj.alpha=.2) ## ----clustering--------------------------------------------------------------- print(gr1$Clustering) ## ----plotGroup2, fig.dim = c(6, 4.5)------------------------------------------ plot(x=gr1, y="time", lwd=4, gcol=c("red", "blue"), gband=TRUE, overall=TRUE, oband=TRUE, ocol="gray", tgt=4.5, tcol="green", tpline=5, tpcol="yellow", name=TRUE, cex.axis=1, cex.lab=1, cex.text=2, cex.main=1, adj.alpha=.3, x.axis=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) ## ----ols---------------------------------------------------------------------- # Model summary summary(assess(hp ~ mpg+wt, data=mtcars, regression="ols")$model) ## ----ols_interpret------------------------------------------------------------ # Model summary interpret(assess(hp ~ mpg+wt, data=mtcars, regression="ols"))$model ## ----lrm---------------------------------------------------------------------- # Model summary summary(assess(am ~ hp+wt, data=mtcars, regression="logistic")$model) ## ----lrm_interpret------------------------------------------------------------ # Specify how many significant digits with digits argument interpret(assess(am ~ hp+wt, data=mtcars, regression="logistic"))$model ## ----pm1---------------------------------------------------------------------- # Model summary summary(assess(HAI ~ Month + offset(log(PatientDays)), data = infections, regression="poisson", link="log")$model) ## ----pm1_interpret------------------------------------------------------------ # Specify how many significant digits with digits argument interpret(assess(HAI ~ Month + offset(log(PatientDays)), data = infections, regression="poisson", link="log"))$model ## ----ols topcoding propensity------------------------------------------------- m1 <- assess(formula=cost ~ month * program, data=hosprog, intervention = "program", regression="ols", topcode=17150, propensity=c("female","age","risk"), newdata=TRUE) ## ----ols topcoding propensity results----------------------------------------- summary(m1$model) ## ----summary topcoding propensity--------------------------------------------- summary(m1$newdata[, c( "cost","top.cost", "pscore")]) ## ----ols ipw weights---------------------------------------------------------- ipw1 <- assess(formula=cost ~ month * program, data=hosprog, intervention = "program", regression="ols", topcode=17150, propensity=c("female","age","risk"), newdata=TRUE, weights="ipw") ## ----ols ipw results---------------------------------------------------------- summary(ipw1$model) ## ----importance--------------------------------------------------------------- importance(m1$model) ## ----plotImportance, fig.dim = c(6, 4.5)-------------------------------------- #Consider using these graphical parameters par(mar=c(4.2, 2, 3.5, 3)) par(oma = c(0, 0, 0, 3)) plot(importance(m1$model)) ## ----plotreview1, fig.dim = c(6, 4.5)----------------------------------------- m01 <- lm(mpg ~ wt+hp+am, data=mtcars) # Most basic code plot(review(m01)) ## ----review model 2 results--------------------------------------------------- m02 <- assess(formula=mpg ~ wt+hp+am, data=mtcars, regression="ols") # Using the assess function, notice 'm02$model' object below print(review(m02$model), digits=4) ## ----plotreview2, fig.dim = c(6, 4.5)----------------------------------------- m02 <- assess(formula=mpg ~ wt+hp+am, data=mtcars, regression="ols") plot(x=review(m02$model, increase=c(hp= 83.5))) ## ----plotpartpred1, fig.dim = c(6, 4.5)--------------------------------------- m02 <- assess(formula=mpg ~ wt+hp+am, data=mtcars, regression="ols") plot(x=m02, y='hp') ## ----plotreview3, fig.dim = c(6, 4.5)----------------------------------------- m03 <- glm(vs ~ wt+hp+am, data=mtcars, family="binomial") # Display options, sorted plot, and printing the 95% CIs to review plot(x=review(m03, increase=c(wt=1.1375, hp= 83.5) ), color="aquamarine", lwd=5, pt.cex= 2, tcol="magenta", sort="coef", cex.axis=2, cex.main=2, xlim=c(-.5, 3), print=TRUE, round.c=6) ## ----plotpartpred2, fig.dim = c(6, 4.5)--------------------------------------- m03 <- assess(vs ~ wt+hp+am, data=mtcars, regression="logistic") plot(x=m03, y=c('wt', 'lp'), col='slategray', main="Logisitic Regression: Engine", adj.alpha=.2, cex.axis=1.5) ## ----plotreview4, fig.dim = c(6, 4.5)----------------------------------------- m04 <- assess(formula=HAI ~ Month+ offset(log(PatientDays)), data=infections, regression="poisson") #Because 1 year is a meaningful period in program evaluation, Month=12 plot(x=review(m04$model, increase=c(Month=12)), xlim=c(0.6, 1.01), lwd=7, color="cyan", pcol="salmon", pt.cex=2) ## ----plotpartpred3, fig.dim = c(6, 4.5)--------------------------------------- m04 <- assess(formula=HAI ~ Month+ offset(log(PatientDays)), data=infections, regression="poisson") plot(x=m04, y=c('Month', 'count'), col='skyblue', tgt=100, tgtcol='magenta', lwd=4) ## ----plotpartpred4, fig.dim = c(6, 4.5)--------------------------------------- infections$rand <- as.factor(sample(1:3, 41, replace = T)) m05 <- assess(HAI ~ Month + rand + offset(log(PatientDays)), data = infections, regression="poisson") plot(x=m05, y=c('rand', 'count'), col='skyblue', tgt=75, tgtcol='orange', lwd=4, cex=2) ## ----plotpartpred140b, fig.dim = c(6, 4.5)------------------------------------ set.seed(1) infections$rand <- as.factor(sample(1:3, 41, replace = T)) m05 <- assess(HAI ~ Month + rand + offset(log(PatientDays)), data = infections, regression="poisson") plot(x=m05, y=c('Month', 'count', 'rand'), col=c('purple', 'gold','green'), add.legend= "topright", lwd=4, cex=2) ## ----did model 1-------------------------------------------------------------- dm1 <- assess(formula= los ~ ., data=hosprog, intervention = "program", int.time="month", treatment= 5, did="two") ## ----did model 1 results------------------------------------------------------ summary(dm1$DID) ## ----interpret did model 1---------------------------------------------------- interpret(dm1)$did ## ----plotDID1, fig.dim = c(6, 4.5)-------------------------------------------- plot(x=dm1, y="DID", add.legend="topleft", xlim=c(-.05, 1.05), ylim=c(2, 9), main="DID: Length of Stay", col=c("cyan","magenta"), lwd=7, cex=2, cex.axis=2, cex.lab=1.5, cex.main=3, arrow=TRUE, xshift=c(.045), cex.text=1.5, coefs=TRUE, round.c=2, cfact=T, conf.int=TRUE, adj.alpha=0.2 ) ## ----did model 3-------------------------------------------------------------- dm3 <- assess(formula= rdm30 ~ ., data=hosprog, intervention = "program", int.time="month", treatment= 5, did="two") ## ----did model 3 results------------------------------------------------------ summary(dm3$DID) ## ----its model 1-------------------------------------------------------------- im11 <- assess(formula=los ~ ., data=hosp1, intervention = "program", int.time="month", interrupt= 5, its="one") ## ----its model 1 results------------------------------------------------------ summary(im11$ITS) ## ----its model 1 interpretations---------------------------------------------- interpret(im11)$its ## ----itsEffect model 1-------------------------------------------------------- itsEffect(model= im11$ITS, type= "sgst", interruptions= 1) ## ----plotITS1, fig.dim = c(6, 4.5)-------------------------------------------- plot(im11, "ITS", add.legend="topright", xlim=c(-1, 14), ylim=c(2, 9), main="ITS: Intervention LOS", col="thistle", lwd=7, cex=3, cex.axis=2, cex.lab=1.5, cex.main=3, arrow=TRUE, xshift=c(.25, .25), cex.text=1.5, coefs=TRUE, round.c=2, cfact=T, conf.int=TRUE, adj.alpha=0.2, pos.text= list("ITS.Time"=3, "Intercept"=4), cex.legend=1.25, add.means=TRUE ) ## ----its model 2-------------------------------------------------------------- im12 <- assess(formula=los ~ ., data=hosp1, intervention = "program", int.time="month", interrupt= c(5, 9), its="one") ## ----its model 2 results------------------------------------------------------ summary(im12$ITS) ## ----its model 2 interpretations---------------------------------------------- interpret(im12)$its ## ----itsEffect model 2-------------------------------------------------------- itsEffect(model= im12$ITS, type= "sgmt", interruptions= 2) ## ----plotITS02, fig.dim = c(6, 4.5)------------------------------------------- plot(im12, "ITS", add.legend="topright", xlim=c(-1, 14), ylim=c(2, 9), main="ITS: Intervention LOS", col=("thistle"), lwd=7, cex=3, cex.axis=2, cex.lab=1.5, cex.main=3, arrow=TRUE, xshift=c(.25, .25), cex.text=1.0, coefs=TRUE, round.c=2, cfact=T, conf.int=TRUE, adj.alpha=0.2, pos.text= list("Intercept"=1, "txp5"=3, "post5"=3, "post9"=3), cex.legend=1.25, add.means=TRUE ) ## ----its model 3-------------------------------------------------------------- im22 <- assess(formula=los ~ ., data=hosprog, intervention = "program", int.time="month", interrupt= c(5, 9), its="two") ## ----plotreview5, fig.dim = c(6, 4.5)----------------------------------------- im22 <- assess(formula=los ~ ., data=hosprog, intervention = "program", int.time="month", interrupt= c(5, 9), its="two") # The intervention group had the biggest change between the baseline and month 5 plot(review(im22$ITS), sort="coef", decreasing=TRUE, main="ITS review", color="red", pcol="green", lwd=2, pt.cex=1.5) ## ----its model 3 results------------------------------------------------------ summary(im22$ITS) ## ----its model 3 interpretations---------------------------------------------- interpret(im22)$its ## ----itsEffect model 3-------------------------------------------------------- itsEffect(model= im22$ITS, type= "mgmt", interruptions= 2) ## ----plotAssess, fig.dim = c(6, 4.5)------------------------------------------ plot(im22, "ITS", add.legend="top", xlim=c(-.75, 13.1), ylim=c(2, 9), main="ITS: Length of Stay", col=c("springgreen","thistle"), lwd=7, cex=2, cex.axis=2, cex.lab=1.5, cex.main=3, cex.legend=1.25, arrow=TRUE, xshift=c(0, .5), cex.text=1, coefs=TRUE, round.c=1, pos.text= list("Intercept"=1, "txp5"=3, "post5"=3, "post9"=3), tcol="dodgerblue", conf.int=TRUE, adj.alpha=0.3, add.means=TRUE) ## ----its model 4-------------------------------------------------------------- id22 <- assess(formula=death30 ~ ., data=hosprog, intervention = "program", int.time="month", interrupt= c(5, 9), its="two") ## ----its model 4 results------------------------------------------------------ summary(id22$ITS) ## ----its model 4 interpretations---------------------------------------------- interpret(id22)$its ## ----itsEffect model 4-------------------------------------------------------- itsEffect(model= id22$ITS, type= "mgmt", interruptions= 2) ## ----plotITS04, fig.dim = c(6, 4.5)------------------------------------------- plot(id22, "ITS", xlim=c(-1, 14), ylim=c(0, .37), main="ITS: 30-day Death", col=c("sienna3","slategray"), lwd=3, cex=2, cex.axis=2, cex.lab=1.5, tcol="green", cex.main=3, arrow=TRUE, xshift=c(.25, .25), cex.text=1.0, coefs=TRUE, round.c=2, cfact=T, conf.int=TRUE, adj.alpha=0.15, pos.text= list("Intercept"=1, "txp5"=3, "txi"=2, "post5"=3, "post9"=4, "ITS.Time"=3), cex.legend=.9, add.means=TRUE ) ## ----its model 5-------------------------------------------------------------- #Key interruption periods key_time <- c(5, 14, 17, 29, 42, 59, 69, 73, 80,92) im10 <- assess(formula=rate ~ ., data=unemployment, intervention = "usa", int.time="year", its="one", interrupt= key_time, newdata=TRUE) ## ----plotITS2, fig.dim = c(7.5, 5.25)----------------------------------------- plot(im10, "ITS", add.means = TRUE, coefs=TRUE, conf.int=TRUE, adj.alpha= .2, lwd=1.75, col="slategray", tcol= "orange", main="US unemployment rate", xlab="Years (1929-2024)", ylab= "Proportion of labor market", cex.main=2, cex.axis = 1.25, cex.lab = 1.25, cex=2, cex.text= .75, pos.text=list("ITS.Time"=4, "post42"=1,"txp42"=3,"txp92"=3), x.axis=unemployment$Year) for(i in 1:length(key_time)) { text(key_time[i], .22-(.01*i), cex=.85, labels = paste0(unemployment[ key_time[i], "Year"], ": ", unemployment[ key_time[i], "event"])) } ## ----cronbachs alpha example-------------------------------------------------- alpha(items=c("i1","i2","i3","i4","i5"), data=cas) ## ----cronbachs alpha interpret------------------------------------------------ interpret(alpha(items=c("i1","i2","i3","i4","i5"), data=cas))