library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly)      # interactive figures and ALL 3-D graphics
library(forecast)
library(tseries)

How this chapter uses graphics

Every two-dimensional figure is drawn with ggplot2 and rendered statically. Immediately after each one, the equivalent plot_ly() code appears in a chunk marked eval=FALSE, echo=TRUE.

Every three-dimensional figure is drawn with plot_ly() and evaluated. Temporal dependence is inherently two-parameter, the autocorrelation is a function of lag and of the model coefficient, prediction intervals of horizon and persistence, covariance structures of lag and $\rho$, and reading these from a fixed viewpoint discards the interaction that is the subject.

1 Learning objectives

After completing this chapter you will be able to:

State the stationarity and invertibility conditions for ARMA models in terms of the roots of the characteristic polynomials.
Interpret the ACF and PACF, with correct sampling bounds, and use them to identify candidate orders.
Test for a unit root using ADF and KPSS together, and read the four possible outcomes.
Explain why smoothing a series before fitting an ARIMA model fabricates autocorrelation and collapses forecast intervals.
Fit seasonal ARIMA models and check residuals with the Ljung–Box test.
Derive prediction intervals from the MA($\infty$) representation and explain why they widen with horizon.
Evaluate forecasts against naive benchmarks using scale-free metrics (MASE, RMSSE) under rolling-origin cross-validation.
Fit linear and generalized linear mixed models, and choose a covariance structure.
Distinguish GEE from GLMM by target of inference, apply the sandwich variance, and use the attenuation relation between them.
Explain which longitudinal methods remain valid under MCAR, MAR, and MNAR dropout.

Estimated time: 14–18 hours including exercises. Prerequisites: Chapter 3 (matrix computing), Chapter 6 (temporal validation, §6.9), Chapter 9 (rolling-origin CV and forecast metrics, §9.17), and Chapter 10 (streaming and drift, §10.13).

2 PART I: TEMPORAL FOUNDATIONS

3 What makes temporal data different

Every method in Chapters 5–11 assumed observations were exchangeable. Time series and repeated measures violate that in two distinct ways, and the two require different machinery.

	Time series	Longitudinal / panel
Structure	One (or few) long sequences	Many short sequences, one per unit
$n$ grows by	More time points	More units
Dependence	Within the sequence, across time	Within unit, across occasions
Asymptotics	$T\to\infty$	$N\to\infty$, $n_i$ fixed
Typical goal	Forecast the future	Estimate an effect, accounting for correlation
Core tools	ARIMA, state space, spectral	Mixed models, GEE, SEM

The distinction matters because the same word, “correlated observations” — points at opposite estimation problems. Forecasting exploits the dependence; longitudinal modeling treats it as a nuisance to be accounted for so that effect estimates and their standard errors are honest.

4 Stationarity

A process $\{X_t\}$ is weakly (second-order) stationary if

\[ \begin{aligned} \mathbb E[X_t]&=\mu &&\text{(constant mean)}\\ \operatorname{Var}(X_t)&=\sigma^2<\infty &&\text{(constant variance)}\\ \operatorname{Cov}(X_t,X_{t+h})&=\gamma(h) &&\text{(covariance depends only on the lag $h$)} \end{aligned} \]

Stationarity is what makes estimation possible at all: it is the assumption that lets a single realization stand in for the ensemble average, so that $\hat\gamma(h)=\frac1T\sum_t(x_t-\bar x)(x_{t+h}-\bar x)$ estimates something fixed.

Common misconception: “stationary means the series has no trend.” A constant mean is only one of three requirements. A series can be trend-free and still non-stationary through changing variance (volatility clustering, common in financial and physiological signals) or through a changing autocorrelation structure (a regime shift).

Conversely, a series can look strongly patterned and be stationary: an AR(1) with $\phi=0.95$ wanders in long excursions that resemble a trend, but its mean, variance, and autocovariance are all constant and it reverts. Distinguishing “trend-stationary” from “difference-stationary” is precisely what the unit-root tests of §12.5 are for.

4.1 Conditions on the polynomials

Write the backshift operator $B x_t = x_{t-1}$. An ARMA($p,q$) process is

\[\phi(B)\,X_t=\theta(B)\,\varepsilon_t, \qquad \phi(z)=1-\phi_1z-\cdots-\phi_pz^p, \quad \theta(z)=1+\theta_1z+\cdots+\theta_qz^q,\]

with $\varepsilon_t\sim\mathrm{WN}(0,\sigma^2)$, white noise of variance $\sigma^2$, a parameter to be estimated, not fixed at 1.

Stationarity and invertibility. \[\boxed{\;\text{stationary}\iff \text{all roots of }\phi(z)\text{ lie strictly \emph{outside} the unit circle}\;}\] \[\boxed{\;\text{invertible}\iff \text{all roots of }\theta(z)\text{ lie strictly \emph{outside} the unit circle}\;}\]

For AR(1), $\phi(z)=1-\phi z$ has root $1/\phi$, so stationarity is $|\phi|<1$. A root on the unit circle ($\phi=1$) is a unit root, a random walk, which is not stationary and does not revert.

Invertibility matters because a non-invertible MA has an observationally equivalent invertible twin with the same autocovariance; requiring invertibility makes the parameterization unique and lets $\varepsilon_t$ be recovered from the observed past.

set.seed(11)
T_s <- 400
sim_ar1 <- function(phi, n = T_s) as.numeric(arima.sim(list(ar = phi), n))

bind_rows(
  data.frame(t = 1:T_s, x = sim_ar1(0.5),  series = "phi = 0.50  (root 2.0)"),
  data.frame(t = 1:T_s, x = sim_ar1(0.95), series = "phi = 0.95  (root 1.05)"),
  data.frame(t = 1:T_s, x = cumsum(rnorm(T_s)), series = "phi = 1.00  (unit root)")) |>
  ggplot(aes(t, x)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_line(linewidth = 0.4, color = "steelblue") +
  facet_wrap(~ series, ncol = 1, scales = "free_y") +
  labs(title = "Approaching the unit circle",
       subtitle = "As the root approaches 1 the excursions lengthen; at the root the process stops reverting",
       x = "Time", y = NULL) +
  theme_dspa(10)

The middle panel is stationary despite looking trended, the root at 1.05 is outside the unit circle, so the process reverts, just slowly. The bottom panel has its root on the circle and does not revert at all. Telling these apart from a finite sample is genuinely hard, which is why unit-root tests have low power.

5 Autocorrelation

\[\rho(h)=\frac{\gamma(h)}{\gamma(0)},\qquad \hat\rho(h)=\frac{\sum_{t=1}^{T-h}(x_t-\bar x)(x_{t+h}-\bar x)}{\sum_{t=1}^{T}(x_t-\bar x)^2}\]

The partial autocorrelation $\alpha(h)$ is the correlation between $X_t$ and $X_{t+h}$ after removing the linear effect of the intervening $X_{t+1},\dots,X_{t+h-1}$.

The pair identifies model orders through a duality:

	ACF	PACF
AR($p$)	tails off (geometric / damped sine)	cuts off after lag $p$
MA($q$)	cuts off after lag $q$	tails off
ARMA($p,q$)	tails off	tails off

Why the PACF cuts off for AR. An AR($p$) is exactly a regression of $X_t$ on its own $p$ previous values. Once those are conditioned on, nothing further back adds information, so $\alpha(h)=0$ for $h>p$ by construction.

The sampling bounds must come from the series being plotted. Under a white-noise null, $\hat\rho(h)\approx N(0,1/T)$, giving the familiar $\pm1.96/\sqrt T$. $T$ is the length of the series whose ACF is shown, a bound computed from a different, longer or shorter series mis-scales every significance judgment on the plot.

For residual ACFs after fitting ARMA($p,q$), even that bound is wrong: the low-lag residual autocorrelations have variance smaller than $1/T$ because $p+q$ parameters were fitted to remove them. The correct omnibus check is the Ljung–Box statistic with degrees of freedom reduced by $p+q$ (§12.10).

# One helper, used everywhere. The bound is always computed from THIS series.
acf_frame <- function(x, lag.max = 36, type = c("correlation", "partial")) {
  type <- match.arg(type)
  a <- stats::acf(x, lag.max = lag.max, plot = FALSE, type = type,
                  na.action = na.pass)
  n <- sum(!is.na(x))
  data.frame(lag = as.numeric(a$lag), value = as.numeric(a$acf),
             bound = stats::qnorm(0.975) / sqrt(n),
             kind = ifelse(type == "correlation", "ACF", "PACF"))
}

plot_acf <- function(x, lag.max = 36, title = "") {
  d <- bind_rows(acf_frame(x, lag.max, "correlation") |> filter(lag > 0),
                 acf_frame(x, lag.max, "partial"))
  ggplot(d, aes(lag, value)) +
    geom_hline(yintercept = 0, color = "grey60") +
    geom_ribbon(aes(ymin = -bound, ymax = bound), fill = "grey88") +
    geom_segment(aes(xend = lag, yend = 0), linewidth = 0.7, color = "steelblue") +
    facet_wrap(~ kind, ncol = 1) +
    labs(title = title,
         subtitle = sprintf("Shaded band: +/- 1.96/sqrt(n) with n = %d, the length of THIS series",
                            sum(!is.na(x))),
         x = "Lag", y = NULL) +
    theme_dspa(10)
}

set.seed(21)
x_ar2 <- as.numeric(arima.sim(list(ar = c(0.6, 0.25)), 600))
x_ma2 <- as.numeric(arima.sim(list(ma = c(0.7, 0.5)), 600))

plot_acf(x_ar2, 24, "AR(2): PACF cuts off after lag 2, ACF tails off") /
  plot_acf(x_ma2, 24, "MA(2): ACF cuts off after lag 2, PACF tails off")

# --- Interactive equivalent ------------------------------------------------
d <- bind_rows(acf_frame(x_ar2, 24, "correlation") |> filter(lag > 0),
               acf_frame(x_ar2, 24, "partial"))
plot_ly(d, x = ~lag, y = ~value, color = ~kind, type = "bar") |>
  add_lines(x = ~range(lag), y = ~c(d$bound[1], d$bound[1]), name = "Upper bound",
            line = list(dash = "dash", color = "grey")) |>
  add_lines(x = ~range(lag), y = ~c(-d$bound[1], -d$bound[1]), name = "Lower bound",
            line = list(dash = "dash", color = "grey")) |>
  layout(title = "ACF and PACF", barmode = "group",
         xaxis = list(title = "Lag"), yaxis = list(title = "Correlation"))

The theoretical ACF of an AR(1) is $\rho(h)=\phi^{|h|}$, a two-parameter family, and worth rotating:

lag_grid <- 0:30
phi_grid <- seq(-0.95, 0.95, length.out = 60)
Zacf <- outer(phi_grid, lag_grid, function(p, h) p^h)

plot_ly(x = lag_grid, y = phi_grid, z = Zacf, type = "surface",
        colorscale = "RdBu", colorbar = list(title = "rho(h)")) |>
  layout(title = "Theoretical AR(1) autocorrelation over lag and phi",
         scene = list(xaxis = list(title = "Lag h"),
                      yaxis = list(title = "AR coefficient phi"),
                      zaxis = list(title = "rho(h)")))

Rotate along the $\phi$ axis. Near $\phi=0$ the surface collapses to a spike at lag 0, white noise. As $|\phi|\to1$ the decay flattens toward a ridge that never reaches zero, which is the geometry of a near-unit root. Negative $\phi$ gives the alternating sheet.

6 The Wold decomposition

Why should ARMA models be the right family at all? Because of a representation theorem.

Wold decomposition. Every zero-mean weakly stationary process $\{X_t\}$ admits a unique decomposition \[X_t=\underbrace{\sum_{j=0}^{\infty}\psi_j\varepsilon_{t-j}}_{\text{stochastic, }\psi_0=1,\ \sum\psi_j^2<\infty}+\underbrace{V_t}_{\text{deterministic}},\] where $\{\varepsilon_t\}$ is white noise uncorrelated with $V_t$.

Every stationary process is a linear filter of white noise. ARMA models are the family of rational such filters, $\psi(B)=\theta(B)/\phi(B)$, and any $\mathrm{MA}(\infty)$ with square-summable weights can be approximated arbitrarily well by a rational one with finitely many parameters. That is the justification for the whole enterprise.

The $\psi$-weights also give the forecast variance directly:

\[\boxed{\;\operatorname{Var}\big(X_{T+h}\mid \mathcal F_T\big)=\sigma^2\sum_{j=0}^{h-1}\psi_j^2\;}\]

which is why prediction intervals widen with horizon, and by exactly how much (§12.11).

# psi-weights of an ARMA, and the forecast variance they imply
psi_of <- function(ar = numeric(0), ma = numeric(0), n = 20)
  c(1, ARMAtoMA(ar = ar, ma = ma, lag.max = n - 1))

data.frame(
  h = 1:8,
  `AR(1) phi=0.5`  = round(cumsum(psi_of(ar = 0.5)[1:8]^2), 4),
  `AR(1) phi=0.9`  = round(cumsum(psi_of(ar = 0.9)[1:8]^2), 4),
  `MA(2)`          = round(cumsum(psi_of(ma = c(0.7, 0.5))[1:8]^2), 4),
  check.names = FALSE)

Read down the columns: the MA(2) forecast variance stops growing after $h=3$, because an MA($q$) is unpredictable beyond $q$ steps and the interval saturates at the unconditional variance. The AR variances keep growing, and faster for larger $\phi$.

7 Unit-root testing

Two tests with opposite nulls, and they are meant to be read together.

\[ \begin{aligned} \textbf{ADF (augmented Dickey--Fuller): }\quad &H_0:\ \text{unit root (non-stationary)}\\ \textbf{KPSS: }\quad &H_0:\ \text{(trend-)stationary} \end{aligned} \]

ADF	KPSS	Conclusion
Reject	Fail to reject	Stationary — both agree
Fail to reject	Reject	Unit root — both agree; difference the series
Fail to reject	Fail to reject	Inconclusive — not enough data to decide
Reject	Reject	Neither — suggests heteroscedasticity or a structural break

Common misconception: “the ADF test rejected, so the series is stationary.” Rejecting a unit-root null is evidence against a unit root, not evidence for stationarity, and the ADF test is well known to have low power against near-unit-root alternatives, so failing to reject frequently means only that $T$ is small.

Running KPSS alongside converts a one-sided statement into a four-way decision, and the two “agree” cells are the only ones that support a confident conclusion.

set.seed(31)
series_list <- list(
  `white noise`          = rnorm(400),
  `AR(1), phi = 0.5`     = sim_ar1(0.5),
  `AR(1), phi = 0.95`    = sim_ar1(0.95),
  `random walk`          = cumsum(rnorm(400)),
  `RW, differenced`      = diff(cumsum(rnorm(400))))

do.call(rbind, lapply(names(series_list), function(nm) {
  x <- series_list[[nm]]
  adf <- suppressWarnings(tseries::adf.test(x, alternative = "stationary"))
  kps <- suppressWarnings(tseries::kpss.test(x, null = "Level"))
  data.frame(series = nm,
             adf_p = round(adf$p.value, 4),
             adf_rejects_unit_root = adf$p.value < 0.05,
             kpss_p = round(kps$p.value, 4),
             kpss_rejects_stationarity = kps$p.value < 0.05,
             ndiffs_suggested = forecast::ndiffs(x))
}))

forecast::ndiffs() automates the decision by applying the test repeatedly until it stops rejecting. Note the fourth row: the random walk needs one difference, and the fifth row confirms that one difference suffices.

8 PART II: ARIMA

9 The Beijing PM2.5 series

The Beijing PM2.5 dataset records hourly concentrations of particulate matter under 2.5 microns, expressed as an Air Quality Index, from 2008 to 2016, roughly 69,000 hourly observations.

pm_raw <- if (ONLINE) {
  dspa_read("https://umich.instructure.com/files/1823138/download?download_frd=1",
            "beijing_pm25.csv")
} else NULL

if (is.null(pm_raw)) {
  # Synthetic substitute with the same structure: daily cycle, annual cycle,
  # AR(1) persistence, right skew, and a realistic missingness rate.
  set.seed(41); n_h <- 69335
  tt <- seq_len(n_h)
  base <- 90 + 25 * sin(2*pi*tt/24 - 1.2) + 40 * sin(2*pi*tt/8766 + 0.5)
  ar <- as.numeric(arima.sim(list(ar = 0.92), n_h, sd = 18))
  val <- pmax(0, round(base + ar))
  val[sample(n_h, round(0.065 * n_h))] <- -999
  pm_raw <- data.frame(
    Date..LST. = format(seq(as.POSIXct("2008-04-08 15:00:00", tz = "UTC"),
                            by = "hour", length.out = n_h), "%m/%d/%Y %H:%M"),
    Value = val)
  message("Note: Beijing file unavailable; using a synthetic substitute.")
}
c(rows = nrow(pm_raw), columns = ncol(pm_raw))

#>    rows columns 
#>   69335      11

library(lubridate)

# The timestamp format uses %H (24-hour), not %h -- which is an alias for the
# abbreviated MONTH name and would silently drop the time component.
pm <- pm_raw |>
  mutate(datetime = as.POSIXct(Date..LST., format = "%m/%d/%Y %H:%M", tz = "UTC"),
         value = ifelse(Value == -999, NA_real_, Value)) |>
  filter(!is.na(datetime)) |>
  arrange(datetime) |>
  select(datetime, value)

c(rows = nrow(pm),
  missing = sum(is.na(pm$value)),
  missing_pct = round(100 * mean(is.na(pm$value)), 2),
  range_start = format(min(pm$datetime)), range_end = format(max(pm$datetime)),
  negative_values_remaining = sum(pm$value < 0, na.rm = TRUE))

#>                      rows                   missing               missing_pct 
#>                   "69335"                    "4408"                    "6.36" 
#>               range_start                 range_end negative_values_remaining 
#>     "2008-04-08 15:00:00"     "2016-04-30 23:00:00"                       "8"

The $-999$ codes are converted to NA once; no replacement sentinel is introduced. Every downstream step treats missingness as missingness.

9.1 Imputing temporal gaps

Common misconception: “fill the gaps with the mean, then proceed.” For cross-sectional data mean imputation is merely inefficient. For a time series it is destructive: it inserts flat runs, breaks the local autocorrelation at every imputed point, and biases the sample ACF toward zero, precisely the quantity the whole analysis is about to estimate.

Temporal imputation must respect temporal structure. Linear interpolation is the minimum; Kalman smoothing (imputeTS::na_kalman()) is the principled choice, since it returns the conditional mean of the missing values under a fitted state-space model and therefore preserves both the level and the dependence.

library(imputeTS)

pm_mean <- pm$value; pm_mean[is.na(pm_mean)] <- mean(pm$value, na.rm = TRUE)
pm_interp <- imputeTS::na_interpolation(pm$value)
pm_kalman <- dspa_try(imputeTS::na_kalman(ts(pm$value, frequency = 24),
                                          model = "auto.arima"),
                      fallback = imputeTS::na_seadec(ts(pm$value, frequency = 24)),
                      label = "Kalman imputation")
pm_kalman <- as.numeric(pm_kalman)

acf_at <- function(x, h) as.numeric(acf(x, lag.max = h, plot = FALSE)$acf[h + 1])
data.frame(
  method = c("mean imputation", "linear interpolation", "Kalman smoothing"),
  acf_lag1  = round(c(acf_at(pm_mean, 1),  acf_at(pm_interp, 1),  acf_at(pm_kalman, 1)), 4),
  acf_lag24 = round(c(acf_at(pm_mean, 24), acf_at(pm_interp, 24), acf_at(pm_kalman, 24)), 4),
  sd        = round(c(sd(pm_mean), sd(pm_interp), sd(pm_kalman)), 2))

Mean imputation attenuates the lag-1 autocorrelation and shrinks the variance, because 6.5% of the series has been replaced by a constant. Kalman smoothing preserves both.

win <- 3000:3400
bind_rows(
  data.frame(t = win, v = pm_mean[win],   method = "Mean imputation"),
  data.frame(t = win, v = pm_kalman[win], method = "Kalman smoothing")) |>
  ggplot(aes(t, v, color = method)) +
  geom_line(linewidth = 0.5) +
  facet_wrap(~ method, ncol = 1) +
  scale_color_manual(values = c(`Mean imputation` = "#D8433B",
                                 `Kalman smoothing` = "#3B7DD8"), guide = "none") +
  labs(title = "The same window under two imputation schemes",
       subtitle = "Mean imputation inserts flat segments at the imputed points; Kalman smoothing does not",
       x = "Hour index", y = "PM2.5 AQI") +
  theme_dspa(10)

pm$value_imp <- pm_kalman

# Modeling uses DAILY means; the hourly series is retained for the daily-cycle
# diagnostics. Aggregating is a modeling choice, stated once, not a filter
# applied silently to the series being modeled.
pm_daily <- pm |>
  mutate(date = as.Date(datetime)) |>
  summarise(value = mean(value_imp), .by = date) |>
  arrange(date)

series <- if (PM25_RESOLUTION == "daily") {
  ts(pm_daily$value, frequency = 7)      # weekly period for the daily series
} else {
  ts(pm$value_imp, frequency = 24)       # daily period for the hourly series
}
c(resolution = PM25_RESOLUTION, n = length(series),
  frequency = frequency(series),
  mean = round(mean(series), 2), sd = round(sd(series), 2))

#> resolution          n  frequency       mean         sd 
#>    "daily"     "2890"        "7"     "93.7"    "73.95"

ggplot(pm_daily, aes(date, value)) +
  geom_line(linewidth = 0.3, color = "grey45") +
  geom_smooth(method = "loess", span = 0.15, se = FALSE,
              color = "#3B7DD8", linewidth = 0.9) +
  labs(title = "Beijing daily mean PM2.5, 2008-2016",
       subtitle = "Grey: daily means. Blue: a loess smooth, shown for ORIENTATION only -- the models below are fitted to the grey series",
       x = NULL, y = "PM2.5 AQI") +
  theme_dspa()

# --- Interactive equivalent ------------------------------------------------
plot_ly(pm_daily, x = ~date, y = ~value, type = "scatter", mode = "lines",
        name = "Daily mean", line = list(width = 1)) |>
  layout(title = "Beijing daily mean PM2.5",
         xaxis = list(title = "Date"), yaxis = list(title = "PM2.5 AQI"),
         hovermode = "x unified")

The loess curve is labelled as orientation, and the models below are fitted to the unsmoothed series. §12.15 shows what happens when that distinction is not maintained.

9.2 The daily cycle

# Hour extracted from the TIMESTAMP, not from row position -- with 6.5% of
# hours missing, a positional modular index would mis-assign every hour after
# the first gap.
hourly_profile <- pm |>
  mutate(hour = lubridate::hour(datetime)) |>
  summarise(mean = mean(value_imp),
            q20 = quantile(value_imp, 0.20),
            q80 = quantile(value_imp, 0.80),
            n = dplyr::n(), .by = hour) |>
  arrange(hour)
head(hourly_profile, 4)

ggplot(hourly_profile, aes(hour, mean)) +
  geom_ribbon(aes(ymin = q20, ymax = q80), fill = "grey86") +
  geom_line(linewidth = 1, color = "steelblue") +
  geom_point(size = 1.8) +
  scale_x_continuous(breaks = seq(0, 23, 3)) +
  labs(title = "Diurnal PM2.5 profile, averaged over eight years",
       subtitle = "Ribbon: 20th to 80th percentile. Hours taken from the timestamp, so gaps do not shift the alignment",
       x = "Hour of day (local)", y = "PM2.5 AQI") +
  theme_dspa()

10 Differencing and integration

The “I” in ARIMA is integration order $d$: the number of differences needed to reach stationarity.

\[\nabla X_t=(1-B)X_t=X_t-X_{t-1}, \qquad \nabla^2X_t=X_t-2X_{t-1}+X_{t-2}\]

An ARIMA($p,d,q$) is an ARMA($p,q$) fitted to $\nabla^dX_t$:

\[\phi(B)\,(1-B)^d X_t=\theta(B)\,\varepsilon_t.\]

Over-differencing has a signature and a cost. Differencing a series that is already stationary introduces a unit root into the MA polynomial: if $X_t$ is white noise, $\nabla X_t$ has $\rho(1)=-0.5$ and is non-invertible. The tell-tale signs are a large negative lag-1 autocorrelation and an inflated residual variance.

Difference the minimum number of times the tests require, usually $d\le2$, almost never more.

set.seed(51)
wn <- rnorm(500)
data.frame(
  series = c("white noise", "differenced once", "differenced twice"),
  acf_lag1 = round(c(acf_at(wn, 1), acf_at(diff(wn), 1), acf_at(diff(wn, differences = 2), 1)), 4),
  variance = round(c(var(wn), var(diff(wn)), var(diff(wn, differences = 2))), 4),
  ndiffs_says = c(ndiffs(wn), ndiffs(diff(wn)), ndiffs(diff(wn, differences = 2))))

Each unnecessary difference drives $\rho(1)$ toward $-0.5$ and roughly doubles the variance. ndiffs() correctly says zero differences are needed in all three cases.

d_needed <- forecast::ndiffs(series)
D_needed <- forecast::nsdiffs(series)
c(regular_differences = d_needed, seasonal_differences = D_needed,
  adf_p = round(suppressWarnings(tseries::adf.test(series)$p.value), 4),
  kpss_p = round(suppressWarnings(tseries::kpss.test(series)$p.value), 4))

#>  regular_differences seasonal_differences                adf_p 
#>               0.0000               0.0000               0.0100 
#>               kpss_p 
#>               0.0624

plot_acf(as.numeric(series), lag.max = 40,
         "Beijing daily PM2.5: ACF and PACF of the unsmoothed series")

11 Seasonal ARIMA

Series with a repeating cycle of period $m$ need seasonal terms. $\text{SARIMA}(p,d,q)(P,D,Q)_m$ is

\[\boxed{\;\Phi(B^m)\,\phi(B)\,(1-B^m)^{D}(1-B)^{d}\,X_t=\Theta(B^m)\,\theta(B)\,\varepsilon_t\;}\]

with $\Phi,\Theta$ polynomials in $B^m$ acting across cycles while $\phi,\theta$ act within them.

Common misconception: “just use a large MA order to capture the seasonality.” Fitting $\text{ARIMA}(1,1,24)$ to hourly data to capture a daily cycle estimates 24 free MA parameters where $\text{SARIMA}(1,1,1)(0,1,1)_{24}$ captures the same structure with three.

The seasonal parameterization is not merely more parsimonious. It encodes the right structure, that lag 24, 48, 72 are related to each other in a way lags 5, 11, 17 are not, and it generalizes to horizons beyond one cycle, which a long non-seasonal MA cannot.

# Fit on a training span; the final year is held out for honest evaluation.
n_ser <- length(series)
h_out <- if (PM25_RESOLUTION == "daily") 180 else 24 * 30
train <- window(series, end = time(series)[n_ser - h_out])
test  <- window(series, start = time(series)[n_ser - h_out + 1])
c(train_length = length(train), test_length = length(test))

#> train_length  test_length 
#>         2710          180

fit_auto <- forecast::auto.arima(train, seasonal = TRUE, stepwise = TRUE,
                                 approximation = TRUE)
fit_auto

#> Series: train 
#> ARIMA(3,1,2)(0,0,2)[7] 
#> 
#> Coefficients:
#>          ar1    ar2     ar3     ma1     ma2    sma1   sma2
#>       -0.273  0.417  -0.129  -0.048  -0.921  -0.033  0.029
#> s.e.   0.038  0.032   0.023   0.033   0.032   0.020  0.020
#> 
#> sigma^2 = 3366:  log likelihood = -14842.4
#> AIC=29700.8   AICc=29700.8   BIC=29748

fit_ns <- forecast::Arima(train, order = arimaorder(fit_auto)[1:3],
                          seasonal = c(0, 0, 0))
data.frame(
  model = c(sprintf("auto: %s", paste(arimaorder(fit_auto), collapse = ",")),
            "same orders, no seasonal terms"),
  n_params = c(length(coef(fit_auto)), length(coef(fit_ns))),
  AIC = round(c(AIC(fit_auto), AIC(fit_ns)), 2),
  BIC = round(c(BIC(fit_auto), BIC(fit_ns)), 2),
  sigma2 = round(c(fit_auto$sigma2, fit_ns$sigma2), 3))

12 Estimation and search cost

ARIMA parameters are estimated by maximum likelihood, and the likelihood is evaluated by putting the model in state-space form and running the Kalman filter:

\[\ell(\phi,\theta,\sigma^2)=-\frac12\sum_{t=1}^{T}\left[\log(2\pi F_t)+\frac{v_t^2}{F_t}\right],\]

where $v_t$ is the one-step prediction error and $F_t$ its variance, both produced by the filter recursion. This costs $O(T)$ per likelihood evaluation for fixed $(p,q)$, and it handles missing values exactly, the filter simply skips the update step, which is why state-space methods are the right tool for gappy series.

Task	Cost
One likelihood evaluation (Kalman)	$O(T(p+q+m)^2)$
One model fit (optimizer, $I$ iterations)	$O(I\,T(p+q+m)^2)$
`auto.arima` stepwise search	$\approx O\big((p_{\max}+q_{\max})\big)$ fits
`auto.arima` exhaustive (`stepwise = FALSE`)	$O\big(p_{\max}q_{\max}P_{\max}Q_{\max}\big)$ fits
`stl` decomposition	$O(T)$ per inner iteration

The stepwise search is the default for a reason. An exhaustive search over $p,q\le5$ and $P,Q\le2$ is $5\times5\times3\times3\times2=450$ fits; the stepwise algorithm typically evaluates a few dozen. It can miss the global optimum, and stepwise = FALSE is worth the cost on a final model.

13 Diagnostics

If the orders are right, the residuals should be indistinguishable from white noise. The omnibus test is Ljung–Box:

\[Q(m)=T(T+2)\sum_{h=1}^{m}\frac{\hat\rho_\varepsilon(h)^2}{T-h}\ \overset{H_0}{\sim}\ \chi^2_{m-p-q}\]

The degrees of freedom must be reduced by the number of fitted ARMA parameters. Testing residuals with $m$ degrees of freedom instead of $m-p-q$ makes the test conservative, it under-rejects, and so fails to flag models that are genuinely inadequate. Box.test(..., fitdf = p + q) supplies the correction.

res <- residuals(fit_auto)
n_par <- length(coef(fit_auto))

lb <- do.call(rbind, lapply(c(10, 20, 30), function(m) {
  bt <- Box.test(res, lag = m, type = "Ljung-Box", fitdf = n_par)
  data.frame(lag = m, statistic = round(bt$statistic, 2),
             df = bt$parameter, p_value = round(bt$p.value, 4),
             residuals_look_white = bt$p.value > 0.05)
}))
lb

plot_acf(as.numeric(res), 36,
         sprintf("Residual ACF/PACF: ARIMA(%s)", paste(arimaorder(fit_auto), collapse = ",")))

rdf <- data.frame(t = seq_along(res), r = as.numeric(res))
p_ts <- ggplot(rdf, aes(t, r)) +
  geom_hline(yintercept = 0, color = "grey60") +
  geom_line(linewidth = 0.3, color = "steelblue") +
  labs(title = "Residuals over time", x = NULL, y = NULL) + theme_dspa(9)
p_qq <- ggplot(rdf, aes(sample = r)) +
  stat_qq(size = 0.5, alpha = 0.5) + stat_qq_line(color = "firebrick") +
  labs(title = "Normal Q-Q", x = NULL, y = NULL) + theme_dspa(9)
p_hi <- ggplot(rdf, aes(r)) +
  geom_histogram(aes(y = after_stat(density)), bins = 50,
                 fill = "steelblue", color = "white") +
  stat_function(fun = dnorm, args = list(mean = 0, sd = sd(rdf$r)),
                color = "firebrick", linewidth = 0.8) +
  labs(title = "Residual distribution vs. Normal", x = NULL, y = NULL) + theme_dspa(9)

(p_ts / (p_qq | p_hi)) +
  patchwork::plot_annotation(
    title = "ARIMA residual diagnostics",
    subtitle = "Heavy tails are common in air-quality data; they widen true prediction intervals beyond the Gaussian ones")

# --- Interactive equivalent ------------------------------------------------
plot_ly(rdf, x = ~t, y = ~r, type = "scatter", mode = "lines",
        name = "Residuals") |>
  add_lines(x = ~range(t), y = rep(1.96 * sd(rdf$r), 2), name = "+1.96 SD",
            line = list(dash = "dash")) |>
  add_lines(x = ~range(t), y = rep(-1.96 * sd(rdf$r), 2), name = "-1.96 SD",
            line = list(dash = "dash")) |>
  layout(title = "ARIMA residuals", hovermode = "x unified",
         xaxis = list(title = "Index"), yaxis = list(title = "Residual"))

The Q-Q plot is the one to read carefully. Air-quality residuals are typically heavy-tailed, which means the Gaussian prediction intervals of §12.11 are too narrow in the tails, worth stating alongside any interval quoted from this model.

14 PART III: FORECASTING AND ITS EVALUATION

15 Prediction intervals

The $h$-step point forecast is the conditional mean $\hat X_{T+h}=\mathbb E[X_{T+h}\mid\mathcal F_T]$, and its variance follows directly from the Wold $\psi$-weights (§12.4):

\[\boxed{\;\operatorname{Var}\big(X_{T+h}\mid\mathcal F_T\big)=\sigma^2\sum_{j=0}^{h-1}\psi_j^2\;}\]

giving the Gaussian interval $\hat X_{T+h}\pm z_{1-\alpha/2}\,\sigma\sqrt{\sum_{j<h}\psi_j^2}$.

Three consequences worth stating explicitly.

Intervals widen with horizon, by exactly $\sqrt{\sum_{j<h}\psi_j^2}$, not arbitrarily, and not linearly.

For a stationary model the width saturates. As $h\to\infty$, $\sum_{j<h}\psi_j^2\to\operatorname{Var}(X_t)/\sigma^2$, so long-horizon intervals converge to the unconditional variance: the forecast becomes the series mean and the interval becomes the series spread.

For a differenced ($d\ge1$) model it does not saturate. The $\psi$-weights do not decay, and the interval grows without bound, which is the correct behavior for a random walk and the reason long-horizon forecasts of integrated series are nearly uninformative.

## OLD:
# psi_var <- function(ar = numeric(0), ma = numeric(0), d = 0, h = 40) {
#   # psi-weights of the ARIMA, obtained by expanding (1-B)^-d into the AR side
#   ar_full <- if (d == 0) ar else {
#     poly_ar <- c(1, -ar); poly_d <- c(1, -1)
#     for (k in seq_len(d)) poly_ar <- convolve(poly_ar, rev(poly_d), type = "open")
#     -poly_ar[-1]
#   }
#   psi <- c(1, ARMAtoMA(ar = ar_full, ma = ma, lag.max = h - 1))
#   cumsum(psi^2)
# }
psi_var <- function(ar = numeric(0), ma = numeric(0), d = 0, h = 40) {
    # FIX: Handle h = 1 explicitly to prevent lag.max = 0 in ARMAtoMA
    if (h == 1) { return(1) }
    # psi-weights of the ARIMA, obtained by expanding (1-B)^-d into the AR side
    ar_full <- if (d == 0) ar else {
        poly_ar <- c(1, -ar); poly_d <- c(1, -1)
        for (k in seq_len(d)) 
          poly_ar <- convolve(poly_ar, rev(poly_d), type = "open") - poly_ar[-1]
    }
    psi <- c(1, ARMAtoMA(ar = ar_full, ma = ma, lag.max = h - 1))
    cumsum(psi^2)
}
bind_rows(
  data.frame(h = 1:40, w = sqrt(psi_var(ar = 0.5)),          model = "AR(1), phi = 0.5"),
  data.frame(h = 1:40, w = sqrt(psi_var(ar = 0.9)),          model = "AR(1), phi = 0.9"),
  data.frame(h = 1:40, w = sqrt(psi_var(ar = 0.5, d = 1)),   model = "ARIMA(1,1,0), phi = 0.5")) |>
  ggplot(aes(h, w, color = model)) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#3B7DD8", "#7FB069", "#D8433B")) +
  labs(title = "Prediction interval width against forecast horizon",
       subtitle = "Stationary models saturate at the unconditional SD; the integrated model does not",
       x = "Horizon h", y = expression(sqrt(sum(psi[j]^2, j<h))), color = NULL) +
  theme_dspa()

h_grid <- 1:36
phi_grid2 <- seq(0.05, 0.98, length.out = 50)
Zint <- outer(phi_grid2, h_grid, Vectorize(function(p, h) sqrt(psi_var(ar = p, h = h)[h])))

plot_ly(x = h_grid, y = phi_grid2, z = Zint, type = "surface",
        colorscale = "Viridis",
        colorbar = list(title = "Interval\nwidth factor")) |>
  layout(title = "Prediction-interval width over horizon and AR persistence",
         scene = list(xaxis = list(title = "Horizon h"),
                      yaxis = list(title = "AR coefficient phi"),
                      zaxis = list(title = "sqrt(sum psi_j^2)")))

Rotate along the horizon axis at low $\phi$: the surface flattens almost immediately, a weakly persistent process becomes unpredictable within a few steps and the interval stops growing. At $\phi\to1$ the surface keeps climbing across the whole horizon range, because a near-unit-root process carries information arbitrarily far forward and its uncertainty accumulates.

16 Benchmarks

Common misconception: “my model has an $R^2$ of 0.95, so it forecasts well.” On a persistent series the naive forecast, predict the last observed value, is extremely hard to beat, and it has no parameters at all. Any forecast accuracy figure quoted without a benchmark is uninterpretable.

\[ \begin{aligned} \textbf{Naive: }\quad &\hat X_{T+h}=X_T\\ \textbf{Seasonal naive: }\quad &\hat X_{T+h}=X_{T+h-m\lceil h/m\rceil}\\ \textbf{Drift: }\quad &\hat X_{T+h}=X_T+h\cdot\frac{X_T-X_1}{T-1}\\ \textbf{Mean: }\quad &\hat X_{T+h}=\bar X \end{aligned} \]

h_test <- length(test)
fc_list <- list(
  `ARIMA (auto)`   = forecast::forecast(fit_auto, h = h_test),
  `Naive`          = forecast::naive(train, h = h_test),
  `Seasonal naive` = forecast::snaive(train, h = h_test),
  `Drift`          = forecast::rwf(train, h = h_test, drift = TRUE),
  `Mean`           = forecast::meanf(train, h = h_test),
  `ETS`            = forecast::forecast(forecast::ets(train), h = h_test))

17 Scale-free accuracy

RMSE and MAE are in the units of the series, so they cannot be compared across series. Two scale-free alternatives, both defined relative to the naive benchmark on the training data:

\[\boxed{\;\mathrm{MASE}=\frac{\frac1h\sum_{t=T+1}^{T+h}\big|X_t-\hat X_t\big|}{\frac{1}{T-m}\sum_{t=m+1}^{T}\big|X_t-X_{t-m}\big|}\;}\]

\[\mathrm{RMSSE}=\sqrt{\frac{\frac1h\sum_t\big(X_t-\hat X_t\big)^2}{\frac{1}{T-m}\sum_t\big(X_t-X_{t-m}\big)^2}}\]

MASE $<1$ means the forecast beats the in-sample naive benchmark; MASE $>1$ means it loses to it. The denominator uses the training data, so the metric is well defined even when the test window is short, and it is finite even when the series contains zeros, unlike MAPE, which is undefined at zero and asymmetric in over- versus under-prediction.

acc_tab <- do.call(rbind, lapply(names(fc_list), function(nm) {
  a <- forecast::accuracy(fc_list[[nm]], test)
  data.frame(model = nm,
             RMSE = round(a["Test set", "RMSE"], 3),
             MAE  = round(a["Test set", "MAE"], 3),
             MASE = round(a["Test set", "MASE"], 4),
             beats_naive = a["Test set", "MASE"] < 1)
}))
acc_tab |> arrange(MASE)

Read the beats_naive column first. A model that cannot beat snaive on a seasonal series is not adding value, whatever its RMSE happens to be.

ggplot(acc_tab, aes(reorder(model, MASE), MASE, fill = beats_naive)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "firebrick") +
  geom_col(width = 0.65) +
  coord_flip() +
  scale_fill_manual(values = c(`TRUE` = "#3B7DD8", `FALSE` = "grey70"),
                    labels = c("loses to naive", "beats naive"), name = NULL) +
  labs(title = "Out-of-sample MASE against the naive benchmark",
       subtitle = "Dashed line at 1: the in-sample naive forecast. Below it is better",
       x = NULL, y = "MASE") +
  theme_dspa()

fc <- forecast::forecast(fit_auto, h = h_test, level = c(80, 95))

recent <- max(1, length(train) - 3 * h_test)
obs_df <- data.frame(
  t = c(time(train)[recent:length(train)], time(test)),
  v = c(as.numeric(train)[recent:length(train)], as.numeric(test)),
  part = c(rep("Training (observed)", length(train) - recent + 1),
           rep("Held out (observed)", length(test))))
fc_df <- data.frame(t = time(test), mean = as.numeric(fc$mean),
                    lo80 = fc$lower[, 1], hi80 = fc$upper[, 1],
                    lo95 = fc$lower[, 2], hi95 = fc$upper[, 2])

ggplot() +
  geom_ribbon(data = fc_df, aes(t, ymin = lo95, ymax = hi95), fill = "grey88") +
  geom_ribbon(data = fc_df, aes(t, ymin = lo80, ymax = hi80), fill = "grey75") +
  geom_line(data = obs_df, aes(t, v, color = part), linewidth = 0.45) +
  geom_line(data = fc_df, aes(t, mean), color = "#D8433B", linewidth = 0.9) +
  scale_color_manual(values = c(`Training (observed)` = "grey35",
                                 `Held out (observed)` = "#3B7DD8")) +
  labs(title = sprintf("Forecast from ARIMA(%s), with genuinely held-out data",
                       paste(arimaorder(fit_auto), collapse = ",")),
       subtitle = "Bands are 80% (dark) and 95% (light), read from the columns matching their labels",
       x = "Time", y = "PM2.5 AQI", color = NULL) +
  theme_dspa()

# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
  add_ribbons(x = fc_df$t, ymin = fc_df$lo95, ymax = fc_df$hi95,
              name = "95% interval", line = list(width = 0), opacity = 0.3) |>
  add_ribbons(x = fc_df$t, ymin = fc_df$lo80, ymax = fc_df$hi80,
              name = "80% interval", line = list(width = 0), opacity = 0.4) |>
  add_lines(x = obs_df$t, y = obs_df$v, name = "Observed",
            line = list(color = "grey40", width = 1)) |>
  add_lines(x = fc_df$t, y = fc_df$mean, name = "Point forecast",
            line = list(color = "red", width = 2)) |>
  layout(title = "ARIMA forecast with 80% and 95% intervals",
         xaxis = list(title = "Time"), yaxis = list(title = "PM2.5 AQI"),
         hovermode = "x unified")

The observed held-out series is plotted, not the fitted values, the point of holding data out is to see the forecast next to what actually happened.

18 Rolling-origin cross-validation

A single train/test split gives one number from one origin. Rolling-origin evaluation (Chapter 9, §9.17.2) repeats the exercise at many origins and averages, producing an estimate with far less variance and a breakdown by horizon.

far_arima <- function(x, h, order, seasonal) {
  fit <- tryCatch(forecast::Arima(x, order = order, seasonal = seasonal),
                  error = function(e) NULL)
  if (is.null(fit)) return(forecast::naive(x, h = h))
  forecast::forecast(fit, h = h)
}
ord <- arimaorder(fit_auto)
ord_ns <- ord[1:3]
ord_s  <- if (length(ord) > 3) ord[4:6] else c(0, 0, 0)

set.seed(61)
e_arima <- forecast::tsCV(series, far_arima, h = TSCV_HORIZONS,
                          initial = floor(0.7 * length(series)),
                          order = ord_ns, seasonal = ord_s)
e_snaive <- forecast::tsCV(series, function(x, h) forecast::snaive(x, h = h),
                           h = TSCV_HORIZONS,
                           initial = floor(0.7 * length(series)))

scale_denom <- mean(abs(diff(as.numeric(series), lag = frequency(series))), na.rm = TRUE)
cv_tab <- data.frame(
  horizon = 1:TSCV_HORIZONS,
  MASE_arima  = round(colMeans(abs(e_arima),  na.rm = TRUE) / scale_denom, 4),
  MASE_snaive = round(colMeans(abs(e_snaive), na.rm = TRUE) / scale_denom, 4))
cv_tab$arima_wins <- cv_tab$MASE_arima < cv_tab$MASE_snaive
cv_tab

cv_tab |>
  select(horizon, ARIMA = MASE_arima, `Seasonal naive` = MASE_snaive) |>
  pivot_longer(-horizon, names_to = "model", values_to = "MASE") |>
  ggplot(aes(horizon, MASE, color = model)) +
  geom_line(linewidth = 1) + geom_point(size = 2.2) +
  scale_color_manual(values = c(ARIMA = "#3B7DD8", `Seasonal naive` = "#D8433B")) +
  scale_x_continuous(breaks = 1:TSCV_HORIZONS) +
  labs(title = "Rolling-origin MASE by forecast horizon",
       subtitle = "Averaged over many origins, so the comparison does not depend on one arbitrary split point",
       x = "Horizon (steps ahead)", y = "MASE", color = NULL) +
  theme_dspa()

The advantage typically shrinks with horizon. At $h=1$ the ARIMA exploits recent dependence; by $h=6$ the forecast is approaching the unconditional mean and the structured model has little left to contribute. Reporting a single horizon hides that.

origins <- seq(floor(0.7 * length(series)), length(series) - TSCV_HORIZONS - 1,
               length.out = 12) |> round()
Zcv <- outer(seq_along(origins), 1:TSCV_HORIZONS, Vectorize(function(oi, h) {
  o <- origins[oi]
  abs(e_arima[o, h]) / scale_denom
}))
Zcv[!is.finite(Zcv)] <- NA

plot_ly(x = 1:TSCV_HORIZONS, y = as.numeric(time(series))[origins], z = Zcv,
        type = "surface", colorscale = "Inferno", reversescale = TRUE,
        colorbar = list(title = "Scaled\nabs. error")) |>
  layout(title = "Rolling-origin scaled error over horizon and forecast origin",
         scene = list(xaxis = list(title = "Horizon"),
                      yaxis = list(title = "Forecast origin (time)"),
                      zaxis = list(title = "Scaled absolute error")))

The ridges running along the origin axis are the informative feature: accuracy is not uniform in time. Periods of high volatility produce large errors at every horizon, and averaging over origins, as the table above does — summarizes a surface that is far from flat.

19 The smoothing trap

Common misconception: “smooth the series first to reveal the signal, then fit ARIMA to the smooth version.” This is the single most damaging preprocessing error in applied time-series work, because it manufactures exactly the structure the model is meant to discover.

Let $m_t=\frac1k\sum_{j=0}^{k-1}x_{t-j}$ be a $k$-point moving average. Even when $\{x_t\}$ is independent white noise, \[\operatorname{Corr}(m_t,m_{t-h})=\Big(1-\frac{|h|}{k}\Big)_{+},\] so the smoothed series is an MA($k-1$) process with autocorrelation declining linearly out to lag $k$. Any ACF/PACF read off it, any unit-root test applied to it, and any order selected from it describe the filter.

The variance consequence is equally severe: $\operatorname{Var}(m_t)=\sigma^2/k$ for independent inputs, so the innovation variance estimated from a smoothed series is too small by roughly $k$, and the prediction intervals too narrow by roughly $\sqrt k$.

Smoothing is a display device. Model the data.

set.seed(71)
pure_noise <- rnorm(2000)
ma_filter <- function(x, k) as.numeric(stats::filter(x, rep(1/k, k), sides = 1))

bind_rows(lapply(c(1, 5, 20, 60), function(k) {
  s <- ma_filter(pure_noise, k)
  a <- acf(s, lag.max = 80, plot = FALSE, na.action = na.pass)
  data.frame(lag = as.numeric(a$lag), acf = as.numeric(a$acf),
             theory = pmax(0, 1 - abs(as.numeric(a$lag))/k),
             window = sprintf("k = %d", k))
})) |>
  mutate(window = factor(window, levels = sprintf("k = %d", c(1, 5, 20, 60)))) |>
  ggplot(aes(lag)) +
  geom_hline(yintercept = 0, color = "grey65") +
  geom_line(aes(y = acf, color = "Observed on smoothed noise"), linewidth = 0.7) +
  geom_line(aes(y = theory, color = "Theory: (1 - |h|/k)+"),
            linetype = "dashed", linewidth = 0.7) +
  facet_wrap(~ window) +
  scale_color_manual(values = c(`Observed on smoothed noise` = "#3B7DD8",
                                 `Theory: (1 - |h|/k)+` = "#D8433B")) +
  labs(title = "Smoothing white noise creates autocorrelation that was never there",
       subtitle = "The input is independent at every k. The apparent structure is entirely the moving-average filter",
       x = "Lag", y = "ACF", color = NULL) +
  theme_dspa(10)

The input is the same independent white noise in all four panels. At $k=1$ (no smoothing) the ACF is flat, correctly. By $k=60$ it declines linearly to lag 60 and would lead any identification procedure to a high-order MA model.

k_grid <- c(1, 2, 3, 5, 8, 12, 20, 35, 60, 100)
lag_grid2 <- 0:60
Zsm <- outer(k_grid, lag_grid2, Vectorize(function(k, h) {
  s <- ma_filter(pure_noise, k)
  as.numeric(acf(s, lag.max = 61, plot = FALSE, na.action = na.pass)$acf[h + 1])
}))

plot_ly(x = lag_grid2, y = k_grid, z = Zsm, type = "surface",
        colorscale = "RdBu", reversescale = TRUE,
        colorbar = list(title = "Induced\nACF")) |>
  layout(title = "Autocorrelation induced in white noise by a k-point moving average",
         scene = list(xaxis = list(title = "Lag h"),
                      yaxis = list(title = "Smoothing window k", type = "log"),
                      zaxis = list(title = "ACF")))

The surface is flat at $k=1$ and rises into a broad plateau as $k$ grows. The input never changes. Every feature visible here was created by the filter.

# What smoothing does to the model that follows it
set.seed(73)
demo <- as.numeric(arima.sim(list(ar = 0.4), 1200))
demo_sm <- ma_filter(demo, 30); demo_sm <- demo_sm[!is.na(demo_sm)]

fit_raw <- forecast::auto.arima(demo, seasonal = FALSE, stepwise = TRUE)
fit_sm  <- forecast::auto.arima(demo_sm, seasonal = FALSE, stepwise = TRUE)

data.frame(
  fitted_to = c("raw series", "30-point moving average"),
  selected_order = c(paste(arimaorder(fit_raw), collapse = ","),
                     paste(arimaorder(fit_sm), collapse = ",")),
  sigma2 = round(c(fit_raw$sigma2, fit_sm$sigma2), 5),
  interval_width_h1 = round(1.96 * sqrt(c(fit_raw$sigma2, fit_sm$sigma2)), 4),
  truth = c("AR(1), phi = 0.4", "same data, smoothed"))

The smoothed fit selects a different, more elaborate model and reports an innovation variance smaller by roughly the window length, so its prediction intervals are narrower by roughly $\sqrt{30}\approx5.5$ while describing the same underlying series. Those intervals would fail their nominal coverage badly.

20 Dynamic regression

When exogenous predictors are available, the right formulation is regression with ARIMA errors:

\[\boxed{\;X_t=\boldsymbol\beta^\top\mathbf z_t+\eta_t,\qquad \phi(B)(1-B)^d\eta_t=\theta(B)\varepsilon_t\;}\]

This is not the same as putting lagged $X$ on the right-hand side. In the “ARIMAX” form $X_t=\boldsymbol\beta^\top\mathbf z_t+\phi_1X_{t-1}+\cdots+\varepsilon_t$, the coefficient $\beta_k$ is not the effect of raising $z_k$ by one unit — because $X_{t-1}$ is itself a function of $z_{k,t-1}$, so the covariate’s influence is spread across the autoregressive terms and the interpretation is conditional on the lagged response.

In the regression-with-ARIMA-errors form the covariates enter the mean and the dependence is confined to the error process, so $\beta_k$ has its ordinary regression interpretation. forecast::Arima(..., xreg = ) fits this parameterization, which is why it is the one to prefer.

Both series must be differenced consistently. If $X_t$ requires differencing, the covariates must be differenced too, or the regression is spurious.

set.seed(81)
weeks <- 100; days <- 7 * weeks

# A generating process whose signal SURVIVES: every term contributes to the
# final value of y, and the holiday indicator marks the days that actually
# receive the holiday effect.
weekday  <- rep(1:7, weeks)
holiday  <- as.integer(seq_len(days) %% 28 == 0)     # period 28, matching the effect
weekend  <- as.integer(weekday %in% c(6, 7))

beta_holiday <- 260
beta_weekend <- -120
trend        <- 0.35

eta <- as.numeric(arima.sim(list(ar = c(0.55, 0.20)), days, sd = 45))
y   <- 1000 + trend * seq_len(days) +
       beta_holiday * holiday + beta_weekend * weekend + eta

sim_df <- data.frame(day = seq_len(days), y = y, holiday = holiday,
                     weekend = weekend, weekday = factor(weekday))
c(true_holiday_effect = beta_holiday, true_weekend_effect = beta_weekend,
  true_trend_per_day = trend,
  holidays_in_series = sum(holiday),
  mean_on_holidays = round(mean(y[holiday == 1])),
  mean_off_holidays = round(mean(y[holiday == 0])))

#> true_holiday_effect true_weekend_effect  true_trend_per_day  holidays_in_series 
#>              260.00             -120.00                0.35               25.00 
#>    mean_on_holidays   mean_off_holidays 
#>             1281.00             1100.00

The last two lines confirm the effect is present in the realized data, the check that a simulation has produced what it intended.

h_sim <- 56
train_idx <- 1:(days - h_sim); test_idx <- (days - h_sim + 1):days

xreg_train <- as.matrix(sim_df[train_idx, c("holiday", "weekend", "day")])
# Future covariates are built for the FORECAST period, not copied from the start
xreg_test  <- as.matrix(sim_df[test_idx,  c("holiday", "weekend", "day")])
c(train_day_range = range(xreg_train[, "day"]),
  test_day_range  = range(xreg_test[, "day"]))

#> train_day_range1 train_day_range2  test_day_range1  test_day_range2 
#>                1              644              645              700

fit_dyn <- forecast::auto.arima(ts(sim_df$y[train_idx], frequency = 7),
                                xreg = xreg_train, seasonal = TRUE)
fit_dyn

#> Series: ts(sim_df$y[train_idx], frequency = 7) 
#> Regression with ARIMA(2,0,0)(1,0,0)[7] errors 
#> 
#> Coefficients:
#>         ar1    ar2    sar1  intercept  holiday   weekend   day
#>       0.493  0.294  -0.069   1007.885  257.860  -124.627  0.37
#> s.e.  0.038  0.038   0.041     14.867    8.421     3.234  0.04
#> 
#> sigma^2 = 1922:  log likelihood = -3345.39
#> AIC=6706.77   AICc=6707   BIC=6742.51

est <- coef(fit_dyn); se <- sqrt(diag(fit_dyn$var.coef))
data.frame(
  term = c("holiday", "weekend", "day"),
  truth = c(beta_holiday, beta_weekend, trend),
  estimate = round(est[c("holiday", "weekend", "day")], 3),
  std_error = round(se[c("holiday", "weekend", "day")], 3),
  covers_truth = abs(est[c("holiday","weekend","day")] -
                       c(beta_holiday, beta_weekend, trend)) <
                 1.96 * se[c("holiday", "weekend", "day")])

Each estimate recovers its true value within two standard errors, which is what a simulation study is for, and what cannot be checked when the generating process has overwritten its own signal.

fc_dyn <- forecast::forecast(fit_dyn, xreg = xreg_test, level = c(80, 95))

fc_dyn_df <- data.frame(day = test_idx, mean = as.numeric(fc_dyn$mean),
                        lo80 = fc_dyn$lower[, 1], hi80 = fc_dyn$upper[, 1],
                        lo95 = fc_dyn$lower[, 2], hi95 = fc_dyn$upper[, 2],
                        actual = sim_df$y[test_idx])

ggplot(fc_dyn_df, aes(day)) +
  geom_ribbon(aes(ymin = lo95, ymax = hi95), fill = "grey88") +
  geom_ribbon(aes(ymin = lo80, ymax = hi80), fill = "grey75") +
  geom_line(data = sim_df[(days - 3*h_sim):(days - h_sim), ],
            aes(day, y), color = "grey40", linewidth = 0.4) +
  geom_line(aes(y = actual), color = "#3B7DD8", linewidth = 0.5) +
  geom_line(aes(y = mean), color = "#D8433B", linewidth = 0.9) +
  labs(title = "Dynamic regression forecast with correctly aligned future covariates",
       subtitle = "Blue: held-out actuals. Red: point forecast. Bands from the matching interval columns",
       x = "Day", y = "Arrivals") +
  theme_dspa()

c(MASE = round(forecast::accuracy(fc_dyn, sim_df$y[test_idx])["Test set", "MASE"], 4),
  MASE_snaive = round(forecast::accuracy(
    forecast::snaive(ts(sim_df$y[train_idx], frequency = 7), h = h_sim),
    sim_df$y[test_idx])["Test set", "MASE"], 4))

#>        MASE MASE_snaive 
#>      0.8302      1.6461

21 Decomposition and Prophet

An alternative to ARIMA is additive decomposition: model trend, seasonality, and holidays as separate components.

\[X_t=g(t)+s(t)+h(t)+\varepsilon_t\]

stl() does this non-parametrically with loess; Prophet (Taylor & Letham, 2018) does it with a piecewise-linear or logistic trend, Fourier-series seasonality, and explicit holiday indicators.

Prophet’s seasonal terms must match the data’s sampling frequency. weekly.seasonality = TRUE on monthly observations asks the model to estimate a within-week pattern from data that contains no within-week variation. And make_future_dataframe(periods = n) defaults to freq = "day", for monthly data, freq = "month" must be given, or the horizon is wrong by a factor of about 30.

dec <- stats::stl(series, s.window = "periodic", robust = TRUE)
dec_df <- data.frame(
  t = rep(as.numeric(time(series)), 4),
  value = c(as.numeric(series), as.numeric(dec$time.series[, "seasonal"]),
            as.numeric(dec$time.series[, "trend"]),
            as.numeric(dec$time.series[, "remainder"])),
  component = factor(rep(c("Observed", "Seasonal", "Trend", "Remainder"),
                         each = length(series)),
                     levels = c("Observed", "Seasonal", "Trend", "Remainder")))

ggplot(dec_df, aes(t, value)) +
  geom_line(linewidth = 0.35, color = "steelblue") +
  facet_wrap(~ component, ncol = 1, scales = "free_y") +
  labs(title = "STL decomposition of the Beijing series",
       subtitle = sprintf("Seasonal period = %d, matching the declared frequency of the series",
                          frequency(series)),
       x = "Time", y = NULL) +
  theme_dspa(9)

var_share <- apply(dec$time.series, 2, var) / var(as.numeric(series))
round(var_share, 4)

#>  seasonal     trend remainder 
#>    0.0001    0.2209    0.7237

# --- Interactive equivalent ------------------------------------------------
p1 <- plot_ly(x = time(series), y = as.numeric(series), type = "scatter",
              mode = "lines", name = "Observed")
p2 <- plot_ly(x = time(series), y = dec$time.series[, "seasonal"],
              type = "scatter", mode = "lines", name = "Seasonal")
p3 <- plot_ly(x = time(series), y = dec$time.series[, "trend"],
              type = "scatter", mode = "lines", name = "Trend")
p4 <- plot_ly(x = time(series), y = dec$time.series[, "remainder"],
              type = "scatter", mode = "lines", name = "Remainder")
subplot(p1, p2, p3, p4, nrows = 4, shareX = TRUE) |>
  layout(title = "STL decomposition", hovermode = "x unified")

# --- Fetching Google Trends data (run once; results cached to disk) ---------
# The API rate-limits aggressively, so retries with exponential backoff are
# needed. Values are a RELATIVE index (0-100) normalized WITHIN each request:
# they are not comparable across keywords, windows, or geographies.
library(gtrendsR)
fetch_trends <- function(keywords, geo = "US", time = "2015-01-01 2025-01-01",
                         attempts = 5) {
  for (i in seq_len(attempts)) {
    out <- tryCatch(gtrends(keywords, geo = geo, gprop = "web", time = time)[[1]],
                    error = function(e) NULL)
    if (!is.null(out)) return(out)
    Sys.sleep(60 * 2^(i - 1))
  }
  NULL
}
trends <- fetch_trends(c("data science", "machine learning"))
write.csv(trends, file.path(dspa_cache_dir(), "gtrends_us.csv"), row.names = FALSE)

PROPHET_OK <- requireNamespace("prophet", quietly = TRUE)

trends <- dspa_read("", "gtrends_us.csv")     # cached from the chunk above
if (is.null(trends)) {
  set.seed(91); n_m <- 120
  d <- seq(as.Date("2015-01-01"), by = "month", length.out = n_m)
  trends <- data.frame(
    date = d, keyword = "data science",
    hits = pmin(100, pmax(0, round(20 + 0.55 * seq_len(n_m) +
      12 * sin(2*pi*seq_len(n_m)/12) + as.numeric(arima.sim(list(ar = .5), n_m, sd = 4))))))
  message("Note: Google Trends cache absent; using a synthetic monthly substitute.")
}

ts_df <- trends |>
  filter(keyword == "data science") |>
  transmute(ds = as.Date(date), y = suppressWarnings(as.numeric(hits))) |>
  filter(!is.na(y)) |> arrange(ds)

# Infer the sampling frequency from the data rather than assuming it
step_days <- as.numeric(median(diff(ts_df$ds)))
freq_label <- if (step_days <= 1) "day" else if (step_days <= 8) "week" else "month"
c(observations = nrow(ts_df), median_step_days = step_days,
  inferred_frequency = freq_label)

#>       observations   median_step_days inferred_frequency 
#>              "120"               "31"            "month"

if (PROPHET_OK) {
  m <- prophet::prophet(ts_df,
                        weekly.seasonality = (freq_label == "day"),
                        daily.seasonality  = FALSE,
                        yearly.seasonality = TRUE)
  fut <- prophet::make_future_dataframe(m, periods = 24, freq = freq_label)
  pp <- predict(m, fut)
  tail(pp[, c("ds", "yhat", "yhat_lower", "yhat_upper")], 3)
}

Weekly seasonality is enabled only when the data are daily, and make_future_dataframe receives the inferred frequency, so 24 periods means 24 months for monthly data, not 24 days.

22 PART IV: REPEATED MEASURES

Time series ask what happens next. Longitudinal studies ask what an effect is, when each subject contributes several correlated observations.

23 The correlated-data problem

Suppose $n_i$ measurements on each of $N$ subjects, and we fit ordinary least squares ignoring the clustering. The point estimates remain unbiased. The standard errors do not.

For a design with $n$ observations per cluster and within-cluster correlation $\rho$, the variance of a cluster-level mean is inflated by the design effect

\[\boxed{\;\mathrm{DEFF}=1+(n-1)\rho\;}\]

so the effective sample size is $N_{\text{eff}}=Nn/\mathrm{DEFF}$.

Common misconception: “more measurements per subject means more data.” Only partly. With $\rho=0.7$ and 10 measurements per subject, the design effect is $1+9(0.7)=7.3$: 100 subjects × 10 visits carry about the information of 137 independent observations, not 1,000.

Ignoring the clustering makes standard errors too small by a factor of $\sqrt{\mathrm{DEFF}}\approx2.7$ here, so nominal 95% intervals achieve far less than 95% coverage and p-values are correspondingly anti-conservative.

set.seed(101)
## OLD: error occurs because out — the result of replicate — does not have row names "ols" and "lmm" as you expect. Consequently, out["lmm", ] fails with subscript out of bounds. This can happen if replicate fails to simplify the list of results into a matrix with the correct row names (e.g., when some replicates return NA or have inconsistent names), or if the lmer model occasionally fails in a way that interrupts the normal return.
# sim_cluster <- function(N = 60, n = 10, rho = 0.7, beta = 0, reps = 400) {
#   sig_b <- sqrt(rho); sig_e <- sqrt(1 - rho)
#   out <- replicate(reps, {
#     b <- rnorm(N, sd = sig_b)
#     id <- rep(seq_len(N), each = n)
#     x  <- rep(rnorm(N), each = n)              # cluster-level covariate
#     y  <- beta * x + b[id] + rnorm(N * n, sd = sig_e)
#     m_ols <- lm(y ~ x)
#     p_ols <- summary(m_ols)$coefficients["x", 4]
#     m_lmm <- suppressMessages(lme4::lmer(y ~ x + (1 | id)))
#     se_l  <- sqrt(diag(vcov(m_lmm)))["x"]
#     p_lmm <- 2 * pnorm(-abs(lme4::fixef(m_lmm)["x"] / se_l))
#     c(ols = p_ols < 0.05, lmm = p_lmm < 0.05)
#   })
#   c(rho = rho, DEFF = 1 + (n - 1) * rho,
#     ols_type_I = mean(out["ols", ]), lmm_type_I = mean(out["lmm", ]))
# }
set.seed(101)
sim_cluster <- function(N = 60, n = 10, rho = 0.7, beta = 0, reps = 400) {
  sig_b <- sqrt(rho)
  sig_e <- sqrt(1 - rho)
  
  # Use vapply to ensure a matrix with 2 rows and reps columns
  out <- vapply(seq_len(reps), function(r) {
    b <- rnorm(N, sd = sig_b)
    id <- rep(seq_len(N), each = n)
    x  <- rep(rnorm(N), each = n)              # cluster-level covariate
    y  <- beta * x + b[id] + rnorm(N * n, sd = sig_e)
    
    # OLS
    m_ols <- lm(y ~ x)
    p_ols <- summary(m_ols)$coefficients["x", 4]
    
    # LMM with error handling
    p_lmm <- tryCatch({
      m_lmm <- suppressMessages(lme4::lmer(y ~ x + (1 | id)))
      se_l  <- sqrt(diag(vcov(m_lmm)))["x"]
      2 * pnorm(-abs(lme4::fixef(m_lmm)["x"] / se_l))
    }, error = function(e) NA)   # return NA if lmer fails
    
    c(ols = p_ols < 0.05, lmm = p_lmm < 0.05)
  }, numeric(2))
  
  # Ensure row names are set (vapply should preserve them, but just in case)
  rownames(out) <- c("ols", "lmm")
  
  c(rho = rho,
    DEFF = 1 + (n - 1) * rho,
    ols_type_I = mean(out[1, ]),   # use index, not name
    lmm_type_I = mean(out[2, ]))
}

library(lme4)
as.data.frame(do.call(rbind, lapply(c(0.1, 0.4, 0.7), function(r) sim_cluster(rho = r)))) |>
  mutate(across(everything(), \(z) round(z, 4)))

With no true effect ($\beta=0$), a valid test rejects 5% of the time. OLS rejects far more often, and the inflation tracks the design effect. The mixed model holds its nominal rate.

24 Linear mixed models

\[\boxed{\;\mathbf y_i=X_i\boldsymbol\beta+Z_i\mathbf b_i+\boldsymbol\varepsilon_i,\qquad \mathbf b_i\sim N(\mathbf 0,G),\quad \boldsymbol\varepsilon_i\sim N(\mathbf 0,R_i)\;}\]

with $\boldsymbol\beta$ the fixed effects (population-level), $\mathbf b_i$ the random effects (subject-level departures), and the implied marginal covariance

\[V_i=\operatorname{Var}(\mathbf y_i)=Z_iGZ_i^\top+R_i.\]

Two properties make this the workhorse for longitudinal data.

Shrinkage. The predicted random effect is $\hat{\mathbf b}_i=GZ_i^\top V_i^{-1}(\mathbf y_i-X_i\hat{\boldsymbol\beta})$, a weighted compromise between the subject’s own data and the population mean. Subjects with fewer or noisier observations are pulled harder toward the average, an automatic, optimal borrowing of strength.

Validity under MAR. Because the model specifies a full likelihood, it uses all available observations and remains valid when dropout depends on observed history (§12.22).

REML, not ML, for the variance components. Maximum likelihood estimates of variance components are biased downward because they do not account for the degrees of freedom consumed by $\hat{\boldsymbol\beta}$, the same correction that makes $s^2$ divide by $n-1$. Restricted maximum likelihood maximizes the likelihood of error contrasts and removes the bias.

The consequence for practice: compare models differing in fixed effects with ML (REML likelihoods for different fixed-effect structures are not comparable), and models differing only in random effects with REML.

ppmi_raw <- if (ONLINE) {
  dspa_read("https://umich.instructure.com/files/330397/download?download_frd=1",
            "ppmi_long.csv")
} else NULL

if (is.null(ppmi_raw)) {
  set.seed(111); N_s <- 200; visits <- c(0, 6, 12, 18, 24, 36)
  b0 <- rnorm(N_s, 0, 6); b1 <- rnorm(N_s, 0, 0.12)
  grp <- sample(c("PD", "Control"), N_s, TRUE, prob = c(0.65, 0.35))
  ppmi_raw <- do.call(rbind, lapply(seq_len(N_s), function(i) data.frame(
    FID_IID = i, time_visit = visits,
    ResearchGroup = grp[i],
    Age = round(rnorm(1, 62, 9)), Sex = rbinom(1, 1, 0.6),
    Weight = round(rnorm(1, 78, 14)),
    UPDRS_part_I   = pmax(0, round(2 + b0[i]*0.2 + (0.03 + b1[i]*0.3)*visits +
                                     ifelse(grp[i]=="PD", 3, 0) + rnorm(6, 0, 1.5))),
    UPDRS_part_II  = pmax(0, round(6 + b0[i]*0.4 + (0.06 + b1[i]*0.6)*visits +
                                     ifelse(grp[i]=="PD", 6, 0) + rnorm(6, 0, 2.5))),
    UPDRS_part_III = pmax(0, round(14 + b0[i] + (0.11 + b1[i])*visits +
                                     ifelse(grp[i]=="PD", 12, 0) + rnorm(6, 0, 4))),
    L_cingulate_gyrus_Volume = rnorm(6, 5800, 400),
    R_cingulate_gyrus_Volume = rnorm(6, 5750, 400))))
  message("Note: PPMI file unavailable; using a synthetic longitudinal substitute.")
}
c(rows = nrow(ppmi_raw), columns = ncol(ppmi_raw))

#>    rows columns 
#>    1764      31

ppmi <- ppmi_raw |>
  filter(!is.na(time_visit), !is.na(UPDRS_part_III)) |>
  mutate(
    subject = factor(FID_IID),
    # time_visit stays NUMERIC: a longitudinal trend analysis needs a SLOPE.
    # The factor version is used only where a per-visit contrast is wanted.
    time = as.numeric(time_visit),
    time_f = factor(time_visit),
    # The reference level is stated explicitly rather than left to alphabetics
    group = factor(ResearchGroup),
    updrs3 = as.numeric(UPDRS_part_III))
levels(ppmi$group)

#> [1] "Control" "PD"      "SWEDD"

ppmi$group <- stats::relevel(ppmi$group, ref = levels(ppmi$group)[1])

c(subjects = nlevels(ppmi$subject),
  observations = nrow(ppmi),
  visits_per_subject = round(nrow(ppmi) / nlevels(ppmi$subject), 2),
  reference_group = levels(ppmi$group)[1])

#>           subjects       observations visits_per_subject    reference_group 
#>              "440"             "1210"             "2.75"          "Control"

table(ppmi$group)

#> 
#> Control      PD   SWEDD 
#>     127     936     147

set.seed(113)
show_ids <- sample(levels(ppmi$subject), min(60, nlevels(ppmi$subject)))

ggplot(filter(ppmi, subject %in% show_ids),
       aes(time, updrs3, group = subject, color = group)) +
  geom_line(alpha = 0.35, linewidth = 0.4) +
  geom_smooth(aes(group = group), method = "lm", se = TRUE, linewidth = 1.1) +
  scale_color_manual(values = c("lightblue", "red", "green")) +
  labs(title = "Individual trajectories with group-level trends",
       subtitle = "Thin lines: subjects. Thick lines: fitted group slopes. The spread between subjects is what the random intercept absorbs",
       x = "Months since baseline", y = "UPDRS Part III", color = NULL) +
  theme_dspa()

# --- Interactive equivalent ------------------------------------------------
plot_ly(filter(ppmi, subject %in% show_ids),
        x = ~time, y = ~updrs3, split = ~subject, color = ~group,
        type = "scatter", mode = "lines", opacity = 0.4,
        line = list(width = 1), showlegend = FALSE) |>
  layout(title = "Individual UPDRS Part III trajectories",
         xaxis = list(title = "Months since baseline"),
         yaxis = list(title = "UPDRS Part III"))

library(lmerTest)      # adds Satterthwaite degrees of freedom to lmer

m_ols <- lm(updrs3 ~ group * time + Age + Sex, data = ppmi)
m_ri  <- lmerTest::lmer(updrs3 ~ group * time + Age + Sex + (1 | subject),
                        data = ppmi, REML = TRUE)
m_rs  <- lmerTest::lmer(updrs3 ~ group * time + Age + Sex + (1 + time | subject),
                        data = ppmi, REML = TRUE)

summary(m_rs)$coefficients |> round(4)

#>                 Estimate Std. Error      df t value Pr(>|t|)
#> (Intercept)      -7.0952     2.9021 512.600 -2.4448   0.0148
#> groupPD          21.2968     1.4076 545.238 15.1298   0.0000
#> groupSWEDD       12.7163     1.8715 471.069  6.7948   0.0000
#> time              0.0009     0.0437 655.242  0.0196   0.9844
#> Age               0.1352     0.0376 491.967  3.6008   0.0003
#> Sex               0.1187     0.8056 474.432  0.1473   0.8829
#> groupPD:time      0.0001     0.0448 691.302  0.0024   0.9981
#> groupSWEDD:time   0.0080     0.0499 748.771  0.1610   0.8721

# Random-effect structures differ only in the RANDOM part, so REML is valid here
data.frame(
  model = c("OLS (ignores clustering)", "Random intercept", "Random intercept + slope"),
  n_params = c(length(coef(m_ols)),
               length(lme4::fixef(m_ri)) + 2, length(lme4::fixef(m_rs)) + 4),
  AIC = round(c(AIC(m_ols), AIC(m_ri), AIC(m_rs)), 1),
  BIC = round(c(BIC(m_ols), BIC(m_ri), BIC(m_rs)), 1),
  se_group_time = round(c(
    summary(m_ols)$coefficients[grep(":time", rownames(summary(m_ols)$coefficients))[1], 2],
    sqrt(diag(vcov(m_ri)))[grep(":time", names(lme4::fixef(m_ri)))[1]],
    sqrt(diag(vcov(m_rs)))[grep(":time", names(lme4::fixef(m_rs)))[1]]), 5))

anova(m_ri, m_rs, refit = FALSE)

The se_group_time column is the practical point: the OLS standard error for the group-by-time interaction is the smallest of the three, and it is the one that is wrong.

24.1 The intraclass correlation

\[\mathrm{ICC}=\frac{\sigma_b^2}{\sigma_b^2+\sigma_\varepsilon^2}\]

— the proportion of total variance attributable to between-subject differences, and equivalently the correlation between two observations on the same subject.

vc <- as.data.frame(lme4::VarCorr(m_ri))
sig_b2 <- vc$vcov[vc$grp == "subject"]
sig_e2 <- vc$vcov[vc$grp == "Residual"]
n_bar <- nrow(ppmi) / nlevels(ppmi$subject)

c(between_subject_var = round(sig_b2, 3),
  residual_var = round(sig_e2, 3),
  ICC = round(sig_b2 / (sig_b2 + sig_e2), 4),
  mean_visits = round(n_bar, 2),
  design_effect = round(1 + (n_bar - 1) * sig_b2 / (sig_b2 + sig_e2), 3),
  effective_n = round(nrow(ppmi) / (1 + (n_bar - 1) * sig_b2/(sig_b2 + sig_e2))))

#> between_subject_var        residual_var                 ICC         mean_visits 
#>             52.8200             24.1090              0.6866              2.7500 
#>       design_effect         effective_n 
#>              2.2020            550.0000

25 Covariance structures

The random-effects specification implies a marginal covariance. When the implied structure is wrong, an explicit one can be modelled instead.

\[ \begin{aligned} \textbf{Independence: }&\ R=\sigma^2I\\ \textbf{Compound symmetry: }&\ R_{jk}=\sigma^2\big[\rho+(1-\rho)\mathbb 1\{j=k\}\big]\\ \textbf{AR(1): }&\ R_{jk}=\sigma^2\rho^{|t_j-t_k|}\\ \textbf{Unstructured: }&\ R_{jk}=\sigma_{jk}\quad (n(n+1)/2\text{ parameters}) \end{aligned} \]

Structure	Parameters	Right when
Independence	1	No within-subject correlation (rare)
Compound symmetry	2	Correlation is constant regardless of spacing
AR(1)	2	Correlation decays with time separation
Toeplitz	$n$	Depends on lag but not geometrically
Unstructured	$n(n+1)/2$	Few time points, plenty of subjects

Compound symmetry and AR(1) make opposite claims about time. Compound symmetry says visits one month apart and thirty months apart are equally correlated; AR(1) says correlation decays geometrically. For anything measured repeatedly over an extended period, AR(1) is usually the more plausible starting point, and a random intercept alone implies compound symmetry, which is worth knowing when the residual ACF shows decay.

library(nlme)
times <- sort(unique(ppmi$time))

fit_cs <- nlme::gls(updrs3 ~ group * time + Age + Sex, data = ppmi,
                    correlation = nlme::corCompSymm(form = ~ time | subject),
                    method = "REML")
fit_ar <- nlme::gls(updrs3 ~ group * time + Age + Sex, data = ppmi,
                    correlation = nlme::corAR1(form = ~ time | subject),
                    method = "REML")
fit_id <- nlme::gls(updrs3 ~ group * time + Age + Sex, data = ppmi,
                    method = "REML")

data.frame(
  structure = c("Independence", "Compound symmetry", "AR(1)"),
  logLik = round(c(logLik(fit_id), logLik(fit_cs), logLik(fit_ar)), 1),
  AIC = round(c(AIC(fit_id), AIC(fit_cs), AIC(fit_ar)), 1),
  BIC = round(c(BIC(fit_id), BIC(fit_cs), BIC(fit_ar)), 1))

lag_grid3 <- 0:12
rho_grid <- seq(0.05, 0.95, length.out = 50)
Z_ar <- outer(rho_grid, lag_grid3, function(r, h) r^h)
Z_cs <- outer(rho_grid, lag_grid3, function(r, h) ifelse(h == 0, 1, r))

plot_ly() |>
  add_surface(x = lag_grid3, y = rho_grid, z = Z_ar, opacity = 0.9,
              showscale = FALSE, colorscale = "Blues", name = "AR(1)") |>
  add_surface(x = lag_grid3, y = rho_grid, z = Z_cs - 1.3, opacity = 0.9,
              showscale = FALSE, colorscale = "Reds",
              name = "Compound symmetry (offset)") |>
  layout(title = "Implied correlation: AR(1) (upper, blue) vs. compound symmetry (lower, red, offset)",
         scene = list(xaxis = list(title = "Lag between visits"),
                      yaxis = list(title = "rho"),
                      zaxis = list(title = "Correlation")))

Rotate along the lag axis. The AR(1) surface decays geometrically; the compound-symmetry surface drops once from 1 to $\rho$ and is then perfectly flat. That flatness is the assumption, visits thirty months apart as correlated as consecutive ones.

26 Generalized linear mixed models

For non-Gaussian outcomes, add a link:

\[g\big(\mathbb E[y_{ij}\mid\mathbf b_i]\big)=\mathbf x_{ij}^\top\boldsymbol\beta+\mathbf z_{ij}^\top\mathbf b_i, \qquad \mathbf b_i\sim N(\mathbf 0,G)\]

The likelihood now requires integrating out $\mathbf b_i$,

\[L(\boldsymbol\beta,G)=\prod_{i=1}^{N}\int\prod_{j=1}^{n_i}f\big(y_{ij}\mid\mathbf b_i\big)\,\phi(\mathbf b_i;G)\,d\mathbf b_i,\]

and the integral has no closed form for non-identity links. Three standard approximations: adaptive Gauss–Hermite quadrature (accurate, cost grows exponentially in the random-effect dimension), the Laplace approximation (quadrature with one point; the lme4 default), and penalized quasi-likelihood (fast, and known to be biased for binary outcomes with few observations per cluster).

set.seed(121)
ppmi <- ppmi |>
  mutate(worsened = as.integer(updrs3 > median(updrs3, na.rm = TRUE)))

m_glmm <- lme4::glmer(worsened ~ group * time + Age + Sex + (1 | subject),
                      data = ppmi, family = binomial(),
                      control = glmerControl(optimizer = "bobyqa"))
round(summary(m_glmm)$coefficients, 4)

#>                 Estimate Std. Error z value Pr(>|z|)
#> (Intercept)     -24.9936  3891.4280 -0.0064   0.9949
#> groupPD          20.8665  3891.4276  0.0054   0.9957
#> groupSWEDD       17.0167  3891.4277  0.0044   0.9965
#> time             -0.0234   161.8976 -0.0001   0.9999
#> Age               0.0731     0.0218  3.3461   0.0008
#> Sex               0.1792     0.4296  0.4172   0.6766
#> groupPD:time      0.0214   161.8976  0.0001   0.9999
#> groupSWEDD:time   0.0487   161.8976  0.0003   0.9998

c(random_intercept_sd = round(sqrt(as.numeric(lme4::VarCorr(m_glmm)$subject)), 4))

#> random_intercept_sd 
#>              3.0868

27 Generalized estimating equations

GEE takes a different route: rather than a full likelihood, specify only the mean and a working correlation, and solve

\[\boxed{\;U(\boldsymbol\beta)=\sum_{i=1}^{N}D_i^\top V_i^{-1}\big(\mathbf y_i-\boldsymbol\mu_i(\boldsymbol\beta)\big)=\mathbf 0\;}\]

with $D_i=\partial\boldsymbol\mu_i/\partial\boldsymbol\beta$ and $V_i=A_i^{1/2}R_i(\alpha)A_i^{1/2}$.

27.1 The sandwich estimator

This is the entire reason GEE exists. The variance of $\hat{\boldsymbol\beta}$ is estimated by \[\boxed{\;\widehat{\operatorname{Var}}(\hat{\boldsymbol\beta})=\underbrace{\mathcal I_0^{-1}}_{\text{bread}}\ \underbrace{\left(\sum_i D_i^\top V_i^{-1}\widehat{\operatorname{Cov}}(\mathbf y_i)V_i^{-1}D_i\right)}_{\text{meat}}\ \underbrace{\mathcal I_0^{-1}}_{\text{bread}},\qquad \mathcal I_0=\sum_i D_i^\top V_i^{-1}D_i\;}\] where the “meat” uses the empirical residual cross-products $(\mathbf y_i-\hat{\boldsymbol\mu}_i)(\mathbf y_i-\hat{\boldsymbol\mu}_i)^\top$ rather than the assumed $V_i$.

The consequence is remarkable: $\hat{\boldsymbol\beta}$ is consistent and its standard errors are valid even when the working correlation is wrong — provided the mean model is right and $N$ is large. Choosing the working correlation well buys efficiency, not validity.

The caveat is the sample size. The sandwich estimator is downward biased when the number of clusters is small; below roughly 40 clusters, apply a small-sample correction (Mancl–DeRouen or Kauermann–Carroll).

library(geepack)
ppmi_g <- ppmi |> arrange(subject, time)

gee_ind <- geeglm(worsened ~ group * time + Age + Sex, id = subject,
                  data = ppmi_g, family = binomial, corstr = "independence")
gee_exc <- geeglm(worsened ~ group * time + Age + Sex, id = subject,
                  data = ppmi_g, family = binomial, corstr = "exchangeable")
gee_ar1 <- geeglm(worsened ~ group * time + Age + Sex, id = subject,
                  data = ppmi_g, family = binomial, corstr = "ar1")

term <- grep(":time", names(coef(gee_exc)), value = TRUE)[1]
data.frame(
  working_correlation = c("independence", "exchangeable", "AR(1)"),
  estimate = round(c(coef(gee_ind)[term], coef(gee_exc)[term], coef(gee_ar1)[term]), 4),
  robust_SE = round(c(summary(gee_ind)$coefficients[term, 2],
                      summary(gee_exc)$coefficients[term, 2],
                      summary(gee_ar1)$coefficients[term, 2]), 4),
  QIC = round(c(geepack::QIC(gee_ind)[1], geepack::QIC(gee_exc)[1],
                geepack::QIC(gee_ar1)[1]), 1))

The estimates are close across all three working structures, that is the robustness property in action. The standard errors differ modestly, which is the efficiency cost of a poor choice. QIC (the quasi-likelihood information criterion) is the GEE analogue of AIC and is the tool for comparing working correlations; ordinary AIC is unavailable because there is no likelihood.

# Model-based ("naive") vs. robust standard errors under a deliberately wrong
# working correlation
naive_se <- sqrt(diag(summary(gee_ind)$cov.unscaled))
robust_se <- summary(gee_ind)$coefficients[, 2]
data.frame(term = names(coef(gee_ind)),
           model_based_SE = round(naive_se, 4),
           robust_sandwich_SE = round(robust_se, 4),
           ratio = round(robust_se / naive_se, 3))

The working correlation here is independence, which is wrong by construction. The model-based standard errors inherit that error; the robust ones do not, and the ratio shows how much the correction matters for each term.

28 GEE versus GLMM

They estimate different quantities, and the difference is quantifiable.

\[ \begin{aligned} \textbf{GLMM (subject-specific): }&\quad \operatorname{logit}P(y_{ij}=1\mid b_i)=\mathbf x_{ij}^\top\boldsymbol\beta^{\mathrm{c}}+b_i\\ \textbf{GEE (population-averaged): }&\quad \operatorname{logit}P(y_{ij}=1)=\mathbf x_{ij}^\top\boldsymbol\beta^{\mathrm{m}} \end{aligned} \]

The GLMM coefficient answers “how does this subject’s odds change?”; the GEE coefficient answers “how do the odds change across the population?”

The attenuation relation. For a logit link with $b_i\sim N(0,\sigma_b^2)$, \[\boxed{\;\beta^{\mathrm{m}}\ \approx\ \frac{\beta^{\mathrm{c}}}{\sqrt{1+0.346\,\sigma_b^2}}\;}\] (Zeger, Liang & Albert, 1988).

The marginal effect is always smaller in magnitude, and the gap widens with the random-effect variance. At $\sigma_b=0$ they coincide; at $\sigma_b=2$ the marginal coefficient is about 65% of the conditional one. For an identity link there is no attenuation at all, the two coincide exactly, which is why the distinction never arises in the linear case.

set.seed(131)
attenuation_sim <- function(sigma_b, beta_c = 1.0, N = 400, n = 6) {
  b <- rnorm(N, sd = sigma_b)
  id <- rep(seq_len(N), each = n)
  x <- rnorm(N * n)
  y <- rbinom(N * n, 1, plogis(-0.5 + beta_c * x + b[id]))
  d <- data.frame(y, x, id = factor(id))
  gl <- suppressMessages(lme4::glmer(y ~ x + (1 | id), data = d, family = binomial,
                                     control = glmerControl(optimizer = "bobyqa")))
  ge <- geepack::geeglm(y ~ x, id = id, data = d[order(d$id), ],
                        family = binomial, corstr = "exchangeable")
  c(sigma_b = sigma_b,
    beta_conditional = unname(lme4::fixef(gl)["x"]),
    beta_marginal = unname(coef(ge)["x"]),
    predicted_marginal = unname(lme4::fixef(gl)["x"]) / sqrt(1 + 0.346 * sigma_b^2))
}
att <- as.data.frame(do.call(rbind, lapply(c(0.25, 0.75, 1.5, 2.5), attenuation_sim)))
att |> mutate(across(everything(), \(z) round(z, 4)))

att |> select(sigma_b, `GLMM (conditional)` = beta_conditional,
              `GEE (marginal)` = beta_marginal,
              `Predicted marginal` = predicted_marginal) |>
  pivot_longer(-sigma_b, names_to = "quantity", values_to = "beta") |>
  ggplot(aes(sigma_b, beta, color = quantity, linetype = quantity)) +
  geom_line(linewidth = 1) + geom_point(size = 2.2) +
  scale_color_manual(values = c("#3B7DD8", "#D8433B", "grey35")) +
  scale_linetype_manual(values = c("solid", "solid", "dashed")) +
  labs(title = "Marginal and conditional logistic effects diverge as the random-effect SD grows",
       subtitle = "Dashed grey: the attenuation formula applied to the fitted GLMM coefficient",
       x = expression(sigma[b]), y = "Coefficient on x", color = NULL, linetype = NULL) +
  theme_dspa()

sb_grid <- seq(0, 3, length.out = 50)
bc_grid <- seq(0, 2, length.out = 50)
Zatt <- outer(sb_grid, bc_grid, function(s, b) b / sqrt(1 + 0.346 * s^2))

plot_ly(x = bc_grid, y = sb_grid, z = Zatt, type = "surface",
        colorscale = "Viridis",
        colorbar = list(title = "Marginal\ncoefficient")) |>
  add_trace(x = bc_grid, y = rep(0, length(bc_grid)), z = bc_grid,
            type = "scatter3d", mode = "lines", name = "No attenuation",
            line = list(width = 7, color = "black")) |>
  layout(title = "Marginal effect implied by a conditional effect and the random-effect SD",
         scene = list(xaxis = list(title = "Conditional beta (GLMM)"),
                      yaxis = list(title = "Random-effect SD"),
                      zaxis = list(title = "Marginal beta (GEE)")))

The black ridge at $\sigma_b=0$ is the identity line, where the two coincide. The surface bends below it everywhere else, and the bend deepens with both the random-effect variance and the size of the effect itself.

29 Missing data

Longitudinal studies lose subjects. Which methods survive depends on why.

\[ \begin{aligned} \textbf{MCAR: }&\quad P(R\mid Y_{\text{obs}},Y_{\text{mis}})=P(R) &&\text{missingness unrelated to anything}\\ \textbf{MAR: }&\quad P(R\mid Y_{\text{obs}},Y_{\text{mis}})=P(R\mid Y_{\text{obs}}) &&\text{explained by observed data}\\ \textbf{MNAR: }&\quad \text{depends on }Y_{\text{mis}} &&\text{depends on the unseen values} \end{aligned} \]

This distinction decides between GEE and GLMM. Likelihood-based methods — LMM, GLMM, are valid under MAR, because the likelihood conditions on the observed data and the missingness mechanism factors out.

GEE requires MCAR. It is not likelihood-based, so the mechanism does not factor out, and dropout that depends on observed history biases the estimates. Weighted GEE (inverse-probability-of-censoring weights) restores validity under MAR, at the cost of having to model the dropout process.

This is often the deciding consideration in a clinical study, where dropout is almost never MCAR, sicker patients leave.

set.seed(141)
missing_sim <- function(mechanism, N = 250, n = 5, reps = 60) {
  est <- replicate(reps, {
    b <- rnorm(N, sd = 2); id <- rep(seq_len(N), each = n)
    tt <- rep(0:(n-1), N); x <- rep(rbinom(N, 1, 0.5), each = n)
    y <- 10 + 2*x + 0.8*tt + 1.2*x*tt + b[id] + rnorm(N*n)
    keep <- rep(TRUE, N*n)
    if (mechanism == "MCAR") keep <- runif(N*n) > 0.3
    if (mechanism == "MAR") {                      # dropout depends on PREVIOUS y
      prev <- c(NA, head(y, -1)); prev[tt == 0] <- NA
      keep <- is.na(prev) | runif(N*n) > plogis(-2 + 0.12*(prev - mean(y)))
    }
    d <- data.frame(y, x, tt, id = factor(id))[keep, ]
    lm_e <- tryCatch(lme4::fixef(suppressMessages(
      lme4::lmer(y ~ x*tt + (1|id), data = d)))["x:tt"], error = function(e) NA)
    ge_e <- tryCatch(coef(geepack::geeglm(y ~ x*tt, id = id,
                          data = d[order(d$id), ], corstr = "exchangeable"))["x:tt"],
                     error = function(e) NA)
    c(lmm = lm_e, gee = ge_e)
  })
  c(mechanism = mechanism, truth = 1.2,
    LMM = mean(est[1, ], na.rm = TRUE), GEE = mean(est[2, ], na.rm = TRUE))
}
mm <- as.data.frame(do.call(rbind, lapply(c("none", "MCAR", "MAR"), missing_sim)))
mm |> mutate(across(c(truth, LMM, GEE), \(z) round(as.numeric(z), 4)),
             LMM_bias = round(as.numeric(LMM) - 1.2, 4),
             GEE_bias = round(as.numeric(GEE) - 1.2, 4))

Under MCAR both methods recover the true interaction of 1.2. Under MAR — dropout driven by the previous observed outcome, the mixed model remains close while GEE drifts. That is the practical content of the distinction.

30 PART V: LATENT VARIABLES

31 Structural equation modeling

Many constructs, disease severity, cognitive function, socioeconomic status — are not measured directly. SEM models them as latent variables with multiple noisy indicators, and simultaneously estimates relationships among the latents.

\[ \begin{aligned} \textbf{Measurement model: }&\quad \mathbf x=\Lambda_x\boldsymbol\xi+\boldsymbol\delta,\qquad \mathbf y=\Lambda_y\boldsymbol\eta+\boldsymbol\epsilon\\ \textbf{Structural model: }&\quad \boldsymbol\eta=B\boldsymbol\eta+\Gamma\boldsymbol\xi+\boldsymbol\zeta \end{aligned} \]

with $\boldsymbol\xi$ exogenous latents, $\boldsymbol\eta$ endogenous latents, $\Lambda$ the loading matrices, $B$ the paths among endogenous latents, and $\Gamma$ the exogenous-to-endogenous paths.

31.1 The discrepancy function

SEM fits a covariance structure. The model implies $\Sigma(\boldsymbol\theta)$; the data give $S$; estimation minimizes their discrepancy. Under multivariate normality,

\[\boxed{\;F_{ML}(\boldsymbol\theta)=\ln\big|\Sigma(\boldsymbol\theta)\big|+\operatorname{tr}\!\big(S\,\Sigma(\boldsymbol\theta)^{-1}\big)-\ln|S|-p\;}\]

where $p$ is the number of observed variables. Two facts follow.

$F_{ML}\ge0$, with equality if and only if $\Sigma(\boldsymbol\theta)=S$. The final $-p$ is $-\operatorname{tr}(SS^{-1})=-\operatorname{tr}(I_p)$, the value the first two terms take at a perfect fit; subtracting it makes the function a proper discrepancy anchored at zero.

$(n-1)\hat F_{ML}\overset{d}{\to}\chi^2_{\mathrm{df}}$ under the null that the model is correct, with $\mathrm{df}=\frac{p(p+1)}{2}-\#\{\text{free parameters}\}$. This is where the SEM $\chi^2$ test and every fit index derived from it come from.

Common misconception: “a non-significant $\chi^2$ means the model is right.” The test’s power grows with $n$, so in large samples trivial misspecifications are rejected and in small samples serious ones are not. The $\chi^2$ is a badness-of-fit statistic that answers “is the discrepancy distinguishable from zero,” which is rarely the question of interest.

Report it alongside indices that behave differently: RMSEA $=\sqrt{\max(0,(\chi^2-\mathrm{df})/(\mathrm{df}(n-1)))}$, which rewards parsimony ($\le0.05$ close, $\le0.08$ reasonable); CFI, comparing against the independence model ($\ge0.95$); and SRMR, the standardized residual covariance ($\le0.08$).

31.2 Identification

A latent variable has no scale of its own, so one must be imposed, either by fixing one loading to 1 (marker variable) or by fixing the latent variance to 1 (standardized). Beyond scaling, the $t$-rule requires

\[\#\{\text{free parameters}\}\ \le\ \frac{p(p+1)}{2},\]

and for the measurement model each latent needs at least two indicators (three for a standalone factor to be identified without cross-model constraints), with uncorrelated measurement errors.

## OLD:
# library(lavaan)
# 
# # Build a wide-format frame: one row per subject, UPDRS parts at each visit
# sem_wide <- ppmi |>
#   filter(time %in% c(0, 12, 24)) |>
#   select(subject, group, Age, Sex, time,
#          UPDRS_part_I, UPDRS_part_II, UPDRS_part_III) |>
#   pivot_wider(names_from = time,
#               values_from = c(UPDRS_part_I, UPDRS_part_II, UPDRS_part_III),
#               names_sep = "_T") |>
#   filter(complete.cases(across(starts_with("UPDRS"))))
# c(subjects_complete = nrow(sem_wide), columns = ncol(sem_wide))
# 
# model_cfa <- '
#   # Measurement model: one severity factor per visit, marker-variable scaling
#   sev0 =~ UPDRS_part_I_T0 + UPDRS_part_II_T0 + UPDRS_part_III_T0
#   sev1 =~ UPDRS_part_I_T12 + UPDRS_part_II_T12 + UPDRS_part_III_T12
#   sev2 =~ UPDRS_part_I_T24 + UPDRS_part_II_T24 + UPDRS_part_III_T24
# 
#   # Structural model: severity carries forward
#   sev1 ~ sev0
#   sev2 ~ sev1
# '
# fit_cfa <- lavaan::sem(model_cfa, data = sem_wide, missing = "fiml")
# lavaan::fitMeasures(fit_cfa,
#   c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "rmsea.ci.upper", "srmr")) |>
#   round(4)
library(lavaan)

# Build wide-format frame without forcing complete cases
sem_wide <- ppmi |>
  filter(time %in% c(0, 12, 24)) |>
  select(subject, group, Age, Sex, time,
         UPDRS_part_I, UPDRS_part_II, UPDRS_part_III) |>
  pivot_wider(names_from = time,
              values_from = c(UPDRS_part_I, UPDRS_part_II, UPDRS_part_III),
              names_sep = "_T")

# Optional: check missingness
# colSums(is.na(sem_wide))

model_cfa <- '
  # Measurement model: one severity factor per visit, marker-variable scaling
  sev0 =~ UPDRS_part_I_T0 + UPDRS_part_II_T0 + UPDRS_part_III_T0
  sev1 =~ UPDRS_part_I_T12 + UPDRS_part_II_T12 + UPDRS_part_III_T12
  sev2 =~ UPDRS_part_I_T24 + UPDRS_part_II_T24 + UPDRS_part_III_T24

  # Structural model: severity carries forward
  sev1 ~ sev0
  sev2 ~ sev1
'

fit_cfa <- lavaan::sem(model_cfa, data = sem_wide, missing = "fiml")
lavaan::fitMeasures(fit_cfa,
  c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "rmsea.ci.upper", "srmr")) |>
  round(4)

#>          chisq             df         pvalue            cfi            tli 
#>        112.449         25.000          0.000          0.606          0.433 
#>          rmsea rmsea.ci.upper           srmr 
#>          0.125          0.149          0.208

# Loadings are EXTRACTED from the fit, never transcribed
pe <- lavaan::parameterEstimates(fit_cfa, standardized = TRUE)
pe |> filter(op %in% c("=~", "~")) |>
  select(lhs, op, rhs, est, se, z, pvalue, std.all) |>
  mutate(across(where(is.numeric), \(z) round(z, 4))) |>
  head(12)

32 Latent growth curves

Common misconception: “build a composite from the estimated loadings, then model it.” Forming $\text{score}=y_1+\hat\lambda_2y_2+\hat\lambda_3y_3$ and passing it to a mixed model as an observed outcome makes three errors at once. It discards the measurement error the latent variable was introduced to represent; it treats estimated loadings as known constants, so the downstream standard errors omit their sampling variability and are too small; and it reuses the same data for both stages.

The latent growth curve model fits the measurement and structural parts jointly: an intercept factor with all loadings fixed at 1, a slope factor with loadings fixed at the observation times, and their variances and covariance freely estimated. Every source of uncertainty is propagated.

model_lgc <- '
  # Intercept: loadings fixed at 1 -- the baseline level
  i =~ 1*UPDRS_part_III_T0 + 1*UPDRS_part_III_T12 + 1*UPDRS_part_III_T24
  # Slope: loadings fixed at the observation times (in years)
  s =~ 0*UPDRS_part_III_T0 + 1*UPDRS_part_III_T12 + 2*UPDRS_part_III_T24

  # Both growth factors regressed on covariates
  i ~ group + Age
  s ~ group + Age

  # Intercept and slope are allowed to covary
  i ~~ s
'
fit_lgc <- dspa_try(lavaan::growth(model_lgc, data = sem_wide, missing = "fiml"),
                    label = "latent growth curve")

if (!is.null(fit_lgc)) {
  lavaan::fitMeasures(fit_lgc, c("chisq", "df", "pvalue", "cfi", "rmsea", "srmr")) |>
    round(4)
  lavaan::parameterEstimates(fit_lgc) |>
    filter(op %in% c("~", "~~"), lhs %in% c("i", "s")) |>
    select(lhs, op, rhs, est, se, z, pvalue) |>
    mutate(across(where(is.numeric), \(z) round(z, 4)))
}

Read the s ~ group row: that coefficient is the difference in rate of change between groups, estimated jointly with the measurement model, the quantity a two-stage composite would estimate with understated uncertainty.

The i ~~ s covariance is substantively interesting in its own right. A negative value means subjects starting higher decline more slowly (or regress toward the mean); a positive one means the gap widens over time.

if (!is.null(fit_lgc)) {
  fs <- lavaan::lavPredict(fit_lgc)
  fs_df <- data.frame(intercept = fs[, "i"], slope = fs[, "s"],
                      group = sem_wide$group[seq_len(nrow(fs))])
  ggplot(fs_df, aes(intercept, slope, color = group)) +
    geom_hline(yintercept = 0, color = "grey70") +
    geom_point(size = 1.8, alpha = 0.7) +
    geom_smooth(method = "lm", se = FALSE, linewidth = 0.8) +
    scale_color_manual(values = c("#3B7DD8", "#D8433B")) +
    labs(title = "Estimated growth factors: baseline level against rate of change",
         subtitle = "Each point is a subject. The slope of the fitted line is the intercept-slope covariance",
         x = "Intercept factor (baseline severity)",
         y = "Slope factor (change per year)", color = NULL) +
    theme_dspa()
}

33 PART VI: SEQUENCE NEURAL NETWORKS

34 Recurrent networks

A recurrent network carries a hidden state across time steps:

\[\mathbf h_t=\sigma\big(W_{hh}\mathbf h_{t-1}+W_{xh}\mathbf x_t+\mathbf b\big), \qquad \hat y_t=W_{hy}\mathbf h_t+\mathbf c\]

The same weights are applied at every step, parameter sharing across time, exactly analogous to the weight sharing across space in a convolutional network (Chapter 6, §6.11), so the model handles sequences of any length with a fixed parameter count.

34.1 Why gradients vanish

Training uses backpropagation through time, and the gradient of a loss at step $T$ with respect to the state at step $t$ is a product of Jacobians:

\[\frac{\partial\mathcal L_T}{\partial\mathbf h_t}=\frac{\partial\mathcal L_T}{\partial\mathbf h_T}\prod_{k=t+1}^{T}\frac{\partial\mathbf h_k}{\partial\mathbf h_{k-1}}, \qquad \frac{\partial\mathbf h_k}{\partial\mathbf h_{k-1}}=W_{hh}^\top\operatorname{diag}\big(\sigma'(\cdot)\big)\]

The product of $T-t$ Jacobians decays or explodes geometrically. Bounding the norm, \[\left\|\prod_{k=t+1}^{T}\frac{\partial\mathbf h_k}{\partial\mathbf h_{k-1}}\right\|\ \le\ \big(\|W_{hh}\|\cdot\gamma\big)^{T-t}, \qquad \gamma=\sup|\sigma'|,\] so if $\|W_{hh}\|\gamma<1$ the gradient vanishes exponentially in the lag and the network cannot learn long-range dependence; if $\|W_{hh}\|\gamma>1$ it explodes.

For $\tanh$, $\gamma=1$; for the logistic sigmoid, $\gamma=1/4$, which makes vanishing nearly certain. Exploding gradients are easily handled by clipping; vanishing gradients are the hard problem, and they are what LSTMs were designed to solve.

set.seed(151)
grad_norm <- function(spectral_radius, T_lag, d = 20, gamma = 1) {
  W <- matrix(rnorm(d * d), d, d)
  W <- W / max(abs(eigen(W, only.values = TRUE)$values)) * spectral_radius
  J <- diag(d)
  vapply(seq_len(T_lag), function(k) { J <<- J %*% (t(W) * gamma); norm(J, "2") },
         numeric(1))
}
lags <- 1:40
bind_rows(lapply(c(0.7, 0.95, 1.0, 1.05), function(r)
  data.frame(lag = lags, norm = grad_norm(r, max(lags)),
             radius = sprintf("spectral radius = %.2f", r)))) |>
  ggplot(aes(lag, norm, color = radius)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "grey45") +
  geom_line(linewidth = 1) +
  scale_y_log10() +
  scale_color_viridis_d(option = "plasma", end = 0.9) +
  labs(title = "Gradient norm through backpropagation over time",
       subtitle = "Below the dashed line the signal vanishes; above it, explodes. The decay is geometric in the lag",
       x = "Lag (time steps back)", y = "Gradient norm (log scale)", color = NULL) +
  theme_dspa()

# --- Interactive equivalent ------------------------------------------------
p <- plot_ly()
for (r in c(0.7, 0.95, 1.0, 1.05))
  p <- add_lines(p, x = lags, y = grad_norm(r, max(lags)),
                 name = sprintf("radius = %.2f", r))
p |> layout(title = "Gradient norm vs. backpropagation lag",
            xaxis = list(title = "Lag"),
            yaxis = list(title = "Gradient norm", type = "log"))

At spectral radius 0.7 the gradient has fallen below $10^{-5}$ by lag 40, the network is effectively blind to anything further back than about 15 steps.

35 LSTM and GRU

The long short-term memory cell (Hochreiter & Schmidhuber, 1997) adds an explicit cell state $\mathbf c_t$ with an additive update:

\[ \begin{aligned} \mathbf f_t&=\sigma\big(W_f[\mathbf h_{t-1},\mathbf x_t]+\mathbf b_f\big) &&\text{forget gate}\\ \mathbf i_t&=\sigma\big(W_i[\mathbf h_{t-1},\mathbf x_t]+\mathbf b_i\big) &&\text{input gate}\\ \tilde{\mathbf c}_t&=\tanh\big(W_c[\mathbf h_{t-1},\mathbf x_t]+\mathbf b_c\big) &&\text{candidate}\\ \mathbf c_t&=\mathbf f_t\odot\mathbf c_{t-1}+\mathbf i_t\odot\tilde{\mathbf c}_t &&\textbf{additive update}\\ \mathbf o_t&=\sigma\big(W_o[\mathbf h_{t-1},\mathbf x_t]+\mathbf b_o\big) &&\text{output gate}\\ \mathbf h_t&=\mathbf o_t\odot\tanh(\mathbf c_t) \end{aligned} \]

The fourth line is the whole idea. Because $\partial\mathbf c_t/\partial\mathbf c_{t-1}=\operatorname{diag}(\mathbf f_t)$, the gradient along the cell state is multiplied by the forget gate rather than by a weight matrix and a derivative. When $\mathbf f_t\approx1$ the gradient passes through essentially unattenuated, a “constant error carousel”, so information can be carried hundreds of steps without decay.

The GRU achieves the same effect with two gates instead of three and about 25% fewer parameters, usually at comparable accuracy.

Architecture	Parameters (hidden size $d$, input $m$)	Cost per step
Vanilla RNN	$d^2+dm+d$	$O(d^2+dm)$
GRU	$3(d^2+dm+d)$	$O(3(d^2+dm))$
LSTM	$4(d^2+dm+d)$	$O(4(d^2+dm))$

Training a sequence of length $T$ for $E$ epochs over $N$ sequences costs $O(E\,N\,T\,d^2)$, linear in sequence length, which is why RNNs scale to long sequences where a Transformer’s $O(T^2)$ attention does not, and why they remain relevant for streaming and embedded settings.

# --- LSTM forecasting with keras3 (set KERAS_EVAL = TRUE to run) -----------
library(keras3)

make_windows <- function(x, lookback, horizon = 1) {
  n <- length(x) - lookback - horizon + 1
  X <- t(vapply(seq_len(n), function(i) x[i:(i + lookback - 1)], numeric(lookback)))
  y <- x[(lookback + horizon):length(x)]
  list(X = array(X, dim = c(nrow(X), lookback, 1)), y = y)
}

z <- as.numeric(scale(train))                       # LSTMs need scaled inputs
w <- make_windows(z, lookback = 28)

model <- keras_model_sequential() |>
  layer_lstm(units = 32, input_shape = c(28, 1)) |>
  layer_dropout(0.2) |>
  layer_dense(units = 1)
model |> compile(optimizer = optimizer_adam(1e-3), loss = "mse")
model |> fit(w$X, w$y, epochs = 30, batch_size = 64,
             validation_split = 0.2, verbose = 0)

# Chronological split throughout: never sample() a temporal index
c(train_windows = dim(w$X)[1], lookback = 28)

The LSTM section is written against keras3 and gated by KERAS_EVAL so the chapter renders without a configured Python environment while the code remains visible.

36 When do neural methods beat ARIMA?

Common misconception: “deep learning beats classical methods on time series.” In the M4 and M5 forecasting competitions, pure neural models did not win; the leading entries were statistical methods, and hybrids that embedded exponential smoothing inside a neural architecture. On a single univariate series of a few hundred points, ARIMA and ETS are usually better and far cheaper.

Neural methods earn their cost under specific conditions: many related series to learn shared structure from, long-range dependence beyond what a parsimonious ARMA captures, many exogenous inputs, or genuine non-linearity. Absent those, the extra capacity buys variance, not accuracy.

Situation	Prefer
One short series, clear seasonality	ARIMA / ETS
Long series, non-linear dynamics	LSTM / GRU
Thousands of related series	Global neural model (DeepAR, N-BEATS)
Many exogenous covariates	Dynamic regression or neural
Interpretability required	ARIMA / dynamic regression
Prediction intervals required	ARIMA gives them analytically; neural needs quantile loss or bootstrapping

# A like-for-like comparison: same origin, same horizon, same scale-free metric,
# with the naive benchmark included so the numbers can be read.
set.seed(161)
comp_h <- min(30, length(test))

# Create list of forecasts
forecast_list <- list(
  `ARIMA (auto)`   = forecast::forecast(fit_auto, h = comp_h),
  `ETS`            = forecast::forecast(forecast::ets(train), h = comp_h),
  `Seasonal naive` = forecast::snaive(train, h = comp_h),
  `TBATS`          = dspa_try(forecast::forecast(forecast::tbats(train), h = comp_h),
                              fallback = forecast::snaive(train, h = comp_h),
                              label = "TBATS")
)

# Compute accuracy for each and combine
bench <- do.call(rbind, lapply(names(forecast_list), function(nm) {
  a <- forecast::accuracy(forecast_list[[nm]], head(test, comp_h))
  data.frame(model = nm,
             MASE = round(a["Test set", "MASE"], 4),
             RMSE = round(a["Test set", "RMSE"], 3))
}))

bench

Every model is scored at the same origin, over the same horizon, with the same scale-free metric, against the naive benchmark. That is what makes the ordering meaningful, and it is the comparison any neural forecaster must also submit to.

37 PART VII: SYNTHESIS

38 Computational complexity summary

$T$ = series length, $N$ = subjects, $n$ = observations per subject, $p,q$ = ARMA orders, $m$ = seasonal period, $d$ = hidden size, $E$ = epochs, $q_r$ = random-effect dimension, $k$ = smoothing window.

Task	Cost	Note
Sample ACF to lag $L$	$O(TL)$, or $O(T\log T)$ via FFT	`acf()` uses FFT
Differencing	$O(T)$	—
ADF / KPSS	$O(T)$	Plus lag-order selection
ARIMA likelihood (Kalman)	$O\big(T(p+q+m)^2\big)$	Handles gaps exactly
ARIMA fit ($I$ optimizer steps)	$O\big(I\,T(p+q+m)^2\big)$	—
`auto.arima` stepwise	$\approx O(p_{\max}+q_{\max})$ fits	Default; can miss the optimum
`auto.arima` exhaustive	$O(p_{\max}q_{\max}P_{\max}Q_{\max})$ fits	450 fits at $5{,}5{,}2{,}2$
STL decomposition	$O(T)$ per inner iteration	—
Rolling-origin CV	$O(\#\text{origins}\times\text{fit})$	The honest cost of honest evaluation
Forecast $h$ steps	$O(h(p+q))$	Intervals need $\psi$-weights: $O(h)$
LMM via REML	$O\big(N n q_r^2+q_r^3\big)$ per iteration	Sparse $Z$ makes this near-linear in $N$
GLMM, Laplace	$O(N n q_r^2)$ per iteration	One quadrature point
GLMM, adaptive GH ($Q$ points)	$O\big(N n Q^{q_r}\big)$	Exponential in $q_r$
GEE	$O\big(N n^3\big)$ per iteration	$V_i^{-1}$ dominates
Sandwich variance	$O(N n^2 p^2)$	One pass after convergence
SEM (ML, $p$ indicators)	$O\big(I p^3\big)$	$\Sigma^{-1}$ and $\|\Sigma\|$ each iteration
RNN / GRU / LSTM training	$O(E\,N\,T\,d^2)$	Linear in $T$
Transformer self-attention	$O(E\,N\,T^2 d)$	Quadratic in $T$

Four consequences.

The Kalman filter is what makes ARIMA practical on gappy data. It evaluates the exact likelihood in $O(T)$ and skips the update step at missing times, no imputation required for estimation, though imputation is still needed for the ACF diagnostics that precede it.

Adaptive quadrature is exponential in the random-effect dimension. A GLMM with a random intercept ($q_r=1$) can use many quadrature points cheaply; with random intercept and slope ($q_r=2$) the cost is $Q^2$, and beyond three random effects the Laplace approximation is effectively the only option.

GEE is cubic in the cluster size. With few observations per subject this is irrelevant; with dense repeated measures ($n$ in the hundreds) it dominates, and a structured working correlation with a closed-form inverse, AR(1), exchangeable, becomes necessary rather than merely convenient.

RNNs are linear in sequence length where attention is quadratic. For a sequence of 10,000 steps that is a factor of 10,000, which is why recurrent architectures remain the practical choice for long streaming signals.

39 Common pitfalls

#	Pitfall	Consequence	Fix
1	Fitting ARIMA to a smoothed series	Fabricates autocorrelation; intervals too narrow by $\sqrt k$	Smooth for display only; model the data
2	Mean-imputing a time series	Flat runs; ACF biased toward zero	`imputeTS::na_kalman()`
3	Introducing a new sentinel after cleaning the first	The sentinel enters the likelihood as data	Convert to `NA` once
4	ACF bounds from a different series’ length	Every significance call mis-scaled	$\pm1.96/\sqrt{T}$ with $T$ from this series
5	Reading residual ACF against $\pm1.96/\sqrt T$	Ignores the $p+q$ fitted parameters	Ljung–Box with `fitdf = p + q`
6	“ADF rejected, so it’s stationary”	Rejecting a unit root is not proving stationarity	Run KPSS too; read the four-way table
7	Over-differencing	$\rho(1)\to-0.5$; variance inflated; non-invertible	`ndiffs()`; $d\le2$
8	Long MA order in place of seasonal terms	24 parameters where 3 suffice; no generalization past one cycle	SARIMA $(p,d,q)(P,D,Q)_m$
9	Plotting fitted values as “observed”	Hides how the forecast compares to reality	Plot the held-out actuals
10	Forecast bands read from the wrong column	Bands labelled 95% that are 80% or 20%	Match columns to `level =`
11	Forecasting without a benchmark	An accuracy number with no scale	Naive, seasonal naive, drift, mean
12	Correlation as a forecast metric	Affine-invariant; $\hat y=100+0.01y$ scores 1	MASE, RMSSE, RMSE
13	One train/test split	High variance; origin-dependent	Rolling-origin CV
14	Random split of temporal data	Trains on the future	Chronological or rolling origin
15	Future covariates copied from the start of the series	Trend and indicators phase-shifted	Build covariates for the forecast window
16	Error bars from predictions rather than `$se`	No statistical content	`predict()$se`, or `forecast()` intervals
17	A simulation whose signal is overwritten	Nothing to recover; the study is vacuous	Verify the effect is present before fitting
18	Prophet weekly seasonality on monthly data	Estimating a pattern the data cannot contain	Match seasonality and `freq` to the sampling
19	Treating a relative index as an absolute measure	Google Trends is normalized per request	State the caveat; avoid cross-request comparison
20	Ignoring clustering in repeated measures	SEs too small by $\sqrt{1+(n-1)\rho}$	Mixed model or GEE
21	Time as a factor in a trend analysis	Estimates contrasts, not a slope	Keep time numeric
22	Composite from estimated loadings, then modelled	Discards measurement error; SEs too small	Latent growth curve
23	Comparing REML fits with different fixed effects	REML likelihoods are not comparable	ML for fixed effects; REML for random
24	GEE under MAR dropout	Requires MCAR; biased under MAR	LMM/GLMM, or weighted GEE

40 Practice problems

40.1 Problem 1: Quantify the smoothing trap

Measure how much a $k$-point moving average distorts the estimated AR(1) coefficient and the innovation variance.

Solution

set.seed(201)
truth_phi <- 0.5; truth_sd <- 1
p1 <- do.call(rbind, lapply(c(1, 2, 5, 10, 30, 60), function(k) {
  est <- replicate(30, {
    x <- as.numeric(arima.sim(list(ar = truth_phi), 1500, sd = truth_sd))
    s <- ma_filter(x, k); s <- s[!is.na(s)]
    f <- tryCatch(forecast::Arima(s, order = c(1, 0, 0)), error = function(e) NULL)
    if (is.null(f)) c(NA, NA) else c(coef(f)["ar1"], sqrt(f$sigma2))
  })
  data.frame(k = k, phi_hat = mean(est[1, ], na.rm = TRUE),
             sigma_hat = mean(est[2, ], na.rm = TRUE))
}))
p1 |> mutate(phi_bias = round(phi_hat - truth_phi, 4),
             sigma_ratio = round(sigma_hat / truth_sd, 4),
             predicted_ratio = round(1/sqrt(k), 4),
             across(c(phi_hat, sigma_hat), \(z) round(z, 4)))

p1 |> ggplot(aes(k, phi_hat)) +
  geom_hline(yintercept = truth_phi, linetype = "dashed", color = "firebrick") +
  geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.4) +
  scale_x_log10(breaks = p1$k) +
  labs(title = "Estimated AR(1) coefficient against smoothing window",
       subtitle = "Dashed: the truth. Smoothing drives the estimate toward 1 regardless of the data",
       x = "Smoothing window k (log scale)", y = expression(hat(phi))) +
  theme_dspa()

The estimated $\phi$ climbs toward 1 as $k$ grows, the smoothed series looks ever more persistent, while $\hat\sigma$ falls roughly as $1/\sqrt k$, matching the predicted_ratio column. Both distortions push prediction intervals in the same dangerous direction: more persistence and less noise means narrower intervals for a series that has not changed.

40.2 Problem 2: ADF and KPSS across the unit circle

Trace both tests’ rejection rates as $\phi$ approaches 1.

Solution

set.seed(203)
p2 <- do.call(rbind, lapply(c(0.5, 0.8, 0.9, 0.95, 0.99, 1.0), function(phi) {
  r <- replicate(60, {
    x <- if (phi == 1) cumsum(rnorm(250)) else as.numeric(arima.sim(list(ar = phi), 250))
    c(adf = suppressWarnings(tseries::adf.test(x)$p.value) < 0.05,
      kpss = suppressWarnings(tseries::kpss.test(x)$p.value) < 0.05)
  })
  data.frame(phi = phi, adf_rejects = mean(r["adf", ]),
             kpss_rejects = mean(r["kpss", ]))
}))
p2 |> mutate(across(-phi, \(z) round(z, 3)),
             both_agree_stationary = adf_rejects > 0.8 & kpss_rejects < 0.2)

p2 |> pivot_longer(-phi, names_to = "test", values_to = "rate") |>
  filter(test != "both_agree_stationary") |>
  ggplot(aes(phi, rate, color = test)) +
  geom_line(linewidth = 1) + geom_point(size = 2.4) +
  scale_color_manual(values = c(adf_rejects = "#3B7DD8", kpss_rejects = "#D8433B"),
                      labels = c("ADF rejects unit root", "KPSS rejects stationarity")) +
  labs(title = "Rejection rates as the AR root approaches the unit circle",
       subtitle = "ADF power collapses well before phi = 1; the two tests disagree in the intermediate zone",
       x = expression(phi), y = "Rejection rate", color = NULL) +
  theme_dspa()

ADF power collapses from $\phi\approx0.9$ onward, at $n=250$ it cannot distinguish $\phi=0.95$ from a unit root. Between roughly 0.9 and 1.0 the two tests disagree, and that disagreement is itself the finding: the data do not determine the answer.

40.3 Problem 3: Verify the interval-width formula

Confirm that empirical forecast-error variance matches $\sigma^2\sum\psi_j^2$.

Solution

set.seed(205)
check_intervals <- function(phi = 0.7, sigma = 1, T0 = 300, h_max = 8, reps = 3000) {
  psi <- c(1, ARMAtoMA(ar = phi, lag.max = h_max - 1))
  theory <- sigma^2 * cumsum(psi^2)
  err <- replicate(reps, {
    x <- as.numeric(arima.sim(list(ar = phi), T0 + h_max, sd = sigma))
    xt <- x[1:T0]; fut <- x[(T0+1):(T0+h_max)]
    f <- forecast::Arima(xt, order = c(1, 0, 0), include.mean = FALSE)
    as.numeric(forecast::forecast(f, h = h_max)$mean) - fut
  })
  data.frame(h = 1:h_max,
             empirical_var = round(apply(err, 1, var), 4),
             theoretical_var = round(theory, 4),
             ratio = round(apply(err, 1, var) / theory, 4))
}
check_intervals()

Empirical and theoretical variances agree to within simulation error at every horizon, confirming $\operatorname{Var}(X_{T+h}\mid\mathcal F_T)=\sigma^2\sum_{j<h}\psi_j^2$. The small excess at large $h$ reflects parameter-estimation uncertainty, which the formula treats as zero.

40.4 Problem 4: Design effect and effective sample size

Measure the Type I error of OLS on clustered data across $\rho$ and $n$.

Solution

set.seed(207)
p4 <- expand.grid(rho = c(0.1, 0.3, 0.5, 0.7), n = c(3, 6, 12))
p4$type_I <- vapply(seq_len(nrow(p4)), function(i) {
  rho <- p4$rho[i]; n <- p4$n[i]; N <- 50
  mean(replicate(250, {
    b <- rnorm(N, sd = sqrt(rho)); id <- rep(seq_len(N), each = n)
    x <- rep(rnorm(N), each = n)
    y <- b[id] + rnorm(N * n, sd = sqrt(1 - rho))
    summary(lm(y ~ x))$coefficients["x", 4] < 0.05
  }))
}, numeric(1))
p4$DEFF <- 1 + (p4$n - 1) * p4$rho
p4 |> mutate(across(c(type_I, DEFF), \(z) round(z, 3))) |> arrange(DEFF)

ggplot(p4, aes(DEFF, type_I, color = factor(n))) +
  geom_hline(yintercept = 0.05, linetype = "dashed", color = "firebrick") +
  geom_line(linewidth = 1) + geom_point(size = 2.4) +
  scale_color_viridis_d(option = "plasma", end = 0.85, name = "Visits per subject") +
  labs(title = "OLS Type I error against the design effect",
       subtitle = "Dashed: the nominal 5% rate. The inflation is a function of DEFF, not of rho or n separately",
       x = expression(DEFF == 1 + (n-1)*rho), y = "Type I error rate") +
  theme_dspa()

The points collapse onto a single curve in DEFF, confirming that $1+(n-1)\rho$ is the right summary. At DEFF near 8 the nominal 5% test rejects several times too often.

40.5 Problem 5: Sandwich robustness

Show that GEE standard errors remain valid under a misspecified working correlation while model-based ones do not.

Solution

set.seed(209)
p5 <- do.call(rbind, lapply(c("independence", "exchangeable", "ar1"), function(cs) {
  cov_rob <- cov_nai <- logical(0)
  for (r in 1:120) {
    N <- 80; n <- 5; rho <- 0.6
    Sig <- rho^abs(outer(1:n, 1:n, "-"))            # truth is AR(1)
    e <- MASS::mvrnorm(N, rep(0, n), Sig)
    id <- rep(seq_len(N), each = n)
    x  <- rep(rnorm(N), each = n)
    y  <- as.vector(t(e)) + 0.5 * x                 # true beta = 0.5
    d  <- data.frame(y, x, id = factor(id))
    g  <- tryCatch(geepack::geeglm(y ~ x, id = id, data = d, corstr = cs),
                   error = function(e) NULL)
    if (is.null(g)) next
    b <- coef(g)["x"]
    se_r <- summary(g)$coefficients["x", 2]
    se_n <- sqrt(diag(summary(g)$cov.unscaled))["x"]
    cov_rob <- c(cov_rob, abs(b - 0.5) < 1.96 * se_r)
    cov_nai <- c(cov_nai, abs(b - 0.5) < 1.96 * se_n)
  }
  data.frame(working_correlation = cs,
             coverage_robust = round(mean(cov_rob), 3),
             coverage_model_based = round(mean(cov_nai), 3), nominal = 0.95)
}))
p5

Robust coverage is near 0.95 for every working correlation, including the deliberately wrong independence, that is the sandwich property. Model-based coverage degrades when the working structure is wrong, because those standard errors take the assumed $V_i$ at face value.

40.6 Problem 6: Attenuation across link functions

Confirm that the marginal/conditional gap exists for the logit link and vanishes for the identity link.

Solution

set.seed(211)
p6 <- do.call(rbind, lapply(c(0.5, 1.5, 2.5), function(sb) {
  N <- 300; n <- 6; bc <- 1.0
  b <- rnorm(N, sd = sb); id <- rep(seq_len(N), each = n); x <- rnorm(N * n)
  # Logit link
  yb <- rbinom(N * n, 1, plogis(bc * x + b[id]))
  db <- data.frame(y = yb, x, id = factor(id))
  gl_b <- suppressMessages(lme4::glmer(y ~ x + (1|id), data = db, family = binomial,
                                       control = glmerControl(optimizer = "bobyqa")))
  ge_b <- geepack::geeglm(y ~ x, id = id, data = db[order(db$id), ],
                          family = binomial, corstr = "exchangeable")
  # Identity link
  yg <- bc * x + b[id] + rnorm(N * n)
  dg <- data.frame(y = yg, x, id = factor(id))
  gl_g <- suppressMessages(lme4::lmer(y ~ x + (1|id), data = dg))
  ge_g <- geepack::geeglm(y ~ x, id = id, data = dg[order(dg$id), ],
                          corstr = "exchangeable")
  data.frame(sigma_b = sb,
             logit_conditional = round(unname(lme4::fixef(gl_b)["x"]), 4),
             logit_marginal = round(unname(coef(ge_b)["x"]), 4),
             identity_conditional = round(unname(lme4::fixef(gl_g)["x"]), 4),
             identity_marginal = round(unname(coef(ge_g)["x"]), 4))
}))
p6 |> mutate(logit_ratio = round(logit_marginal / logit_conditional, 3),
             identity_ratio = round(identity_marginal / identity_conditional, 3),
             predicted_logit_ratio = round(1/sqrt(1 + 0.346*sigma_b^2), 3))

For the logit link the ratio tracks the predicted $1/\sqrt{1+0.346\sigma_b^2}$ and falls as $\sigma_b$ grows. For the identity link the ratio is 1 throughout: the linear model’s marginal and conditional effects coincide exactly, because the expectation of a linear function of the random effect is the function of its expectation. Non-collapsibility is a property of the link, not of the clustering.

40.7 Problem 7: Two-stage composite versus latent growth

Compare standard errors from a transcribed-loadings composite against a joint latent growth model.

Solution

set.seed(213)
N7 <- 250
b_i <- rnorm(N7, 0, 2); b_s <- rnorm(N7, 0, 0.5)
grp <- rbinom(N7, 1, 0.5)
mk <- function(t) {
  eta <- 10 + b_i + 2*grp + (0.5 + b_s + 0.8*grp) * t
  cbind(1.0*eta + rnorm(N7, 0, 2), 0.8*eta + rnorm(N7, 0, 2), 1.2*eta + rnorm(N7, 0, 2))
}
w0 <- mk(0); w1 <- mk(1); w2 <- mk(2)
d7 <- data.frame(group = grp,
                 y1_T0 = w0[,1], y2_T0 = w0[,2], y3_T0 = w0[,3],
                 y1_T1 = w1[,1], y2_T1 = w1[,2], y3_T1 = w1[,3],
                 y1_T2 = w2[,1], y2_T2 = w2[,2], y3_T2 = w2[,3])

# Stage 1: estimate loadings; Stage 2: build a composite and model it
cfa1 <- lavaan::cfa('f =~ y1_T0 + y2_T0 + y3_T0', data = d7)
ld <- lavaan::parameterEstimates(cfa1) |> filter(op == "=~") |> pull(est)
comp <- function(a, b, c_) ld[1]*a + ld[2]*b + ld[3]*c_
long7 <- data.frame(
  id = rep(seq_len(N7), 3), group = rep(grp, 3), t = rep(0:2, each = N7),
  y = c(comp(d7$y1_T0, d7$y2_T0, d7$y3_T0),
        comp(d7$y1_T1, d7$y2_T1, d7$y3_T1),
        comp(d7$y1_T2, d7$y2_T2, d7$y3_T2)))
m_two <- suppressMessages(lme4::lmer(y ~ group * t + (1 + t | id), data = long7))
se_two <- sqrt(diag(vcov(m_two)))["group:t"]

# Joint latent growth model
lgc7 <- '
  i =~ 1*y1_T0 + 1*y1_T1 + 1*y1_T2
  s =~ 0*y1_T0 + 1*y1_T1 + 2*y1_T2
  i ~ group
  s ~ group
  i ~~ s
'
f7 <- lavaan::growth(lgc7, data = d7)
pe7 <- lavaan::parameterEstimates(f7) |> filter(lhs == "s", op == "~", rhs == "group")

data.frame(
  approach = c("Two-stage composite (LMM)", "Joint latent growth curve"),
  estimate = round(c(unname(lme4::fixef(m_two)["group:t"]), pe7$est), 4),
  std_error = round(c(unname(se_two), pe7$se), 4),
  note = c("loadings treated as known constants",
           "measurement and structure fitted jointly"))

The two-stage estimate is on a different scale (the composite absorbs the loadings), so compare the relative precision. The joint model propagates the uncertainty in the loadings; the two-stage approach treats them as known, and its standard error omits that source entirely.

40.8 Problem 8: A fair neural-versus-classical comparison

Compare a simple neural forecaster against ARIMA and the naive benchmark at matched origins.

Solution

set.seed(215)
# A deliberately NON-LINEAR series: a threshold autoregression that no linear
# ARMA can represent
T8 <- 900
x8 <- numeric(T8)
for (t in 3:T8)
  x8[t] <- if (x8[t-1] > 0) 0.6*x8[t-1] - 0.3*x8[t-2] + rnorm(1) else
                            -0.5*x8[t-1] + 0.2*x8[t-2] + rnorm(1)
s8 <- ts(x8)
tr8 <- window(s8, end = 800); te8 <- window(s8, start = 801)

# Nonlinear series (first part)
# A small feed-forward autoregressive network -- forecast's nnetar, no Python
fit_nn <- forecast::nnetar(tr8, p = 5, size = 6, repeats = 15)
fit_ar8 <- forecast::auto.arima(tr8, seasonal = FALSE)

h8 <- length(te8)

forecast_list <- list(
  `ARIMA` = forecast::forecast(fit_ar8, h = h8),
  `Neural AR (nnetar)` = forecast::forecast(fit_nn, h = h8),
  `Naive` = forecast::naive(tr8, h = h8)
)

p8 <- do.call(rbind, lapply(names(forecast_list), function(nm) {
  a <- forecast::accuracy(forecast_list[[nm]], te8)
  data.frame(model = nm,
             MASE = round(a["Test set", "MASE"], 4),
             RMSE = round(a["Test set", "RMSE"], 4))
}))

p8 |> arrange(MASE)

# The same comparison on a LINEAR series
lin8 <- ts(as.numeric(arima.sim(list(ar = c(0.6, -0.3)), T8)))
trl <- window(lin8, end = 800)
tel <- window(lin8, start = 801)

forecast_list_linear <- list(
  `ARIMA` = forecast::forecast(forecast::auto.arima(trl, seasonal = FALSE), h = 100),
  `Neural AR (nnetar)` = forecast::forecast(
    forecast::nnetar(trl, p = 5, size = 6, repeats = 15), h = 100),
  `Naive` = forecast::naive(trl, h = 100)
)

results_linear <- do.call(rbind, lapply(names(forecast_list_linear), function(nm) {
  a <- forecast::accuracy(forecast_list_linear[[nm]], tel)
  data.frame(series = "linear AR(2)",
             model = nm,
             MASE = round(a["Test set", "MASE"], 4))
}))

results_linear

On the threshold series the neural model has non-linearity to exploit and can beat ARIMA. On the linear series it does not, and the extra capacity costs accuracy. That contrast is the entire answer to “should I use a neural forecaster”, and it is only visible because both are scored against the same benchmark at the same origin.

41 Checkpoint

A colleague smooths a noisy series with a 50-point moving average “to reveal the signal,” then fits ARIMA and reports very tight forecast intervals. What are the two separate problems?
Your ADF test gives $p = 0.31$. Is the series non-stationary?
You have 80 patients with 12 visits each and an ICC of 0.6. How much information do you actually have?
GEE and GLMM give logistic coefficients of 0.62 and 0.95 for the same predictor. Which is wrong?
Your study has substantial dropout, and sicker patients leave earlier. Which analysis method is safe?
Your LSTM achieves RMSE 4.2 on held-out data where ARIMA achieves 4.6. Has the LSTM won?

Answers

Fabricated autocorrelation, collapsed intervals, and their compounding. A $k$-point moving average of white noise has $\operatorname{Corr}(m_t, m_{t-h}) = (1-|h|/k)_+$, so the ACF, the PACF, the unit-root test, and the selected orders all describe the filter rather than the data. Separately, $\operatorname{Var}(m_t) = \sigma^2/k$ for independent inputs, so the estimated innovation variance is too small by roughly $k$ and the intervals too narrow by roughly $\sqrt{50} \approx 7$. Those intervals will fail their nominal coverage badly. Smooth for display; model the data.
Not established. Failing to reject a unit-root null is not evidence for the null, and the ADF test is known to have low power against near-unit-root alternatives, so $p = 0.31$ at moderate $n$ is consistent with $\phi = 0.9$, $\phi = 0.98$, or a genuine unit root. Run KPSS, whose null is stationarity, and read the pair: both agreeing supports a conclusion, and the “fail to reject / fail to reject” cell means the data cannot decide. Consider also whether a structural break is producing the appearance of a unit root.
About 148 independent observations, not 960. The design effect is $1 + (n-1)\rho = 1 + 11(0.6) = 7.6$, so $N_{\text{eff}} = 960/7.6 \approx 126$, and since the covariate of interest is often between-subject, the effective sample size for that effect is closer to the 80 patients. The practical implication: adding visits buys much less than adding patients, and any analysis ignoring the clustering will report standard errors too small by roughly $\sqrt{7.6} \approx 2.8$.
Neither. They estimate different quantities. The GLMM coefficient is subject-specific, the change in log-odds for a given patient, while the GEE coefficient is population-averaged. For a logit link they are related by $\beta^{\mathrm{m}} \approx \beta^{\mathrm{c}}/\sqrt{1+0.346\sigma_b^2}$, so $0.62/0.95 = 0.65$ implies $\sigma_b \approx 1.9$, a substantial between-subject variance, entirely consistent. The marginal effect is always the smaller in magnitude. For an identity link the two coincide exactly, which is why this never arises in the linear case.
A likelihood-based method: LMM or GLMM. Dropout driven by observed severity is MAR, and likelihood methods remain valid under MAR because the missingness mechanism factors out of the likelihood. GEE requires MCAR and will be biased here. If a marginal interpretation is essential, use weighted GEE with inverse-probability-of-censoring weights, which restores validity under MAR at the cost of having to model the dropout process. And if dropout depends on the unobserved outcome (MNAR), no method is safe without a sensitivity analysis over assumed mechanisms.
Not established. Three things are missing. There is no benchmark, if the seasonal naive forecast achieves 4.5, both models are near-worthless. There is no scale-free metric: report MASE so the number means something beyond this series. And it is one origin, so the difference may be origin-specific noise; rolling-origin CV would show whether the ordering holds. Also check whether the LSTM produces calibrated prediction intervals, which ARIMA gives analytically and a point-forecasting network does not.

42 Summary

Temporal foundations

Stationarity requires constant mean, constant variance, and lag-only covariance. AR stationarity means the roots of $\phi(z)$ lie outside the unit circle; MA invertibility means the same of $\theta(z)$.
The Wold decomposition says every stationary process is a linear filter of white noise, which is why ARMA is the right family, and which gives the forecast variance $\sigma^2\sum_{j<h}\psi_j^2$ directly.
ACF bounds must come from the series being plotted; residual ACFs need Ljung–Box with df reduced by $p+q$.
ADF and KPSS have opposite nulls and should be read together.

ARIMA and forecasting

Difference the minimum number of times; over-differencing drives $\rho(1)\to-0.5$ and inflates the variance.
Use seasonal terms for periodic structure, three parameters where a long MA needs twenty-four, and generalizable beyond one cycle.
Never fit ARIMA to a smoothed series. A $k$-point moving average of white noise has $\operatorname{Corr}(h)=(1-|h|/k)_+$ and variance $\sigma^2/k$: it fabricates the dependence and collapses the intervals.
Always report a benchmark. Naive and seasonal naive are hard to beat; MASE below 1 is the minimum bar.
Correlation is not a forecast metric. Use MASE or RMSSE, evaluated by rolling origin, broken down by horizon.
Prediction intervals widen as $\sqrt{\sum_{j<h}\psi_j^2}$, saturating for stationary models, growing without bound for integrated ones.

Repeated measures

Clustering inflates variance by $\mathrm{DEFF}=1+(n-1)\rho$; ignoring it makes standard errors too small by $\sqrt{\mathrm{DEFF}}$.
Mixed models borrow strength through shrinkage and remain valid under MAR. Use ML to compare fixed effects, REML for variance components.
A random intercept implies compound symmetry; AR(1) is usually more plausible over extended follow-up.
GEE’s sandwich variance is valid even when the working correlation is wrong, the working structure buys efficiency, not validity, but GEE requires MCAR.
GEE and GLMM estimate different quantities, related by $\beta^{\mathrm m}\approx\beta^{\mathrm c}/\sqrt{1+0.346\sigma_b^2}$ for the logit link, and identical for the identity link.

Latent variables

$F_{ML}=\ln|\Sigma|+\operatorname{tr}(S\Sigma^{-1})-\ln|S|-p\ge0$, and $(n-1)\hat F_{ML}\to\chi^2$, the source of the SEM test and its fit indices.
Never build a composite from estimated loadings and model it as observed. A latent growth curve fits measurement and structure jointly and propagates every source of uncertainty.

Sequence networks

Vanilla RNN gradients decay geometrically as $(\|W\|\gamma)^{T-t}$; the LSTM’s additive cell update replaces that product with the forget gate, so gradients survive.
RNNs cost $O(T d^2)$, linear in sequence length, against attention’s $O(T^2 d)$.
Neural methods earn their cost with many related series, long-range non-linearity, or many covariates. On one short seasonal series, ARIMA and ETS usually win.

Where these threads continue

Thread	Continues in
Optimization behind mixed-model and SEM estimation	Function optimization
Convolutional, recurrent, and attention architectures	Deep learning

43 Chapter roadmap

Chapter 1: Foundations. R toolchain, reproducibility, simulation.
Chapter 2: Data quality and visual analytics. Missingness, graphical perception.
Chapter 3: Linear algebra and regression. SVD, conditioning, the hat matrix.
Chapter 4: Dimensionality reduction. PCA, spurious structure when $p\gg n$.
Chapter 5: Supervised classification. Metrics, the leakage taxonomy.
Chapter 6: Black-box methods. Kernels, ensembles, temporal validation.
Chapter 7: Text mining and association rules. Multiplicity over exponential search.
Chapter 8: Unsupervised clustering. Internal and external validation.
Chapter 9: Model assessment. Optimism, calibration, rolling-origin CV.
Chapter 10: Formats, streams, performance. Sketches, drift, Amdahl’s law.
Chapter 11: Feature selection. LASSO, stability selection, knockoffs.
Function optimization. Gradient methods, proximal algorithms, EM.
Deep learning. Convolutional networks, sequence models, attention, transfer learning.

44 Session information

sessionInfo()

#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] lavaan_0.6-18   geepack_1.3.12  nlme_3.1-165    lmerTest_3.1-3 
#>  [5] lme4_1.1-35.5   Matrix_1.6-5    imputeTS_3.3    lubridate_1.9.3
#>  [9] tseries_0.10-55 forecast_8.22.0 plotly_4.12.1   patchwork_1.3.0
#> [13] tidyr_1.3.1     dplyr_1.1.4     ggplot2_4.0.1  
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1    viridisLite_0.4.2   timeDate_4032.109  
#>  [4] farver_2.1.2        loo_2.8.0           S7_0.2.1           
#>  [7] fastmap_1.2.0       digest_0.6.37       timechange_0.3.0   
#> [10] lifecycle_1.0.5     StanHeaders_2.32.9  magrittr_2.0.3     
#> [13] compiler_4.3.3      rlang_1.1.5         sass_0.4.9         
#> [16] tools_4.3.3         yaml_2.3.10         data.table_1.16.4  
#> [19] knitr_1.51          labeling_0.4.3      htmlwidgets_1.6.4  
#> [22] mnormt_2.1.1        pkgbuild_1.4.4      curl_6.2.0         
#> [25] xml2_1.3.6          TTR_0.24.4          RColorBrewer_1.1-3 
#> [28] numDeriv_2016.8-1.1 withr_3.0.2         purrr_1.0.2        
#> [31] nnet_7.3-19         grid_4.3.3          stats4_4.3.3       
#> [34] xts_0.13.2          colorspace_2.1-1    inline_0.3.19      
#> [37] scales_1.4.0        MASS_7.3-60.0.1     cli_3.6.3          
#> [40] rmarkdown_2.31      generics_0.1.3      otel_0.2.0         
#> [43] RcppParallel_5.1.7  rstudioapi_0.18.0   httr_1.4.7         
#> [46] minqa_1.2.7         cachem_1.1.0        rstan_2.32.6       
#> [49] splines_4.3.3       prophet_1.0         parallel_4.3.3     
#> [52] urca_1.3-3          matrixStats_1.3.0   vctrs_0.6.5        
#> [55] V8_6.0.3            boot_1.3-30         jsonlite_1.8.9     
#> [58] stinepack_1.5       crosstalk_1.2.1     jquerylib_0.1.4    
#> [61] quantmod_0.4.26     glue_1.8.0          nloptr_2.1.1       
#> [64] codetools_0.2-20    ggtext_0.1.2        gtable_0.3.6       
#> [67] QuickJSR_1.2.2      quadprog_1.5-8      lmtest_0.9-40      
#> [70] tibble_3.2.1        pillar_1.10.1       htmltools_0.5.8.1  
#> [73] R6_2.6.1            pbivnorm_0.6.0      evaluate_1.0.3     
#> [76] lattice_0.22-6      extraDistr_1.10.0   backports_1.5.0    
#> [79] gridtext_0.1.5      broom_1.0.6         fracdiff_1.5-3     
#> [82] bslib_0.9.0         Rcpp_1.0.14         gridExtra_2.3      
#> [85] mgcv_1.9-1          xfun_0.52           zoo_1.8-12         
#> [88] pkgconfig_2.0.3

Architecture	Parameters (hidden size \(d\), input \(m\))	Cost per step
Vanilla RNN	\(d^2+dm+d\)	\(O(d^2+dm)\)
GRU	\(3(d^2+dm+d)\)	\(O(3(d^2+dm))\)
LSTM	\(4(d^2+dm+d)\)	\(O(4(d^2+dm))\)

Task	Cost	Note
Sample ACF to lag \(L\)	\(O(TL)\), or \(O(T\log T)\) via FFT	`acf()` uses FFT
Differencing	\(O(T)\)	—
ADF / KPSS	\(O(T)\)	Plus lag-order selection
ARIMA likelihood (Kalman)	\(O\big(T(p+q+m)^2\big)\)	Handles gaps exactly
ARIMA fit (\(I\) optimizer steps)	\(O\big(I\,T(p+q+m)^2\big)\)	—
`auto.arima` stepwise	\(\approx O(p_{\max}+q_{\max})\) fits	Default; can miss the optimum
`auto.arima` exhaustive	\(O(p_{\max}q_{\max}P_{\max}Q_{\max})\) fits	450 fits at \(5{,}5{,}2{,}2\)
STL decomposition	\(O(T)\) per inner iteration	—
Rolling-origin CV	\(O(\#\text{origins}\times\text{fit})\)	The honest cost of honest evaluation
Forecast \(h\) steps	\(O(h(p+q))\)	Intervals need \(\psi\)-weights: \(O(h)\)
LMM via REML	\(O\big(N n q_r^2+q_r^3\big)\) per iteration	Sparse \(Z\) makes this near-linear in \(N\)
GLMM, Laplace	\(O(N n q_r^2)\) per iteration	One quadrature point
GLMM, adaptive GH (\(Q\) points)	\(O\big(N n Q^{q_r}\big)\)	Exponential in \(q_r\)
GEE	\(O\big(N n^3\big)\) per iteration	\(V_i^{-1}\) dominates
Sandwich variance	\(O(N n^2 p^2)\)	One pass after convergence
SEM (ML, \(p\) indicators)	\(O\big(I p^3\big)\)	\(\Sigma^{-1}\) and \(\|\Sigma\|\) each iteration
RNN / GRU / LSTM training	\(O(E\,N\,T\,d^2)\)	Linear in \(T\)
Transformer self-attention	\(O(E\,N\,T^2 d)\)	Quadratic in \(T\)

	ACF	PACF
AR(\(p\))	tails off (geometric / damped sine)	cuts off after lag \(p\)
MA(\(q\))	cuts off after lag \(q\)	tails off
ARMA(\(p,q\))	tails off	tails off

Task	Cost
One likelihood evaluation (Kalman)	\(O(T(p+q+m)^2)\)
One model fit (optimizer, \(I\) iterations)	\(O(I\,T(p+q+m)^2)\)
`auto.arima` stepwise search	\(\approx O\big((p_{\max}+q_{\max})\big)\) fits
`auto.arima` exhaustive (`stepwise = FALSE`)	\(O\big(p_{\max}q_{\max}P_{\max}Q_{\max}\big)\) fits
`stl` decomposition	\(O(T)\) per inner iteration

DSPA3: Data Science and Predictive Analytics (UMich HS650)

Chapter 12: Longitudinal Data, Time Series, and Forecasting

SOCR/MCAIM (Ivo Dinov)

August 28, 2026

How this chapter uses graphics

1 Learning objectives

2 PART I: TEMPORAL FOUNDATIONS

3 What makes temporal data different

4 Stationarity

4.1 Conditions on the polynomials

5 Autocorrelation

6 The Wold decomposition

7 Unit-root testing

8 PART II: ARIMA

9 The Beijing PM2.5 series

9.1 Imputing temporal gaps

9.2 The daily cycle

10 Differencing and integration

11 Seasonal ARIMA

12 Estimation and search cost

13 Diagnostics

14 PART III: FORECASTING AND ITS EVALUATION

15 Prediction intervals

16 Benchmarks

17 Scale-free accuracy

18 Rolling-origin cross-validation

19 The smoothing trap

20 Dynamic regression

21 Decomposition and Prophet

22 PART IV: REPEATED MEASURES

23 The correlated-data problem

24 Linear mixed models

24.1 The intraclass correlation

25 Covariance structures

26 Generalized linear mixed models

27 Generalized estimating equations

27.1 The sandwich estimator

28 GEE versus GLMM

29 Missing data

30 PART V: LATENT VARIABLES

31 Structural equation modeling

31.1 The discrepancy function

31.2 Identification

32 Latent growth curves

33 PART VI: SEQUENCE NEURAL NETWORKS

34 Recurrent networks

34.1 Why gradients vanish

35 LSTM and GRU

36 When do neural methods beat ARIMA?

37 PART VII: SYNTHESIS

38 Computational complexity summary

39 Common pitfalls

40 Practice problems

40.1 Problem 1: Quantify the smoothing trap

40.2 Problem 2: ADF and KPSS across the unit circle

40.3 Problem 3: Verify the interval-width formula

40.4 Problem 4: Design effect and effective sample size

40.5 Problem 5: Sandwich robustness

40.6 Problem 6: Attenuation across link functions

40.7 Problem 7: Two-stage composite versus latent growth

40.8 Problem 8: A fair neural-versus-classical comparison

41 Checkpoint

42 Summary

43 Chapter roadmap

44 Session information