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", dataf=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.571974For 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", dataf=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 12There 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.30240495Graph 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)
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)
#>
#> 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-06ham provides coefficient interpretations with interpret()
# Model summary
interpret(assess(hp ~ mpg+wt, data=mtcars, regression="ols"))$model
#> Interpretations: 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.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.
model results
summary(m1$model)
#>
#> Call:
#> stats::lm(formula = primary_formula, data = combined_df)
#>
#> 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-16Descriptive 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.5609importance
We can examine variable importance based on partial chi-square (i.e., which variables explain the outcome the most).
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))
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()
Code to view DID model results
summary(dm1$DID)
#>
#> Call:
#> stats::lm(formula = DID_formula, data = combined_df)
#>
#> 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-16Hint: 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.494.
#>
#> Post.All is the change in the control group's los
#> value in the 2nd time period (Time 2). There was a
#> significant increase for the control group
#> at time 2: 1.563.
#>
#> 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.066.
#>
#> 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.545.
#> In other words, there was a significant decrease in the
#> mean los trend by -3.545 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.
View DID 30-day readmissions results
summary(dm3$DID)
#>
#> Call:
#> stats::lm(formula = DID_formula, data = combined_df)
#>
#> 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-12Significant 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
View ITS model 1 results
summary(im11$ITS)
#>
#> Call:
#> stats::lm(formula = ITS_formula, data = combined_df)
#>
#> 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-16Hint: 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.003 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 was a non-significant change in the los trend. This gives
#> a total post-intervention trend in the los of -0.087
#> over time (i.e., the total combined value of change 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
View ITS model 2 results
summary(im12$ITS)
#>
#> Call:
#> stats::lm(formula = ITS_formula, data = combined_df)
#>
#> 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-16Now 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.
View ITS model 3 results
summary(im22$ITS)
#>
#> Call:
#> stats::lm(formula = ITS_formula, data = combined_df)
#>
#> 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-16Hint: 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 the 1st intervention start. The coefficient is
#> non-significant.
#>
#> txp5 is -0.038 and the difference between pre- and post-intervention
#> control group slopes (e.g., change in the pre-intervention
#> slope). The coefficient is non-significant.
#>
#> ixp5 is -1.584 and the difference between the intervention and
#> control groups (intervention - control) in the period immediately
#> after the 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 pre-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 the 2nd intervention start. The coefficient is
#> non-significant.
#>
#> txp9 is 0.28 and the difference between 1st and 2nd intervention
#> control group slopes (e.g., change in the 1st 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 the 2nd intervention started (e.g., 1st year of intervention 2).
#> The coefficient is non-significant.
#>
#> txip9 is -0.05 and non-significant. This is the difference in both group's
#> slope changes since the 1st intervention (1st intervention slope compared
#> to the 2nd). For example, both have 1st intervention slopes of 2, the control
#> group's slope remained the same, therefore the 2nd intervention slope is 0. And
#> the intervention group's slope increased by 2, then txip2 = 2 (= 2 - 0).
#>
#> Summary 1: For the 1st intervention period, the results show that the
#> intervention group's non-significant shift in los,
#> post 1st intervention was -0.037. The control group's non-significant
#> shift in los, post 1st intervention was 0.126. The non-significant
#> difference between both groups is -0.163.
#>
#> Summary 2: For the 2nd intervention period, the results show that the
#> intervention group's non-significant shift in los,
#> post 2nd intervention was 0.192. The control group's significant
#> shift in los, post 2nd intervention was 0.406. The non-significant
#> difference between both groups 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 on the actual number in a lower (or higher) LOS over time for the intervention and/or control groups.
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("txp5"=3, "post9"=4), 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
View ITS model 4 results
summary(id22$ITS)
#>
#> Call:
#> stats::lm(formula = ITS_formula, data = combined_df)
#>
#> 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.00012086. 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 = 0Interpret 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.92. This is
#> generally considered as being in the 'excellent' range.
#>
#> The scale mean is 3.77 and has a standard deviation of 0.67.
#>
#> 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.