Introduction
Interrupted time series analysis (ITSA) is a study design used to study the effects of an intervention across time. An important feature of the ITSA is the time when the intervention occurs. The time before and after the intervention are of interest because we want to visualize if the trends are similar or different. Additionally, we want to visualize the change immediately after the intervention is implemented. I call this period the index date.
In this article, I’ll review the single-group ITSA and multiple
groups ITSA. Then I’ll review how to perform an ITSA in
Stata.
Single-group ITSA
Let’s look at Figure 1 and visualize an ITSA (this
is a modified version based on Linden’s
paper; I highly encourage you to read his paper to learn more about
ITSA for Stata). Figure 1 illustrates the
basic features of a single-group ITSA. The single group has an
increasing slope in the period before the intervention and then there is
a jump in the level immediately after the intervention at the index date
\(X\) followed by a steeper slope in
the period after the intervention.
The structure of the ITSA is written as:
\(Y_i = \beta_0 + \beta_1(T_i) + \beta_2(X_i) + \beta_3(T_i * X_i) + \epsilon_i\),
where \(Y_i\) denotes the dependent variable or outcome of interest, \(\beta_0\) denotes the intercept, \(\beta_1\) denotes the change outcome due to a one-unit change in \(T_i\), \(\beta_2\) denotes the change in outcome due to a change in \(X_i\) (which is another way of saying that this is the chance in the outcome when a subject goes from the before \(X\) period to the after \(X\) period), \(\beta_3\) denotes the differences in the trends (e.g., slopes) after and before the intervention \(X\), \(T_i\) denotes the time for the \(i\)-th subject, \(X_i\) denotes the index date when the intervention was implemented, and \(\epsilon_i\) denotes the error term. The dash line denotes the counterfactual.
Multiple groups ITSA
Things get a little more complicated when we have more than one group in the ITSA. Let’s take a look at Figure 2, which provides an ITSA for two groups. The first group is the same one from Figure 1, which we’ll call the control. The new group is the Treatment group, which we denote as \(G\). We can see that the Treatment group starts out lower than the Control group and has a less steep slope. We also can see that immediately after implementation of the intervention at the index date \(X\), there is a jump in the level and the slope increases.
The regression model can be expanded to include the Treatment group:
\(Y_i = \beta_0 + \beta_1(T_i) + \beta_2(X_i) + \beta_3(T_i * X_i) + \beta_4(G_i) + \beta_5(G_i * T_i) + \beta_6(G_i * X_i) + \beta_7(G_i * T_i * X_i) + \epsilon_i\),
where \(G_i\) is a new variable that denotes the groups, \(\beta_4\) denotes the difference in the outcomes between the groups at beginning of the study (\(T = 0\)), \(\beta_5\) denotes the difference in the slopes between the groups before the intervention, \(\beta_6\) denotes the difference in the level changes between the groups immediately after implementation of the intervention \(X\), and \(\beta_7\) denotes the difference in the slopes after and before the intervention between the groups (e.g., difference-in-differences).
ITSA with Stata
We will use Stata to explore how to perform an ITSA.
First, we will need to download the data. You can download the data from my GitHub site.
The sleep_cat3pw.dta comes from the excellent text book
by Michael N. Mitchell, whcih you can find on the Stata
website (link).
(Note: I have the first edition where this data came from. The link is
for the second edition of the text.)
Once you’ve downloaded the data, you can load it onto
Stata. Make sure you use the correct path. Here is an
example of how I loaded the sleep_cat3pw.dta onto
Stata.
clear all
cd "C:\Users\carn3_a5ada\Dropbox\STATA references\Mitchels book - Visualing in Stata\ivrm2"
use sleep_cat3pw.dta, replaceC:\Users\carn3_a5ada\Dropbox\STATA references\Mitchels book - Visualing in Stat
> a\ivrm2Once you’ve loaded the data, let’s take a look using the
summarize command.
summarizeRunning C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
id | 600 38 21.66677 1 75
group | 600 2 .8171778 1 3
obsday | 600 23.815 15.02877 1 54
sleep | 600 362.4817 39.72333 207 531The data sleep_cat3pw.dta contains information on the
amount of sleep for subjects in a study. The amount of sleep was
captured once per week for a total of eight weeks. This means that a
subject can have their sleep captured any day between Sunday and
Saturday for each week. For example Subject A can have their sleep
measured on Monday of Week 1, and Subject B can have their sleep
measured on Wednesday of Week 1. This corresponds to Day 2 for Subject A
and Day 4 for Subject B. You will notice that the days correspond to
when the subject’s sleep was measured for that week.
The data contains four variables:
id: This unique identifier of the subjectgroup: Intervention group (Control, Medication, and Education)obsday: The observation day when the sleep measurement was capturedsleep: Duration of sleep (minutes)
There are there groups in the sleep_cat3pw.dta data:
Control, Medication, and Education.
The Control denotes the subjects who did not receive an intervention.
The Medication denotes the subjects who received a medication for their sleep.
The Education denotes the subjects who received education to help them with their sleep.
You can view the number of observations for each group using the
tabulate command in Stata.
tabulate group, mRunning C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Treatment |
group | Freq. Percent Cum.
------------+-----------------------------------
Control | 200 33.33 33.33
Medication | 200 33.33 66.67
Education | 200 33.33 100.00
------------+-----------------------------------
Total | 600 100.00The time variable is obsday, which is the day when the
measurement (sleep duration) was captured.
The outcome variable is sleep, which is the duration
(minutes) of sleep that the subject experienced during that week of the
study.
We can view the data structure for a subject in each of the three
groups using the list command:
list if id == 12 | id == 45 | id == 68Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
+----------------------------------+
| id group obsday sleep |
|----------------------------------|
89. | 12 Control 1 360 |
90. | 12 Control 8 356 |
91. | 12 Control 14 342 |
92. | 12 Control 20 360 |
93. | 12 Control 26 366 |
|----------------------------------|
94. | 12 Control 32 350 |
95. | 12 Control 38 383 |
96. | 12 Control 45 385 |
353. | 45 Medication 1 332 |
354. | 45 Medication 7 338 |
|----------------------------------|
355. | 45 Medication 13 349 |
356. | 45 Medication 18 357 |
357. | 45 Medication 26 370 |
358. | 45 Medication 32 415 |
359. | 45 Medication 37 402 |
|----------------------------------|
360. | 45 Medication 44 413 |
537. | 68 Education 1 360 |
538. | 68 Education 6 358 |
539. | 68 Education 13 367 |
540. | 68 Education 19 358 |
|----------------------------------|
541. | 68 Education 25 368 |
542. | 68 Education 33 389 |
543. | 68 Education 40 424 |
544. | 68 Education 46 419 |
+----------------------------------+Notice that the days when the sleep measurement was captured varies between the subjects. However, they all occurred within the same corresponding week.
Conducting an ITSA in Stata
I know two methods to perform an ITSA in Stata. The
first method uses the mkplines command and the second uses
several interaction terms.
The first ITSA example will be a single-group ITSA.
The second ITSA example will be a multiple groups ITSA.
Single-group ITSA
To start off, we’ll go over how to construct an ITSA model with a single group. This will get us comfortable with understanding the parameters that estimates the slopes before and after the intervention, the level change immediately after the intervention, and the differences in the slopes after and before the intervention.
Method 1: Splines
The splines method takes advantage of the mkspline
command which will generate several new variables that take into
consideration the time after the intervention. This makes interpretation
of the coefficients easier than using a lot of interaction terms.
To begin, we need to figure out where on the X-axis the intervention
occurred. In the sleep_cat3pw.dta data, the intervention
occurred at day 31 (obsday = 31). We will use the
mkspline command to create a knot at
obsday = 31. These knots are the same as time sequences or
intervals. They help to keep track of when to start the clock with
respect to the index date. For example, knot1 should
correspond to the main time variable obsday, but
knot2 should restart the clock to start after the
implementation of the intervention.
But before we do, we only want to keep the Control group for the single-ITSA example. Let’s put all these into a code chunk:
// Keep only the control group
keep if group == 1
// Use mkspline to create knots. We will create a knot at obsday = 31 (index date)
/* The marginal option => estimate the change in slope b/w knot2 v. knot1 */
mkspline knot1 31 knot2 = obsday, marginal Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
(400 observations deleted)Once the knots have been created, we want to create a variable that partitions time periods to before and after the implementation of the intervention. Since there are only two time periods (before the intervention and after the intervention), we will create a binary indicator.
// Create a variable to separate the treatment periods.
gen period = .
replace period = 0 if obsday < 31 /* Baseline phase */
replace period = 1 if obsday >= 31 & !missing(obsday) /* Treatment phase */
tab period, m Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
(200 missing values generated)
(126 real changes made)
(74 real changes made)
period | Freq. Percent Cum.
------------+-----------------------------------
0 | 126 63.00 63.00
1 | 74 37.00 100.00
------------+-----------------------------------
Total | 200 100.00Let’s take a look at the date. Notice the sequence of the numbers for
knot1, knot2 and period.
knot1 should reflect the obsday variable
exactly, but knot2 restarts the time to after the
implementation of the intervention, and period should be
0 during the period before the intervention and
1 after the intervention.
list if id == 1 | id == 3 | id == 6Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
+--------------------------------------------------------+
| id group obsday sleep knot1 knot2 period |
|--------------------------------------------------------|
1. | 1 Control 1 353 1 0 0 |
2. | 1 Control 9 345 9 0 0 |
3. | 1 Control 15 340 15 0 0 |
4. | 1 Control 21 337 21 0 0 |
5. | 1 Control 26 324 26 0 0 |
|--------------------------------------------------------|
6. | 1 Control 33 332 33 2 1 |
7. | 1 Control 41 326 41 10 1 |
8. | 1 Control 49 329 49 18 1 |
17. | 3 Control 1 345 1 0 0 |
18. | 3 Control 9 360 9 0 0 |
|--------------------------------------------------------|
19. | 3 Control 16 347 16 0 0 |
20. | 3 Control 24 348 24 0 0 |
21. | 3 Control 29 351 29 0 0 |
22. | 3 Control 37 348 37 6 1 |
23. | 3 Control 44 334 44 13 1 |
|--------------------------------------------------------|
24. | 3 Control 50 335 50 19 1 |
41. | 6 Control 1 353 1 0 0 |
42. | 6 Control 8 340 8 0 0 |
43. | 6 Control 15 334 15 0 0 |
44. | 6 Control 22 330 22 0 0 |
|--------------------------------------------------------|
45. | 6 Control 28 327 28 0 0 |
46. | 6 Control 35 306 35 4 1 |
47. | 6 Control 41 321 41 10 1 |
48. | 6 Control 49 319 49 18 1 |
+--------------------------------------------------------+Linear mixed effects model - Using splines
Now that we have our necessary variables, we can begin to construct our ITSA model. We will use the linear mixed effects model for this ITSA. You can use other models such as the generalized estimating equation (GEE) model or a simple linear regression. I chose to use the linear mixed effects model because this was the one that Mitchell used in his textbook, but also because I was working on another study that used a similar model using triple interaction terms, and I wanted to see if I can replicate this result. (Note: In Method 2, we will see that we generate the same answers.)
After the model completes, it is a good idea to use the
predict command to estimate the outcomes for each subject,
which will generate a new column called linear_spline. We
will use this to plot our results.
// Linear mixed effects model for the single-group ITSA
xtmixed sleep i.period c.knot1 c.knot2 || id: period knot1 knot2, covariance(unstructured)
predict linear_splineRunning C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Performing EM optimization:
Performing gradient-based optimization:
Iteration 0: Log likelihood = -815.91977
Iteration 1: Log likelihood = -815.91239
Iteration 2: Log likelihood = -815.91238
Computing standard errors:
Mixed-effects ML regression Number of obs = 200
Group variable: id Number of groups = 25
Obs per group:
min = 8
avg = 8.0
max = 8
Wald chi2(3) = 2.16
Log likelihood = -815.91238 Prob > chi2 = 0.5404
------------------------------------------------------------------------------
sleep | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.period | -4.061266 3.17595 -1.28 0.201 -10.28601 2.163481
knot1 | -.0393757 .3208933 -0.12 0.902 -.668315 .5895635
knot2 | -.1248309 .4277671 -0.29 0.770 -.963239 .7135771
_cons | 351.3585 2.465165 142.53 0.000 346.5268 356.1901
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
id: Unstructured |
sd(period) | 10.22294 3.868926 4.868934 21.46434
sd(knot1) | 1.553309 .2347396 1.155108 2.088784
sd(knot2) | 1.896657 .3445179 1.328533 2.707728
sd(_cons) | 10.30786 2.098042 6.916981 15.36104
corr(period,knot1) | -.389369 .2899772 -.7935565 .2532219
corr(period,knot2) | .5930494 .3983103 -.4791216 .9550717
corr(period,_cons) | .1307331 .3606115 -.5281599 .6913624
corr(knot1,knot2) | -.5681853 .1626752 -.8060611 -.1722627
corr(knot1,_cons) | .077858 .2518023 -.3956807 .5186924
corr(knot2,_cons) | -.1536023 .2735296 -.6068313 .3750053
-----------------------------+------------------------------------------------
sd(Residual) | 8.254314 .6082031 7.144336 9.536744
------------------------------------------------------------------------------
LR test vs. linear model: chi2(10) = 421.77 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.
(option xb assumed)We can plot the predicted lines. Notice that the slope in the period before the intervention is slightly negative. Then, at the index date (immediately after the intervention) there is a level decrease or drop followed by a steeper negative slope in the period after the intervention.
*** Plot single-group ITSA
graph twoway (lfit linear_spline obsday if group == 1 & period == 0, lcol(navy)) ///
(lfit linear_spline obsday if group == 1 & period == 1, lcol(navy) ///
xline(31) ///
legend(order(1 "Before intervention" 2 "After intervention")))Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Slope before the intervention
To estimate the slope before the intervention, we can use the
margins command:
margins, dydx(knot1) /* Slope before the intervention */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Average marginal effects Number of obs = 200
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: knot1
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
knot1 | -.0393757 .3208933 -0.12 0.902 -.668315 .5895635
------------------------------------------------------------------------------The slope in the period before the intervention was -0.04 minutes per week (95% CI: -0.67, 0.59), which was not statistically significant. Since the slope is negative, this means that the predicted fit should also decrease. We can confirm this with the figure above.
Level change immediately after the intervention
To estimate the level change immediately after the intervention, we
can use the margins command:
margins, dydx(period) /* Level change immediately after the intervention */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Average marginal effects Number of obs = 200
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: 1.period
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.period | -4.061266 3.17595 -1.28 0.201 -10.28601 2.163481
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.The sleep duration dropped by -4.06 minutes (95% CI: -10.29, 2.16) at
obsday = 31, but this was not statistically
significant.
Comparing the Slopes after the intervention
The slope after the intervention can be estimated in two ways using
(1) margins command and (2) lincom
command:
/* Slope after the intervention */
margins, expression(_b[knot1] + _b[knot2]) /* Method 1 */
lincom knot1+knot2 /* Method 2 */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
warning: option expression() does not contain option predict() or xb().
warning: prediction constant over observations.
Predictive margins Number of obs = 200
Expression: _b[knot1] + _b[knot2]
------------------------------------------------------------------------------
| Delta-method
| Margin std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_cons | -.1642067 .3733492 -0.44 0.660 -.8959576 .5675442
------------------------------------------------------------------------------
( 1) [sleep]knot1 + [sleep]knot2 = 0
------------------------------------------------------------------------------
sleep | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
(1) | -.1642067 .3733492 -0.44 0.660 -.8959576 .5675442
------------------------------------------------------------------------------The slope after the intervention was decreasing by -0.16 minutes per week (95% CI: -0.90, 0.57), which was not statistically significant.
Difference in slopes after and before the intervention
The differences in the slopes after and before the intervention can
be estimated using the margins command with the
pwcompare(effects) option:
margins, dydx(knot2) /* Diff in slopes after and before the intervention */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Average marginal effects Number of obs = 200
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: knot2
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
knot2 | -.1248309 .4277671 -0.29 0.770 -.963239 .7135771
------------------------------------------------------------------------------The slope in the period after the intervention has a greater decrease compared to the slope in the period before the intervention by -0.12 minutes per week (95% CI: -0.96, 0.71), which is not statistically significant. In other words, the sleep duration was decreasing at a faster rate after the intervention compared to before.
Method 2: Using an interaction term
The mkspline method makes it easy to perform an ITSA in
Stata. However, we can do a lot of the
mkspline work ourselves using interaction terms,
particularly the interaction between the group and
obsday variables. We also want to restrict our sample to
the Control group.
Let’s start over and reload the data:
clear all
cd "C:\Users\carn3_a5ada\Dropbox\STATA references\Mitchels book - Visualing in Stata\ivrm2"
use sleep_cat3pw.dta, replace
// Keep only the control group
keep if group == 1Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
C:\Users\carn3_a5ada\Dropbox\STATA references\Mitchels book - Visualing in Stat
> a\ivrm2
(400 observations deleted)Next, we want to create a variable that partitions time periods to before and after the implementation of the intervention. Since there are only two time periods (before the intervention and after the intervention), we will create a binary indicator.
// Create a variable to separate the treatment periods.
gen period = .
replace period = 0 if obsday < 31 /* Baseline phase */
replace period = 1 if obsday >= 31 & !missing(obsday) /* Treatment phase */
tab period, m Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
(200 missing values generated)
(126 real changes made)
(74 real changes made)
period | Freq. Percent Cum.
------------+-----------------------------------
0 | 126 63.00 63.00
1 | 74 37.00 100.00
------------+-----------------------------------
Total | 200 100.00Linear mixed effects model - Using an interaction term
Now that we have our necessary variables, we can begin to construct
our ITSA model. We’ll construct the same type of model: linear mixed
effects model. We’ll also created fitted values using the
predict linear_int command. Not how the interaction terms
are constructed:
// Linear mixed effects model for the single-group ITSA
xtmixed sleep i.period c.obsday i.period#c.obsday || id: period obsday i.period#c.obsday, covariance(unstructured)
predict linear_intRunning C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Performing EM optimization:
Performing gradient-based optimization:
Iteration 0: Log likelihood = -819.10509
Iteration 1: Log likelihood = -816.20946 (not concave)
Iteration 2: Log likelihood = -815.98347
Iteration 3: Log likelihood = -815.91256
Iteration 4: Log likelihood = -815.91238
Iteration 5: Log likelihood = -815.91238
Computing standard errors:
Mixed-effects ML regression Number of obs = 200
Group variable: id Number of groups = 25
Obs per group:
min = 8
avg = 8.0
max = 8
Wald chi2(3) = 2.16
Log likelihood = -815.91238 Prob > chi2 = 0.5404
------------------------------------------------------------------------------
sleep | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.period | -.1915202 12.97583 -0.01 0.988 -25.62367 25.24063
obsday | -.0393758 .3208933 -0.12 0.902 -.6683152 .5895636
|
period#|
c.obsday |
1 | -.1248306 .427767 -0.29 0.770 -.9632385 .7135773
|
_cons | 351.3585 2.465166 142.53 0.000 346.5268 356.1901
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
id: Unstructured |
sd(period) | 53.37212 11.43859 35.06605 81.23479
sd(obsday) | 1.55331 .2347451 1.1551 2.088798
sd(1.period#obsday) | 1.896656 .3445258 1.328522 2.70775
sd(_cons) | 10.30787 2.098045 6.916985 15.36105
corr(period,obsday) | .5513518 .1901583 .0846336 .8196685
corr(period,1.period#obsday) | -.9880365 .0144753 -.9988931 -.8772388
corr(period,_cons) | .1942535 .2899708 -.3746878 .6569181
corr(obsday,1.period#obsday) | -.5681866 .1627069 -.8060944 -.1721743
corr(obsday,_cons) | .0778578 .2518162 -.395704 .5187123
corr(1.period#obsday,_cons) | -.1536021 .2735914 -.6069096 .3751121
-----------------------------+------------------------------------------------
sd(Residual) | 8.254318 .6082069 7.144333 9.536756
------------------------------------------------------------------------------
LR test vs. linear model: chi2(10) = 421.77 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.
(option xb assumed)*** Plot single-group ITSA
graph twoway (lfit linear_int obsday if group == 1 & period == 0, lcol(navy)) ///
(lfit linear_int obsday if group == 1 & period == 1, lcol(navy) ///
xline(31) ///
legend(order(1 "Before intervention" 2 "After intervention")))Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
The figure looks exactly the same as the one created with the
mkspline method.
Slope before the intervention
To estimate the slope before the intervention, we can use the
margins command. We have to use the if
statement to limit the estimation to the period before the intervention
\(obsday < 31\):
margins if obsday < 31, dydx(obsday) /* Slope before the intervention */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Average marginal effects Number of obs = 126
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday | -.0393758 .3208933 -0.12 0.902 -.6683152 .5895636
------------------------------------------------------------------------------This is exactly the same as the slope we estimate with the
mkpline approach. The slope in the period before the
intervention was -0.04 minutes per week (95% CI: -0.67, 0.59), which was
not statistically significant. Since the slope is negative, this means
that the predicted fit should also decrease. We can confirm this with
the figure above.
Level change immediately after the intervention
To estimate the level change immediately after the intervention, we
can use the margins command. We have to use the
at operator to estimate the level change at
obsday = 31:
margins, dydx(period) at(obsday = 31) /* Estimate the level change at obsday == 31 */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Conditional marginal effects Number of obs = 200
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: 1.period
At: obsday = 31
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
1.period | -4.061269 3.175945 -1.28 0.201 -10.28601 2.163468
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.This is exactly the same as the slope we estimate with the
mkpline approach. The sleep duration dropped by -4.06
minutes (95% CI: -10.29, 2.16) at obsday = 31, but this was
not statistically significant.
Slope after the intervention
The slope after the intervention can be estimated in two ways using
(1) margins command and (2) lincom
command:
/* Slope after the intervention */
margins, expression(_b[obsday] + _b[1.period#c.obsday]) /* Method 1 */
lincom obsday + 1.period#c.obsday /* Method 2 */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
warning: option expression() does not contain option predict() or xb().
warning: prediction constant over observations.
Predictive margins Number of obs = 200
Expression: _b[obsday] + _b[1.period#c.obsday]
------------------------------------------------------------------------------
| Delta-method
| Margin std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_cons | -.1642064 .3733487 -0.44 0.660 -.8959565 .5675437
------------------------------------------------------------------------------
( 1) [sleep]obsday + [sleep]1.period#c.obsday = 0
------------------------------------------------------------------------------
sleep | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
(1) | -.1642064 .3733487 -0.44 0.660 -.8959565 .5675437
------------------------------------------------------------------------------This is exactly the same as the slope we estimate with the
mkpline approach. The slope after the intervention was
decreasing by -0.16 minutes per week (95% CI: -0.90, 0.57), which was
not statistically significant.
Note: You can estimate the both the slopes before and after the
intervention with the following margins command:
margins, dydx(obsday) over(period) /* Est slopes before and after intervention */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Average marginal effects Number of obs = 200
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
Over: period
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday |
period |
0 | -.0393758 .3208933 -0.12 0.902 -.6683152 .5895636
1 | -.1642064 .3733487 -0.44 0.660 -.8959565 .5675437
------------------------------------------------------------------------------Difference in slopes after and before the intervention
The differences in the slopes after and before the intervention can
be estimated using the margins command with the
pwcompare(effects) option:
margins, dydx(obsday) over(period) pwcompare(effects) /* Diff between slopes after and before the intervention */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Pairwise comparisons of average marginal effects
Number of obs = 200
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
Over: period
------------------------------------------------------------------------------
| Contrast Delta-method Unadjusted Unadjusted
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday |
period |
1 vs 0 | -.1248306 .427767 -0.29 0.770 -.9632385 .7135773
------------------------------------------------------------------------------The answer is exactly the same with the mkspline method.
The slope in the period after the intervention has a greater decrease
compared to the slope in the period before the intervention by -0.12
minutes per week (95% CI: -0.96, 0.71), which is not statistically
significant. In other words, the sleep duration was decreasing at a
faster rate after the intervention compared to before.
Multiple groups ITSA
In this section, we’ll examine how to construct ITSA models with
multiple groups. This will be a little more complicated because of the
interpretations of the coefficients due to the introduction of the
group variable. The example will include three groups:
Control, Medication, and
Education.
Method 1: Splines
Let’s start over and reload the data:
clear all
cd "C:\Users\carn3_a5ada\Dropbox\STATA references\Mitchels book - Visualing in Stata\ivrm2"
use sleep_cat3pw.dta, replaceRunning C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
C:\Users\carn3_a5ada\Dropbox\STATA references\Mitchels book - Visualing in Stat
> a\ivrm2The splines method takes advantage of the mkspline
command which will generate several new variables that take into
consideration the time after the intervention. This makes interpreting
the difference-in-differences estimation much easier for an ITSA. I like
this method because it takes the guess-work out of trying to make sure
your interaction terms are properly set up because the
mkspline command takes care that for you.
To begin, we need to figure out where on the X-axis the intervention
occurred. In the sleep_cat3pw.dta data, the intervention
occurred at day 31 (obsday = 31). We will use the
mkspline command to create a knot at
obsday = 31. These knots are the same as time sequences or
intervals. They help to keep track of when to start the clock with
respect to the index date. For example, knot1 should
correspond to the main time variable obsday, but
knot2 should restart the clock to start after the
implementation of the intervention.
We use the marginal option because we want to make
comparisons to the slopes before and after the intervention.
// Use mkspline to create knots. We will create a knot at obsday = 31 (index date)
/* The marginal option => estimate the change in slope b/w knot2 v. knot1 */
mkspline knot1 31 knot2 = obsday, marginal Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...Once the knots have been created, let’s take a look at the data
structure. You should notice that two new variables have been added
(knot1 and knot2). knot1 should
reflect the obsday variable exactly, but knot2
restarts the time to after the implementation of the intervention:
list if id == 12 | id == 45 | id == 68Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
+--------------------------------------------------+
| id group obsday sleep knot1 knot2 |
|--------------------------------------------------|
89. | 12 Control 1 360 1 0 |
90. | 12 Control 8 356 8 0 |
91. | 12 Control 14 342 14 0 |
92. | 12 Control 20 360 20 0 |
93. | 12 Control 26 366 26 0 |
|--------------------------------------------------|
94. | 12 Control 32 350 32 1 |
95. | 12 Control 38 383 38 7 |
96. | 12 Control 45 385 45 14 |
353. | 45 Medication 1 332 1 0 |
354. | 45 Medication 7 338 7 0 |
|--------------------------------------------------|
355. | 45 Medication 13 349 13 0 |
356. | 45 Medication 18 357 18 0 |
357. | 45 Medication 26 370 26 0 |
358. | 45 Medication 32 415 32 1 |
359. | 45 Medication 37 402 37 6 |
|--------------------------------------------------|
360. | 45 Medication 44 413 44 13 |
537. | 68 Education 1 360 1 0 |
538. | 68 Education 6 358 6 0 |
539. | 68 Education 13 367 13 0 |
540. | 68 Education 19 358 19 0 |
|--------------------------------------------------|
541. | 68 Education 25 368 25 0 |
542. | 68 Education 33 389 33 2 |
543. | 68 Education 40 424 40 9 |
544. | 68 Education 46 419 46 15 |
+--------------------------------------------------+Once the knots have been created, we want to create a variable that partitions time periods to before and after the implementation of the intervention. Since there are only two time periods (before the intervention and after the intervention), we will create a binary indicator.
// Create a variable to separate the treatment periods.
gen period = .
replace period = 0 if obsday < 31 /* Baseline phase */
replace period = 1 if obsday >= 31 & !missing(obsday) /* Treatment phase */
tab period, m Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
(600 missing values generated)
(378 real changes made)
(222 real changes made)
period | Freq. Percent Cum.
------------+-----------------------------------
0 | 378 63.00 63.00
1 | 222 37.00 100.00
------------+-----------------------------------
Total | 600 100.00Let’s re-examine our data. You should notice that a new variable
called period has been added:
list if id == 12 | id == 45 | id == 68Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
+-----------------------------------------------------------+
| id group obsday sleep knot1 knot2 period |
|-----------------------------------------------------------|
89. | 12 Control 1 360 1 0 0 |
90. | 12 Control 8 356 8 0 0 |
91. | 12 Control 14 342 14 0 0 |
92. | 12 Control 20 360 20 0 0 |
93. | 12 Control 26 366 26 0 0 |
|-----------------------------------------------------------|
94. | 12 Control 32 350 32 1 1 |
95. | 12 Control 38 383 38 7 1 |
96. | 12 Control 45 385 45 14 1 |
353. | 45 Medication 1 332 1 0 0 |
354. | 45 Medication 7 338 7 0 0 |
|-----------------------------------------------------------|
355. | 45 Medication 13 349 13 0 0 |
356. | 45 Medication 18 357 18 0 0 |
357. | 45 Medication 26 370 26 0 0 |
358. | 45 Medication 32 415 32 1 1 |
359. | 45 Medication 37 402 37 6 1 |
|-----------------------------------------------------------|
360. | 45 Medication 44 413 44 13 1 |
537. | 68 Education 1 360 1 0 0 |
538. | 68 Education 6 358 6 0 0 |
539. | 68 Education 13 367 13 0 0 |
540. | 68 Education 19 358 19 0 0 |
|-----------------------------------------------------------|
541. | 68 Education 25 368 25 0 0 |
542. | 68 Education 33 389 33 2 1 |
543. | 68 Education 40 424 40 9 1 |
544. | 68 Education 46 419 46 15 1 |
+-----------------------------------------------------------+Linear mixed effects model - Using splines
Now that we have our necessary variables, we can begin to construct our ITSA model. We will use the linear mixed effects model for this ITSA. You can use other models such as the generalized estimating equation (GEE) model or a simple linear regression. I chose to use the linear mixed effects model because this was the one that Mitchell used in his textbook, but also because I was working on another study that used a similar model using triple interaction terms, and I wanted to see if I can replicate this result. (Note: In Method 2, we will see that we generate the same answers.)
After the model completes, it is a good idea to use the
predict command to estimate the outcomes for each subject,
which will generate a new column called linear1. We will
use this to plot our results.
// Linear mixed effects regression model
xtmixed sleep i.group i.period c.knot1 c.knot2 i.group#i.period i.group#c.knot1 i.group#c.knot2 || id: period knot1 knot2, covariance(unstructured)
predict linear1 /* Predict the linear estimates */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Performing EM optimization:
Performing gradient-based optimization:
Iteration 0: Log likelihood = -2427.9356
Iteration 1: Log likelihood = -2427.8869
Iteration 2: Log likelihood = -2427.8865
Iteration 3: Log likelihood = -2427.8865
Computing standard errors:
Mixed-effects ML regression Number of obs = 600
Group variable: id Number of groups = 75
Obs per group:
min = 8
avg = 8.0
max = 8
Wald chi2(11) = 194.53
Log likelihood = -2427.8865 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
sleep | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
group |
Medication | -2.072872 3.411818 -0.61 0.543 -8.759912 4.614168
Education | 16.31454 3.410371 4.78 0.000 9.630336 22.99875
|
1.period | -4.203853 2.793214 -1.51 0.132 -9.678451 1.270745
knot1 | -.0415024 .2998389 -0.14 0.890 -.6291758 .5461711
knot2 | -.0966276 .5077213 -0.19 0.849 -1.091743 .8984879
|
group#period |
Medication #|
1 | 34.28463 3.997248 8.58 0.000 26.45016 42.11909
Education#1 | 1.647258 3.973666 0.41 0.678 -6.140983 9.4355
|
group#|
c.knot1 |
Medication | .4016104 .4238518 0.95 0.343 -.4291239 1.232345
Education | -.1437658 .4240338 -0.34 0.735 -.9748567 .6873251
|
group#|
c.knot2 |
Medication | -.2913187 .7196931 -0.40 0.686 -1.701891 1.119254
Education | 2.478031 .7176651 3.45 0.001 1.071433 3.884629
|
_cons | 351.3643 2.414271 145.54 0.000 346.6324 356.0962
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
id: Unstructured |
sd(period) | 7.622454 2.261546 4.261351 13.6346
sd(knot1) | 1.449165 .1262265 1.221732 1.718937
sd(knot2) | 2.352807 .2239956 1.95231 2.835461
sd(_cons) | 10.19626 1.173629 8.136996 12.77666
corr(period,knot1) | -.1047645 .2214208 -.495967 .3217899
corr(period,knot2) | .1578074 .2444491 -.3205053 .5719962
corr(period,_cons) | -.0429446 .2586461 -.501153 .4340715
corr(knot1,knot2) | -.6348598 .079245 -.7656375 -.4536861
corr(knot1,_cons) | -.1454689 .1370307 -.3976797 .1271808
corr(knot2,_cons) | .1491897 .1433531 -.1361996 .411715
-----------------------------+------------------------------------------------
sd(Residual) | 7.89317 .3226308 7.285491 8.551534
------------------------------------------------------------------------------
LR test vs. linear model: chi2(10) = 1169.63 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.
(option xb assumed)We can plot the results of the linear mixed effects model. We can see that the Control trend before the intervention is slightly negative with a small drop immediately after the intervention, and a similar negative trend after the intervention. The Medication trend is increasing in the period before the intervention, it has a large jump in the level immediately after the intervention, but the slope is neutral after the intervention. The Education tend is slightly negative in the period before the intervention, there is a very small drop immediately after the intervention, and a steep increase in the slope after the intervention.
*** Plot linear estimates
graph twoway (lfit linear1 obsday if group == 1 & period == 0, lcol(navy) ) ///
(lfit linear1 obsday if group == 1 & period == 1, lcol(navy)) ///
(lfit linear1 obsday if group == 2 & period == 0, lcol(green) lpattern(dash)) ///
(lfit linear1 obsday if group == 2 & period == 1, lcol(green) lpattern(dash)) ///
(lfit linear1 obsday if group == 3 & period == 0, lcol(cranberry) lpattern(shortdash)) ///
(lfit linear1 obsday if group == 3 & period == 1, lcol(cranberry) lpattern(shortdash) ///
xline(31) ///
legend(title("Treatment groups", size(small)) size(vsmall) position(3) ///
order(1 "Control" 3 "Medication" 5 "Education")))Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Slopes before the intervention
We can estimate the marginal effects at different points along the ITSA.
To get the slopes for the three groups at the period before the intervention, we can write the following code:
*** Estimate the slopes for the three groups in the period == 0 (baseline phase)
margins, dydx(knot1) over(group)Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Average marginal effects Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: knot1
Over: group
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
knot1 |
group |
Control | -.0415024 .2998389 -0.14 0.890 -.6291758 .5461711
Medication | .360108 .299578 1.20 0.229 -.2270542 .9472702
Education | -.1852682 .2998354 -0.62 0.537 -.7729347 .4023984
------------------------------------------------------------------------------The slope for the Control group was -0.04 minutes per week (95% CI: -0.63, 0.55) before the intervention.
The slope for the Medication group was 0.36 minutes per week (95% CI: -0.23, 0.95) before the intervention.
The slope for the Education group was -0.19 minutes per week (95% CI: -0.77, 0.40) before the intervention.
None of these slopes were significant.
To test the difference in the slopes between the groups, we can add
the pwcompare(effects) option:
*** Estimate the slopes for the three groups in the period == 0 (baseline phase)
margins, dydx(knot1) over(group) pwcompareRunning C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Pairwise comparisons of average marginal effects
Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: knot1
Over: group
--------------------------------------------------------------------------
| Contrast Delta-method Unadjusted
| dy/dx std. err. [95% conf. interval]
-------------------------+------------------------------------------------
knot1 |
group |
Medication vs Control | .4016104 .4238518 -.4291239 1.232345
Education vs Control | -.1437658 .4240338 -.9748567 .6873251
Education vs Medication | -.5453762 .4238493 -1.376106 .2853533
--------------------------------------------------------------------------The P-values are all above 0.05; hence, there were no statistically significant differences in the slopes between the groups prior to the intervention.
Level change immediately after the intervention
Next, we can compare the level change immediately after the
implementation of the intervention. This analysis provides use with
information on how much the levels jumped or dropped immediately after
the implementation of the intervention at obsday = 31.
There are two ways to write this code (Note: I prefer the
margins command):
*** Estimate the amount of level change for each group
/* Method 1: using the contrast command */
contrast period@group, pveffects nowald
/* Method 2: using the margins command */
margins, dydx(period) over(group)Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Contrasts of marginal linear predictions
Margins: asbalanced
----------------------------------------------------------------
| Contrast Std. err. z P>|z|
------------------------+---------------------------------------
sleep |
period@group |
(1 vs base) Control | -4.203853 2.793214 -1.51 0.132
(1 vs base) Medication | 30.08077 2.859362 10.52 0.000
(1 vs base) Education | -2.556595 2.826301 -0.90 0.366
----------------------------------------------------------------
Average marginal effects Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: 1.period
Over: group
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
0.period | (base outcome)
-------------+----------------------------------------------------------------
1.period |
group |
Control | -4.203853 2.793214 -1.51 0.132 -9.678451 1.270745
Medication | 30.08077 2.859362 10.52 0.000 24.47653 35.68502
Education | -2.556595 2.826301 -0.90 0.366 -8.096042 2.982853
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.In the Control group, the level changed was -4.20 minutes (95% CI:
-9.68, 1.27) immediately after implementation of the intervention;
however, this was not statistically significant (P=0.132). We visualize
this in the Figure above by the decrease in the level at
obsday = 31.
In the Medication group, the level change was +30.08 minutes (95% CI:
24.48, 35.69) immediately after implementation of the intervention; this
was statistically significant (P<0.001). We visualize this in the
Figure above by the large increase in the level at
obsday = 31.
In the Education group, the level change was -2.56 minutes (95% CI:
-8.10, 2.98) immediately after implementation of the intervention;
however, this was not statistically significant (P=0.366). We visualize
this in the Figure above by the very small decrease in the level at
obsday = 31.
Comparing the Slopes after the intervention
Lastly, we can measure the differences in the slopes after and before the intervention between the groups, which is also called the difference-in-differences. This will yield information about the slope changes between the groups and the differences in the slopes after and before the intervention.
Let’s see what the differences in the slopes after and before the intervention for each group without comparing them to each other:
*** Estimate the change in slopes for each group
margins, dydx(knot2) over(group)Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Average marginal effects Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: knot2
Over: group
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
knot2 |
group |
Control | -.0966276 .5077213 -0.19 0.849 -1.091743 .8984879
Medication | -.3879464 .5100757 -0.76 0.447 -1.387676 .6117837
Education | 2.381403 .5072103 4.70 0.000 1.38729 3.375517
------------------------------------------------------------------------------In the Control group, the slope in the period after the intervention was lower than before the intervention by -0.10 minutes per week (95% CI: -1.09, 0.90); however, this was not statistically significant (P=0.849).
In the Medication group, the slope in the period after the intervention was lower than before the intervention by -0.39 minutes per week (95% CI: -1.39, 0.61); however, this was not statistically significant (P=0.447).
In the Eduation group, the slope in the period after the intervention was greater than before the intervention by 2.38 minute per week (95% CI: 1.39, 3.38); this was statistically significant (P<0.001).
We can compare these difference to each other using the
margins command with the pwcompare(effectS)
option:
*** Estimate the change in slopes for each group
margins, dydx(knot2) over(group) pwcompare(effects)Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Pairwise comparisons of average marginal effects
Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: knot2
Over: group
------------------------------------------------------------------------------
| Contrast Delta-method Unadjusted Unadjusted
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
knot2 |
group |
Medication |
vs |
Control | -.2913187 .7196931 -0.40 0.686 -1.701891 1.119254
Education |
vs |
Control | 2.478031 .7176651 3.45 0.001 1.071433 3.884629
Education |
vs |
Medication | 2.76935 .7193327 3.85 0.000 1.359484 4.179216
------------------------------------------------------------------------------In the comparison between the Medication and Control groups, the difference in the slope change was -0.29 minute per week (95% CI: -1.70, 1.12); however, this was not statistically signficiant (P=0.686).
In the comaprison between the Eduation and Control groups, the difference in the slope change was 2.48 minutes per week (95% CI: 1.07, 3.88); this was statistically significant (P=0.001).
In the comparison between the Education versus and Medication groups, the difference in teh slope change was 2.77 mintues per week (95% CI: 1.36, 4.18); this was statistically significant (P<0.001).
Method 2: Triple Interaction Term
The mkspline method makes it easy to perform an ITSA in
Stata. However, we can do a lot of the
mkspline work ourselves using interaction terms,
particularly the triple interaction between the group,
obsday, and period variables.
clear all
cd "C:\Users\carn3_a5ada\Dropbox\STATA references\Mitchels book - Visualing in Stata\ivrm2"
use sleep_cat3pw.dta, replaceRunning C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
C:\Users\carn3_a5ada\Dropbox\STATA references\Mitchels book - Visualing in Stat
> a\ivrm2Recall the following formula:
\(Y_i = \beta_0 + \beta_1(T_i) + \beta_2(X_i) + \beta_3(T_i * X_i) + \beta_4(G_i) + \beta_5(G_i * T_i) + \beta_6(G_i * X_i) + \beta_7(G_i * T_i * X_i) + \epsilon_i\)
We will use this to our advantage to create the \(\beta_7\) coefficient in our model.
First, we’ll need to generate the binary period variable
that partitions the time period before and after the intervention.
// Generate a variable for the pre- and post-index date (obsday == 31)
gen period = .
replace period = 0 if obsday < 31 /* Baseline phase */
replace period = 1 if obsday >= 31 & !missing(obsday) /* Treatment phase */
tab period, m Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
(600 missing values generated)
(378 real changes made)
(222 real changes made)
period | Freq. Percent Cum.
------------+-----------------------------------
0 | 378 63.00 63.00
1 | 222 37.00 100.00
------------+-----------------------------------
Total | 600 100.00Second, we will need to create lots of interaction terms. Pay attention to how the terms are created. There are a total of three two-way interaction terms and one three-way interaction term (triple interaction).
After the model completes, it is a good idea to use the
predict command to estimate the outcomes for each subject,
which will generate a new column called linear2. We will
use this to plot our results.
// Linear mixed effects regression model
xtmixed sleep i.group c.obsday i.period i.group#c.obsday i.group#i.period c.obsday#i.period i.group#c.obsday#i.period || id: i.period obsday c.obsday#i.period, covariance(unstructured)
predict linear2 /* Predict the linear estimates */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Performing EM optimization:
Performing gradient-based optimization:
Iteration 0: Log likelihood = -2440.1403
Iteration 1: Log likelihood = -2429.4305
Iteration 2: Log likelihood = -2428.2399
Iteration 3: Log likelihood = -2427.8876
Iteration 4: Log likelihood = -2427.8865
Iteration 5: Log likelihood = -2427.8865
Computing standard errors:
Mixed-effects ML regression Number of obs = 600
Group variable: id Number of groups = 75
Obs per group:
min = 8
avg = 8.0
max = 8
Wald chi2(11) = 194.53
Log likelihood = -2427.8865 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
sleep | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
group |
Medication | -2.072872 3.41182 -0.61 0.543 -8.759916 4.614172
Education | 16.31454 3.410373 4.78 0.000 9.630332 22.99875
|
obsday | -.0415024 .2998387 -0.14 0.890 -.6291755 .5461707
1.period | -1.208401 16.09127 -0.08 0.940 -32.7467 30.3299
|
group#|
c.obsday |
Medication | .4016104 .4238516 0.95 0.343 -.4291235 1.232344
Education | -.1437658 .4240335 -0.34 0.735 -.9748562 .6873247
|
group#period |
Medication #|
1 | 43.31551 22.86307 1.89 0.058 -1.495271 88.1263
Education#1 | -75.1717 22.75371 -3.30 0.001 -119.7682 -30.57525
|
period#|
c.obsday |
1 | -.0966275 .5077213 -0.19 0.849 -1.091743 .8984879
|
group#|
period#|
c.obsday |
Medication #|
1 | -.2913189 .7196931 -0.40 0.686 -1.701891 1.119254
Education#1 | 2.478031 .7176651 3.45 0.001 1.071433 3.884629
|
_cons | 351.3643 2.414273 145.54 0.000 346.6324 356.0962
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
id: Unstructured |
sd(1.period) | 72.12794 7.416293 58.9633 88.23181
sd(obsday) | 1.449164 .1262334 1.221719 1.718952
sd(1.period#obsday) | 2.352807 .224007 1.952292 2.835488
sd(_cons) | 10.19626 1.173635 8.136995 12.77668
corr(1.period,obsday) | .6309106 .0869351 .4299745 .7722948
corr(1.period, |
1.period# |
c.obsday) | -.9945401 .003436 -.9984118 -.9813177
corr(1.period,_cons) | -.1554014 .1485937 -.4260955 .1408325
corr(obsday,1.period#obsday) | -.6348605 .0792723 -.7656752 -.4536155
corr(obsday,_cons) | -.1454687 .1370455 -.3977045 .1272102
corr(1.period#obsday,_cons) | .1491895 .1434002 -.1362924 .4117932
-----------------------------+------------------------------------------------
sd(Residual) | 7.893171 .3226315 7.285491 8.551537
------------------------------------------------------------------------------
LR test vs. linear model: chi2(10) = 1169.63 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.
(option xb assumed)Once the model run is complete, we can plot the predicted values.
Notice that these plots are exactly the same as the previous plot based
on the mkspline method.
*** Plot linear estimates
graph twoway (lfit linear2 obsday if group == 1 & period == 0, lcol(navy) ) ///
(lfit linear2 obsday if group == 1 & period == 1, lcol(navy)) ///
(lfit linear2 obsday if group == 2 & period == 0, lcol(green) lpattern(dash)) ///
(lfit linear2 obsday if group == 2 & period == 1, lcol(green) lpattern(dash)) ///
(lfit linear2 obsday if group == 3 & period == 0, lcol(cranberry) lpattern(shortdash)) ///
(lfit linear2 obsday if group == 3 & period == 1, lcol(cranberry) lpattern(shortdash) ///
xline(31) ///
legend(title("Treatment groups", size(small)) size(vsmall) position(3) ///
order(1 "Control" 3 "Medication" 5 "Education")))Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Slopes before the intervention
We can estimate the marginal effects at different points along the ITSA.
There are two methods to estimate the slopes for the three groups at the period before the intervention:
*** Estimate the slopes for the three groups in the period == 0 (baseline phase)
margins if period == 0, dydx(obsday) over(group) /* Method 1 */
margins group, dydx(obsday) at(period = (0 1)) /* Method 2 */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Average marginal effects Number of obs = 378
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
Over: group
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday |
group |
Control | -.0415024 .2998387 -0.14 0.890 -.6291755 .5461707
Medication | .360108 .2995779 1.20 0.229 -.2270538 .9472699
Education | -.1852682 .2998352 -0.62 0.537 -.7729344 .4023981
------------------------------------------------------------------------------
Average marginal effects Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
1._at: period = 0
2._at: period = 1
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday |
_at#group |
1#Control | -.0415024 .2998387 -0.14 0.890 -.6291755 .5461707
1 #|
Medication | .360108 .2995779 1.20 0.229 -.2270538 .9472699
1#Education | -.1852682 .2998352 -0.62 0.537 -.7729344 .4023981
2#Control | -.1381299 .4038305 -0.34 0.732 -.9296232 .6533634
2 #|
Medication | -.0278384 .4069828 -0.07 0.945 -.82551 .7698332
2#Education | 2.196135 .4031956 5.45 0.000 1.405886 2.986384
------------------------------------------------------------------------------Compare these results to the mkspline method. They are
exactly the same.
We can also compare the slopes between the groups using the
pwcompare(effects) option.
margins if period == 0, dydx(obsday) over(group) pwcompare(effects) /* Method 1 */
margins group, dydx(obsday) at(period = 0) pwcompare(effects) /* Method 2 */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Pairwise comparisons of average marginal effects
Number of obs = 378
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
Over: group
------------------------------------------------------------------------------
| Contrast Delta-method Unadjusted Unadjusted
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday |
group |
Medication |
vs |
Control | .4016104 .4238516 0.95 0.343 -.4291235 1.232344
Education |
vs |
Control | -.1437658 .4240335 -0.34 0.735 -.9748562 .6873247
Education |
vs |
Medication | -.5453762 .4238491 -1.29 0.198 -1.376105 .2853528
------------------------------------------------------------------------------
Pairwise comparisons of average marginal effects
Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
At: period = 0
------------------------------------------------------------------------------
| Contrast Delta-method Unadjusted Unadjusted
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday |
group |
Medication |
vs |
Control | .4016104 .4238516 0.95 0.343 -.4291235 1.232344
Education |
vs |
Control | -.1437658 .4240335 -0.34 0.735 -.9748562 .6873247
Education |
vs |
Medication | -.5453762 .4238491 -1.29 0.198 -1.376105 .2853528
------------------------------------------------------------------------------Once again, these results are exactly the same as the
mkspline method.
Level change immediately after the intervention
We can perform this analysis using the margins commnand
in Stata.
margins if obsday == 31, dydx(period) over(group)Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Average marginal effects Number of obs = 7
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: 1.period
Over: group
------------------------------------------------------------------------------
| Delta-method
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
0.period | (base outcome)
-------------+----------------------------------------------------------------
1.period |
group |
Control | -4.203854 2.793213 -1.51 0.132 -9.67845 1.270742
Medication | 30.08077 2.859361 10.52 0.000 24.47653 35.68502
Education | -2.556595 2.8263 -0.90 0.366 -8.096041 2.982851
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.When compared to the mkspline method, these results are
exactly the same.
Comparing the Slopes after the intervention
Lastly, we can measure the differences in the slopes after and before the intervention between the groups, which is also called the difference-in-differences. This will yield information about the slope changes between the groups and the differences in the slopes after and before the intervention.
There is no easy way to do this like the mkspline
approach due to the numerous interactions that are occurring. Rather
than give you the comparisons of the changes in the slopes between the
groups, the margins command will give you multiple
combinations of comparisons. For example, this will yield comparisons
between the Medication group in the period before the intervention with
the Control group in the period before the intervention, which we
already acquired in the previous analysis. However, the
difference-in-differences that we are looking for is in the output, we
just need to search for it.
There are three ways to evaluate the triple interaction:
*** Estimate the change in slopes for each group
margins group, dydx(obsday) at(period = (0 1)) pwcompare(effects) /* Method 1 */
margins, dydx(obsday) over(group period) pwcompare(effects) /* Method 2 */
margins group, dydx(obsday) over(period) pwcompare(effects) /* Method 3 */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Pairwise comparisons of average marginal effects
Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
1._at: period = 0
2._at: period = 1
------------------------------------------------------------------------------
| Contrast Delta-method Unadjusted Unadjusted
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday |
_at#group |
(1 #|
Medication) |
vs |
(1 #|
Control) | .4016104 .4238516 0.95 0.343 -.4291235 1.232344
(1 #|
Education) |
vs |
(1 #|
Control) | -.1437658 .4240335 -0.34 0.735 -.9748562 .6873247
(2#Control) |
vs |
(1#Control) | -.0966275 .5077213 -0.19 0.849 -1.091743 .8984879
(2 #|
Medication) |
vs |
(1 #|
Control) | .013664 .5055079 0.03 0.978 -.9771133 1.004441
(2 #|
Education) |
vs |
(1 #|
Control) | 2.237638 .5024639 4.45 0.000 1.252826 3.222449
(1 #|
Education) |
vs |
(1 #|
Medication) | -.5453762 .4238491 -1.29 0.198 -1.376105 .2853528
(2 #|
Control) |
vs |
(1 #|
Medication) | -.4982379 .5028181 -0.99 0.322 -1.483743 .4872674
(2 #|
Medication) |
vs |
(1 #|
Medication) | -.3879464 .5100757 -0.76 0.447 -1.387676 .6117836
(2 #|
Education) |
vs |
(1 #|
Medication) | 1.836027 .5023083 3.66 0.000 .8515211 2.820533
(2 #|
Control) |
vs |
(1 #|
Education) | .0471382 .5029714 0.09 0.925 -.9386677 1.032944
(2 #|
Medication) |
vs |
(1 #|
Education) | .1574297 .5055058 0.31 0.755 -.8333435 1.148203
(2 #|
Education) |
vs |
(1 #|
Education) | 2.381403 .5072102 4.70 0.000 1.38729 3.375517
(2 #|
Medication) |
vs |
(2 #|
Control) | .1102915 .5733359 0.19 0.847 -1.013426 1.234009
(2 #|
Education) |
vs |
(2 #|
Control) | 2.334265 .5706538 4.09 0.000 1.215804 3.452726
(2 #|
Education) |
vs |
(2 #|
Medication) | 2.223974 .5728889 3.88 0.000 1.101132 3.346815
------------------------------------------------------------------------------
Pairwise comparisons of average marginal effects
Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
Over: group period
------------------------------------------------------------------------------
| Contrast Delta-method Unadjusted Unadjusted
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday |
group#period |
(Control#1) |
vs |
(Control#0) | -.0966275 .5077213 -0.19 0.849 -1.091743 .8984879
(Medication #|
0) |
vs |
(Control #|
0) | .4016104 .4238516 0.95 0.343 -.4291235 1.232344
(Medication #|
1) |
vs |
(Control #|
0) | .013664 .5055079 0.03 0.978 -.9771133 1.004441
(Education #|
0) |
vs |
(Control #|
0) | -.1437658 .4240335 -0.34 0.735 -.9748562 .6873247
(Education #|
1) |
vs |
(Control #|
0) | 2.237638 .5024639 4.45 0.000 1.252826 3.222449
(Medication #|
0) |
vs |
(Control #|
1) | .4982379 .5028181 0.99 0.322 -.4872674 1.483743
(Medication #|
1) |
vs |
(Control #|
1) | .1102915 .5733359 0.19 0.847 -1.013426 1.234009
(Education #|
0) |
vs |
(Control #|
1) | -.0471382 .5029714 -0.09 0.925 -1.032944 .9386677
(Education #|
1) |
vs |
(Control #|
1) | 2.334265 .5706538 4.09 0.000 1.215804 3.452726
(Medication #|
1) |
vs |
(Medication #|
0) | -.3879464 .5100757 -0.76 0.447 -1.387676 .6117836
(Education #|
0) |
vs |
(Medication #|
0) | -.5453762 .4238491 -1.29 0.198 -1.376105 .2853528
(Education #|
1) |
vs |
(Medication #|
0) | 1.836027 .5023083 3.66 0.000 .8515211 2.820533
(Education #|
0) |
vs |
(Medication #|
1) | -.1574297 .5055058 -0.31 0.755 -1.148203 .8333435
(Education #|
1) |
vs |
(Medication #|
1) | 2.223974 .5728889 3.88 0.000 1.101132 3.346815
(Education #|
1) |
vs |
(Education #|
0) | 2.381403 .5072102 4.70 0.000 1.38729 3.375517
------------------------------------------------------------------------------
Pairwise comparisons of average marginal effects
Number of obs = 600
Expression: Linear prediction, fixed portion, predict()
dy/dx wrt: obsday
Over: period
------------------------------------------------------------------------------
| Contrast Delta-method Unadjusted Unadjusted
| dy/dx std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
obsday |
period#group |
(0 #|
Medication) |
vs |
(0 #|
Control) | .4016104 .4238516 0.95 0.343 -.4291235 1.232344
(0 #|
Education) |
vs |
(0 #|
Control) | -.1437658 .4240335 -0.34 0.735 -.9748562 .6873247
(1#Control) |
vs |
(0#Control) | -.0966275 .5077213 -0.19 0.849 -1.091743 .8984879
(1 #|
Medication) |
vs |
(0 #|
Control) | .013664 .5055079 0.03 0.978 -.9771133 1.004441
(1 #|
Education) |
vs |
(0 #|
Control) | 2.237638 .5024639 4.45 0.000 1.252826 3.222449
(0 #|
Education) |
vs |
(0 #|
Medication) | -.5453762 .4238491 -1.29 0.198 -1.376105 .2853528
(1 #|
Control) |
vs |
(0 #|
Medication) | -.4982379 .5028181 -0.99 0.322 -1.483743 .4872674
(1 #|
Medication) |
vs |
(0 #|
Medication) | -.3879464 .5100757 -0.76 0.447 -1.387676 .6117836
(1 #|
Education) |
vs |
(0 #|
Medication) | 1.836027 .5023083 3.66 0.000 .8515211 2.820533
(1 #|
Control) |
vs |
(0 #|
Education) | .0471382 .5029714 0.09 0.925 -.9386677 1.032944
(1 #|
Medication) |
vs |
(0 #|
Education) | .1574297 .5055058 0.31 0.755 -.8333435 1.148203
(1 #|
Education) |
vs |
(0 #|
Education) | 2.381403 .5072102 4.70 0.000 1.38729 3.375517
(1 #|
Medication) |
vs |
(1 #|
Control) | .1102915 .5733359 0.19 0.847 -1.013426 1.234009
(1 #|
Education) |
vs |
(1 #|
Control) | 2.334265 .5706538 4.09 0.000 1.215804 3.452726
(1 #|
Education) |
vs |
(1 #|
Medication) | 2.223974 .5728889 3.88 0.000 1.101132 3.346815
------------------------------------------------------------------------------The output yields lots of comparisons, but the ones we want are:
(1#Control) vs (0#Control), which yields the differences
in the slopes after and before the intervention for the Control
group.
(1#Medication) vs (0#Medication), which yields the
differences in the slopes after and before the intervention for the
Medication group.
(1#Education) vs (0#Education), which yields the
differences in the slopes after and before the intervention for the
Education group.
When compared to the mkspline results, these are exactly
the same.
Finally, to get the difference-in-differences estimates, we need to
use the constrast command in Stata. I could
not get the contrast command to give me all the pair-wise
comparisons between the three groups, so I had to write another
contrast command using the ar.factor operator,
which means that the difference will be made with the previous level; in
our scenario, this will make Medication as the reference for the
Education comparison.
contrast i.group#c.obsday#i.period, effects /* difference-in-differences */
contrast ar.group#c.obsday#i.period, effects /* difference-in-differences */Running C:\Users\carn3_a5ada\Dropbox\Marks blog\ITSA\R Markdown code\profile.do
> ...
Contrasts of marginal linear predictions
Margins: asbalanced
---------------------------------------------------------
| df chi2 P>chi2
----------------------+----------------------------------
sleep |
group#period#c.obsday | 2 17.95 0.0001
---------------------------------------------------------
------------------------------------------------------------------------------
| Contrast Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
sleep |
group#|
period#|
c.obsday |
(Medication |
vs |
base) |
(1 vs base) | -.2913189 .7196931 -0.40 0.686 -1.701891 1.119254
(Education |
vs |
base) |
(1 vs base) | 2.478031 .7176651 3.45 0.001 1.071433 3.884629
------------------------------------------------------------------------------
Contrasts of marginal linear predictions
Margins: asbalanced
----------------------------------------------------------------------
| df chi2 P>chi2
-----------------------------------+----------------------------------
sleep |
group#period#c.obsday |
(Medication vs Control) (joint) | 1 0.16 0.6856
(Education vs Medication) (joint) | 1 14.82 0.0001
Joint | 2 17.95 0.0001
----------------------------------------------------------------------
------------------------------------------------------------------------------
| Contrast Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
sleep |
group#|
period#|
c.obsday |
(Medication |
vs |
Control) |
(1 vs base) | -.2913189 .7196931 -0.40 0.686 -1.701891 1.119254
(Education |
vs |
Medication) |
(1 vs base) | 2.76935 .7193326 3.85 0.000 1.359484 4.179216
------------------------------------------------------------------------------Based on these findings, we can see that difference in the changes in the slopes between Medication and Control is -0.29 minutes per week (95% CI: -1.70, 1.12), which is not statistically significant (P=0.686).
However, there was a greater difference in the changes in the slope between the Education versus Control, which was +2.47 minutes per week 995% CI: 1.07, 3.89); this is statistically significant (P=0.001).
Similarly, the difference in the changes in the slope was greater for the Medication group versus the Education group by +2.77 minutes per week (95% CI: 1.36, 4.18), which was statistically significant (P<0.001).
Once again, these results are exactly the same as the
mkspline approach.
Conclusions
ITSA is a power tool that can help us quantify the trends before and
after the implementation of an intervention. This can be used for a
single group design or a multiple groups design. The ease of the
mkspline methods facilitates our ability to create and
interpret the linear mixed effects model. However, the triple
interaction term method is also viable despite requiring a little bit of
additional coding with the margins and
constrast commands.
Acknowledgements
I have been exposed to many different approaches to performing ITSA, but the two that I found to be the most helpful were by Ariel Linden and Michael N. Mitchell.
Linden’s paper, “Conducting interrupted time-series analysis for single- and multiple-group comparisons,” is available here. This was a great paper on the theory and principles of ITSA; and it was written in a way that was easy to absorb and understand. The citation is Linden A. Conducting interrupted time-series analysis for single- and multiple-group comparisons. The Stata Journal. 2015;15(2):480-500.
Micheal N. Mitchell’s book, “Interpreting and Visualizing Regression
Models Using Stata, Second Edition,” is available on the Stata
website. I recommend reading Chapters 12 and 16; I found these to be
the most helpful in writing code for ITSA.
Another useful guide is the Stata documentation on
contrast options, which is available here.
My code for this exercise is available on my GitHub site.
Work in progress
This is a work in progress, and I plan to make updates in the future. So stay tuned.
Lastly, this is for education purposes only.