Title:	Healthcare Analysis Methods
Version:	1.2.0
Description:	Conducts analyses for healthcare program evaluations or intervention studies. Calculates regression analyses for standard ordinary least squares (OLS or linear) or logistic models. Performs regression models used for causal modeling such as differences-in-differences (DID) and interrupted time series (ITS) models. Provides limited interpretations of model results and a ranking of variable importance in models. Performs propensity score models, top-coding of model outcome variables, and can return new data with the newly formed variables. Conducts Bayesian analysis summaries and graphs, decision curve analysis, and produces some Shewhart control charts. Also performs Cronbach's alpha for various scale items (e.g., survey questions). See Github URL for examples in the README file. For more details on the statistical methods, see Allen & Yen (1979, ISBN:0-8185-0283-5), Angrist & Pischke (2009, ISBN:9780691120355), Cohen (1988, ISBN:0-8058-0283-5), Gebski (2012) <doi:10.1017/S0950268812000179>, Gelman & Goodrich (2019) <doi:10.1080/00031305.2018.1549100>, Harrell (2016, ISBN:978-3-319-19424-0), Kline (1999, ISBN:9780415211581), Kruschke (2014, ISBN:9780124058880), Linden (2015) <doi:10.1177/1536867X1501500208>, Merlo (2006) <doi:10.1136/jech.2004.029454>, Muthen & Satorra (1995) <doi:10.2307/271070>, Rabe-Hesketh & Skrondal (2008, ISBN:978-1-59718-040-5), Ryan (2011, ISBN:978-0-470-59074-4), and Vickers & Elkin (2006) <doi:10.1177/0272989X06295361>.
License:	MIT + file LICENSE
Encoding:	UTF-8
RoxygenNote:	7.3.2
Depends:	R (≥ 3.5.0)
LazyData:	true
URL:	https://github.com/szuniga07/ham
BugReports:	https://github.com/szuniga07/ham/issues
Imports:	methods
Suggests:	knitr, rmarkdown, testthat (≥ 3.0.0)
Config/testthat/edition:	3
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2026-03-19 23:49:39 UTC; szuni
Author:	Stephen Zuniga [aut, cre, cph]
Maintainer:	Stephen Zuniga <rms.shiny@gmail.com>
Repository:	CRAN
Date/Publication:	2026-03-20 00:20:02 UTC

Summarize Bayesian Markov Chain Monte Carlo (MCMC) object

Description

Convert a list of Bayesian analysis chains (e.g., coda package mcmc.list objects) into a data frame for analysis and creating plots. For example, models from JAGS or Stan that were converted into coda class objects can be used to create the data frames. Calculates a set of descriptive statistics that summarize MCMC parameters. And values can be calculated for use in descriptive graphs such as values associated with specific percentiles and vice-versa to help set targets, summaries on hierarchical or multilevel models (up to 3 levels), and the R2 for Bayesian regression models with metric level predictors.

Usage

Bayes(
  x,
  y = "mcmc",
  parameter = NULL,
  mass = 0.95,
  compare = NULL,
  rope = NULL,
  newdata = FALSE,
  type = NULL,
  center = "mode",
  data = NULL,
  dv = NULL,
  iv = NULL,
  expand = NULL,
  targets = NULL
)

Arguments

Value

data frame of summary statistics for the MCMC parameter's distribution and/or MCMC data frame. Statistics include highest density interval, effective sample size, proportion of distribution within and outside of a ROPE, distribution compared with a set value, and the parameter's mean, median, and mode. And distribution summaries for multilevel models, target summaries, and regression model R2.

References

Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences, Second Edition. Hillsdale, NJ: Lawrence Erlbaum Associates, Publishers. ISBN 0-8058-0283-5

Gelman, A., Goodrich, B., Gabry, J., & Vehtari, A. (2019). R-squared for Bayesian Regression Models. The American Statistician, 73, 3, 307–309. https://doi.org/10.1080/00031305.2018.1549100

Kruschke, J. (2014). Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan, Second Edition. New York: Academic Press. ISBN: 9780124058880

Examples

# Posterior estimates of length of stay (LOS) #
# What is the 'average' LOS and are we beyond our goal of it being <= 3 days?
blos1 <- Bayes(losmcmc, y="post", parameter="muOfY", newdata=TRUE,
compare=4.5, rope=c(1,3))   #Where are we in relation to 4.5 days?
print(blos1$Posterior.Summary) #looks like LOS is substantially higher than 3 days

# Multilevel or hierarchical model summaries #
# Code below does not run, no 'mcmc_sample' object
# bmulti0 <- Bayes(x=mcmc_sample, parameter=c("theta", "omega","omegaO"),
# y="multi", type="bern", data=mydf, dv="upbin", iv= c("Plant", "Group"))

# Targets for length of stay (LOS) #
# Our administrators ask how far are we from our goals, they ask about targets
# in increments of 5 points of probability or specific days. We answer both.
btarget1 <- Bayes(x=losmcmc, y="target", type="n", parameter=c("muOfY","sigmaOfY"),
newdata=TRUE, target=list(p=c(.35,.4,.45, .5, .55),  y=c(3,4))) # 'newdata' for plots
print(btarget1$Target$Est.Quantile.P)  # 5-point increments
print(btarget1$Target$Est.Prob.GT.Y)   # specific day values, 4 days is plausible

# Gelamn's R^2 #
# the regression model using Base R data, CO2: update ~ conc
bR2 <- Bayes(x=co2mcmc, y='r2', data=CO2, iv="uptake",
parameter=c("b0", "b1", "sigma"), newdata=TRUE)
# bR2 returns various information
print(bR2$R2.Summary$R2)                   # R^2
print(bR2$R2.Summary$Variance.Pred.Y)      # Variance of predicted outcome
print(bR2$R2.Summary$Variance.Residuals)   # Variance of residuals
print(head(bR2$R2.Summary$yPRED))          # Predicted outcomes

Calculates Cronbach's alpha on scale items

Description

Performs Cronbach's alpha of specified items from a data frame. Cronbach's Alpha is a formula for estimating the internal consistency reliability of a measurement instrument such as survey items (see Allen & Yang, 1979; Kline, 1999). Survey items can have 2 or more categories such as 5-point scales and contain 2 or more items.

Usage

alpha(items, data)

Arguments

Value

A list object with Cronbach's alpha summary statistics.

References

Allen, M. J., & Yen, W. M. (1979). Introduction to Measurement Theory. Brooks/Cole. ISBN: 0-8185-0283-5. Kline, Paul (1999). Handbook of Psychological Testing (2nd ed). Routledge, New York. ISBN: 9780415211581.

Examples

alpha(items=c("i1","i2","i3","i4","i5"), data=cas)

# remove i1 as suggested in the previous example, returns higher alpha
alpha(items=c("i2","i3","i4","i5"), data=cas)

Assess models with regression

Description

Fit ordinary least squares (OLS) and logistic models. And fit models for causal inference such as differences-in-differences and interrupted time series. Run these models to evaluate program performance or test intervention effects (e.g., healthcare programs). Options are available for top coding the outcome variable as well as propensity scores. New data can optionally be returned that has these additional variables and constructed variables that are used for DID and ITS models.

Usage

assess(
  formula,
  data,
  regression = "none",
  did = "none",
  its = "none",
  intervention = NULL,
  int.time = NULL,
  treatment = NULL,
  interrupt = NULL,
  subset = NULL,
  stagger = NULL,
  topcode = NULL,
  propensity = NULL,
  newdata = FALSE
)

Arguments

Value

a list of results from selected regression models. Will return new data if selected. And returns relevant model information such as variable names, type of analysis, formula, study information, and summary of ITS effects if analyzed.

References

Angrist, J. D., & Pischke, J. S. (2009). Mostly Harmless Econometrics: An Empiricist's Companion. Princeton University Press. ISBN: 9780691120355.

Gebski, V., et al. (2012). Modelling Interrupted Time Series to Evaluate Prevention and Control of Infection in Healthcare. Epidemiology & Infections, 140, 2131–2141. https://doi.org/10.1017/S0950268812000179

Linden, A. (2015). Conducting Interrupted Time-series Analysis for Single- and Multiple-group Comparisons. The Stata Journal, 15, 2, 480-500. https://doi.org/10.1177/1536867X1501500208

Examples

# ordinary least squares R^2
summary(assess(hp ~ mpg+wt, data=mtcars, regression="ols")$model)

# logistic
summary(assess(formula=vs~mpg+wt+hp, data=mtcars, regression="logistic")$model)

# OLS with a propensity score
summary(assess(formula=los ~ month+program, data=hosprog, intervention = "program",
regression="ols", propensity=c("female","age","risk"))$model)

# OLS: top coding los at 8.27 and propensity score means (top.los and pscore)
summary(assess(formula=los ~ month+program, data=hosprog, intervention = "program",
regression="ols", topcode=8.27, propensity=c("female","age","risk"),
newdata=TRUE)$newdata[, c("los", "top.los", "pscore")])

# differences-in-differences model: using 2 time periods, pre- and post-intervention
summary(assess(formula=los ~ ., data=hosprog, intervention = "program",
int.time="month", treatment = 5, did="two")$DID)

# DID model: using time points
summary(assess(formula=los ~ ., data=hosprog, intervention = "program",
int.time="month", treatment = 5, did="many")$DID)

#interrupted time series model: two groups and 1 interruption (interrupt= 5)
summary(assess(formula=los ~ ., data=hosprog, intervention = "program",
int.time="month", its="two", interrupt = 5)$ITS)

#interrupted time series model: two groups and 2 interruptions (interrupt= c(5,9))
summary(assess(formula=los ~ ., data=hosprog, intervention = "program",
int.time="month", its="two", interrupt = c(5,9))$ITS)

Patient survey data

Description

Artificial data of a 5 item hospital satisfaction survey for a Cronbach's alpha scale (cas).

Usage

cas

Format

cas

An artificial data frame with 100 rows and 5 columns:

i1 - i5

5 survey items

...

Source

Artificial dataset created with rbinom for 5 items. For example, rbinom(100, 5, .9) generates 1 item. The prob argument is modified to give more or less consistent ratings per item.

Markov Chain Monte Carlo linear regression estimates of plant's CO2 uptake regressed on ambient carbon dioxide concentrations

Description

Markov Chain Monte Carlo linear regression estimates of plant's CO2 uptake regressed on ambient carbon dioxide concentrations

Usage

co2mcmc

Format

co2mcmc

MCMC list of 3 chains with 2000 rows and 4 columns in each of the list elements:

b0

Regression intercept.

b1

Regression coefficient of the carbon dioxide concentration effect.

sigma

Square root of model residual variance.

lp__

Total log unnormalized joint posterior density of the model parameters.

...

Source

co2mcmc is a list of 3 MCMC simulations. It estimates a linear regression using the CO2 data frame from base R. A regression model was specified using a normal distribution maximum likelihood and normal and uniform distribution priors for the beta coefficients and sigma parameters.

Bayes class object with summarized Markov Chain Monte Carlo estimates of plant's CO2 uptake as a binary variable that is above or below the median level

Description

Bayes class object with summarized Markov Chain Monte Carlo estimates of plant's CO2 uptake as a binary variable that is above or below the median level

Usage

co2multi

Format

co2multi

Bayes class list of 7 elements with varying rows and columns in each of the list elements:

Posterior.Summary

Standard element that is NA in this object because it is not applicable to this analysis.

MCMC

Dataframe of 4 combined MCMC chain simulations with 8000 rows and 18 columns. Column 1 is CHAIN indicates the 4 chains, columns 2-13 are 'theta' and indicates each plant's rate of having above the median uptake level, columns 14-17 are 'omega' group rates of 'Type' by 'Treatment', and column 18 is the overall group rate.

Multilevel

Posterior estimate summaries in alphabetical order of plant or group names and numerical order plant or group rates.

Target

Standard element that is NA in this object because it is not applicable to this analysis.

targets

Standard element that is NA in this object because it is not applicable to this analysis.

R2.Summary

Standard element that is NA in this object because it is not applicable to this analysis.

parameter

3 parameter names in the multilevel analysis, 'theta', 'omega', and 'omegaO'.

...

Source

co2multi is a list of summarized Markov Chain Monte Carlo estimates of plant's CO2 uptake as a binary variable that is above or below the median level. It summarizes measurements nested within plants nested within the 4 levels from a combination of 'Treatment' by 'Type'. An estimation used the CO2 data frame from base R. A hierarchical estimation model was specified using a gamma distribution maximum likelihood and normal and uniform distribution priors at various group levels.

Statistics for Shewhart control charts

Description

Calculate statistics that can be used to produce X-bar charts, p-charts, and u-charts. This includes producing means for center lines, 3-sigma upper and lower control limits. Users can also calculate values before and after an intervention to see if a change in the control process happened. Values are returned in a data frame.

Usage

control(
  x,
  y = NULL,
  time,
  data,
  type = "x",
  subset = NULL,
  n.equal = TRUE,
  intervention = NULL
)

Arguments

Value

data frame of control chart statistics for X-bar charts, p-charts, and u-charts. Includes means, standard deviations, and 3-sigma upper and lower control limit values.

References

Ryan, T. (2011). Statistical Methods for Quality Improvement, Third Edition. New Jersey: John Wiley & Sons, Inc. ISBN: 978-0-470-59074-4.

Examples

## Hospital LOS and readmissions ##
# X-bar chart statistics
spc_x <- control(x="los", time="month", data=hosprog, type="x", n.equal=TRUE)
print(spc_x) # get data frame output

# X-bar chart statistics not assuming equal sample sizes, subsetting for females
spc_x <- control(x="los", time="month", data=hosprog, type="x", n.equal=FALSE,
  subset=hosprog$female==1)
print(spc_x) # get data frame output

# p-chart statistics, using only the numerator (i.e., y=NULL). Specify unequal sample sizes
spc_p <- control(x="rdm30", time="month", data=hosprog, type="p", n.equal=FALSE)
print(spc_p) # get data frame output

# u-chart for infection rates with an intervention at the 22nd month
spc_u <- control(x="HAI", y="PatientDays", time="Month", data=infections,
type="u", n.equal=FALSE, intervention=22)
print(spc_u) # get data frame output

Statistics for Decision Curve Analysis and logistic regression model classification

Description

Calculate statistics such as sensitivity, specificity, positive and negative predictive values, and net benefit and interventions saved from a decision curve analysis.

Usage

decide(x, threshold)

Arguments

Value

summary of model classification based on the selected threshold/cutoff, area under the curve (AUC)/c-statistic, predicted outcomes (transformed if applicable), decision curve analysis values at various percentiles, and sensitivity and specificity related statistics for the regression model at the specified threshold.

References

Vickers, A. & Elkin, E. (2006). Decision Curve Analysis: A Novel Method for Evaluating Prediction Models. Society for Medical Decision Making, 26, 6, 565-574. https://doi.org/10.1177/0272989X06295361

Examples

## Predicting car engine shape type, v or straight  ##
# run the model
car_m1 <- assess(formula=vs ~ hp + am, data=mtcars, regression="logistic")
# create a decide object, enter the model name and a threshold on the logit scale
d1 <- decide(x=car_m1, threshold= -0.767)
# View model classification related statistics
print(d1$Model.Summary$Classification)

# View decision curve analysis results like 'net benefit' at various thresholds
print(d1$DCA)

Group level confidence intervals and between-group variation

Description

Group level confidence intervals and between-group variation

Usage

group(
  x,
  y,
  z = NULL,
  data,
  dist = "t",
  conf.int = 0.95,
  increment = 1,
  rolling = NULL,
  quarts = FALSE,
  cluster = FALSE,
  subset = NULL,
  asis = FALSE,
  dataf = NULL
)

Arguments

Value

list of confidence intervals for outcomes by groups, over time, and clustering measures. Some values returned in alphabetical and numerical order based on the group.

References

Merlo, J. (2006). A brief conceptual tutorial of multilevel analysis in social epidemiology: using measures of clustering in multilevel logistic regression to investigate contextual phenomena. Journal of Epidemiological Health, 60, 4, 290-297. https://doi.org/10.1136/jech.2004.029454.

Muthen, B. & Satorra, A. (1995). Complex Sample Data in Structural Equation Modeling. Sociological Methodology, 25, 267-316. https://doi.org/10.2307/271070.

Rabe-Hesketh, S. & Skrondal, A. (2008). Multilevel and Longitudinal Modeling Using Stata, Second Edition. ISBN: 978-1-59718-040-5.

Examples

#default t distribution results
group(x="program", y="los", data=hosprog)
#Rounding LOS to integers
hp2 <- hosprog; hp2$los2 <- round(hp2$los, 0)
#Exact Poisson confidence intervals
group(x="program", y="los2", data=hp2, dist="p")
#Rolling 6-months of data
group(x="program", y="los", z="month", data=hosprog, dist="t", rolling=6)
#Data returned separately for rolling 6-months of data and 3-month increments (e.g., quarters)
group(x="program", y="los", z="month", data=hosprog, dist="t", increment=3, rolling=6)
#Quartile groups for continuous risk score and returned clustering info
group(x="risk", y="los", data=hosprog, quarts=TRUE, cluster=TRUE)
#Binomial distribution with less conservative 90% confidence intervals
group(x="risk", y="rdm30", data=hosprog, quarts=TRUE, dist="b", conf.int=0.90)

Patient hospital program/intervention data, intervention group only

Description

Patient hospital program/intervention data, intervention group only

Usage

hosp1

Format

hosp1

An artificial data frame with 352 rows and 10 columns, intervention patients only:

survey

Patient satisfaction survey mean score.

los

Hospital length of stay (los)

cost

Hospital stay cost

rdm30

Patient readmission within 30 days of discharge

death30

Patient death within 30 days of discharge

female

Patient sex, 1 indicates female, 0 otherwise

age

Patient age

risk

Patient health risk score ranging from 0 to 1

month

12 month indicator (1 to 12)

program

Indicates patient program participation. 1='yes', 0='no'

...

Source

hosp1 is a subset of the artificial dataset hosprog. It is the intervention group's data used for single group interrupted time series.

Patient hospital program/intervention data

Description

Patient hospital program/intervention data

Usage

hosprog

Format

hosprog

An artificial data frame with 720 rows and 10 columns:

survey

Patient satisfaction survey mean score.

los

Hospital length of stay (los)

cost

Hospital stay cost

rdm30

Patient readmission within 30 days of discharge

death30

Patient death within 30 days of discharge

female

Patient sex, 1 indicates female, 0 otherwise

age

Patient age

risk

Patient health risk score ranging from 0 to 1

month

12 month indicator (1 to 12)

program

Indicates patient program participation. 1='yes', 0='no'

...

Source

Artificial dataset created by using runif. The strength in the association between each variable is weighted by multiplying each subsequent predictor in increments of 1. For example, Y equals runif(720) multiplied by 1 plus runif(720) multiplied by 2 and so on. This allows some predictors to have stronger correlations with Y.

Importance of variables based on partial chi-square statistic

Description

Calculates partial chi-square (Wald chi-square for individual coefficients) from assess class objects. The importance is the partial chi-square minus its degrees of freedom based on the regression coefficients (Harrell, 2015). A higher chi-square indicates a larger effect by the predictors. Therefore, the rank of the chi-square can indicate which predictors can contribute more in explaining the variation in the outcome variable.

Usage

importance(model)

Arguments

Value

a data.frame object with partial X^2 summary statistics.

References

Harrell, F. E., Jr. (2016). Regression Modeling Strategies. Springer International Publishing. ISBN: 978-3-319-19424-0.

Examples

# OLS regression
importance(assess(mpg ~ hp + wt + cyl, data=mtcars, regression= "ols")$model)

# logistic regression
importance(assess(vs~mpg+wt+hp, data=mtcars, regression= "logistic")$model)

Hospital acquired infections (HAI) during 41 months.

Description

Hospital acquired infections (HAI) during 41 months.

Usage

infections

Format

infections

An artificial data frame with 41 rows and 3 columns:

Month

Calendar month as an integer, ranging from 1 to 41.

HAI

Count of hospital acquired infections (HAI) within 1 month.

PatientDays

Count of hospital patient days within 1 month (i.e., the sum of days in the hospital for all patients).

...

Source

infections is an artificial data frame of reasonable hospital acquired infections (HAI) estimates made entirely for the purpose of demonstrating control charts.

Interpret model output

Description

Provides simple interpretations of regression coefficients and Cronbach's alpha from assess and alpha function classes. The interpretations describe coefficients and significance values as well as modifying item scales. The interpretations are text comments associated with specific parameters of the various analyses.

Usage

interpret(object)

Arguments

Value

a list with interpretations of Cronbach's alpha scales or regression model results.

Examples

# Interpret Cronbach's alpha
interpret(alpha(items=c("i1","i2","i3","i4","i5"), data=cas))

# interpret a standard linear (OLS) regression
hos1 <- assess(formula=survey ~ program + month, data=hosprog, regression= "ols")
interpret(hos1)$model

# interpret a differences-in-differences model
hos2 <- assess(formula=survey ~ ., data=hosprog, intervention = "program",
int.time="month", treatment = 5, did="two", newdata=TRUE)
interpret(hos2)$did  #interpret(hos2) also runs, returns ITS results if present

# interpret an interrupted time series model
hos3 <- assess(formula=survey ~ ., data=hosprog, intervention = "program",
int.time="month", its="two", interrupt = 5)
interpret(hos3)$its

Interrupted time series analysis effects

Description

Calculates effects for intervention and control groups based on interrupted time series models from an assess class object. Within a period (or interruption), the effect that represents the trend during the period is calculated for both groups as well as the difference between the groups. Summary statistics are provided that include the effect sizes, t-statistic, standard errors, p-values, and 95% confidence intervals of the effect sizes. These values are provided for the intervention group, control group (when applicable), and the differences between the two groups (Linden, 2015). These values are automatically generated while running a model in assess.

Usage

itsEffect(model, type, interruptions)

Arguments

Value

a data.frame object of ITS effects and summary statistics. Row names are indicated with first, second, and additional numbers when there are multiple interruptions (e.g., 'Intervention.2'). Generally run within the assess function.

References

Linden, Ariel. (2015). Conducting Interrupted Time-series Analysis for Single- and Multiple-group Comparisons. The Stata Journal, 2015, 15(2), 480-500, https://doi.org/10.1177/1536867X1501500208

Examples

i21 <- assess(formula=survey ~ ., data=hosprog, intervention = "program",topcode =NULL,
int.time="month", regression="none", interrupt=5, its="two", newdata=TRUE, propensity=NULL)
itsEffect(model= i21$ITS, type= "mgst", interruptions= 1)

Markov Chain Monte Carlo estimates of hospital length of stay from the hosprog data frame

Description

Markov Chain Monte Carlo estimates of hospital length of stay from the hosprog data frame

Usage

losmcmc

Format

losmcmc

MCMC list of 4 chains with 5000 rows and 2 columns in each of the list elements:

muOfY

Estimate of the mean parameter.

sigmaOfY

Estimate of the square root of variance parameter.

...

Source

losmcmc is a list of 4 MCMC simulations. It estimates the mean and standard deviation using the artificial hosprog data frame from ham. An estimation was specified using a normal distribution maximum likelihood and normal and uniform distribution priors for the mean and sigma parameters.

Bayesian plots for various analyses

Description

Graph Bayesian diagnostic traceplots, density plots on convergence, autocorrelation factor, calculate Monte Carlo Standard Errors, and effective sample sizes. And plot summaries of posterior distributions, posterior predictive checks, summaries on hierarchical or multilevel models (up to 3 levels), and summary graphs of values associated with specific percentiles (and vice-versa) that can be used to help set targets. Plots are developed from Bayes class objects that converted multiple chains into data frames.

Usage

## S3 method for class 'Bayes'
plot(
  x,
  y = NULL,
  type = "n",
  parameter = NULL,
  center = "mode",
  mass = 0.95,
  compare = NULL,
  rope = NULL,
  data = NULL,
  dv = NULL,
  iv = NULL,
  add.data = "n",
  group = NULL,
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  xlim = NULL,
  ylim = NULL,
  vlim = NULL,
  curve = FALSE,
  lwd = NULL,
  breaks = 15,
  bcol = NULL,
  lcol = NULL,
  pcol = NULL,
  xpt = NULL,
  tgt = NULL,
  tgtcol = "gray",
  tpline = NULL,
  tpcol = NULL,
  pline = 20,
  pct = 95,
  add.legend = NULL,
  legend = NULL,
  cex = 1,
  cex.lab = NULL,
  cex.axis = NULL,
  cex.main = NULL,
  cex.text = NULL,
  cex.legend = NULL,
  HDItext = 0.7,
  math = "n",
  es = "n",
  subset = NULL,
  level = NULL,
  aorder = TRUE,
  round.c = 2,
  ...
)

Arguments

Value

plots assessing model quality, posterior distributions, posterior predictive checks, hierarchical or multilevel model summaries, and target summaries.

References

Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences, Second Edition. Hillsdale, NJ: Lawrence Erlbaum Associates, Publishers. ISBN 0-8058-0283-5 Gelman, A., Goodrich, B., Gabry, J., & Vehtari, A. (2019). R-squared for Bayesian Regression Models. The American Statistician, 73, 3, 307–309. https://doi.org/10.1080/00031305.2018.1549100

Kruschke, J. (2014). Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan, Second Edition. New York: Academic Press. ISBN: 9780124058880

Examples

# Posterior estimates of length of stay #
# Set up the Bayes object first, convert the chains
blos1 <- Bayes(losmcmc, newdata=TRUE)

#' # Examine Bayesian model diagnostics for estimated length of stay (LOS) #
# Review traceplots, density plots, autocorrelation factor, and Gelman-Rubin statistic
plot(blos1, y="dxt", parameter="muOfY")
plot(blos1, y="dxd", parameter="muOfY")
plot(blos1, y="dxa", parameter="muOfY")
plot(blos1, y="dxg", parameter="muOfY")

# And we can calculate statistics that combine multiple parameters (e.g., calculate
# the mean difference between intervention and control groups). Here we get the
# coefficient of variation by dividing the standard deviation by the mean. We use
# the 'math' argument and arrange our parameters in the proper order.
# We only need the first 4 arguments but we'll add a little more.
plot(x=blos1, y="post", parameter=list("sigmaOfY", "muOfY" ),math="divide",
bcol="cyan", HDItext=.3, main= "Coefficient of Variation")

# Posterior Predictive Checks on how well our model fits the data #
# Estimating center and spread for hospital length of stay
# Generally, we only need the first 6 arguments but we'll modify the plot.
# A model with a gamma likelihood would fit better but this is ok for now.
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")

# Estimating the regression trend line
# Now lets look at the trend of conc on CO2 uptake from the CO2 data.
# Using a quadratic model with conc^2 would help and an option in ham.
# First, create the Bayes object
bco2 <- Bayes(x=co2mcmc, newdata=TRUE )
# We generally only need the first 7 arguments below.
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

# Hierarchical or Multilevel Model Summary #
# We generally only need the first 3 arguments below. But we'll subset
# on 8 of 12 plants in the level 2 model (observations nested in plants)
# and modify other settings.
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=1, bcol="yellow",
xlim=c(-.1, 1), legend=NULL, add.legend="topright", lwd=3, cex.lab=1.2,
cex= 2, cex.main=1.5, cex.axis=.75, cex.legend=1.5, X.Lab=NULL)
# And now the level 3 plot (observation in plants in Treatment by type groups)
plot(x=co2multi, y="multi", level=3, aorder=FALSE, lcol="blue", pcol= c("green", "pink"),
round.c=1, bcol="lavender", xlim=c(-.1, 1), legend=NULL, add.legend="right", lwd=3,
cex.lab =1.2, cex= 2, cex.main=1.5, cex.axis=.75, cex.legend=1.5, X.Lab=NULL)

# Targets for length of stay (LOS) #
# Our administrators ask how far are we from our goals, they ask about targets
# in increments of 5 points of probability or specific days. We answer both.
btarget1 <- Bayes(x=losmcmc, y="target", type="n", parameter=c("muOfY","sigmaOfY"),
newdata=TRUE, target=list(p=c(.35,.4,.45, .5, .55),  y=c(3,4))) # 'newdata' for plots
# Target graph using the standard option, overlaying estimate of mode parameter
plot(x=btarget1, y="target", type="n", lcol="purple", tgtcol="blue", xlim=c(3.5, 5))
# Target graph using a posterior predictive check, more intuitive
plot(x=btarget1, y="target", type="n", data=hosprog, dv="los", breaks=30,
cex.axis=1.3, lwd=3, pline=20, vlim=c(-1, 12), xlim=c(-1, 10),
parameter=c("muOfY","sigmaOfY"), add.legend="topright", main="Length of Stay",
cex.main=1.5, xpt=5, pcol="black", lcol="cyan", tgtcol="blue", bcol="orange",
cex.legend=1.5, cex.text = 2)

Prediction plot of treatment and control groups for DID and ITS models

Description

Provides partial prediction plots for treatment and control groups from difference-in-difference (DID) and interrupted time series (ITS) models. The graph will produce lines for treatment/intervention and control groups to gain understanding through a visual representation of the regression coefficients. By default, the treatment/intervention group is represented with a blue line, the control group is represented with a red line, and the counterfactual line, when available, is a dashed line. There are many options to change the plot.

Usage

## S3 method for class 'assess'
plot(
  x,
  y,
  xlim = NULL,
  ylim = NULL,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  lwd = NULL,
  col = NULL,
  tcol = NULL,
  cfact = FALSE,
  conf.int = FALSE,
  adj.alpha = NULL,
  add.means = FALSE,
  add.legend = NULL,
  legend = NULL,
  tgt = NULL,
  tgtcol = "gray",
  cex = 1,
  cex.axis = NULL,
  cex.lab = NULL,
  cex.main = NULL,
  cex.text = NULL,
  cex.legend = NULL,
  name = FALSE,
  coefs = FALSE,
  round.c = NULL,
  pos.text = NULL,
  arrow = FALSE,
  xshift = NULL,
  x.axis = NULL,
  ...
)

Arguments

Value

plot of partial predictions for treatment and control groups.

Examples

am2 <- assess(formula= los ~ ., data=hosprog, intervention = "program",
topcode =NULL, int.time="month", regression="none", treatment= 5,
interrupt=c(5,9), did="two", its="two", newdata=TRUE, propensity=NULL)
plot(am2, "DID", add.legend="bottomleft", ylim=c(2, 8))  #DID model, basic plot
plot(am2, "ITS", add.legend="top", ylim=c(2, 8))         #ITS model, basic plot
plot(am2, "DID", add.legend="topleft", main="DID study", col=c("dodgerblue","magenta"),
ylim=c(2, 8), lwd=6, cex=3, cex.axis=2, cex.lab=1.5, cex.main=3, cex.text=2,
arrow=TRUE, xshift=0.02, coefs=TRUE, round.c=2 )
plot(am2, "ITS", add.legend="top", xlim=c(-.5, 13), ylim=c(2, 8), main="ITS study",
col=c("cyan","hotpink"), tcol="springgreen", lwd=7, cex=2, cex.axis=2, cex.lab=2,
cex.main=3, cex.text=1.2, cex.legend=1.25, name=FALSE, coefs=TRUE, round.c=1,
pos.text= list("txp5"=3, "post9"=4), arrow=TRUE, xshift=c(.5, 1.5),
cfact=TRUE, conf.int=TRUE, adj.alpha=0.2)
# ITS: USA's unemployment rate with a single group, 10 interruption model
key_time <- c(5, 14, 17, 29, 42, 59, 69, 73, 80,92)# key interruption periods
u110 <- assess(formula=rate ~ ., data=unemployment, intervention="usa",# model
               int.time="year", its="one", interrupt= key_time)
# Graph the results, note shift in intercept and slope during 1933's New Deal
plot(u110, "ITS", add.means = TRUE, coefs=TRUE, conf.int=TRUE, adj.alpha = .2,
     lwd=2, col="slategray", tcol= "orange", main="US unemployment rate",
     xlab="Years (1929-2024)", ylab="Proportion of labor market",
     cex.main=2, cex.axis = 1.5, cex.lab = 1.5, cex=2, cex.text = 1.25,
     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=1.25, labels =
         paste0(unemployment[ key_time[i], "Year"], ": ", unemployment[ key_time[i], "event"]))
}

Shewhart control charts

Description

Graph X-bar charts, p-charts, and u-charts. This includes producing means center lines, 3-sigma upper and lower control limits. Users can also calculate values before and after an intervention to see if a change in the control process happened. Values are returned in a data frame.

Usage

## S3 method for class 'control'
plot(
  x,
  y = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  lwd = NULL,
  col = NULL,
  iname = NULL,
  icol = "black",
  trend = FALSE,
  trcol = "gray",
  tgt = NULL,
  tgtcol = "gray",
  tpline = NULL,
  tpcol = NULL,
  cex = 1,
  cex.lab = NULL,
  cex.axis = NULL,
  cex.main = NULL,
  cex.text = NULL,
  x.axis = NULL,
  y.axis = NULL,
  round.c = NULL,
  ...
)

Arguments

Value

plot of Shewhart control charts: X-bar charts, p-charts, and u-charts with 3-sigma control limits.

References

Ryan, T. (2011). Statistical Methods for Quality Improvement, Third Edition. John Wiley & Sons, Inc. ISBN: 978-0-470-59074-4.

Examples

## Hospital LOS and readmissions ##
# X-bar chart
spc_x <- control(x="los", time="month", data=hosprog, type="x", n.equal=TRUE)
# Basic X-bar chart
plot(spc_x)

# p-chart, using only the numerator (i.e., y=NULL). Specify unequal sample sizes
spc_p <- control(x="rdm30", time="month", data=hosprog, type="p", n.equal=FALSE)
# p-chart, adding target and time point lines
plot(spc_p, tgt=c(0,.25), tgtcol="green", ylim=c(0,0.4), tpline=c(4,8),
tpcol= c("yellow","black"))

# u-chart for infection rates with an intervention
spc_u <- control(x="HAI", y="PatientDays", time="Month", data=infections,
type="u", n.equal=FALSE, intervention=22)
# u-chart with trend lines, various graphing options, x.axis start at 2nd year
# and y.axis changed to show HAIs per 1,000 patient days
plot(spc_u, main="u-Chart: HAI per 1,000 Patient Days Pre/Post Intervention",
col=c("green","dodgerblue"), trend=TRUE, trcol="red", x.axis=c((1:41+12)), round.c=1,
y.axis=seq(min(spc_u$HAI)*1000, max(spc_u$HAI)*1000, length.out=nrow(spc_u)),
xlab="Months (starting at year 2)", icol="gray", lwd=2, cex=2,
cex.axis=1.1, cex.main=1.25, cex.text=1.25)

Decision Curve Analysis plots and regression model classification graphs on sensitivity and specificity

Description

Graph decide class object's model classification results as well as Decision Curve analysis output for net benefit and interventions saved.

Usage

## S3 method for class 'decide'
plot(
  x,
  y = NULL,
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  xlim = NULL,
  ylim = NULL,
  lwd = NULL,
  bcol = NULL,
  lcol = NULL,
  add.legend = NULL,
  legend = NULL,
  cex = 1,
  cex.lab = NULL,
  cex.axis = NULL,
  cex.main = NULL,
  cex.legend = NULL,
  round.c = 2,
  ...
)

Arguments

Value

plots of model classification, net benefit, and interventions saved.

References

Vickers, A. & Elkin, E. (2006). Decision Curve Analysis: A Novel Method for Evaluating Prediction Models. Society for Medical Decision Making, 26, 6, 565-574. https://doi.org/10.1177/0272989X06295361

Examples

## Predicting car engine shape type, v or straight  ##
# run the model
car_m1 <- assess(formula=vs ~ hp + am, data=mtcars, regression="logistic")
# create a decide object, enter the model name and a threshold on the logit scale
d1 <- decide(x=car_m1, threshold= -0.767)
#Plot the classification results
plot(x=d1, y= "cl")

#Plot the net benefit
plot(x=d1, y= "nb")

#Plot the interventions saved
plot(x=d1, y= "is")

Confidence interval graphs for group class objects

Description

Confidence interval graphs for group class objects

Usage

## S3 method for class 'group'
plot(
  x,
  y = "group",
  order = "alpha",
  gcol = "blue",
  gband = FALSE,
  pcol = "red",
  overall = FALSE,
  ocol = "gray",
  oband = FALSE,
  tgt = NULL,
  tcol = "gray",
  tpline = NULL,
  tpcol = "gray",
  xlim = NULL,
  ylim = NULL,
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  lwd = 1,
  adj.alpha = 0.4,
  cex = 1,
  cex.axis = 1,
  cex.lab = 1,
  cex.main = 1,
  cex.text = 1,
  round.c = 2,
  name = FALSE,
  abbrv = 5,
  roll.x = "Stop",
  ...
)

Arguments

Value

plot of group level confidence intervals, including estimates over time periods.

Examples

#Simple graph for confidence intervals using the t-distribution
gr1 <- group(x="program", y="los", z="month", data=hosprog, dist="t",
increment=3, rolling=6)
# Group level confidence intervals
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=1, cex.axis=1, cex.lab=1, cex.text=2,
cex.main=1.25, name=TRUE, adj.alpha=.2)
#Trend plots over time in the 3 month increments (i.e., quarters)
plot(x=gr1, y="time", lwd=4, gcol=c("red", "blue"), gband=TRUE, overall=TRUE,
  oband=TRUE, ocol="gray", tcol="green", tgt=4, tpline=3,
  tpcol="yellow", name=TRUE, cex.axis=1, cex.lab=1, cex.text=2,
  cex.main=1.25, adj.alpha=.3)
#Plot for rolling 6-month averages
plot(x=gr1, y="roll", lwd=4, gcol=c("red", "blue"), gband=TRUE, overall=TRUE,
  oband=TRUE, ocol="gray", tcol="green", tgt=4, tpline=c(4,6),
  tpcol="yellow", name=TRUE, cex.axis=1, cex.lab=1, cex.text=2,
  cex.main=1.25, adj.alpha=.3)

  # Helpful note: You can plot aggregated data if you used the argument asis=TRUE in group(),

Plot of variable importance ranked by partial chi-square statistic

Description

Plots an importance class object. Produces a dot chart that places the predictor variable with the highest partial chi-square (Wald chi-square for individual coefficients) at the bottom. It is a metric of the partial chi-square minus its degrees of freedom (Harrell, 2015). Predictor variables with significant p-values at the 0.05 alpha are highlighted red. Consider graphical parameters of mar=c(4.2, 2, 3.5, 3) and oma = c(0, 0, 0, 3).

Usage

## S3 method for class 'importance'
plot(
  x,
  y,
  main = NULL,
  cex = NULL,
  pt.cex = NULL,
  pch = NULL,
  color = NULL,
  lcolor = NULL,
  ...
)

Arguments

Value

plot of variable importance, significant variables highlighted in red.

References

Harrell, F. E., Jr. (2016). Regression Modeling Strategies. Springer International Publishing. ISBN: 978-3-319-19424-0.

Examples

# OLS regression
plot(importance(assess(mpg ~ hp + wt + cyl, data=mtcars, regression= "ols")$model))

# logistic regression
plot(importance(assess(vs~mpg+wt+hp, data=mtcars, regression= "logistic")$model))

Print alpha results

Description

Formats alpha class results to display summary statistics of scale information. These include the overall alpha, scale mean and standard deviation, item statistics, scale statistics if an item is removed from the scale, and the total sample size.

Usage

## S3 method for class 'alpha'
print(x, ...)

Arguments

Value

formatted alpha results.

Examples

print(alpha(items=c("i1","i2","i3","i4","i5"), data=cas))

Print interpret object

Description

Formats interpretations from interpret class objects. Provides simple interpretations of regression coefficients and Cronbach's alpha. Print specific model interpretations (or run all), returned in sentence and paragraph formats.

Usage

## S3 method for class 'interpret'
print(x, ...)

Arguments

Value

formatted interpret object results.

Examples

#Cronbach's alpha
print(interpret(alpha(items=c("i1","i2","i3","i4","i5"), data=cas)))

#' # interpret a standard linear (OLS) regression
hos1 <- assess(formula=survey ~ program + month, data=hosprog, regression= "ols")
print(interpret(hos1)$model)

# interpret a differences-in-differences model
hos2 <- assess(formula=survey ~ ., data=hosprog, intervention = "program",
int.time="month", treatment = 5, did="two", newdata=TRUE)
interpret(hos2)$did

# interpret an interrupted time series model
hos3 <- assess(formula=survey ~ ., data=hosprog, intervention = "program",
int.time="month", its="two", interrupt = 5)
interpret(hos3)$its

USA's unemployment rate between 1929 and 2024

Description

USA's unemployment rate between 1929 and 2024

Usage

unemployment

Format

unemployment

An artificial data frame with 96 rows and 5 columns:

Year

Calendar year as an integer from 1929 to 2024.

rate

Unemployment rate defined as the proportion of the US Labor Market that were unemployed.

event

Notable events in the United States that may have impacted the unemployment rate. Other than notable events (or unspecified) listed as 'other'.

year

Integer for the number of years for 1929 to 2024, ranging from 1 to 96.

usa

Value of 1 to indicate data for the USA.

...

Source

unemployment is an artificial data frame of reasonable estimates made entirely for the purpose of demonstrating interrupted time series with many interruptions. For precise rates, please see a reliable source.