# ham-package

    library(ham)

# Healthcare Analysis Methods (ham)

The goal of ham is to provide different modeling approaches to evaluating healthcare programs (or interventions) with descriptive statistics, point estimates and confidence intervals, and regression analyses. This includes standard regression methods like linear (OLS) and logistic regression. And ham adds options for differences-in-differences models as well as interrupted time series models. DID and ITS models offer options for causal modeling. What is unique about ham is that it creates datasets with constructed variables for DID and ITS models, optionally it can add top coded outcome variables, propensity scores, and provides some interpretation of model results. Additionally, Cronbach’s alpha can be calculated for such things as patient surveys. As the logo of Dr. Ham suggests, the ham package can help illuminate your results.

## 1. Introduction

Did our hospital intervention make a difference? Does our work improve healthcare? Did our efforts impact health outcomes? Should we spread our efforts to our other healthcare facilities?

These are questions we are constantly asked each year in healthcare. Another important question we should ask is how do we answer these questions? This vignette shows how we can assess program effects by building upon a set of analyses and the associated ham functions.

This is the ham package for healthcare analysis methods. This package can help when performing program evaluations or intervention studies in healthcare. Or simply to test if a program has an impact on an outcome of interest.

ham can help with research or evaluation studies. When working in healthcare systems, multi-site evaluations can strongly resemble research studies and commonly use regression methods. For an introduction to evaluation, please see Michael Patton’s book (1997).

What is unique about ham, is that it provides options for running standard linear or ordinary least squares (OLS) and logistic regression as well as methods used in causal modeling such as differences-in-differences (DID) and interrupted time series analysis (ITS). It also optionally makes data with the newly created variables (i.e., this saves you time).

This vignette will introduce ham’s features in the following functions:

- alpha: Conduct Cronbach’s alpha on scale items (e.g., survey questions).

- assess: Perform various regression methods (OLS, logistic, differences-in-differences, and interrupted time series).

- group: Confidence Intervals for descriptive and inferential statistics. Compare groups in cross-sectional or repeated cross-sectional data.

- importance: Rank variable importance from regression coefficients using the partial chi-square statistic.

- interpret: Provides simple coefficient interpretations. This is a helpful reminder, especially as models have increased coefficients (e.g., ITS).

- There are also printing and plotting options to help review your results.

Below will cover 4 sections using an intervention program example and reduced length of stay in the hospital.

## 2. Descriptive and inferential statistics

This example shows group level point estimates and confidence intervals.

### Group level 95% confidence intervals over 12 months (t distribution)

    gr1 <- group(x="program", y="los", z="month", data=hosprog, dist="t", cluster=TRUE)
    print(gr1$Group.CI)
    #> $adf_alpha
    #>   Group PointEst    Lower    Upper
    #> 1     1 4.247012 4.034842 4.459181
    #> 0     0 4.585967 4.391132 4.780802
    #>
    #> $adf_numeric
    #>   Group PointEst    Lower    Upper
    #> 0     0 4.585967 4.391132 4.780802
    #> 1     1 4.247012 4.034842 4.459181
    #>
    #> $adf_all
    #>   PointEst    Lower    Upper
    #> 1 4.428259 4.284543 4.571974

For many health outcomes, we may track data with high variation across months. This code that uses the ‘increment’ and ‘rolling’ arguments can help make better sense in those situations with: 1) 3-month increments, and 2) 6-month rolling averages.

And we may have a continuous variable that we would want to explore by forming it into quartiles. You can turn a continuous variable into quartiles using the quarts=TRUE argument. Below, we get 6-month rolling average LOS confidence intervals by quartiles of age.

## Rolling 6-month averages by quartiles of age (Q1 - Q4)

    gr2 <- group(x="age", y="los", z="month", data=hosprog, dist="t", increment=3, rolling=6, quarts=TRUE)
    print(gr2$Roll.CI$Rolling)
    #>     Group PointEst    Lower    Upper Rolling Start Stop
    #> Q4     Q4 3.970417 3.611306 4.329528       6     1    6
    #> Q3     Q3 4.444200 3.993038 4.895363       6     1    6
    #> Q2     Q2 4.677006 4.239698 5.114313       6     1    6
    #> Q1     Q1 4.625708 4.198619 5.052797       6     1    6
    #> Q41    Q4 3.973603 3.625839 4.321367       6     2    7
    #> Q31    Q3 4.481611 4.024713 4.938510       6     2    7
    #> Q21    Q2 4.501353 4.114604 4.888103       6     2    7
    #> Q11    Q1 4.554040 4.167738 4.940343       6     2    7
    #> Q42    Q4 3.936571 3.583429 4.289712       6     3    8
    #> Q32    Q3 4.377289 3.923197 4.831380       6     3    8
    #> Q22    Q2 4.590305 4.201249 4.979361       6     3    8
    #> Q12    Q1 4.679806 4.297483 5.062129       6     3    8
    #> Q43    Q4 3.918199 3.594263 4.242136       6     4    9
    #> Q33    Q3 4.212241 3.786031 4.638452       6     4    9
    #> Q23    Q2 4.579703 4.228804 4.930601       6     4    9
    #> Q13    Q1 4.522246 4.184928 4.859565       6     4    9
    #> Q44    Q4 3.849373 3.550891 4.147855       6     5   10
    #> Q34    Q3 4.051160 3.636456 4.465863       6     5   10
    #> Q24    Q2 4.601965 4.260039 4.943891       6     5   10
    #> Q14    Q1 4.618284 4.265461 4.971107       6     5   10
    #> Q45    Q4 4.047580 3.747916 4.347243       6     6   11
    #> Q35    Q3 4.079230 3.679460 4.478999       6     6   11
    #> Q25    Q2 4.547797 4.220430 4.875164       6     6   11
    #> Q15    Q1 4.617848 4.262364 4.973333       6     6   11
    #> Q46    Q4 3.985610 3.687051 4.284169       6     7   12
    #> Q36    Q3 4.262495 3.819285 4.705705       6     7   12
    #> Q26    Q2 4.646603 4.230228 5.062979       6     7   12
    #> Q16    Q1 4.786002 4.385044 5.186960       6     7   12

There is also an option to retrieve those estimates over time periods. And there are graphing options to help see how the data looks for both.

Results can be returned for each unit of time, in various increments such as each 3 months (i.e., quarters) or for rolling averages such as a rolling 12-month period. Using increments of multiple months or as rolling averages can help visualize results when there is high variation over time.

Hint: The rolling average can add time to to the function, for example using N= 760,000 rows over 11 years took 35 seconds on my moderately fast computer built 5 years ago. In contrast, without rolling averages, it took 5 seconds.

### Graph the results

There are 3 main graphing options that go into the y argument:

plot(x=your_group_object, y=type_of_plot)

- y= “group” is the group confidence intervals. This is the default graph.
- y= “time” is the standard confidence intervals over time including the increment of multiple time periods (e.g., 3-month quarters).
- y= “roll” is the ‘rolling’ average (e.g., rolling 12-month average).

    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)

It appears there is a significant difference between the groups and the intervention group is at a significantly lower level than the target set by hospital executives (vertical green line, tgt=4.5).

The results show up from “low-to-high” LOS values because the code uses order=‘numeric’. It would be in the opposite order when using order=‘alpha’ which indicates alphabetical ordering based on group names. The gray overall confidence band is included with: overall=TRUE, oband=TRUE, ocol=“gray”.

Hint: Need something fast? Just use 1 argument to get the simplified graph:

plot(gr1)

## Between-group variance and within-group clustering

We can assess the variation between-groups with statistics such as the Intra-class Correlation Coefficient. This reflects the proportion of total variance that is between-group variance and also the within-group clustering. We can have this and similar values such as the Median Odds Ratio (binary outcomes only) and the Design Effect. These will give us an idea of whether we should consider the variation between groups more in our study. For example, Design Effects >= 2 are considered important. Please see the documentation for more information. This may be more helpful with many groups (e.g., 10-20 hospitals in a healthcare system).

    print(gr1$Clustering)
    #> Intraclass.Correlation.Coefficient                  Median.Odds.Ratio
    #>                         0.01198441                                 NA
    #>                     Dessign.Effect
    #>                         5.30240495

### Graph the confidence interval trend results

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

The intervention group’s data makes up the blue line while the control is indicated with the red line. This looks interesting, the intervention and control group swap LOS ranks at the same time as the start of the intervention. This graph is made for descriptive purposes but may shed light on trends for intervention and control groups. It looks like something happened, we should continue exploring.

We see group names when name=TRUE. 95% confidence bands are provided using gband=TRUE. We added the overall gray line and confidence band using overall=TRUE, oband=TRUE, and ocol=“gray”. adj.alpha=0.3 provides the level of shading, higher values are darker. The 4.5 LOS target and time point lines (5) and colors are provided with tgt=4.5, tcol=“green”, tpline=5, tpcol=“yellow”.

Hint: Need something fast? Get a quick graph using 2 arguments:

plot(gr1, “time”)

## 3. Linear and logistic regression

Our regression models are performed with assess(). The example dataset has common variables found in program evaluation or intervention studies (I’ll refer to both as a ‘study’), there are various outcome and predictor variables (or response and explanatory variables or dependent and independent variables or other names common in your field). In these studies, we try to assess the impact of the predictors on the outcome. We often use treatment and control groups to asses a healthcare program or intervention’s impact on our key outcome variable of interest. Because of multiple stakeholders, we often conduct these studies with multiple outcomes to help answer the multiple stakeholder’s questions.

A common approach to answering study questions is using a regression to test a treatment effect while controlling for other covariates or other important variable that may explain the variation in the regression outcome.

Here are OLS and logistic regression examples using assess() on the mtcars data. This ogistic regression model code not run: assess(formula=vs~mpg + wt + hp, data=mtcars, regression=“logistic”)$model

### OLS or linear regression

    # Model summary
    summary(assess(hp ~ mpg+wt, data=mtcars, regression="ols")$model)
    #>
    #> Call:
    #> stats::lm(formula = primary_formula, data = combined_df, weights = wght_obj,
    #>     offset = offst_obj)
    #>
    #> Residuals:
    #>    Min     1Q Median     3Q    Max
    #> -59.42 -30.75 -12.07  24.82 141.84
    #>
    #> Coefficients:
    #>             Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)  349.287    103.509   3.374  0.00212 **
    #> mpg           -9.417      2.676  -3.519  0.00145 **
    #> wt            -4.168     16.485  -0.253  0.80217
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 44.65 on 29 degrees of freedom
    #> Multiple R-squared:  0.6033, Adjusted R-squared:  0.576
    #> F-statistic: 22.05 on 2 and 29 DF,  p-value: 1.505e-06

### ham provides coefficient interpretations with interpret()

    # Model summary
    interpret(assess(hp ~ mpg+wt, data=mtcars, regression="ols"))$model
    #> Interpretations: Regression model
    #> ---------------------------------
    #> These estimates tell you about the relationship between the
    #> independent variables and the dependent variable. These estimates
    #> tell the amount of change in outcome scores that would be
    #> predicted by a 1 unit increase in the predictor.
    #>
    #> The following predictor variable(s) have coefficient(s)
    #> significantly different from 0 using an alpha of 0.05:
    #> mpg
    #>
    #> For every 1 unit increase in these predictor variables,
    #> hp is predicted to increase by the value of the
    #> coefficient, holding all other variables constant. The following
    #> predictor variable(s) have positive coefficient(s) that
    #> increase the predicted value of the outcome:
    #> No positive coefficients in your model were significant.
    #>
    #> For every 1 unit increase in these predictor variables,
    #> hp is predicted to decrease by the value of the
    #> coefficient, holding all other variables constant. The following
    #> predictor variable(s) have negative coefficient(s) that
    #> decrease the predicted value of the outcome:
    #> mpg
    #>
    #> R-Squared (R2) is the proportion of variance in the dependent
    #> variable which can be predicted from the independent
    #> variable(s). For example, if R2 = 0.50, 50% of the variance
    #> in test scores can be predicted from the 5 variables. R2 >= 0.80
    #> may be at a level to reliably make individual predictions.
    #> Lower R2 may be helpful in group level predictions. And low R2 can
    #> still be adequate for hypothesis testing. This model has a R2 of 0.603.

In this example, we predict the binary outcome of Automatic or Manual car transmission.

### Logistic regression

    # Model summary
    summary(assess(am ~ hp+wt, data=mtcars, regression="logistic")$model)
    #>
    #> Call:
    #> stats::glm(formula = primary_formula, family = binomial(link = family_link),
    #>     data = combined_df, weights = wght_obj, offset = offst_obj)
    #>
    #> Coefficients:
    #>             Estimate Std. Error z value Pr(>|z|)
    #> (Intercept) 18.86630    7.44356   2.535  0.01126 *
    #> hp           0.03626    0.01773   2.044  0.04091 *
    #> wt          -8.08348    3.06868  -2.634  0.00843 **
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> (Dispersion parameter for binomial family taken to be 1)
    #>
    #>     Null deviance: 43.230  on 31  degrees of freedom
    #> Residual deviance: 10.059  on 29  degrees of freedom
    #> AIC: 16.059
    #>
    #> Number of Fisher Scoring iterations: 8

### Logistic regression interpretations

    # Specify how many significant digits with digits argument
    interpret(assess(am ~ hp+wt, data=mtcars, regression="logistic"))$model
    #> Interpretations: Regression model
    #> ---------------------------------
    #> These estimates tell you about the relationship between the
    #> independent variables and the dependent variable. These estimates
    #> tell the amount of change in outcome scores that would be
    #> predicted by a 1 unit increase in the predictor.
    #>
    #> The following predictor variable(s) have coefficient(s)
    #> significantly different from 0 using an alpha of 0.05:
    #> hp, wt
    #>
    #> For every 1 unit increase in these predictor variables,
    #> am is predicted to increase by the value of the
    #> coefficient, holding all other variables constant. The following
    #> predictor variable(s) have positive coefficient(s) that
    #> increase the predicted value and odds of the outcome:
    #> hp (3.69% increased odds or odds are 1.04 times higher)
    #>
    #> For every 1 unit increase in these predictor variables,
    #> am is predicted to decrease by the value of the
    #> coefficient, holding all other variables constant. The following
    #> predictor variable(s) have negative coefficient(s) that
    #> decrease the predicted value and odds of the outcome:
    #> wt (100% decreased odds)
    #>
    #> There is no R2 or C-statistic (AUC) information provided.

### Poisson regression

In this example, we examine the Hospital Acquired Infection rate over 3+ years. Notice the offset(log(PatientDays)) part of the formula, offsets are common in health care analyses and allow to take into account varying exposures like patient days or can be used for observed/expected ratios, e.g., offset(log(Expected)). Another option is to log transform your exposure variable as a variable in your data frame and enter it into the offset argument. Please note that in assess, the offset argument uses a character name, thus offset= ‘logPatientDays’.

    # Model summary
    summary(assess(HAI ~ Month + offset(log(PatientDays)),
                    data = infections, regression="poisson", link="log")$model)
    #>
    #> Call:
    #> stats::glm(formula = primary_formula, family = poisson(link = family_link),
    #>     data = combined_df, weights = wght_obj, offset = offst_obj)
    #>
    #> Coefficients:
    #>              Estimate Std. Error z value Pr(>|z|)
    #> (Intercept) -6.388783   0.032283 -197.90   <2e-16 ***
    #> Month       -0.032955   0.001558  -21.16   <2e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> (Dispersion parameter for poisson family taken to be 1)
    #>
    #>     Null deviance: 838.37  on 40  degrees of freedom
    #> Residual deviance: 374.75  on 39  degrees of freedom
    #> AIC: 628.36
    #>
    #> Number of Fisher Scoring iterations: 4

We see that there is a significant Incidence Rate by month as seen in the interpretations below.

### Poisson regression interpretations

    # Specify how many significant digits with digits argument
      interpret(assess(HAI ~ Month + offset(log(PatientDays)),
                    data = infections, regression="poisson", link="log"))$model
    #> Interpretations: Regression model
    #> ---------------------------------
    #> These estimates tell you about the relationship between the
    #> independent variables and the dependent variable. These estimates
    #> tell the amount of change in outcome scores that would be
    #> predicted by a 1 unit increase in the predictor.
    #>
    #> The following predictor variable(s) have coefficient(s)
    #> significantly different from 0 using an alpha of 0.05:
    #> Month
    #>
    #> For every 1 unit increase in these predictor variables,
    #> HAI is predicted to increase by the value of the
    #> coefficient, holding all other variables constant. The following
    #> predictor variable(s) have positive coefficient(s) that
    #> increase the predicted value and incidence rate of the outcome:
    #> No positive coefficients in your model were significant.
    #>
    #> For every 1 unit increase in these predictor variables,
    #> HAI is predicted to decrease by the value of the
    #> coefficient, holding all other variables constant. The following
    #> predictor variable(s) have negative coefficient(s) that
    #> decrease the predicted value and incidence rate of the outcome:
    #> Month (3.24% decreased incidence rate)
    #>
    #> There is no R2 or C-statistic (AUC) information provided.

### OLS top-coding and propensity scores

ham can topcode the outcome and create a propensity score variable. top-coding cost at $17,150 and propensity score based on age, female indicator, and a health risk probability score.

Top-coding provides more accurate results because it lessens the impact of outliers on regression results. And propensity scores offer more balance in data by controlling for factors that are important in understanding how patients are more or less likely to be in the intervention group.

Topcoding can be applied to any model. Propensity scores can be created for any model except the single group interrupted time series because there is no control group (i.e. intervention group only).

We included an interaction so that we can examine both the intervention and control group’s trends.

    m1 <- assess(formula=cost ~ month * program, data=hosprog, intervention = "program",
    regression="ols", topcode=17150, propensity=c("female","age","risk"),
    newdata=TRUE)

### model results

    summary(m1$model)
    #>
    #> Call:
    #> stats::lm(formula = primary_formula, data = combined_df, weights = wght_obj,
    #>     offset = offst_obj)
    #>
    #> Residuals:
    #>     Min      1Q  Median      3Q     Max
    #> -6570.9 -2158.3  -598.1  1962.0 10521.1
    #>
    #> Coefficients:
    #>                Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)    14006.64    1363.65  10.271  < 2e-16 ***
    #> month            424.81      46.57   9.122  < 2e-16 ***
    #> program         5524.74     514.18  10.745  < 2e-16 ***
    #> pscore        -15989.00    2838.57  -5.633 2.55e-08 ***
    #> month:program   -908.68      68.64 -13.239  < 2e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 3165 on 715 degrees of freedom
    #> Multiple R-squared:  0.233,  Adjusted R-squared:  0.2287
    #> F-statistic:  54.3 on 4 and 715 DF,  p-value: < 2.2e-16

Descriptive statistics on newly created variables and the original cost as a comparison. top.cost is the topcoded cost variable. R^2 is not bad.

    summary(m1$newdata[, c( "cost","top.cost", "pscore")])
    #>       cost          top.cost         pscore
    #>  Min.   : 1483   Min.   : 1483   Min.   :0.3708
    #>  1st Qu.: 6410   1st Qu.: 6410   1st Qu.:0.4364
    #>  Median : 8639   Median : 8639   Median :0.4655
    #>  Mean   : 9348   Mean   : 9215   Mean   :0.4653
    #>  3rd Qu.:11487   3rd Qu.:11487   3rd Qu.:0.4942
    #>  Max.   :27540   Max.   :17150   Max.   :0.5609

In the previous example, we used covariate adjustment with the propensity score but now we will use inverse probability of treatment weighting using the propensity score (IPW). There is a weights argument that we can include a character vector for 3 types of IPW weights: “ipw” (Horvitz-Thompson estimator), “nipw” (Normalized IPW by Hajek estimator), and “att”. Trimming of propensity scores is also available with ‘trim’ argument. See assess() documentation to read more.

### Inverse Probability of Treatment Weighting using the propensity score (IPW)

    ipw1 <- assess(formula=cost ~ month * program, data=hosprog, intervention = "program",
    regression="ols", topcode=17150, propensity=c("female","age","risk"),
    newdata=TRUE, weights="ipw")

### Weighted least squares model results

    summary(ipw1$model)
    #>
    #> Call:
    #> stats::lm(formula = primary_formula, data = combined_df, weights = wght_obj,
    #>     offset = offst_obj)
    #>
    #> Weighted Residuals:
    #>     Min      1Q  Median      3Q     Max
    #> -9998.5 -3023.7  -780.7  2639.6 16772.8
    #>
    #> Coefficients:
    #>               Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)    6636.33     376.41  17.631   <2e-16 ***
    #> month           423.80      49.29   8.598   <2e-16 ***
    #> program        5452.85     524.19  10.402   <2e-16 ***
    #> month:program  -917.16      69.98 -13.105   <2e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 4579 on 716 degrees of freedom
    #> Multiple R-squared:  0.2004, Adjusted R-squared:  0.1971
    #> F-statistic: 59.83 on 3 and 716 DF,  p-value: < 2.2e-16

### importance

We can examine variable importance based on partial chi-square (i.e., which variables explain the outcome the most).

    importance(m1$model)
    #>               X    Chi.Sq d.f.      p.value
    #> 1         month 175.30630    2 8.564878e-39
    #> 2       program 179.92844    2 8.492485e-40
    #> 3        pscore  31.72813    1 1.773352e-08
    #> 4 month:program 175.27159    1 5.222734e-40

### Graph importance

We can examine variable importance to see a ranking of variables with a graph. The hospital program has the highest rank, variables highlighted in red are statistically significant. In other words, we get the biggest bang for the buck with our intervention program, it is associated with reduced costs. As saving resources is part of healthcare, this graph can help identify the most important factors to focus on (or identify which things are less important).

    #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))

### Graph a summary of the coefficients

This graph will be added in version 1.3 and it is a way to get a review of the coefficients, showing the direction of the effect with point estimates and 95% confidence intervals. This was inspired by Frank Harrell’s plot.summary in the rms package. Instead ham uses plot.review with the almost identical summary but adds some class information to help with plotting. This will work on assess regression objects for OLS linear, logistic, Poisson and models coming from Base R’s lm() and glm() as well as the survival package and possibly other glm() models. It can work with Cox Proportional Hazards from the survival package. And this will work with the differences-in-differences and interrupted time series models created with assess().

OLS regression–most basic format with Base R lm()

    m01 <- lm(mpg ~ wt+hp+am, data=mtcars)
    # Most basic code
    plot(review(m01))

Take a look at the effect when hp is increased by 83.5

    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)
    #> =================================================
    #>                   Model Review
    #> =================================================
    #>    PointEst    Lower   Upper         P
    #> wt -2.87900 -4.73200 -1.0250 0.0035740
    #> hp -0.03748 -0.05715 -0.0178 0.0005464
    #> am  2.08400 -0.73580  4.9030 0.1413000
    #> =================================================

let’s see the impactful 83.5 hp increase in going from the 1st to 3rd quartiles

    m02 <- assess(formula=mpg ~ wt+hp+am, data=mtcars, regression="ols")
    plot(x=review(m02$model, increase=c(hp= 83.5)))

Logistic model

    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)
    #>    PointEst    Lower     Upper        P
    #> hp 0.000084 0.000000  0.805363 0.044832
    #> am 0.004062 0.000010  1.697294 0.073752
    #> wt 0.136915 0.001783 10.513523 0.369319

Poisson model with an offset using assess()

    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)

## 4. Differences-in-Differences

Let us continue with our length of stay example. By now we have begun to compare results for the intervention and control groups. We’ve spoken with various stakeholders and primary intended users to learn more about the program and its potential effects. Staff have listed various outcomes that would benefit so the following analysis can be applied to those outcomes as well.

And analyzing multiple outcomes would give us a “good story” on program impact and increase the potential utilization of our evaluation’s results.

We have the option of doing the analysis with an OLS or multiple linear regression, controlling for other important variables that can explain the outcome. We could create an interaction between the variables that represent time and intervention/control group participation to test if there is a different trend for each group.

This would work in many situations but we may want to consider what is available in causal modeling, especially if our control and intervention groups begin at different levels. For example, the intervention group may be the hospital (or department, etc.) with the most need of an intervention. They may have outcome rates that are at a worse level than the other hospitals. In this case, we nay want to assess program impact as the difference in trends for the intervention and control groups.

Another reason we would want to consider a DID model is if we are considering spreading this program to other medical centers, we would then want to be highly certain of the program’s efficacy.

DID can be used on binary or continuous outcome variables. Below is an example using the hosprog data with length of stay as the outcome and the created DID variables as the predictors, use ‘.’ on the right-hand side of the formula to indicate only the created variables will be used. Replace ‘.’ with any additional selected variables. The newly created DID variables will be added in all DID models.

This model has a pre/post design (i.e., there are only 2 distinct time points) by selecting: did=“two”, with post starting at month 5. Select did=“many” for a repeated cross-sectional design that has many time points (i.e., more than the 2 time periods in a pre/post test).

### DID model for length of stay using assess()

    dm1 <- assess(formula= los ~ ., data=hosprog, intervention = "program",
                  int.time="month", treatment= 5, did="two")

### Code to view DID model results

    summary(dm1$DID)
    #>
    #> Call:
    #> stats::lm(formula = DID_formula, data = combined_df, weights = wght_obj,
    #>     offset = offst_obj)
    #>
    #> Residuals:
    #>     Min      1Q  Median      3Q     Max
    #> -3.6247 -1.2003 -0.3145  0.8564  8.8642
    #>
    #> Coefficients:
    #>             Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)   3.4940     0.1650  21.175  < 2e-16 ***
    #> Post.All      1.5629     0.1974   7.917 9.25e-15 ***
    #> Int.Var       2.0664     0.2349   8.797  < 2e-16 ***
    #> DID          -3.5448     0.2849 -12.444  < 2e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 1.777 on 716 degrees of freedom
    #> Multiple R-squared:  0.1848, Adjusted R-squared:  0.1814
    #> F-statistic: 54.11 on 3 and 716 DF,  p-value: < 2.2e-16

Hint: The primary variables of the intervention are the ones that begin with ‘DID’. In this model, DID represents the intervention effect and is the differences-in-differences. In a model with many or more than 2 time periods, “DID.Trend” is the primary variable of interest while the “DID” variable is less so because it represents the difference before-and-after the intervention for both intervention and control groups (i.e., intervention group after the intervention vs. control group + intervention before the intervention start) so it represents something other than when DID=“two”. To run with “many” time periods, use similar code to this:

- assess(formula= los ~ ., data=hosprog, intervention = “program”, int.time=“month”, treatment= 5, did=“many”)

### Code to get the coefficient interpretations

    interpret(dm1)$did
    #> Interpretations: DID
    #> --------------------
    #> The intercept represents the mean los value of the
    #> control group at the baseline period (Time 1): 3.49.
    #>
    #> Post.All is the change in the control group's los
    #> value in the 2nd time period (Time 2). There is a
    #> significant increase for the control group
    #> at time 2: 1.56.
    #>
    #> Int.Var is the difference between the intervention
    #> and control group at the baseline period (Time 1). The
    #> intervention group had a significant increase in the
    #> mean los value compared to the control group: 2.07.
    #>
    #> DID estimates the average treatment effect on the
    #> treated group (ATET). This interaction represents the
    #> difference in the trend differences for the intervention and
    #> control groups:
    #> (Int. Time 2 - Int. Time 1) - (Ctl. Time 2 - Ctl. Time 1) = -3.54.
    #> In other words, there is a significant decrease in the
    #> mean los trend by -3.54 for the intervention group.
    #>
    #> If there are additional variables in the model then the coefficients
    #> above represent the effects after controlling for the other variables.

Now lets graph the results with the counterfactual (what could have happened without the intervention program).

The ham options include arrows to highlight the coefficient as the difference between counterfactual and post-intervention value. There are options to present the coefficient names only, coefficient values next to the coefficient names with asterisks (*) to indicate statistical significance, or none of the above. And we can add confidence bands to the expected lines.

All of the above mentioned options are shown below, including lighter confidence bands, with the use of these arguments:

- x=dm1 indicates the model
- y=“DID” indicates the type of model
- arrow=TRUE
- xshift=c(.045) shifts the arrow .045 to the right
- coefs=TRUE add coefficient values, names, and significance
- round.c=2 rounds numbers to the 2nd digit
- cfact=T adds a counterfactual line
- conf.int=TRUE adds confidence interval bands
- adj.alpha=0.2 adjust the shading to a lower level

Hint: Need something fast? Just use a barebones graph to get a quick view (blue line= Intervention, red line= Controls) using:

plot(dm1, “DID”)

### View a partial prediction plot

    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 )

We see above that there is a significant DID result over 3.5 fewer days length of stay. That is quite a high number. We can continue with the analysis to explain that result more but lets first examine the graph to help understand the results.

When we observe the solid blue trend line from the model, we see how much it differs with the blue dashed line, the counterfactual as denoted with the right-shifted blue arrow (xshift= .045). The counterfactual line represents what could have happened and we see that we may have expected to see an increase in LOS like the control group. But we instead see a decrease in LOS for the intervention group.

We can also use DID on binary outcomes like hospital re-admission within 30-days.

    dm3 <- assess(formula= rdm30 ~ ., data=hosprog, intervention = "program",
                  int.time="month", treatment= 5, did="two")

### View DID 30-day readmissions results

    summary(dm3$DID)
    #>
    #> Call:
    #> stats::lm(formula = DID_formula, data = combined_df, weights = wght_obj,
    #>     offset = offst_obj)
    #>
    #> Residuals:
    #>      Min       1Q   Median       3Q      Max
    #> -0.31858 -0.25651 -0.07207 -0.06034  0.93966
    #>
    #> Coefficients:
    #>             Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)  0.06034    0.03422   1.763   0.0782 .
    #> Post.All     0.19616    0.04094   4.792 2.01e-06 ***
    #> Int.Var      0.25824    0.04871   5.301 1.53e-07 ***
    #> DID         -0.44267    0.05907  -7.493 1.98e-13 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 0.3686 on 716 degrees of freedom
    #> Multiple R-squared:  0.0759, Adjusted R-squared:  0.07203
    #> F-statistic:  19.6 on 3 and 716 DF,  p-value: 3.197e-12

Significant DID effect showing reduced re-admissions

## 5. Interrupted Time Series

ITS lets us look at trends for 1 or 2 groups such as an intervention/treatment group without a control group or both a treatment and control group. And we have the option of one or more treatment periods (or interruptions). This gives us 4 options that can be specified using the interrupt and its arguments.

Below are examples using the hosprog data for the patient length of stay and death within 30-days. The dataset hosp1 will be used for the single group examples.

We’ll continue with our length of stay example because our clinical staff have expressed interest in understanding how the trend varies over time. Some hope for an immediate shift in the intercept followed by a ‘continuous’ effect in the slope while others fear a ‘discontinuous’ effect.

In the spirit of “formative” evaluation, we model these questions with ITS. It may be possible that we have a single, linear trend line in which differences-in-differences would provide more than adequate results. However, we want to examine our primary intended user’s questions about changes in the trend with ITS.

We can begin by looking at a single group (intervention group) with a single interruption/treatment period and assessing their LOS scores. We specify it with: interrupt= 5 and its=“one”.

### ITS model 1

    im11 <- assess(formula=los ~ ., data=hosp1, intervention = "program",
                   int.time="month", interrupt= 5, its="one")

### View ITS model 1 results

    summary(im11$ITS)
    #>
    #> Call:
    #> stats::lm(formula = ITS_formula, data = combined_df, weights = wght_obj,
    #>     offset = offst_obj)
    #>
    #> Residuals:
    #>     Min      1Q  Median      3Q     Max
    #> -3.7673 -1.1321 -0.3755  0.5758  9.1718
    #>
    #> Coefficients:
    #>             Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)   6.3299     0.4170  15.178   <2e-16 ***
    #> ITS.Time     -0.2889     0.1442  -2.003   0.0460 *
    #> post5        -1.0026     0.4264  -2.351   0.0193 *
    #> txp5          0.2014     0.1522   1.324   0.1865
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 1.726 on 331 degrees of freedom
    #> Multiple R-squared:  0.2426, Adjusted R-squared:  0.2357
    #> F-statistic: 35.34 on 3 and 331 DF,  p-value: < 2.2e-16

Hint: The primary variables of interest when examining ITS models with only the intervention group are the ones that begin with ‘post’ and ‘txp’ as they represent changes in the intercept and slope of the trend line due to the intervention or program. Please note that the main variable names will change in the models that include a control group.

### interpret coefficients

    interpret(im11)$its
    #> Interpretations: ITS
    #> --------------------
    #> Note: Some variable names below based on time points (or 'interruptions').
    #> This analysis is for a one-group, single intervention period (interruption).
    #>
    #> Intercept is 6.33 and the starting value of the trend
    #> for the intervention group.
    #>
    #> ITS.Time is -0.289 and the slope prior to intervention.
    #> The coefficient is significant.
    #>
    #> post5 is -1 and the immediate shift in the trend line
    #> after the intervention start (e.g., 1st year of intervention).
    #> The coefficient is significant.
    #>
    #> txp5 is 0.201 and the difference between pre- and
    #> post-intervention slopes (e.g., change in the pre-intervention
    #> slope). The coefficient is non-significant.
    #>
    #> Summary: The results show that after the start of the intervention,
    #> there is a non-significant change in the los trend. This gives
    #> a post-intervention trend change per time unit in the los of -0.0874
    #> (i.e., the value of change per-unit-of-time, such as month or year, in
    #> the intervention period; not the change relative to pre-intervention).
    #>
    #> If there are additional variables in the model then the coefficients
    #> above represent effects after controlling for the other variables.

There is a second key period of interest at month 9 which is specified with interrupt= c(5, 9) and its=“one”. Changing its=“one” to its=“two” would allow me to compare two groups (i.e., intervention and control groups). We’ll examine those options later but for now here is how we specify a model for the intervention group only with 1 interruption (or intervention) period.

Hint: In assess(), DID=“two” and ITS=“two” mean two different things. For DID models, DID=“two” refers to 2 time periods while ITS=“two” refers to 2 groups.

### Graph ITS model 1

    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 )

You’ve probably noticed only 1 trend line because we only have 1 group, the intervention group. In the 1 group models, our counterfactual lines appear to continue from the end of the baseline. However, we could continue the line from the the second period when there are 2 interruptions.

### ITS model 2

    im12 <- assess(formula=los ~ ., data=hosp1, intervention = "program",
                   int.time="month", interrupt= c(5, 9), its="one")

### View ITS model 2 results

    summary(im12$ITS)
    #>
    #> Call:
    #> stats::lm(formula = ITS_formula, data = combined_df, weights = wght_obj,
    #>     offset = offst_obj)
    #>
    #> Residuals:
    #>     Min      1Q  Median      3Q     Max
    #> -3.7673 -1.1571 -0.2641  0.5841  8.8966
    #>
    #> Coefficients:
    #>             Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)   6.3299     0.4157  15.226   <2e-16 ***
    #> ITS.Time     -0.2889     0.1438  -2.009   0.0453 *
    #> post5        -0.9877     0.4569  -2.161   0.0314 *
    #> txp5          0.2517     0.2051   1.227   0.2207
    #> post9        -0.7806     0.5100  -1.530   0.1269
    #> txp9          0.2296     0.2114   1.086   0.2781
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 1.72 on 329 degrees of freedom
    #> Multiple R-squared:  0.2519, Adjusted R-squared:  0.2405
    #> F-statistic: 22.16 on 5 and 329 DF,  p-value: < 2.2e-16

Let’s take a look at the interpretations.

### interpret ITS 2 coefficients

    interpret(im12)$its
    #> Interpretations: ITS
    #> --------------------
    #> Note: Some variable names below based on time points (or 'interruptions').
    #> This analysis is for a two-group, single intervention period (interruption).
    #> Positive values indicate higher intervention group values and vice-versa for:
    #> post1, txp1, ixp1, txip1, post2, txp2, ixp2, txip2.
    #>
    #> Intercept is 6.33 and the starting value of the trend for the
    #> control group.
    #>
    #> ITS.Time is -0.289 and the group's slope prior to intervention.
    #> The coefficient is significant.
    #>
    #> post5 is -0.988 and the immediate shift in the group trend
    #> line after this intervention time starts. The coefficient is
    #> significant.
    #>
    #> txp5 is 0.252 and the difference between current and prior intervention
    #> group slopes (e.g., change in the pre-intervention slope).
    #> The coefficient is non-significant.
    #>
    #> post9 is -0.781 and the immediate shift in the group trend
    #> line after this intervention time starts. The coefficient is
    #> non-significant.
    #>
    #> txp9 is 0.23 and the difference between current and prior intervention
    #> group slopes (e.g., change in the pre-intervention slope).
    #> The coefficient is non-significant.
    #>
    #> Summary 1: For this intervention period 1, the results show that the
    #> group's non-significant change in los, post-intervention is -0.0372
    #> (i.e., the value of change per-unit-of-time, such as month or year, in
    #> the intervention period; not the change relative to the prior period).
    #>
    #> Summary 2: For this intervention period 2, the results show that the
    #> group's non-significant change in los, post-intervention is 0.192
    #> (i.e., the value of change per-unit-of-time, such as month or year, in
    #> the intervention period; not the change relative to the prior period).
    #>
    #> If there are additional variables in the model then the coefficients
    #> above represent effects after controlling for the other variables.

Let’s take a look at the graph.

### Graph ITS model 2

    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 )

Now that we see how our intervention group is performing, let’s see how it compares with the control group. When looking at the confidence intervals from the group() results, we see what looks like a change in trends. Let’s test that!

We have an interest in a comparison of the intervention and control group with 2 potential interruptions at months 5 and 9, which is specified with interrupt= c(5, 9) and its=“two”. We also have the option of comparing 2 groups with 1 interruption, if you would like to try that, use 1 time point instead (e.g., its=“two”, interrupt= 5).

### ITS model 3

    im22 <- assess(formula=los ~ ., data=hosprog, intervention = "program",
                   int.time="month", interrupt= c(5, 9), its="two")

Hint: Use top-coding or propensity scores in any DID or ITS model to improve accuracy of results.

### Plot the summary coefficients

plot(review(x)) also works for ITS and DID causal models through ham

    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)

### View ITS model 3 results

    summary(im22$ITS)
    #>
    #> Call:
    #> stats::lm(formula = ITS_formula, data = combined_df, weights = wght_obj,
    #>     offset = offst_obj)
    #>
    #> Residuals:
    #>     Min      1Q  Median      3Q     Max
    #> -3.7673 -1.1755 -0.2888  0.9020  8.8966
    #>
    #> Coefficients:
    #>             Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)  3.08956    0.39418   7.838 1.68e-14 ***
    #> ITS.Time     0.16346    0.14519   1.126   0.2606
    #> ITS.Int      3.24031    0.57773   5.609 2.92e-08 ***
    #> txi         -0.45232    0.20594  -2.196   0.0284 *
    #> post5        0.59672    0.47463   1.257   0.2091
    #> txp5        -0.03780    0.20378  -0.186   0.8529
    #> ixp5        -1.58438    0.66393  -2.386   0.0173 *
    #> txip5        0.28951    0.29147   0.993   0.3209
    #> post9       -0.21293    0.46895  -0.454   0.6499
    #> txp9         0.28005    0.19285   1.452   0.1469
    #> ixp9        -0.56762    0.69886  -0.812   0.4169
    #> txip9       -0.05043    0.28863  -0.175   0.8613
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 1.748 on 708 degrees of freedom
    #> Multiple R-squared:  0.2203, Adjusted R-squared:  0.2081
    #> F-statistic: 18.18 on 11 and 708 DF,  p-value: < 2.2e-16

Hint: The primary variables of interest when examining ITS models with both intervention and control groups are the ones that begin with ‘ixp’ and ‘txip’ as they represent changes in the intercept and slope of the trend line due to the intervention or program.

### interpret coefficients

    interpret(im22)$its
    #> Interpretations: ITS
    #> --------------------
    #> Note: Some variable names below based on time points (or 'interruptions').
    #> This analysis is for a two-group, single intervention period (interruption).
    #> Positive values indicate higher intervention group values and vice-versa for:
    #> post1, txp1, ixp1, txip1, post2, txp2, ixp2, txip2.
    #>
    #> Intercept is 3.09 and the starting value of the trend for the
    #> control group.
    #>
    #> ITS.Time is 0.163 and the control group's slope prior to intervention.
    #> The coefficient is non-significant.
    #>
    #> ITS.Int is 3.24 and the difference in the level between intervention
    #> and control group prior to intervention 1 (intervention - control).
    #> The coefficient is significant.
    #>
    #> txi is -0.452 and the difference between the intervention and
    #> control group's pre-intervention slopes (intervention - control).
    #> The coefficient is significant.
    #>
    #> post5 is 0.597 and the immediate shift in the control group trend
    #> line after this intervention time starts. The coefficient is
    #> non-significant.
    #>
    #> txp5 is -0.0378 and the difference between current and prior intervention
    #> control group slopes (e.g., change in the pre-intervention slope).
    #> The coefficient is non-significant.
    #>
    #> ixp5 is -1.58 and the difference between the intervention and
    #> control groups (intervention - control) in the period immediately
    #> after this intervention started (e.g., 1st year of intervention 1).
    #> The coefficient is significant.
    #>
    #> txip5 is 0.29 and non-significant. This is the difference in both
    #> group's slope changes since the prior intervention (pre-slopes compared
    #> to post-slopes). For example, both have pre-intervention slopes
    #> of 2, the control group's slope remained the same, therefore the
    #> post 1st intervention slope is 0. And the intervention group's slope
    #> increased by 2, then txip1 = 2 (= 2 - 0).
    #>
    #> post9 is -0.213 and the immediate shift in the control group trend
    #> line after this intervention time starts. The coefficient is
    #> non-significant.
    #>
    #> txp9 is 0.28 and the difference between current and prior intervention
    #> control group slopes (e.g., change in the pre-intervention slope).
    #> The coefficient is non-significant.
    #>
    #> ixp9 is -0.568 and the difference between the intervention and
    #> control groups (intervention - control) in the period immediately
    #> after this intervention started (e.g., 1st year of intervention 1).
    #> The coefficient is non-significant.
    #>
    #> txip9 is -0.0504 and non-significant. This is the difference in both
    #> group's slope changes since the prior intervention (pre-slopes compared
    #> to post-slopes). For example, both have pre-intervention slopes
    #> of 2, the control group's slope remained the same, therefore the
    #> post 1st intervention slope is 0. And the intervention group's slope
    #> increased by 2, then txip1 = 2 (= 2 - 0).
    #>
    #> Summary 1: For this intervention period 1, the results show that
    #> the intervention group's non-significant change in los,
    #> post-intervention is -0.0372 (i.e., value of change per-unit-of-time,
    #> such as month or year, in the intervention period; not the change
    #> relative to the prior period). The control group's non-significant
    #> change in los, post-intervention is 0.126. The non-significant
    #> difference between both groups, per-unit-of-time, is -0.163.
    #>
    #> Summary 2: For this intervention period 2, the results show that
    #> the intervention group's non-significant change in los,
    #> post-intervention is 0.192 (i.e., value of change per-unit-of-time,
    #> such as month or year, in the intervention period; not the change
    #> relative to the prior period). The control group's significant
    #> change in los, post-intervention is 0.406. The non-significant
    #> difference between both groups, per-unit-of-time, is -0.213.
    #>
    #> If there are additional variables in the model then the coefficients
    #> above represent effects after controlling for the other variables.

Hint: The 2 summaries at the bottom help us to understand the overall effect on the LOS reduction up to the end of those time periods. In other words, the coefficient interpretations explain what each means while the summaries give us an idea that each value represents the change per unit of time in that period, the degree of lower (or higher) LOS over time for the intervention and/or control groups. For example, the intervention group saw a slight decrease of 0.037 per month after the first interruption but saw an increase of about 0.192 per month after the second interruption.

We will now graph the results based on our coefficients. We use many of the same arguments as the DID model. In this example, I adjusted the position of the text to read better with the pos.text argument. And you’ll notice that within the confidence bands there are dots that represent the means for each time period from the model data. This is done with the argument add.means=TRUE. Although not provided in ITS models with control groups, the ITS models that only have the intervention groups do have counterfactual lines. Counterfactual lines may be added to all ITS models in future versions if can do it in the way I prefer.

### View a partial prediction plot

    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)

The ITS model provides another layer of information. The graph is helpful at showing how the results changed. The control group appears to show a steady increase in LOS, perhaps this reflects a change in practice, patient population or some other issue that impacted results. However, the intervention group may have helped. It appears that their was an “immediate” change in the intercept.

We see a significant decrease of over 1 1/2 days (ixp5). There was a slight decrease at the 2nd interruption denoted by the blue dashed line on the right. This Program+ intervention at month 9 may not have worked as well as the staff would have hoped. The decrease in the intercept (ixp9) isn’t significant and neither is the change in the slope (txip9). Also, we see almost similar trends for the intervention and control groups during the last period of the study. But we do see some possible intervention effect after period 5.

Hint: One option that the literature recommends in graphing DID and ITS results is with lines connecting intervention and control group means. That is a good alternative to graphing the fitted lines for DID and ITS models as above.

We can also perform an ITS on binary outcomes like death within 30-days. We will examine an intervention and control group at months 5 and 9, which is specified with interrupt= c(5, 9) and its =“two”.

### ITS model 4: 30-day death (mortality)

    id22 <- assess(formula=death30 ~ ., data=hosprog, intervention = "program", int.time="month", interrupt= c(5, 9), its="two")

### View ITS model 4 results

    summary(id22$ITS)
    #>
    #> Call:
    #> stats::lm(formula = ITS_formula, data = combined_df, weights = wght_obj,
    #>     offset = offst_obj)
    #>
    #> Residuals:
    #>      Min       1Q   Median       3Q      Max
    #> -0.26332 -0.19587 -0.08391 -0.04646  0.95377
    #>
    #> Coefficients:
    #>              Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)  0.082862   0.075334   1.100   0.2717
    #> ITS.Time    -0.009101   0.027748  -0.328   0.7430
    #> ITS.Int      0.181411   0.110412   1.643   0.1008
    #> txi         -0.013699   0.039359  -0.348   0.7279
    #> post5        0.198078   0.090709   2.184   0.0293 *
    #> txp5        -0.026453   0.038946  -0.679   0.4972
    #> ixp5        -0.302124   0.126886  -2.381   0.0175 *
    #> txip5        0.071580   0.055704   1.285   0.1992
    #> post9        0.170102   0.089622   1.898   0.0581 .
    #> txp9        -0.014163   0.036857  -0.384   0.7009
    #> ixp9        -0.306371   0.133562  -2.294   0.0221 *
    #> txip9        0.020050   0.055162   0.363   0.7164
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 0.334 on 708 degrees of freedom
    #> Multiple R-squared:  0.05044,    Adjusted R-squared:  0.03569
    #> F-statistic: 3.419 on 11 and 708 DF,  p-value: 0.0001208

Let’s take a look at the interpretations.

### interpret ITS 4 coefficients

    interpret(id22)$its
    #> Interpretations: ITS
    #> --------------------
    #> Note: Some variable names below based on time points (or 'interruptions').
    #> This analysis is for a two-group, single intervention period (interruption).
    #> Positive values indicate higher intervention group values and vice-versa for:
    #> post1, txp1, ixp1, txip1, post2, txp2, ixp2, txip2.
    #>
    #> Intercept is 0.0829 and the starting value of the trend for the
    #> control group.
    #>
    #> ITS.Time is -0.0091 and the control group's slope prior to intervention.
    #> The coefficient is non-significant.
    #>
    #> ITS.Int is 0.181 and the difference in the level between intervention
    #> and control group prior to intervention 1 (intervention - control).
    #> The coefficient is non-significant.
    #>
    #> txi is -0.0137 and the difference between the intervention and
    #> control group's pre-intervention slopes (intervention - control).
    #> The coefficient is non-significant.
    #>
    #> post5 is 0.198 and the immediate shift in the control group trend
    #> line after this intervention time starts. The coefficient is
    #> significant.
    #>
    #> txp5 is -0.0265 and the difference between current and prior intervention
    #> control group slopes (e.g., change in the pre-intervention slope).
    #> The coefficient is non-significant.
    #>
    #> ixp5 is -0.302 and the difference between the intervention and
    #> control groups (intervention - control) in the period immediately
    #> after this intervention started (e.g., 1st year of intervention 1).
    #> The coefficient is significant.
    #>
    #> txip5 is 0.0716 and non-significant. This is the difference in both
    #> group's slope changes since the prior intervention (pre-slopes compared
    #> to post-slopes). For example, both have pre-intervention slopes
    #> of 2, the control group's slope remained the same, therefore the
    #> post 1st intervention slope is 0. And the intervention group's slope
    #> increased by 2, then txip1 = 2 (= 2 - 0).
    #>
    #> post9 is 0.17 and the immediate shift in the control group trend
    #> line after this intervention time starts. The coefficient is
    #> non-significant.
    #>
    #> txp9 is -0.0142 and the difference between current and prior intervention
    #> control group slopes (e.g., change in the pre-intervention slope).
    #> The coefficient is non-significant.
    #>
    #> ixp9 is -0.306 and the difference between the intervention and
    #> control groups (intervention - control) in the period immediately
    #> after this intervention started (e.g., 1st year of intervention 1).
    #> The coefficient is significant.
    #>
    #> txip9 is 0.0201 and non-significant. This is the difference in both
    #> group's slope changes since the prior intervention (pre-slopes compared
    #> to post-slopes). For example, both have pre-intervention slopes
    #> of 2, the control group's slope remained the same, therefore the
    #> post 1st intervention slope is 0. And the intervention group's slope
    #> increased by 2, then txip1 = 2 (= 2 - 0).
    #>
    #> Summary 1: For this intervention period 1, the results show that
    #> the intervention group's non-significant change in death30,
    #> post-intervention is 0.0223 (i.e., value of change per-unit-of-time,
    #> such as month or year, in the intervention period; not the change
    #> relative to the prior period). The control group's non-significant
    #> change in death30, post-intervention is -0.0356. The non-significant
    #> difference between both groups, per-unit-of-time, is 0.0579.
    #>
    #> Summary 2: For this intervention period 2, the results show that
    #> the intervention group's non-significant change in death30,
    #> post-intervention is 0.0282 (i.e., value of change per-unit-of-time,
    #> such as month or year, in the intervention period; not the change
    #> relative to the prior period). The control group's significant
    #> change in death30, post-intervention is -0.0497. The significant
    #> difference between both groups, per-unit-of-time, is 0.0779.
    #>
    #> If there are additional variables in the model then the coefficients
    #> above represent effects after controlling for the other variables.

Let’s take a look at the graph (intervention is brown, control is gray).

### Graph ITS model 4

    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 )

In the latest version of ham, we now have the option of an interrupted time series that tests as many interruptions as the data permits. This example demonstrates the US unemployment rate since 1929, highlighting key moments during that period.

### 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)

We see that many things increased or decreased unemployment rates. The biggest shift after the stock market crash in unemployment seems to be FDR’s New Deal.

### Graph ITS model 5 results

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

## 6. Cronbach’s alpha example

A very commonly used method of collecting information on patient experiences, behaviors, beliefs, or attitudes is with a survey or questionnaire. Many states, provinces, and countries use surveys to ask important questions of their constituents.

A widely used survey is the HCAHPS survey for hospitals in the United States. Medicare provides the results for the overall score and subscales on individual dimensions to help inform people when selecting which hospital they would like to receive care.

The developers of the HCAHPS survey use careful measurement and assessment methods to build reliable and valid survey tools. We will only cover parts of reliability here but Cronbach’s alpha should be a good instrument in our pursuit of developing an accurate measure of whichever dimension we are interested in.

In this example, we’ll check how reliably we measure “patient satisfaction” with our 5 survey items (i1 to i5). The function alpha() and interpret() will provide us with helpful information.

### An example of calculating Cronbach’s alpha:

    alpha(items=c("i1","i2","i3","i4","i5"), data=cas)
    #> Scale statistics
    #> Cronbach's alpha   = 0.919
    #> Mean               = 3.772
    #> Variance           = 0.453
    #> Standard Deviation = 0.673
    #> Items              = 5
    #>
    #> Item statistics
    #>    Mean Variance Std. Dev.
    #> i1 3.50    0.434     0.659
    #> i2 3.81    0.681     0.825
    #> i3 3.88    0.652     0.808
    #> i4 3.82    0.594     0.770
    #> i5 3.85    0.634     0.796
    #>
    #> Scale statistics if item deleted
    #>    Alpha  Mean Variance Std. Dev.
    #> i1 0.930 3.840    0.528     0.727
    #> i2 0.891 3.763    0.436     0.660
    #> i3 0.883 3.745    0.433     0.658
    #> i4 0.910 3.760    0.472     0.687
    #> i5 0.885 3.752    0.439     0.662
    #>
    #> Sample
    #> Total    = 100
    #> Valid    = 100
    #> Excluded = 0

### Interpret Cronbach’s alpha:

    interpret(alpha(items=c("i1","i2","i3","i4","i5"), data=cas))
    #> Interpretations: Alpha
    #> ----------------------
    #> Your 5 item scale has a Cronbach's alpha of 0.919. This is
    #> generally considered as being in the 'excellent' range.
    #>
    #> The scale mean is 3.77 and has a standard deviation of 0.673.
    #>
    #> Removing one of these item(s): i1, can improve the Cronbach's
    #> alpha in a new scale to a higher level than the current alpha
    #> based on all items.
    #>
    #> 0 row(s) of data excluded from the analysis because of missing
    #> data.

## 7. Discussion

The ham package can be helpful in identifying differences between intervention and control groups. And can help spot trends in the data with the group() function.

We can then use that information in guiding our regression analyses, especially with differences-in-differences and interrupted time series models.

A model can benefit by using top-coding and/or propensity scores. Both can be added in the assess() function. We can determine which variables are most important in explaining our outcome with importance(), therefore providing an opportunity to make more efficient use of our resources.

If we select a measure based on multiple survey questions, alpha() can help us identify the most reliable set of items that represent the dimension we are interested in.

And we can get an interpretation of all the numbers we are looking at from assess() and alpha() by using interpret().

Hopefully, the ham package can help you from start to finish in your study.

## 8. References

See Description for author list.