Recursos para análise de séries temporais interrompidas em R

Eu sou bastante novo para R. Tentei ler sobre a análise de séries temporais e já terminei

1. Análise de séries temporais de Shumway e Stoffer e suas aplicações 3rd Edition ,
2. Excelente previsão de Hyndman : princípios e prática
3. Avril Coghlan está usando R para análise de séries temporais
4. A. Ian McLeod e cols. Análise de séries temporais com R
5. Análise Aplicada de Séries Temporais do Dr. Marcel Dettling

Edit: Não tenho certeza de como lidar com isso, mas encontrei um recurso útil fora do Cross Validated. Eu queria incluí-lo aqui, caso alguém se deparasse com essa pergunta.

Análise de regressão segmentada de estudos de séries temporais interrompidas em pesquisas sobre uso de medicamentos

Tenho uma série temporal univariada do número de itens consumidos (dados de contagem) medidos diariamente por 7 anos. Uma intervenção foi aplicada à população do estudo aproximadamente no meio da série temporal. Não se espera que esta intervenção produza um efeito imediato e o momento do início do efeito seja essencialmente incognoscível.

Usando o forecastpacote de Hyndman, adaptei um modelo ARIMA aos dados pré-intervenção usando auto.arima(). Mas não tenho certeza de como usar esse ajuste para responder se houve uma mudança estatisticamente significativa na tendência e quantificar o valor.

# for simplification I will aggregate to monthly counts
# I can later generalize any teachings the community supplies
count <- c(2464, 2683, 2426, 2258, 1950, 1548, 1108,  991, 1616, 1809, 1688, 2168, 2226, 2379, 2211, 1925, 1998, 1740, 1305,  924, 1487, 1792, 1485, 1701, 1962, 2896, 2862, 2051, 1776, 1358, 1110,  939, 1446, 1550, 1809, 2370, 2401, 2641, 2301, 1902, 2056, 1798, 1198,  994, 1507, 1604, 1761, 2080, 2069, 2279, 2290, 1758, 1850, 1598, 1032,  916, 1428, 1708, 2067, 2626, 2194, 2046, 1905, 1712, 1672, 1473, 1052,  874, 1358, 1694, 1875, 2220, 2141, 2129, 1920, 1595, 1445, 1308, 1039,  828, 1724, 2045, 1715, 1840)
# for explanatory purposes
# month <- rep(month.name, 7)
# year <- 1999:2005
ts <- ts(count, start(1999, 1))
train_month <- window(ts, start=c(1999,1), end = c(2001,1))
require(forecast)
arima_train <- auto.arima(train_month)
fit_month <- Arima(train_month, order = c(2,0,0), seasonal = c(1,1,0), lambda = 0)
plot(forecast(fit_month, 36)); lines(ts, col="red")

Existem recursos especificamente relacionados à análise de séries temporais interrompidas no R? Eu encontrei isso lidando com ITS no SPSS, mas não consegui traduzir isso para R.

Isso é conhecido como análise de ponto de mudança. O pacote R changepointpode fazer isso por você: consulte a documentação aqui (incluindo referências à literatura): http://www.lancs.ac.uk/~killick/Pub/KillickEckley2011.pdf

Eu sugeriria um modelo hierárquico de medidas repetidas. Esse método deve fornecer resultados robustos, pois cada indivíduo atuará como seu próprio controle. Tente verificar este link da UCLA.

Para uma abordagem bayesiana, você pode usar mcppara ajustar um modelo de Poisson ou Binomial (porque você tem contagens de períodos de intervalo fixo) com regressão automática aplicada aos resíduos (no espaço de log). Em seguida, compare um modelo de dois segmentos com um modelo de um segmento usando validação cruzada.

Antes de começarmos, observe que, para este conjunto de dados, este modelo não se encaixa bem e a validação cruzada parece instável. Portanto, evitaria usar o seguinte em cenários de alto risco, mas ilustra uma abordagem geral:

# Fit the change point model
library(mcp)
model_full = list(
  count ~ 1 + ar(1),  # intercept and AR(1)
  ~ 1  # New intercept
)
fit_full = mcp(model_full, data = df, family = poisson(), par_x = "year")


# Fit the null model
model_null = list(
  count ~ 1 + ar(1)  # just a stable AR(1)
)
fit_null = mcp(model_null, data = df, family = poisson(), par_x = "year")

# Compare predictive performance using LOO cross-validation
fit_full$loo = loo(fit_full)
fit_null$loo = loo(fit_null)
loo::loo_compare(fit_full$loo, fit_null$loo)

Para o presente conjunto de dados, isso resulta em

       elpd_diff se_diff
model2    0.0       0.0 
model1 -459.1      64.3 

Ou seja, uma elpd_diff/se_diffproporção de cerca de 7 a favor do modelo nulo (sem alteração). As possíveis melhorias incluem:

• modelando a tendência periódica usando sin()ou cos().
• adicionando informações anteriores sobre a provável localização da mudança, por exemplo prior = list(cp_1 = dnorm(1999.8, 0.5),.

Leia mais sobre modelagem de regressão automática, comparação de modelos e configurações anteriores ao mcpsite . Divulgação: Eu sou o desenvolvedor de mcp.