## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  message  = TRUE,
  warning  = FALSE
)

## ----load---------------------------------------------------------------------
library(idionomics)

set.seed(42)
panel <- do.call(rbind, lapply(1:9, function(id) {
  a <- rnorm(50)
  b <- 0.4 * a + rnorm(50)
  c <- 0.4 * b + rnorm(50)
  data.frame(id = as.character(id), time = seq_len(50),
             a = a, b = b, c = c,
             y = 0.5 * a + rnorm(50),
             stringsAsFactors = FALSE)
}))

# Subject 10: near-constant "a" — will be caught by i_screener(min_sd = 0.5)
s10 <- data.frame(id = "10", time = seq_len(50),
                  a = rep(3, 50), b = rnorm(50), c = rnorm(50),
                  y = rnorm(50), stringsAsFactors = FALSE)

# Subject 11: full positive loop (a -> b -> c -> a)
a11 <- rnorm(50)
b11 <- 0.6 * a11 + rnorm(50, sd = 0.5)
c11 <- 0.6 * b11 + rnorm(50, sd = 0.5)
s11 <- data.frame(id = "11", time = seq_len(50),
                  a = a11 + 0.4 * c11, b = b11, c = c11,
                  y = 0.5 * a11 + rnorm(50), stringsAsFactors = FALSE)

# Subject 12: negative a -> y effect
a12 <- rnorm(50)
s12 <- data.frame(id = "12", time = seq_len(50),
                  a = a12, b = 0.4 * a12 + rnorm(50),
                  c = rnorm(50), y = -0.5 * a12 + rnorm(50),
                  stringsAsFactors = FALSE)

panel <- rbind(panel, s10, s11, s12)

## ----screener-----------------------------------------------------------------
panel_clean <- i_screener(panel, cols = c("a", "b", "c", "y"), id_var = "id",
                          min_n_subject = 20, min_sd = 0.5, verbose = TRUE)

## ----screener-report----------------------------------------------------------
report <- i_screener(panel, cols = c("a", "b", "c", "y"), id_var = "id",
                     min_n_subject = 20, min_sd = 0.5, max_mode_pct = 0.80,
                     mode = "report", verbose = TRUE)
print(report)

## ----pmstandardize------------------------------------------------------------
panel_std <- pmstandardize(panel_clean, cols = c("a", "b", "c", "y"), id_var = "id",
                           verbose = TRUE)
head(panel_std[, c("id", "time", "a_psd", "b_psd", "c_psd", "y_psd")])

## ----detrender----------------------------------------------------------------
panel_dt <- i_detrender(panel_std, cols = c("a_psd", "b_psd", "c_psd", "y_psd"),
                        id_var = "id", timevar = "time", verbose = TRUE)
head(panel_dt[, c("id", "time", "a_psd_dt", "b_psd_dt", "c_psd_dt", "y_psd_dt")])

## ----iarimax------------------------------------------------------------------
result <- iarimax(panel_dt,
                  y_series  = "y_psd_dt",
                  x_series  = "a_psd_dt",
                  id_var    = "id",
                  timevar   = "time",
                  verbose   = TRUE)

## ----iarimax-fixed-d, eval = FALSE--------------------------------------------
#  # Fix d = 0 for all subjects (coefficients describe effects on levels)
#  result_d0 <- iarimax(panel_dt, y_series = "y_psd_dt", x_series = "a_psd_dt",
#                       id_var = "id", timevar = "time", fixed_d = 0)

## ----iarimax-summary----------------------------------------------------------
summary(result)

## ----iarimax-plot, fig.width = 6, fig.height = 5------------------------------
plot(result,
     y_series_name = "Outcome (y)",
     x_series_name = "Predictor (a)")

## ----results-df---------------------------------------------------------------
head(result$results_df[, c("id", "nAR", "nI", "nMA",
                           "estimate_a_psd_dt", "std.error_a_psd_dt",
                           "n_valid", "n_params", "raw_cor")])

## ----looping-machine----------------------------------------------------------
loop_result <- looping_machine(panel_dt,
                               a_series = "a_psd_dt",
                               b_series = "b_psd_dt",
                               c_series = "c_psd_dt",
                               id_var   = "id",
                               timevar  = "time",
                               verbose  = TRUE)

## ----looping-inspect----------------------------------------------------------
# Proportion of subjects with a detected positive directed loop
mean(loop_result$loop_df$Loop_positive_directed, na.rm = TRUE)

# Per-leg I-ARIMAX results are also accessible
summary(loop_result$iarimax_a_to_b)

## ----i-pval-------------------------------------------------------------------
result_pval <- i_pval(result)
result_pval$results_df[, c("id", "estimate_a_psd_dt", "pval_a_psd_dt")]

## ----sden---------------------------------------------------------------------
sden <- sden_test(result_pval)
summary(sden)

## ----sden-force---------------------------------------------------------------
sden_ent <- sden_test(result, test = "ENT")
summary(sden_ent)