Interrupted Time Series

[jrmorse]

Hi there! I'm trying to do a simple interrupted time series and I am running into a lot of trouble as I have never done this before.

I want to see if renewals of specific medications decreased in March 2020 by comparing what the rate of growth had been at a particular community clinic between 2018-March2020 to what the clinic experienced between March-September 2020.

From my research, I think the itsa.model function should be used from the its.analysis library. Perhaps I should be using a different approach?

Using itsa.model, the independent variable should be any date after March 2020 (this was identified through the creation of a new column and assigning "1" to the values after March and 0 before-- "renewals_itsa") and the dependent variable is the sum of prescription refills-- "renewals". Time is based off prescription fill year and month-- "start_month_year".

Here's what I have so far, but it's not working. Any advice or pointers of if I am on the right track would be very helpful!

itsa.model(data = renewals_itsa, time = "start_month_year", dependent_variable = "renewals",
  interrupt_var = "inter_val", alpha=0.05, bootstrap=TRUE, Reps = 250)
#> Error in itsa.model(data = renewals_itsa, time = "start_month_year", dependent_variable = "renewals", : could not find function "itsa.model"

^{Created on 2020-12-26 by the reprex package (v0.3.0)}

[technocrat]

@jrmorse

could not find function "itsa.model"

just needs

library(its.analysis)

I'm not entirely sure that this is the right framework, however, because medication renewals are typically time dependent (every 30 days for refills, say or some multiple of 30 days for new prescriptions to continue existing medications). That leads to auto-correlation with a 30-day lag and also messes with monthly data, collected at an average scale of 30.41667 or 30.5 days, depending on leap years. For any given medication, the results are discrete—37 or 38 renewals, not 37.43 renewals. That suggests a Poisson distribution but for the fact that successive months are more than weakly dependent.

[jrmorse]

@technocrat

Poisson distribution

Thank you for your response! I copied the reprex incorrectly-- I do have the its.analysis loaded into the rmd, the problem seems to be with the dependent variable. The reprex should read as follows:

library(its.analysis)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo

itsa.model(data = "renewals_itsa", time = "start_month_year", depvar = "renewals",
  interrupt_var = "inter_val", alpha=0.05, bootstrap=FALSE, Reps = 250)
#> Error in data[, depvar]: incorrect number of dimensions

^{Created on 2020-12-28 by the reprex package (v0.3.0)}

Regarding the significance of the Poisson distribution... I'm not sure I completely understand, but because I have segmented the data into months (rather than daily), I think I should be addressing this concern. In other words, if the intervention (COVID lockdown) was implemented in March, then I would expect to see an impact on the data come April prescriptions. Do you think that there may be a better way to conduct an ITS on the data than what I have done above? I'm all ears if so

[technocrat]

This outline discusses the reasons that methods assuming continuous variables are unsuited for discrete data.

In assessing differences between March and April (H_0 = no difference, H_1 = difference) the task is to determine whether they merely due to random variation. In doing so, consideration must be given to seasonality and trend. The reprex below illustrates this with discretized monthly energy usage serving as a proxy for prescriptions.

suppressPackageStartupMessages({
  library(fpp3)
})

series <- tsibble(
  Month = rep(yearmonth("2018,1")) + 0:24,
  Scripts = c(
    7432, 6330, 7086, 6886, 7290, 6976,
    7334, 6539, 5941, 7034, 6814, 5509, 6530, 5126, 5694, 5572, 6027,
    5246, 5922, 5608, 6240, 7070, 6257, 6257, 6433
  )
)
#> Using `Month` as index variable.

series
#> # A tsibble: 25 x 2 [1M]
#>       Month Scripts
#>       <mth>   <dbl>
#>  1 2018 Jan    7432
#>  2 2018 Feb    6330
#>  3 2018 Mar    7086
#>  4 2018 Apr    6886
#>  5 2018 May    7290
#>  6 2018 Jun    6976
#>  7 2018 Jul    7334
#>  8 2018 Aug    6539
#>  9 2018 Sep    5941
#> 10 2018 Oct    7034
#> # … with 15 more rows


decomp <- stl(series, s.window = "periodic")
autoplot(decomp) + theme_minimal()
#> Error: Objects of type stl not supported by autoplot.

## NOTE: Error above is due to `reprex` environment; plot attached separately

^{Created on 2020-12-28 by the reprex package (v0.3.0.9001)}