Recursos para el análisis de series temporales interrumpidas en R

Question

Soy bastante nuevo en R. He intentado leer sobre el análisis de series temporales y ya he terminado

1. Shumway y Stoffer's Análisis de series temporales y sus aplicaciones 3ª Edición ,
2. Hyndman es excelente Previsión: principios y práctica
3. Avril Coghlan's Uso de R para el análisis de series temporales
4. A. Ian McLeod y otros Análisis de series temporales con R
5. El Dr. Marcel Dettling Análisis de Series Temporales Aplicadas

Edición: No estoy seguro de cómo manejar esto, pero he encontrado un recurso útil fuera de Cross Validated. Quería incluirlo aquí en caso de que alguien se tropiece con esta pregunta.

Análisis de regresión segmentada de estudios de series temporales interrumpidas en la investigación del uso de medicamentos

Tengo una serie temporal univariante del número de artículos consumidos (datos de recuento) medidos diariamente durante 7 años. Se aplicó una intervención a la población del estudio aproximadamente a la mitad de la serie temporal. No se espera que esta intervención produzca un efecto inmediato y el momento en el que se produce el efecto es esencialmente desconocido.

Usando la fórmula de Hyndman forecast he ajustado un modelo ARIMA a los datos anteriores a la intervención utilizando auto.arima() . Pero no estoy seguro de cómo utilizar este ajuste para responder si ha habido un cambio de tendencia estadísticamente significativo y cuantificar la cantidad.

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

¿Existen recursos que traten específicamente del análisis de series temporales interrumpidas en R? He encontrado este tratando con ITS en SPSS pero no he podido traducirlo a R.

Answer 1

Esto se conoce como análisis del punto de cambio. El paquete R changepoint puede hacer esto por usted: vea la documentación aquí (incluyendo referencias a la literatura): http://www.lancs.ac.uk/~killick/Pub/KillickEckley2011.pdf

Answer 2

Yo sugeriría un modelo jerárquico de medidas repetidas. Este método debería proporcionar resultados sólidos, ya que cada individuo actuará como su propio control. Pruebe a comprobar este enlace de la UCLA.

Answer 3

Para un enfoque bayesiano, puede utilizar mcp para ajustar un modelo Poisson o Binomial (porque tiene recuentos de períodos de intervalo fijo) con autoregresión aplicada a los residuos (en el espacio logarítmico). A continuación, compare un modelo de dos segmentos con un modelo de un segmento utilizando la validación cruzada.

Antes de empezar, tenga en cuenta que para este conjunto de datos, este modelo no se ajusta bien y la validación cruzada parece inestable. Así que me abstendría de usar lo siguiente en escenarios de alto riesgo, pero ilustra un enfoque general:

# 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 el presente conjunto de datos, esto da como resultado

       elpd_diff se_diff
model2    0.0       0.0 
model1 -459.1      64.3 

Es decir, un elpd_diff/se_diff de alrededor de 7 a favor del modelo nulo (sin cambios). Las posibles mejoras incluyen:

• modelización de la tendencia periódica mediante sin() o cos() .
• añadiendo información previa sobre la ubicación probable del cambio, por ejemplo prior = list(cp_1 = dnorm(1999.8, 0.5) .

Más información sobre la modelización de la autorregresión, la comparación de modelos y el establecimiento de prioridades el mcp sitio web . Divulgación: Soy el desarrollador de mcp .