--- title: "Choosing how a record is read" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Choosing how a record is read} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4.5, dev.args = list(pointsize = 14)) set.seed(1) library(timesift) ``` A sensor records every hour for years. Before any model sees it, the record is reduced: to monthly means, to growing-degree-days, to whatever the analyst settled on once. `timesift` makes that reduction an argument, fits at each setting of it, and shows how much of the record the prediction actually needed. This vignette runs the whole path on a small simulated record, from two tables to a scored comparison and the prediction that follows from it. ## Two tables `targets` is one row per thing to predict. `series` is the long, time-stamped record belonging to those rows, and the two are linked by an identifier both carry. ```{r data} hours <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 120) ids <- sprintf("p%03d", 1:60) warmth <- rnorm(60) series <- data.frame( plot = rep(ids, each = length(hours)), t = rep(hours, times = 60), temp = as.numeric(vapply(warmth, function(w) { 1.5 * w + 6 * sin(seq_along(hours) / (24 * 40)) + rnorm(length(hours), sd = 12) }, numeric(length(hours)))) ) sign <- rep(c(1, -1), length.out = 6) targets <- data.frame(plot = ids, elevation = 2000 + 300 * rnorm(60)) targets[paste0("sp", 1:6)] <- lapply(sign, function(s) rbinom(60, 1, plogis(3 * s * warmth))) str(targets[1:4], give.attr = FALSE) ``` Each plot has a level of its own, and that level is buried in hour-to-hour noise an order of magnitude larger. ## One call ```{r fit} fit <- timesift( targets, series, y = starts_with("sp"), id = plot, time = t, static = elevation, sift = grains("day", "week", "month"), resampling = cv(v = 5), verbose = FALSE ) fit ``` Every representation named in `sift` was built, and `models` defaulting to `list(elasticnet())`, a penalised logistic regression was fitted on each of them over the same five folds. The report has two parts. The candidates are the comparison. Each is scored on the five outer folds, on the same cells, so their means can be read against each other across grains. `won` is how many responses a candidate scored highest on, and `responses` says whether one fitted model covered them all or one was fitted per response. The highest of those means was picked out on the folds it is scored on, so it is not the number to report. The procedure rows are. Inside each outer training fold the candidates were cross-validated again on five inner folds, one was chosen on its inner score and the stack's weights were fitted on the inner out-of-fold predictions; the choice and the weights then predicted the outer test fold once. `selected` and `ensemble` are those held-out scores, choosing and weighting included. ```{r estimate} fit$estimate[c("arm", "metric", "score", "se", "lower", "upper")] fit$selected[c("fold", "candidate", "inner_score")] ``` The interval is across the six responses of this dataset, all fitted and scored on the same plots and folds, so it does not carry the error those share. A column of `targets` reaches the model only where `static` names it: `elevation` is a predictor here because it was asked for, and a column of notes sitting beside it would not be. ```{r candidates} fit$candidates[c("candidate", "representation", "bins", "channels", "status")] ``` Two channels at every grain: the binned temperature, and elevation held constant across the bins. ```{r plot-run, fig.alt = "Mean AUC against representation, with the level the combination reached drawn across it."} plot(fit) ``` `predict()` rebuilds each member's representation for new rows from the settings its own arm was built with, and combines them through the ensemble refitted on every target; `candidate = "selected"` predicts with the candidate the rule chose on every target instead. ```{r predict} p <- predict(fit, targets, series) round(p[1:3, 1:4], 3) ``` ## Representations A representation carries the settings and nothing else, so the same object describes an arm before any record has been read and rebuilds itself for new targets afterwards. ```{r representations} native() grain("week", stats = c("cold_day", "mean", "warm_day")) multigrain(c("week", "month")) lookback("30 days", bins = 3) ``` `grains()` and `lookbacks()` are sets of them, and `grains("auto")` reads off the record every named grain it gives at least two bins. `multigrain()` flattens several grains side by side into one block of features; `lookback()` is a fixed span ending at each target's own instant, which is what a plot carrying several targets through time needs, given as `target_time`. A learner runs across the whole set, or at one representation it is pinned to: ```{r pinned} elasticnet(data = grain("month")) ``` ## What a learner may be handed A learner declares whether the bins reach it as a block of predictors or as a sequence whose order in time is what it reads. A tabular learner given `native()` is refused before any record is touched, since building the array it would have been handed is the expensive half of the call; a sequence learner is refused a representation that turns out to hold one bin, once the array says how many it has. ```{r compatibility, error = TRUE} timesift(targets, series, y = starts_with("sp"), id = plot, time = t, models = elasticnet(), sift = native()) ``` Inside a set such a pair is skipped and listed as `not applicable` by `summary()`; named through a learner's own `data =` it is an error, because a representation named by hand is a decision. ## The split, and the cells a score is defined on One fold map is read by everything that scores, so every candidate is fitted and scored on identical splits. `cv()` deals units into folds balanced on a stratifying value and `grouped_cv()` keeps every target sharing a group value on one side of each split; `resampling` also takes a fold vector or a `fold_map()` result directly. ```{r folds} fit$folds ``` A per-response score needs both classes among the held-out units, and a per-response model needs both classes among the units it was fitted on. The mask says which cells those are, from the response and the fold map alone. ```{r cells} fit$cells ``` Because it involves no model, every candidate is restricted to the same cells: their means share one denominator and every paired difference runs on matched cells. ## Learners, and how they are trained `elasticnet()` and `stepwise()` read a block of features, `forest()` grows a probability forest over one, and the `torch` encoders `mlp()`, `cnn()` and `rescnn()` read a sequence with a joint multi-label head, so every response is predicted together from a shared embedding. Pooling strength across responses is what makes the rarer ones learnable at these sample sizes. `elasticnet()` is the arm a network is measured against, so it is fitted by the same compiled core the Python package calls rather than by a fitter on either side: the penalty path, the standardisation and the cross-validated choice of penalty are one implementation, and the two languages return the same coefficients for the same design. `s` reads that fit at `"lambda.min"`, at `"lambda.1se"`, or at a penalty of your own, and `thresh` trades how close the descent settles to the optimum against what it costs. Architecture belongs to the constructor and training belongs to `train_control()`, which every neural learner of a run reads. A learner given a control of its own overrides the run's on the settings it names and takes the rest from it. ```{r control} train_control(epochs = 200L, device = "cpu") cnn(channels = c(16L, 32L), epochs = 300L) ``` A learner of your own is a fit and a predict pair, and it goes through the same folds, the same cells and the same scoring as the ones that ship. ```{r custom} flat <- function(x) matrix(as.numeric(x), nrow = dim(x)[1]) nearest_neighbour <- learner( "1nn", fit = function(x, y, ...) list(x = flat(x), y = y), predict = function(model, x) { d <- as.matrix(dist(rbind(flat(x), model$x)))[seq_len(dim(x)[1]), -seq_len(dim(x)[1])] model$y[apply(d, 1, which.min), , drop = FALSE] } ) both <- timesift(targets, series, y = starts_with("sp"), id = plot, time = t, models = c(elasticnet(), nearest_neighbour), sift = grains("week", "month"), resampling = cv(v = 5), verbose = FALSE) summary(both) ``` ## The combination The combiner is handed out-of-fold predictions, the response, the mask and the fold map, and never a model. `"stack"` fits non-negative weights summing to one by minimising the response head's own loss over the scorable cells; `"mean"`, `"median"` and `"weighted"` combine without fitting anything. The weights below are fitted on the outer out-of-fold predictions of every target and are what `predict()` uses. The ensemble's score is not read with them, because they were fitted to the responses it would be scored against; each outer fold's weights are in `fit$fold_weights`. ```{r weights} ensemble_weights(fit) ``` The weights say how much of the combination each candidate carries, and a weight at zero is one the combination reached past. Where the ensemble line of the plot above sits over every curve the candidates are carrying different parts of the signal, and where it sits on the best curve they are not. ## Reading a level honestly The true skill statistic is read at the threshold that maximises it, chosen on the same units the score is then read on. That inflates the level, and by more the fewer presences a cell holds. `tss_inflation()` measures the inflation for the presence counts of the design in hand. ```{r inflation} tss_inflation(fit$y, fit$folds, skill = c(0.6, 0.9), replicates = 100) ``` The inflation is an average over the planted model's predictions: a level is optimistic in expectation, and `implied_skill()` inverts the map to say which population skills a level read is consistent with. Its size depends on how a model's predictions are distributed as well as on the presence counts, so two candidates of equal skill scored on the same cells can be inflated by different amounts, and a paired difference in TSS is not free of it. That is why a run is scored by AUC unless told otherwise. ## The arrays on their own `grain_matrix()` is the representation without the fitting layer around it. It bins the readings by the calendar and summarises every bin. ```{r representation} x <- grain_matrix(series, plot, t, temp, grain = "week", stats = c("cold_day", "mean", "warm_day")) x ``` Bins follow the calendar. A month is 28, 30 or 31 days, and a week starts on a Monday, so a bin is a real month or a real week and not a drifting block of 730 or 168 hours. ```{r calendar} attr(grain_matrix(series, plot, t, temp, grain = "month"), "bin_n")[1, ] ``` ### An extreme day is not an extreme reading `min` and `max` take the coldest and warmest single reading of a bin. `cold_day` and `warm_day` reduce each day to its own mean first and then take the extreme over days. `mean_daily_min` and `mean_daily_max` take the mean of the daily extremes, which is the exposure a typical day of the bin brought. One hour at -50 sets `min` to -50 outright; it reaches the day-level statistics only through its twenty-fourth of that day's mean. ```{r statistics} week <- grain_matrix(series, plot, t, temp, grain = "week", stats = c("min", "mean_daily_min", "cold_day", "mean", "warm_day", "mean_daily_max", "max")) round(week[1, 1, ], 2) ``` An encoder that ends in global pooling discards when a thermal event happened, so the position of a bin in the year is given to it as input. ```{r channels} week_mean <- grain_matrix(series, plot, t, temp, grain = "week") dimnames(bind_channels(week_mean, calendar_channels(week_mean)))[[3]] ``` ## One grain at a time Where the arrays are already built, `grain_ladder()` fits every learner at every grain of a set on one split and one mask. It is the ladder a run reports as a curve, reachable on its own. ```{r ladder} set <- grain_matrix(series, plot, t, temp, grain = c("day", "week", "month")) lad <- grain_ladder(set, fit$y, elasticnet(), folds = fit$folds, verbose = FALSE) summary(lad) ``` A claim about one step of that curve rests on the paired contrast. The difference is taken inside each cell both arms scored, averaged within a response, and summarised across the responses, which are the independent replicates. Six responses is few, so the interval is wide and the rank test behind `p_value` has few values to work with. ```{r contrast} paired_contrast(lad, "month|elasticnet", "day|elasticnet") ``` Where the whole curve is the question rather than one step of it, `grain_contrasts()` fits a mixed model on the per-cell scores and compares every grain against the best one, correcting for the comparisons made and no others. ```{r grain-contrasts, eval = all(vapply(c("lme4", "lmerTest", "emmeans"), requireNamespace, logical(1), quietly = TRUE))} grain_contrasts(lad) ``` ## What was read With the per-fold fits kept, `occlusion()` holds each bin of the record back in turn, rescores the held-out units, and records the fall in score as that bin's weight. Nothing is refitted. ```{r occlusion} kept <- timesift(targets, series, y = starts_with("sp"), id = plot, time = t, sift = grains("month"), resampling = cv(v = 5), ensemble = FALSE, keep_fits = TRUE, verbose = FALSE) weight <- occlusion(kept, "elasticnet / month", permutations = 5) head(aggregate(weight ~ part, weight, mean), 4) ``` Holding a channel back instead asks what each statistic of a grain carries, which is the question behind keeping a bin's extremes at all.