---
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}
resource_files:
- timesift_files/py-figure-1.svg
---
```{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)
if (nzchar(Sys.getenv("TIMESIFT_DOCS_CACHE"))) source(Sys.getenv("TIMESIFT_DOCS_CACHE"))
```
Sensor measurements can be summarised as daily, weekly or monthly values. `timesift` fits models using each of these summaries and compares how well they predict observations excluded from fitting. The time scale of a summary is called its *grain*.
This vignette uses a small simulated dataset to show how to prepare the two input tables, compare grains and predict new observations.
## Installation
::: {.lang-pair}
```{r, eval = FALSE}
# Install from CRAN
install.packages("timesift")
# Or install the development version from GitHub
# install.packages("pak")
pak::pak("gcol33/timesift")
```
```bash
# Install from PyPI
pip install timesift
# with the torch encoders, the contrasts and the plots
pip install "timesift[torch,contrasts,plot]"
# Or install the development version from GitHub
pip install git+https://github.com/gcol33/timesift
```
:::
## Two tables
`targets` contains one row per sampling unit, such as a study plot, and the values to predict, called *responses*. Here the responses record whether each species occurs on a plot. `series` contains the time-stamped sensor readings for those plots. A shared plot identifier links the two tables.
::: {.lang-py-only}
In Python either table can be a mapping of column names to arrays, including a pandas or polars data frame.
:::
::: {.lang-pair}
```{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)
```
```python
import numpy as np
import pandas as pd
import timesift as ts
rng = np.random.default_rng(1)
hours = pd.date_range("2021-09-01", periods=24 * 120, freq="h")
ids = [f"p{i:03d}" for i in range(1, 61)]
warmth = rng.normal(size=60)
season = 6 * np.sin(np.arange(1, len(hours) + 1) / (24 * 40))
series = pd.DataFrame({
"plot": np.repeat(ids, len(hours)),
"t": np.tile(hours, 60),
"temp": np.concatenate([1.5 * w + season + rng.normal(scale=12, size=len(hours))
for w in warmth]),
})
targets = pd.DataFrame({"plot": ids, "elevation": 2000 + 300 * rng.normal(size=60)})
for k, sign in enumerate([1, -1, 1, -1, 1, -1], start=1):
targets[f"sp{k}"] = rng.binomial(1, 1 / (1 + np.exp(-3 * sign * warmth)))
targets.iloc[:, :4].head()
#> plot elevation sp1 sp2
#> 0 p001 2159.960722 1 0
#> 1 p002 1516.179723 1 0
#> 2 p003 1818.289129 1 0
#> 3 p004 2233.373010 0 1
#> 4 p005 1779.544877 1 0
```
:::
Each plot has its own temperature offset. The hourly noise is an order of magnitude larger than the differences between these offsets.
## One call
::: {.lang-pair}
```{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
```
```python
fit = ts.timesift(
targets, series,
y="sp*",
id="plot",
time="t",
static="elevation",
sift=ts.grains("day", "week", "month"),
resampling=ts.cv(v=5),
verbose=False,
)
fit
#> timesift 60 targets, 6 responses, 5-fold random CV, roc_auc
#>
#> candidates, scored on the outer folds
#> candidate mean won responses
#> elasticnet / day 0.870 0 separate
#> elasticnet / week 0.899 3 separate
#> elasticnet / month 0.911 3 separate
#>
#> procedure, chosen and weighted inside each outer training fold
#> selected 0.911 se 0.011
#> ensemble 0.907 se 0.012
#> selected elasticnet / month in 5 of 5 folds
#>
#> choice on every target elasticnet / month
#> weights on every target elasticnet / month 0.96 elasticnet / week 0.04
```
:::
The call builds a table of predictors at each time scale listed in `sift`. By default, it fits a penalised logistic regression, which estimates presence probabilities while limiting the size of its coefficients. To evaluate predictions, the plots are divided into five groups, called *folds*. Each fold is predicted by a model fitted on the other four. This is five-fold cross-validation.
Each model and time-scale combination is a *candidate*. The candidate rows compare scores on the same folds for the same responses. A response in one fold is called a *cell* in the scoring output. `won` counts the responses on which a candidate scored highest. `responses` indicates whether the model fits responses jointly or separately. The default score is the area under the receiver operating characteristic curve (AUC), which measures how well predictions rank presences above absences.
The `selected` row evaluates the process of choosing a candidate. The `ensemble` row evaluates a weighted combination of candidates. Both use an additional round of cross-validation within each training set: these *inner folds* choose the candidate and estimate the weights, and the *outer test fold* evaluates their predictions. Use these rows to report predictive performance after selection or combination.
::: {.lang-pair}
```{r estimate}
fit$estimate[c("arm", "metric", "score", "se", "lower", "upper")]
fit$selected[c("fold", "candidate", "inner_score")]
```
```python
pd.DataFrame(fit.estimate)[["arm", "metric", "score", "se", "lower", "upper"]]
#> arm metric score se lower upper
#> 0 selected accuracy 0.902778 0.011719 0.872652 0.932904
#> 1 selected average_precision 0.908116 0.023171 0.848552 0.967679
#> 2 selected bias 1.032778 0.038854 0.932901 1.132654
#> 3 selected boyce 0.701781 0.015503 0.661930 0.741633
#> 4 selected csi 0.819987 0.026700 0.751353 0.888621
#> 5 selected ets 0.678483 0.034219 0.590519 0.766447
#> 6 selected far 0.102672 0.026525 0.034488 0.170856
#> 7 selected kappa 0.660423 0.024039 0.598629 0.722217
#> 8 selected kappa_youden 0.796165 0.023657 0.735353 0.856978
#> 9 selected neg_mae -0.371196 0.013208 -0.405149 -0.337243
#> 10 selected neg_max_error -0.635774 0.007559 -0.655206 -0.616342
#> 11 selected neg_mse -0.162203 0.006823 -0.179742 -0.144664
#> 12 selected neg_poisson_deviance -0.513075 0.018045 -0.559461 -0.466689
#> 13 selected neg_rmse -0.401536 0.008560 -0.423540 -0.379532
#> 14 selected or 22.750000 1.030776 19.469609 26.030391
#> 15 selected ordinal_accuracy 0.808333 0.009379 0.784224 0.832442
#> 16 selected ordinal_f1 0.814673 0.010448 0.787816 0.841529
#> 17 selected ordinal_precision 0.807500 0.009108 0.784086 0.830914
#> 18 selected ordinal_recall 0.822665 0.012320 0.790995 0.854336
#> 19 selected orss 0.986042 0.005010 0.973164 0.998920
#> 20 selected pearson 0.675293 0.014535 0.637930 0.712656
#> 21 selected pod 0.910992 0.015011 0.872405 0.949579
#> 22 selected pofd 0.093545 0.019192 0.044211 0.142879
#> 23 selected r_squared 0.285434 0.029940 0.208471 0.362397
#> 24 selected roc_auc 0.910976 0.011209 0.882161 0.939790
#> 25 selected sr 0.897328 0.026525 0.829144 0.965512
#> 26 selected tss 0.817447 0.019358 0.767686 0.867208
#> 27 ensemble accuracy 0.897222 0.012485 0.865130 0.929315
#> 28 ensemble average_precision 0.905069 0.021508 0.849781 0.960356
#> 29 ensemble bias 1.013267 0.045904 0.895268 1.131267
#> 30 ensemble boyce 0.692614 0.027483 0.621968 0.763260
#> 31 ensemble csi 0.810000 0.026807 0.741092 0.878908
#> 32 ensemble ets 0.665135 0.036505 0.571295 0.758975
#> 33 ensemble far 0.099101 0.028427 0.026027 0.172174
#> 34 ensemble kappa 0.663810 0.029921 0.586894 0.740725
#> 35 ensemble kappa_youden 0.785917 0.025575 0.720174 0.851659
#> 36 ensemble neg_mae -0.370921 0.013551 -0.405754 -0.336088
#> 37 ensemble neg_max_error -0.647332 0.004716 -0.659455 -0.635209
#> 38 ensemble neg_mse -0.163398 0.007436 -0.182512 -0.144283
#> 39 ensemble neg_poisson_deviance -0.516055 0.018174 -0.562774 -0.469336
#> 40 ensemble neg_rmse -0.403008 0.009305 -0.426928 -0.379088
#> 41 ensemble or 23.125000 0.875000 20.340359 25.909641
#> 42 ensemble ordinal_accuracy 0.811111 0.015316 0.771741 0.850481
#> 43 ensemble ordinal_f1 0.817597 0.016364 0.775531 0.859662
#> 44 ensemble ordinal_precision 0.810853 0.015969 0.769804 0.851902
#> 45 ensemble ordinal_recall 0.825112 0.017178 0.780955 0.869270
#> 46 ensemble orss 0.983375 0.006059 0.967800 0.998951
#> 47 ensemble pearson 0.672901 0.015095 0.634098 0.711704
#> 48 ensemble pod 0.896243 0.019159 0.846994 0.945492
#> 49 ensemble pofd 0.086878 0.020548 0.034058 0.139699
#> 50 ensemble r_squared 0.281085 0.030892 0.201674 0.360497
#> 51 ensemble roc_auc 0.906884 0.011914 0.876257 0.937510
#> 52 ensemble sr 0.900899 0.028427 0.827826 0.973973
#> 53 ensemble tss 0.809365 0.022716 0.750971 0.867759
pd.DataFrame(fit.selected)[["fold", "candidate", "inner_score"]]
#> fold candidate inner_score
#> 0 1 elasticnet / month 0.883939
#> 1 2 elasticnet / month 0.899061
#> 2 3 elasticnet / month 0.897351
#> 3 4 elasticnet / month 0.908608
#> 4 5 elasticnet / month 0.896122
```
:::
The interval describes variation across the six responses. All responses use the same plots and folds, so the interval excludes uncertainty shared through that sampling design.
Use `static` to name columns of `targets` that should enter the model as predictors. Here it selects `elevation`; additional columns, such as notes, are ignored.
::: {.lang-pair}
```{r candidates}
fit$candidates[c("candidate", "representation", "bins", "channels", "status")]
```
```python
pd.DataFrame(fit.candidates)[["candidate", "representation", "bins", "channels"]]
#> candidate representation bins channels
#> 0 elasticnet / day day 120 2
#> 1 elasticnet / week week 18 2
#> 2 elasticnet / month month 4 2
```
:::
The model now receives temperature and elevation. Each is stored as a separate *channel*, a set of values over the time bins. A bin is one interval, such as a calendar week. Elevation has the same value in every bin for a given plot.
::: {.lang-pair}
```{r plot-run, fig.alt = "Mean AUC against representation, with the level the combination reached drawn across it."}
plot(fit)
```
```{.python fig-alt="Mean AUC against representation, with the level the combination reached drawn across it."}
ts.plot(fit);
```
:::
`predict()` prepares the new data using each candidate's saved settings, applies models refitted on all training plots, and combines their predictions using the fitted ensemble weights. Set `candidate = "selected"` to use the candidate chosen from all training plots.
::: {.lang-pair}
```{r predict}
p <- predict(fit, targets, series)
round(p[1:3, 1:4], 3)
```
```python
p = fit.predict(targets, series)
p[:3, :4].round(3)
#> array([[0.535, 0.429, 0.629, 0.337],
#> [0.644, 0.213, 0.772, 0.221],
#> [0.588, 0.312, 0.63 , 0.381]])
```
:::
`type = "binary"` cuts each response into presence and absence. The cut is learned from the same
candidate's out-of-fold predictions of the fit's own targets, so it comes from predictions of
plots the model had not been fitted on. `decision_threshold()` on the fit returns the cuts
themselves, and `rule` picks how they are learned: `"youden"` maximises sensitivity plus
specificity, `"kappa"` maximises Cohen's kappa, and `"prevalence"` predicts as many presences as
were observed.
::: {.lang-pair}
```{r binary}
decision_threshold(fit, rule = "prevalence")
predict(fit, targets, series, type = "binary", rule = "prevalence")[1:3, 1:4]
```
```python
pd.Series(ts.decision_threshold(fit, rule="prevalence"))
#> sp1 0.520268
#> sp2 0.502498
#> sp3 0.523192
#> sp4 0.485687
#> sp5 0.513939
#> sp6 0.507207
#> dtype: float64
fit.predict(targets, series, type="binary", rule="prevalence")[:3, :4]
#> array([[1., 0., 1., 0.],
#> [1., 0., 1., 0.],
#> [1., 0., 1., 0.]])
```
:::
A binary map is this prediction with one target per map cell, each carrying a series of its own
at the grain the fit reads.
## Representations
A representation stores the settings used to summarise a record. The same settings are used to build model inputs during fitting and prediction.
::: {.lang-pair}
```{r representations}
native()
grain("week", stats = c("cold_day", "mean", "warm_day"))
multigrain(c("week", "month"))
lookback("30 days", bins = 3)
```
```python
ts.native()
#> native
#> kind : grain (a sequence)
#> grain : native
#> stats : mean
ts.grain("week", stats=("cold_day", "mean", "warm_day"))
#> week
#> kind : grain (a sequence)
#> grain : week
#> stats : cold_day, mean, warm_day
ts.multigrain(("week", "month"))
#> multigrain(week+month)
#> kind : multigrain (a block of features)
#> grains : week, month
#> stats : mean
ts.lookback("30 days", bins=3)
#> 30 days x3
#> kind : lookback (a sequence)
#> span : 30 days in 3 bins ending 0 days before the target
#> stats : mean
```
:::
`grains()` defines several time scales to compare. `grains("auto")` includes each named scale that gives at least two time bins. `multigrain()` puts summaries at several scales into one table. A *lookback window* covers a fixed period before a prediction date, such as the four weeks before a survey. `lookback()` defines one window and `lookbacks()` defines several to compare; `target_time` identifies the date column.
By default, each learner runs on all compatible representations in the set. Use `data` to assign a specific representation:
::: {.lang-pair}
```{r pinned}
elasticnet(data = grain("month"))
```
```python
ts.elasticnet(data=ts.grain("month"))
#> elasticnet
#> reads : tabular ; one model per response: yes, separate
#> data : month
#> settings: alpha = 0.5, n_inner = 10, squares = True, s = lambda.1se, n_lambda = 100, tol = 1e-08, threads = 1, seed = 1
```
:::
## Compatible model inputs
Learners accept either a table of predictors or a sequence of time bins. Tabular learners reject `native()` before the array is built. Sequence learners require more than one bin; this is checked after building the array.
::: {.lang-pair}
```{r compatibility, error = TRUE}
timesift(targets, series, y = starts_with("sp"), id = plot, time = t,
learners = elasticnet(), sift = native())
```
```{.python error=true}
ts.timesift(targets, series, y="sp*", id="plot", time="t",
learners=[ts.elasticnet()], sift=ts.native(), verbose=False)
#> ValueError: no learner can read any of the representations:
#> elasticnet() reads a tabular representation; native gives it one column per reading. Use grain(), multigrain() or lookback().
```
:::
In a set of representations, incompatible learner-representation pairs are skipped and reported as `not applicable`. An incompatible representation assigned explicitly through a learner's `data` argument raises an error.
## The split, and the cells a score is defined on
A *fold map* records which test fold each plot belongs to. Every candidate uses the same map. `cv()` distributes plots across folds with similar values of a chosen variable. `grouped_cv()` keeps observations from the same group together, such as repeated surveys of one plot. You can also supply fold assignments directly through `resampling`.
::: {.lang-pair}
```{r folds}
fit$folds
```
```python
fit.folds
#> 60 units in 5 folds
#> 1: 12 2: 12 3: 12 4: 12 5: 12
```
:::
For a response to be scored in a fold, both presences and absences must occur in the held-out data. A model fitted separately for that response also needs both classes in its training data. The scorable-cell mask identifies the response-fold pairs that meet these conditions.
::: {.lang-pair}
```{r cells}
fit$cells
```
```python
fit.cells
#> 30 cells over 6 variables
#> scorable: 30 (100.0%); variables with at least one scorable fold: 6 of 6
```
:::
The mask is calculated from the responses and folds before fitting. All candidates are evaluated on the same cells, so mean scores and paired differences use matching observations.
## Learners, and how they are trained
A *learner* specifies how to fit and predict with a model. `elasticnet()`, `linear()` and `forest()` use a table of predictors. The neural networks `mlp()`, `cnn()` and `rescnn()` model all responses together, sharing information between them. `mlp()` uses a table; `cnn()` and `rescnn()` use the order of the time bins to detect patterns across neighbouring periods. These three networks require the torch library.
The R and Python versions of `elasticnet()` use the same compiled implementation to standardise predictors, fit the models and choose the penalty by cross-validation. The `s` argument selects the penalty: `"lambda.min"` gives the lowest cross-validation error; `"lambda.1se"` chooses a stronger penalty within one standard error of that minimum. A numeric value specifies the penalty directly. `tol` controls when the fitting algorithm stops.
Set the network architecture in the learner constructor and shared training settings in `train_control()`. Settings supplied to an individual learner override the corresponding run-level settings.
::: {.lang-pair}
```{r control}
train_control(epochs = 200L, device = "cpu")
cnn(channels = c(16L, 32L), epochs = 300L)
```
```python
ts.train_control(epochs=200, device="cpu")
#>
#> epochs 200
#> batch_size 64 (default)
#> learning_rate 0.001 (default)
#> weight_decay 0.0001 (default)
#> optimizer adamw (default)
#> penalty 0.0 (default)
#> alpha 0.5 (default)
#> schedule cosine (default)
#> plateau_factor 0.1 (default)
#> plateau_patience 10 (default)
#> early_stopping inf (default)
#> val_frac 0.0 (default)
#> device cpu
#> seed 1 (default)
#> swa False (default)
#> swa_start 0.7 (default)
ts.cnn(channels=(16, 32), epochs=300)
#> cnn
#> reads : sequence ; one model per response: no, joint
#> data : every representation of the run
#> settings: channels = 16/32, kernel = 7, dropout = 0.3
#> training: epochs = 300
#> needs : torch
```
:::
To add a model, supply functions for fitting and prediction. It will use the same folds and scoring rules as the built-in models. [`as.matrix()`]{.lang-r-only}[`flatten()`]{.lang-py-only} converts the input to a table with one predictor per bin and channel. A fitting function with a `weights` argument receives weights for the observations; the default presence-absence weights give more weight to rare presences.
::: {.lang-pair}
```{r custom}
nearest_neighbour <- learner(
"1nn",
fit = function(x, y, ...) list(x = as.matrix(x), y = y),
predict = function(model, x) {
m <- as.matrix(x)
d <- as.matrix(dist(rbind(m, model$x)))[seq_len(nrow(m)), -seq_len(nrow(m))]
model$y[apply(d, 1, which.min), , drop = FALSE]
}
)
weighted_glm <- learner(
"weighted_glm",
fit = function(x, y, weights, ...) {
f <- as.data.frame(as.matrix(x))
lapply(seq_len(ncol(y)), function(j) {
glm(y[, j] ~ ., data = f, weights = weights[, j], family = quasibinomial(),
control = glm.control(maxit = 100))
})
},
predict = function(model, x) {
f <- as.data.frame(as.matrix(x))
vapply(model, function(m) predict(m, f, type = "response"), numeric(nrow(f)))
},
data = grain("month")
)
both <- timesift(targets, series, y = starts_with("sp"), id = plot, time = t,
learners = c(elasticnet(), nearest_neighbour, weighted_glm),
sift = grains("week", "month"), resampling = cv(v = 5), verbose = FALSE)
summary(both)
```
```python
def nn_fit(x, y, **_):
return {"x": ts.flatten(x), "y": y}
def nn_predict(model, x):
m = ts.flatten(x)
d = ((m[:, None, :] - model["x"][None, :, :]) ** 2).sum(axis=2)
return model["y"][d.argmin(axis=1)].astype(float)
nearest_neighbour = ts.Learner("1nn", fit=nn_fit, predict=nn_predict, multi="joint")
def irls(design, y, w):
beta = np.zeros(design.shape[1])
for _ in range(100):
p = 1 / (1 + np.exp(-design @ beta))
step = np.linalg.solve(design.T @ ((w * p * (1 - p))[:, None] * design),
design.T @ (w * (y - p)))
beta += step
if np.abs(step).max() < 1e-8:
break
return beta
def logistic_fit(x, y, weights, **_):
m = ts.flatten(x)
centre, scale = m.mean(axis=0), m.std(axis=0)
design = np.column_stack([np.ones(len(m)), (m - centre) / scale])
betas = np.column_stack([irls(design, y[:, j], weights[:, j]) for j in range(y.shape[1])])
return {"centre": centre, "scale": scale, "betas": betas}
def logistic_predict(model, x):
m = (ts.flatten(x) - model["centre"]) / model["scale"]
design = np.column_stack([np.ones(len(m)), m])
return 1 / (1 + np.exp(-design @ model["betas"]))
weighted_logistic = ts.Learner("weighted_logistic", fit=logistic_fit,
predict=logistic_predict, data=ts.grain("month"))
both = ts.timesift(targets, series, y="sp*", id="plot", time="t",
learners=[ts.elasticnet(), nearest_neighbour, weighted_logistic],
sift=ts.grains("week", "month"), resampling=ts.cv(v=5), verbose=False)
print(ts.summary(both))
#> timesift 60 targets, 6 responses, 5-fold random CV, roc_auc
#>
#> candidates, scored on the outer folds
#> candidate mean won responses
#> 1nn / week 0.744 0 joint
#> 1nn / month 0.757 0 joint
#> elasticnet / week 0.899 3 separate
#> weighted_logistic / month 0.900 0 separate
#> elasticnet / month 0.912 3 separate
#>
#> procedure, chosen and weighted inside each outer training fold
#> selected 0.912 se 0.011
#> ensemble 0.901 se 0.009
#> selected elasticnet / month in 5 of 5 folds
#>
#> choice on every target elasticnet / month
#> weights on every target weighted_logistic / month 0.71 elasticnet / week 0.24 1nn / week 0.05
```
:::
## The combination
The ensemble combines *out-of-fold predictions*: predictions for observations excluded from the model's training data. The `"stack"` method chooses non-negative weights that sum to one to minimise prediction error, using the loss function defined for the response type. Other methods take the mean, median, score-weighted mean or the share of models predicting presence (`"committee"`).
The weights below use out-of-fold predictions for all targets and are used by `predict()`. The reported ensemble score uses weights estimated separately within each outer training fold, stored in [`fit$fold_weights`]{.lang-r-only}[`fit.fold_weights`]{.lang-py-only}.
::: {.lang-pair}
```{r weights}
ensemble_weights(fit)
```
```python
pd.Series(ts.ensemble_weights(fit))
#> elasticnet / day 5.967019e-269
#> elasticnet / week 3.517667e-02
#> elasticnet / month 9.648233e-01
#> dtype: float64
```
:::
The weights determine each candidate's contribution to the ensemble prediction. A candidate with zero weight contributes nothing. Compare the ensemble score with the individual candidate scores to see whether combining their predictions improved performance in this run.
## TSS and threshold selection
The true skill statistic (TSS) is overestimated when the threshold is chosen to maximise TSS on the same held-out observations used to report the score. The fewer presences there are in the held-out data, the more TSS tends to overestimate population skill. `tss_inflation()` estimates the TSS a given sampling design would report for a known population skill, here the presence counts in your response-fold cells.
::: {.lang-pair}
```{r inflation}
tss_inflation(fit$y, fit$folds, skill = c(0.6, 0.9), replicates = 100)
```
```python
pd.DataFrame(ts.tss_inflation(fit.y, fit.folds, skill=(0.6, 0.9), replicates=100))
#> skill reported inflation lower upper replicates
#> 0 0.6 0.769709 0.169709 0.718532 0.835235 100
#> 1 0.9 0.966596 0.066596 0.942541 0.987280 100
```
:::
`tss_inflation()` averages the bias over simulated predictions, and `implied_skill()` estimates the population skill consistent with an observed score under that simulation. The bias depends on the distribution of predictions as well as the presence counts. It can therefore differ between equally skilled candidates evaluated on the same cells and affect their paired TSS difference. The default scoring metric is area under the receiver operating characteristic curve (AUC).
## The arrays on their own
`grain_matrix()` groups readings into calendar intervals and summarises each interval. It returns the model inputs so you can inspect them or fit a model separately.
::: {.lang-pair}
```{r representation}
x <- grain_matrix(series, plot, t, temp, grain = "week",
stats = c("cold_day", "mean", "warm_day"))
x
```
```python
x = ts.grain_matrix(series, "plot", "t", "temp", grain="week",
stats=("cold_day", "mean", "warm_day"))
x
#> 60 units x 18 bins x 3 channels
#> grain: week stats: cold_day, mean, warm_day
#> from : 2021-08-30T00:00:00Z to 2021-12-27T00:00:00Z
```
:::
Bins follow the calendar: months have 28 to 31 days, and weeks begin on Monday.
::: {.lang-pair}
```{r calendar}
attr(grain_matrix(series, plot, t, temp, grain = "month"), "bin_n")[1, ]
```
```python
ts.grain_matrix(series, "plot", "t", "temp", grain="month").bin_n[0]
#> array([720, 744, 720, 696], dtype=int32)
```
:::
### Daily and individual extremes
`min` and `max` return the lowest and highest individual reading in each bin. `cold_day` and `warm_day` return the lowest and highest daily mean. `mean_daily_min` and `mean_daily_max` average the daily minima and maxima within the bin.
::: {.lang-pair}
```{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)
```
```python
week = ts.grain_matrix(series, "plot", "t", "temp", grain="week",
stats=("min", "mean_daily_min", "cold_day", "mean",
"warm_day", "mean_daily_max", "max"))
pd.Series(week.values[0, 0], index=week.stats).round(2)
#> min -28.68
#> mean_daily_min -22.30
#> cold_day -1.73
#> mean -0.13
#> warm_day 0.94
#> mean_daily_max 21.54
#> max 31.50
#> dtype: float64
```
:::
`calendar_channels()` adds values describing when each bin occurs within the year. This lets a neural network use the season of an event even when it averages information across time later in the network.
::: {.lang-pair}
```{r channels}
week_mean <- grain_matrix(series, plot, t, temp, grain = "week")
dimnames(bind_channels(week_mean, calendar_channels(week_mean)))[[3]]
```
```python
week_mean = ts.grain_matrix(series, "plot", "t", "temp", grain="week")
ts.bind_channels(week_mean, ts.calendar_channels(week_mean)).stats
#> ('mean', 'year_sin', 'year_cos')
```
:::
## One grain at a time
`grain_ladder()` fits and evaluates each learner at each grain using arrays you have already built. All combinations use the same folds and eligible response-fold pairs.
::: {.lang-pair}
```{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)
```
```python
grain_set = ts.grain_matrix(series, "plot", "t", "temp", grain=("day", "week", "month"))
lad = ts.grain_ladder(grain_set, fit.y, [ts.elasticnet()], folds=fit.folds, verbose=False)
pd.DataFrame(lad.summary())
#> learner grain score n_variable best
#> 0 elasticnet day 0.869715 6 False
#> 1 elasticnet week 0.899052 6 False
#> 2 elasticnet month 0.912017 6 True
```
:::
`paired_contrast()` calculates the score difference between two candidates in each shared cell, averages across folds within each response, and summarises those response-level differences. The interval and rank test use responses as replicates. With six responses, the interval is wide and the test has limited p-value resolution.
::: {.lang-pair}
```{r contrast}
paired_contrast(lad, "month|elasticnet", "day|elasticnet")
```
```python
pd.Series(ts.paired_contrast(lad, "month|elasticnet", "day|elasticnet"))
#> a month|elasticnet
#> b day|elasticnet
#> diff 0.042303
#> center 0.042303
#> lower 0.003589
#> upper 0.081016
#> n_variable 6
#> n_cell 30
#> n_favour 6
#> p_value 0.03125
#> p_method exact
#> interval variables
#> dtype: object
```
:::
`grain_contrasts()` compares each grain with the best-scoring grain using a mixed model, which accounts for differences among responses and folds. It adjusts the comparisons for testing several grains.[ It needs scipy, available through the `contrasts` extra.]{.lang-py-only}
::: {.lang-pair}
```{r grain-contrasts, eval = all(vapply(c("lme4", "lmerTest", "emmeans"), requireNamespace, logical(1), quietly = TRUE))}
grain_contrasts(lad)
```
```python
pd.DataFrame(ts.grain_contrasts(lad))
#> learner grain reference diff lower upper p_value
#> 0 elasticnet day month -0.042303 -0.088449 0.003844 0.076834
#> 1 elasticnet week month -0.012965 -0.059111 0.033181 0.752021
```
:::
## Identifying influential bins
`occlusion()` measures how much a score changes when the information in one time bin is disrupted. It uses the saved models for each fold and evaluates the modified inputs on the held-out observations.
::: {.lang-pair}
```{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)
```
```python
kept = ts.timesift(targets, series, y="sp*", id="plot", time="t",
sift=ts.grains("month"), resampling=ts.cv(v=5), ensemble=False,
keep_fits=True, verbose=False)
profile = ts.occlusion(kept, "elasticnet / month", permutations=5)
weight = pd.DataFrame(profile["weight"], index=profile["part"], columns=profile["variable"])
weight.mean(axis=1).head(4)
#> 2021-09-01T00:00:00Z 0.092622
#> 2021-10-01T00:00:00Z 0.078312
#> 2021-11-01T00:00:00Z 0.022786
#> 2021-12-01T00:00:00Z 0.032392
#> dtype: float64
```
:::
Occluding a whole channel measures how predictions depend on a particular statistic, such as the mean or the coldest day.
## Where next
Each part of this path has an article of its own on the package site:
- [Representing a record](https://gillescolling.com/timesift/articles/representations.html): every
grain and statistic, calendars of one's own, partial bins, coverage and lookbacks.
- [Learners](https://gillescolling.com/timesift/articles/learners.html): the shipped learners,
training the encoders, tuning, and registering a learner, a response or a metric.
- [Finding the grain](https://gillescolling.com/timesift/articles/comparing-grains.html): the
ladder, the contrasts and the nested selection on a simulated record whose grain is known.
- [Ensembles, thresholds and maps](https://gillescolling.com/timesift/articles/prediction.html):
spatial and environmental folds, the combiners, thresholds and projection onto rasters.
- [Coming from biomod2](https://gillescolling.com/timesift/articles/biomod2.html) and
[Troubleshooting](https://gillescolling.com/timesift/articles/troubleshooting.html).