Rolling regression slope down every column of a matrix

Fits y ~ 1 + t, with t from 1 to width, to each window of width rows in every column. It is the kernel of the ensemble method of short_term_trend().

Usage

rolling_slope_matrix(
  Y,
  width,
  family = c("identity", "quasipoisson", "binomial"),
  prior_weights = NULL,
  tol = 1e-10,
  maxit = 25L,
  rate_denominator = NULL,
  per = 100
)

Arguments

Y: A numeric matrix, rows = weeks in time order, columns = draws.
width: The window width in rows, at least 2.
family: "identity" is ordinary least squares. "quasipoisson" is a log link and "binomial" a logit link, and for those beta1 is on the link scale.
prior_weights: The binomial denominators, a matrix the shape of Y or with one column. Only family = "binomial" takes it, and it needs it.
tol: The iteration stops when no coefficient moves more than tol.
maxit: The most iterations. A window that still moves more than tol did not converge.
rate_denominator: The denominators of a rate, a matrix the shape of Y or with one column. Only family = "identity" takes it. It turns on the variance floor, see the return value. For a count, give a denominator of 1 and per = 1. Each value MUST be 0 or more and finite, or NA. A 0 counts as NA. An NA week gives NA in the windows that hold it, and in no other window.
per: The scale of the rate in Y, such as 100 for a percentage. Only used with rate_denominator. It MUST be one number above 0 and finite.

Value

A list of four matrices, each the shape of Y: beta0, beta1, se and converged. se is the standard error of beta1, a Wald standard error under IRLS. Under "identity" with rate_denominator, the residual variance of se has a floor at the Poisson variance, mean(Y) * per / mean(rate_denominator) over the window. Without rate_denominator, se is the OLS standard error, with no floor. The logical converged is NA on the first width - 1 rows.

Details

Row i holds the fit of the window that ends on row i, so the first width - 1 rows are NA. The identity family has a closed form. The other two families run iteratively reweighted least squares (IRLS) on every window at once. A window with complete or quasi-complete separation cannot converge: its beta0, beta1 and se are NA, and one warning gives the count. In every family, an NA in Y gives NA in the windows that hold it, and in no other window. Under "identity", Inf and -Inf count as NA.

One deliberate difference from glm

The function agrees with stats::glm() except on an all-zero window. There it returns NA, and stats::glm() returns a slope of 0. stats::glm() stops on the relative change in deviance, which is 0 from the first step. This function stops on the coefficient step, which never settles, because the intercept goes to -Inf. An all-zero window identifies no slope, and NA says so.

Examples

# 10 weeks x 4 draws, all rising at a true slope of 2 per week
set.seed(1)
Y <- matrix(rep(1:10, 4) * 2 + rnorm(40), nrow = 10)

rs <- rolling_slope_matrix(Y, width = 4)

# the first three rows have no complete window, so they are NA
head(rs$beta1, 3)
#>      [,1] [,2] [,3] [,4]
#> [1,]   NA   NA   NA   NA
#> [2,]   NA   NA   NA   NA
#> [3,]   NA   NA   NA   NA

# the later rows find the slope, one estimate per draw
round(rs$beta1[8:10, ], 2)
#>      [,1] [,2] [,3] [,4]
#> [1,] 2.25 1.95 1.36 2.40
#> [2,] 2.44 2.36 1.74 2.49
#> [3,] 1.75 2.17 2.27 2.46

# `se` is the OLS standard error of the slope
round(rs$se[10, ], 2)
#> [1] 0.18 0.20 0.39 0.20

# a log link: counts that rise by 0.2 per week on the log scale
N <- matrix(rpois(40, rep(exp(1 + 0.2 * 1:10), 4)), nrow = 10)
round(rolling_slope_matrix(N, width = 4, family = "quasipoisson")$beta1[10, ], 2)
#> [1] 0.24 0.45 0.19 0.35