Package 'timesift'

Description

One additive model per response, a smooth function of every bin-by-channel column of the representation, fitted as mgcv's gam(y ~ s(x1) + s(x2) + ..., method = "GCV.Cp") fits it and as biomod2 fits GAM. Each column enters as a thin plate regression spline of k basis functions (Wood 2003): a cubic radial function centred on each of the column's distinct values, reduced to its k - 2 directions of greatest eigenvalue, together with the linear function, which the penalty on the spline's squared second derivative leaves free. Each smooth sums to zero over the units, beside one intercept.

Usage

additive(data = NULL, k = 10L, gamma = 1, max_knots = 2000L, threads = 1L)

Arguments

Details

The coefficients maximise the penalised likelihood, by penalised iteratively reweighted least squares, and the smoothing parameters, one per column, minimise the unbiased risk estimator under presence-absence and a count head and the generalised cross-validation score under a squared-error head, by Newton's method with the exact derivatives (Wood 2008). gamma multiplies the charge each effective degree of freedom adds to that criterion, as mgcv's gamma does; above one it gives smoother fits. Over columns as alike as neighbouring weeks the criterion can have more than one local minimum; the one the search settles in then depends on where it starts, and mgcv, starting from a rule of its own parametrisation, can settle in another.

The defaults are mgcv's, which biomod2 passes unchanged: k = 10 and a criterion charge of one. Above max_knots distinct values, the radial functions are centred on that many of them, drawn as mgcv draws them, so that the basis is the one mgcv fits. A column of fewer than k distinct values takes as many basis functions as it holds values, a column of two enters linearly, and a column of one is left out; mgcv refuses the first and the last of those. A column whose linear part the columns before it already span keeps its penalised part and loses the linear one.

The case weights are the response head's, positive_weights() under presence-absence, and enter the likelihood as mgcv's prior weights. The model holds at most as many coefficients as there are units, as mgcv's does: one for the intercept and k - 1 per column. Under the shipped presence-absence head those weights are on, so a default additive() is mgcv's specification fitted under them; a head registered without weights fits it unweighted.

The basis, the fit and the choice of the smoothing parameters live in the core the Python package calls, pinned against mgcv in the fixtures, so the two languages fit the same model. A response holding one value is predicted its mean and named in unfitted; one whose smoothing parameter search stopped short of its tolerance is named in stopped.

Value

A learner().

Examples

additive()
additive(data = grain("season"), k = 5)

Description

The representation is built twice, once per language, from one shared core. The response matrix, the fold map and the mask of scorable cells are not: they are built once and read wherever they are needed, because a fold map drawn from a seed in R and one drawn from the same seed in Python are different maps, and aligning the two random streams would be the wrong fix.

Usage

write_folds(x, file)

read_folds(file, units = NULL)

write_response(y, file)

read_response(file, units = NULL)

write_cells(cells, file)

read_cells(file)

Arguments

Details

These six functions are the handover. The Python side carries write_folds() and its siblings under the same names; they write the same bytes from the same artifact and read the same files, so a fold map built in either language is usable in the other without the caller knowing the format.

The format is normative and is given in inst/spec/representation.md: CSV, UTF-8, a header row, no quoting, LF line endings on every platform, numbers at twelve significant digits, and rows ordered by the identifier under C collation.

Value

The writers return file, invisibly. read_folds() returns a named integer vector of class timesift_folds, read_response() a numeric matrix with the units in its row names, and read_cells() a timesift_cells data frame.

Examples

set.seed(1)
y <- matrix(rbinom(120, 1, 0.3), nrow = 20,
            dimnames = list(sprintf("p%02d", 1:20), paste0("sp", 1:6)))
f <- fold_map(y, v = 4)

path <- tempfile(fileext = ".csv")
write_folds(f, path)
back <- read_folds(path, names(f))
identical(names(back), names(f)) && identical(as.integer(back), as.integer(f))
unlink(path)

Description

Lays a ⁠[unit, bin, channel]⁠ array out as one row per unit and one column per bin of each channel, which is what a model reading a table of predictors takes. It is what the learners that read a tabular block fit on, so a learner of your own that calls it reads the same columns they do.

Usage

## S3 method for class 'timesift_matrix'
as.matrix(x, ...)

Arguments

Details

Columns are named channel@bin, channel by channel and bin by bin within each. A channel that holds the same number in every bin, as a static predictor does, is one column named by the channel alone.

Value

A numeric ⁠[unit, predictor]⁠ matrix, the units as row names.

Examples

t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 60)
d <- data.frame(plot = rep(c("a", "b"), each = length(t)), t = rep(t, 2),
                temp = sin(seq_len(2 * length(t)) / 24))
x <- grain_matrix(d, plot, t, temp, grain = "month", stats = c("min", "max"))
as.matrix(x)

Description

The area under the precision-recall curve, as the step sum over the distinct predictions: at each, the precision of calling every unit at or above it a presence, weighted by the share of presences that cut adds. Units sharing a prediction enter together, so the order they arrived in does not move the score. Its floor is the prevalence rather than one half, which makes it the reading of how well presences are ranked above absences where presences are rare, and it is reported beside roc_auc() for that reason.

Usage

average_precision(y, p)

Arguments

Value

One number, or NA where the cell defines none.

Examples

average_precision(c(0, 0, 1, 1), c(0.1, 0.2, 0.8, 0.9))
average_precision(c(0, 1, 0, 1), c(0.1, 0.2, 0.8, 0.9))

Description

Joins representations of the same units and bins into one array, in the order given. It is how a temperature reading, an external product such as snow cover, and the calendar position of each bin reach a model as one input.

Usage

bind_channels(...)

Arguments

Value

One array carrying every channel, in the order the arguments are given and, inside each argument, in its own channel order. The attributes are the first argument's, with stats the joined names, except static and position, which name channels and so name those of every argument.

Examples

t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 60)
d <- data.frame(plot = rep(c("a", "b"), each = length(t)), t = rep(t, 2),
                temp = rnorm(2 * length(t)))
x <- grain_matrix(d, plot, t, temp, grain = "week")
dimnames(bind_channels(x, calendar_channels(x)))[[3]]

Description

One boosted model per response, over every bin-by-channel column of the representation: a logistic model under a presence-absence head, a squared-error one under a head with a squared-error loss and a Poisson one with a log link under a count head. The score starts at the log-odds of the weighted share of presences, the weighted mean, or the log of the weighted mean count, and each tree is fitted to the loss's gradient at the current score and added to it scaled by shrinkage. Each tree is grown on a subsample of the units drawn without replacement, and reads a subsample of the columns.

Usage

boosting(
  data = NULL,
  method = c("gbm", "xgboost"),
  trees = NULL,
  depth = NULL,
  shrinkage = NULL,
  min_leaf = NULL,
  subsample = NULL,
  colsample = NULL,
  lambda = NULL,
  gamma = NULL,
  n_inner = NULL,
  preset = c("default", "bigboss"),
  seed = 1L,
  threads = 1L
)

Arguments

Details

method picks the trees. Under "gbm" they are Friedman's gradient boosting machine as the gbm package grows it, which is what biomod2 fits as GBM: depth splits grown best first, each the one that most reduces the weighted squared error of the working response with at least min_leaf units on each side, and a leaf that takes one Newton step on the loss. Under "xgboost" they are XGBoost's exact greedy trees (Chen and Guestrin 2016), which biomod2 fits as XGBOOST: grown level by level to depth, each split chosen by the second-order gain under the L2 penalty lambda with at least min_leaf of hessian on each side, pruned where a split gains less than gamma, and a leaf the step -G / (H + lambda). Either way depth is the order of interaction a tree can hold. With subsample = 1 the first-order fit is gbm's own to rounding, and the second-order one xgboost's to its single-precision storage; the model is grown by the core the Python package calls, so the two languages fit the same model.

n_inner folds, when above zero, choose how many trees are kept: the fit is repeated on each fold's complement, and the number of trees of least held-out deviance, summed over the folds and weighted by how many units each holds, is kept, as gbm's cv.folds chooses it. The folds are dealt for each response and stratified on it, as the elastic net's are.

preset says whose defaults the settings left NULL take. "default" is the fitting package's own, which is what biomod2's default option set fits: under gbm 100 trees of one split, shrinkage = 0.1, min_leaf = 10 and subsample = 0.5; under xgboost 100 trees of depth 6, shrinkage = 0.3, min_leaf = 1, lambda = 1 and every unit and column. "bigboss" is biomod2's tuned option set: under gbm 2500 trees of seven splits, shrinkage = 0.001, min_leaf = 5, subsample = 0.5 and three inner folds; under xgboost four trees of depth 2 at shrinkage = 1. A setting given explicitly beats either.

The case weights are the response head's, positive_weights() under presence-absence, and weigh the gradient and every sum a tree is grown on; min_leaf counts units under gbm, as n.minobsinnode does. Under the shipped presence-absence head those weights are on, so a default boosting() is gbm's specification fitted under them; a head registered without weights fits it unweighted.

Value

A learner().

Examples

boosting()
boosting(preset = "bigboss")
boosting(method = "xgboost", depth = 3L)

Description

How consistently the predictions rise with the density of presences, without a cut and without absences (Hirzel et al. 2006). A window of width times the range of the predictions slides over that range in resolution + 1 equal steps. In each, the share of presences falling inside divided by the share of all units falling inside is the predicted-to-expected ratio, undefined where no unit falls inside, and the index is the Spearman correlation of that ratio with the window's midpoint, between -1 and 1. Predictions that rank no better than chance read near 0.

Usage

boyce_index(y, p, resolution = 100L, width = 0.1)

Arguments

Details

The units of p are the background, which is how biomod2 reads BOYCE on a fit. The windows are closed at both ends and the ratio is not thinned of repeated values.

Value

One number, or NA where the cell defines none, where the predictions are all equal, or where the ratio is defined in fewer than three windows or takes one value.

Examples

boyce_index(c(0, 0, 0, 1, 1, 1, 0, 1), c(0.10, 0.20, 0.35, 0.40, 0.60, 0.90, 0.55, 0.70))

Description

Turns a representation and the two tables into the ⁠[target, bin, channel]⁠ array a learner is fitted on. It is the one place the fitting layer builds an array, so timesift() and predict.timesift() reach a record the same way and a candidate refitted on new targets is built from the settings its own arm was.

Usage

build_representation(rep, series, targets, spec)

Arguments

Details

The rows are the targets, in the order the fitting layer keeps them: sorted by identifier where one target row belongs to each unit, and in the targets' own order where target_time anchors them. Columns named in static are appended as channels holding one value per target, constant across the bins.

Value

A timesift_matrix.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 90)
units <- sprintf("p%02d", 1:30)
warmth <- rnorm(30)
logger <- data.frame(
  plot = rep(units, each = length(t)), datetime = rep(t, 30),
  temp = as.numeric(vapply(warmth, function(w) w + sin(seq_along(t) / 300), numeric(length(t)))))
plots <- data.frame(plot = units, sp_a = rbinom(30, 1, plogis(2 * warmth)))
fit <- timesift(plots, logger, y = sp_a, id = plot, time = datetime,
                sift = grains("month"), resampling = cv(v = 3L), ensemble = FALSE,
                verbose = FALSE)
x <- build_representation(grain("week"), logger, plots, fit$spec)
dim(x)

Description

An encoder that ends in global pooling discards when a thermal event happened, so the position of a bin in the year has to be given to it as input if it is to be used at all. These channels carry that position as the sine and cosine of the bin's fractional place in the year, which is continuous across the turn of the year where the fraction itself is not. On a record read finer than a day, the same pair for the place in the day carries where a reading sits in the daily cycle.

Usage

calendar_channels(x, cycles = "year")

Arguments

Details

They are the time index of each bin, not a summary of the readings, so adding them introduces no hand-built thermal feature: whatever a model does with them it could have done with a calendar.

The position is read at the midpoint of the record each bin holds, in UTC, so a bin the record only partly covers sits at the phase it was actually measured over. A site's longitude and a zone's offset move the day's phase by the same amount for every bin, which a model absorbs, where a local clock's summer time would move it twice a year. inst/spec/representation.md is the normative description.

Value

An array of the same units and bins with two channels per cycle, year_sin and year_cos, day_sin and day_cos, identical across units. Combine it with the readings using bind_channels(). The channels are recorded as positions in the position attribute, and the encoders of torch_learners read them at their own amplitude rather than standardise them.

Examples

t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 400)
d <- data.frame(plot = "a", t = t, temp = sin(seq_along(t) / 500))
x <- grain_matrix(d, plot, t, temp, grain = "month")
round(calendar_channels(x)[1, 1:4, ], 3)

hourly <- grain_matrix(d, plot, t, temp, grain = "native")
round(calendar_channels(hourly, cycles = c("year", "day"))[1, 1:4, ], 3)

Description

c() on learners is the set of them, and on representations the set of those. A set handed to c() again splices, so a set can be added to rather than rewritten, which is what list() cannot do: c(base, cnn()) where base is already a set.

Usage

## S3 method for class 'timesift_learner'
c(...)

## S3 method for class 'timesift_models'
c(...)

## S3 method for class 'timesift_representation'
c(...)

## S3 method for class 'timesift_sift'
c(...)

Arguments

Details

learners and sift take either form. A length-one string is the name of a registered learner.

Value

A timesift_models for learners and a timesift_sift for representations.

Examples

base <- c(elasticnet(), forest())
base
c(base, linear())

c(grains("day", "week"), lookback("30 days"))

Description

A representation needs every unit in every bin, and grain_matrix() refuses a record where one is missing rather than pad it. This is the same binning laid out so the gaps can be read: how many readings each unit has in each bin, over every bin the calendar tiles the record with from the first bin any unit touches to the last. A logger that started late, stopped early or lost a month is a row with zeros in it; a bin the whole record skips is a column of zeros.

Usage

coverage(data, id, time, grain = "day", year_start = "01-01")

Arguments

Details

What to do about a gap is the analyst's decision, and this is the table it is made on: drop the units that do not span the record, cut the record to the span every unit covers, or move to a grain the gap does not reach. Nothing here fills a cell.

Value

An integer matrix of reading counts, one row per unit and one column per bin, of class timesift_coverage, with the units and the ISO-8601 bin starts as dimnames and the grain and the bin_start instants as attributes.

Examples

t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 40)
d <- data.frame(plot = rep(c("a", "b"), each = length(t)), t = rep(t, 2),
                temp = rnorm(2 * length(t)))
# Unit b loses the calendar week beginning Monday 6 September.
lost <- d$plot == "b" & d$t >= as.POSIXct("2021-09-06", tz = "UTC") &
  d$t < as.POSIXct("2021-09-13", tz = "UTC")
coverage(d[!lost, ], plot, t, grain = "week")

Description

The resampling timesift() scores on. cv() deals units into folds balanced on the stratifying value; grouped_cv() deals whole groups, keeping every target sharing a group value on one side of each split. block_cv() and env_cv() hold out a block of units whole: a region of space for the first, a region of predictor space for the second.

Usage

cv(v = 10L, seed = 1L, strata = 5L, by = NULL, repeats = 1L)

block_cv(by, v = 4L)

env_cv(by, v = 4L)

grouped_cv(group, v = 10L, seed = 1L, repeats = 1L)

Arguments

Details

The blocks are cut by halving. The units are split at the median of the column with the widest range (the first on a tie), the lower share going to the left and the units tied on that column ordered as they arrived, and each part is cut again until there are v of them, so a block holds as many units as another to within one and a map of v blocks has v folds. Nothing is drawn, so both languages return the same map. env_cv() centres and scales every column first, which a coordinate system does not need. Under the nested selection the inner folds keep blocks whole as well, so n_inner is at most v - 1.

resampling also accepts a fold vector or a fold_map() result directly, which is how a split the package has no constructor for – a spatial block, a season held out whole – reaches the same fitting path.

Value

A timesift_resampling.

Examples

cv(v = 5L)
cv(v = 5L, by = "elevation")
grouped_cv("site")
block_cv(c("x", "y"), v = 4L)
env_cv(c("elevation", "slope"), v = 5L)

Description

The MD5 the contract in inst/spec/representation.md defines, and what the fixtures in ⁠inst/spec/fixtures/⁠ are checked against. timesift.digest_array() on the Python side is the same function: the same array gives the same string in either language, which is what makes a digest a statement about the representation rather than about the machine.

Usage

digest_array(x)

Arguments

Details

Values are traversed in the array's own order, unit fastest, then bin, then channel; each is formatted to twelve decimal places, joined with a line feed and terminated with one, and the UTF-8 bytes of that are hashed. The line ending is a line feed on every platform. It is written explicitly because writeLines() emits CRLF on Windows, which would make a digest depend on the machine that produced it.

The array has to be finite. ⁠%.12f⁠ writes an infinity as Inf here and as inf in Python, and R tells a missing value apart from a not-a-number where Python has one spelling for both, so a digest over such an array would say something about the language rather than about the representation. A representation never holds one: a reading that is not a finite number is refused where the record is read.

Value

The digest, a single string of 32 hexadecimal characters.

Examples

digest_array(array(c(1, -0.5), dim = c(2L, 1L, 1L)))

Description

One discriminant per response, over every bin-by-channel column of the representation, fitted as mda's fda(method = mars) fits it and as biomod2 fits FDA. Optimal scoring gives presence and absence one score each and regresses the scored response on a MARS basis of the columns; the fitted score is the one canonical variate, and a unit's prediction is its posterior probability of presence under two normal classes around the class centroids on that variate, with the classes' shares among the fitting units as priors.

Usage

discriminant(
  data = NULL,
  degree = 1L,
  penalty = NULL,
  max_terms = NULL,
  min_gain = 0.001,
  prune = TRUE,
  calibrate = TRUE,
  threads = 1L
)

Arguments

Details

The basis is mda's own MARS, which biomod2 reaches through fda() and which differs from earth's that mars() reproduces. Each forward step adds a column linearly or a pair of hinges on it, max(0, x - t) and max(0, t - x), choosing by Friedman's running updates, and the pass stops when a step lowers the residuals by less than min_gain of them, when they fall to min_gain of the null model's, when the generalised cross-validation passes ten times the null model's, or at max_terms terms. The pruning pass drops the term of least t statistic, one at a time, and keeps the subset of least generalised cross-validation, which counts each term beyond the intercept as 1 + penalty / 2 degrees of freedom. The defaults are mda's, which biomod2 passes unchanged under its default option set and under "bigboss": degree one, penalty = 2, min_gain = 0.001 and ⁠max_terms = max(21, 2 p + 1)⁠ for p columns.

The case weights are the response head's, positive_weights() under presence-absence. They set the classes' scores and the variate, as fda() reads them, and not the basis: mda's forward pass sets them to one. Under the shipped presence-absence head those weights are on, so a default discriminant() is mda's specification fitted under them; a head registered without weights fits it unweighted.

biomod2 always recalibrates an FDA posterior by a probit regression of the response on it, under the case weights, and predicts through that regression. calibrate = TRUE does the same, on the units the discriminant was fitted on; FALSE predicts the posterior itself. biomod2 rounds the posterior to three decimals before recalibrating it, which is not reproduced.

The forward pass, the pruning, the scoring and the recalibration live in the core the Python package calls, pinned against mda and R's glm() in the fixtures, so the two languages keep the same terms and predict the same probabilities. A response holding one value is predicted its mean, and so is one whose scored response the basis does not reach; both are named in unfitted. The learner needs a presence-absence response, and a head whose loss is not the binary cross-entropy is refused.

Value

A learner().

Examples

discriminant()
discriminant(data = grain("season"), calibrate = FALSE)

Description

One elastic net per variable, over every bin-by-channel column of the representation and, by default, their squares. There is no discrete selection step: the penalty path uses every column and shrinks, and the penalty itself is chosen by an inner cross-validation on the fitting units, so nothing about the model is decided outside the fold it is fitted in.

Usage

elasticnet(
  data = NULL,
  alpha = 0.5,
  n_inner = 10L,
  squares = TRUE,
  s = "lambda.1se",
  n_lambda = 100L,
  tol = 1e-08,
  threads = 1L,
  seed = 1L
)

Arguments

Details

The family is the response head's: a binary cross-entropy loss fits a logistic model, a squared-error loss a linear one and a Poisson-deviance loss a log-linear one, so the learner is the same under a presence-absence head, a continuous one and a count one. So are the case weights: the head's weights, positive_weights() for presence-absence, are what every learner that ships fits under.

The inner folds are dealt for each response and stratified on it, so a rare outcome is spread over them as evenly as its count allows. A presence-absence response whose inner training sets cannot each hold two of each outcome, the fewest a logistic path is fitted to, has too few of one outcome to choose a penalty on. It is predicted its share among the fitting units, as a response holding one outcome is, and the fit names every such response in unfitted.

A reweighted fit that does not settle at a penalty ends its path there and keeps the points before it, which is what glmnet does with the same event, and the penalty is chosen over the points fitted. It happens where a rare outcome is nearly separable at the small end of the path. The fit names every response whose path on the fitting units, or on any inner fold, ended that way in stopped.

This is the aggregate-feature side of the comparison the package was built for, and it is the fair opponent for a network: a per-fold discrete selector pays selection variance a network never pays, so beating that one is not a matched result.

The path is fitted by the same core the Python package calls, so the two return the same coefficients for the same input. Its conventions are glmnet's, which is what the arm is measured against: weights normalised to sum to one, columns centred and scaled by their weighted mean and weighted standard deviation, a hundred penalties down from the smallest that leaves every coefficient at zero, and the held-out deviance read fold by fold.

Value

A learner().

Examples

elasticnet(alpha = 0.5)

Description

Every candidate emits an out-of-fold prediction for every scorable cell over the same folds, so the combination is arithmetic on those predictions and nothing else. ensemble() says which arithmetic.

Usage

ensemble(
  method = c("stack", "mean", "median", "weighted", "committee"),
  scope = c("all", "learners", "representations"),
  metric = NULL,
  response = NULL,
  min_score = NULL,
  decay = NULL,
  rule = NULL
)

Arguments

Details

"stack" fits non-negative weights summing to one on the out-of-fold predictions alone, never on in-sample ones, minimising the response head's loss over the scorable cells: binomial deviance for presence-absence. One weight vector covers every response, because per-response weights would be fitted on the handful of cells a rare response has. "mean" and "median" combine without fitting anything. "weighted" takes each candidate's own mean score, keeps its non-negative part and rescales those to sum to one, so a candidate scoring at or below zero is left out and the rest are weighted by how well they scored; where no candidate scores above zero the scores order nothing worth weighting by, and every candidate weighs the same.

decay changes how "weighted" turns scores into weights, as biomod2's EMwmean.decay does. Given a number d, the candidates scoring above zero are ranked from the best down, the one ranked r of K takes d^(K - r + 1), candidates on the same score share the mean of their ranks' weights, and the weights are rescaled to sum to one. Each rank down weighs 1 / d of the one above it, whatever the gap in score between them.

"committee" is biomod2's committee averaging. Each member cuts each response at its own threshold, learned by decision_threshold() under rule from that member's out-of-fold predictions of every target, and the combined prediction is the share of members voting presence. A member whose predictions of a response give no cut does not vote on it. Under timesift() the thresholds of each outer fold are learned from the inner out-of-fold predictions of its training targets, as a stack's weights are, so the estimate never reads a cut chosen on the targets it scores.

min_score is biomod2's metric.select.thresh: a candidate whose mean score is below it is not eligible, and the filter is applied before scope picks among what is left.

scope says which candidates are eligible. "all" is every candidate. "learners" keeps the several learners that read the representation of the best-scoring candidate, and "representations" keeps the one learner of the best-scoring candidate across the representations it ran on; both are read off the same mean scores the report shows.

Value

A timesift_ensemble.

Examples

ensemble()
ensemble("weighted", scope = "learners")
ensemble("weighted", decay = 1.6, min_score = 0.4)
ensemble("committee", rule = "kappa")

Description

Combine one prediction per member into one prediction

Usage

ensemble_combine(stack, preds)

Arguments

Value

One ⁠[target, response]⁠ matrix.

Examples

set.seed(1)
y <- matrix(rbinom(200, 1, 0.4), nrow = 50,
            dimnames = list(sprintf("p%02d", 1:50), paste0("sp", 1:4)))
folds <- fold_map(y, v = 5)
truth <- matrix(runif(200), nrow = 50, dimnames = dimnames(y))
oof <- list(good = 0.8 * y + 0.2 * truth, noise = truth)
st <- ensemble_fit(oof, y, scorable_cells(y, folds), folds)
dim(ensemble_combine(st, oof))

Description

The combiner sees the out-of-fold predictions, the response, the mask of scorable cells and the fold map, and never a model, so it cannot read anything a candidate fitted in-sample.

Usage

ensemble_fit(oof, y, cells, folds, spec = ensemble(), scores = NULL)

Arguments

Details

The weights are fitted to the response on those predictions, so the combined prediction scored against the same response is scored on the data its weights were fitted to, and that score is optimistic. timesift() evaluates the stack the other way: each outer fold's weights are fitted on inner out-of-fold predictions of its training targets and applied to the outer test fold.

Value

A timesift_stack: method, the named weights, and what they were fitted on.

Examples

set.seed(1)
y <- matrix(rbinom(200, 1, 0.4), nrow = 50,
            dimnames = list(sprintf("p%02d", 1:50), paste0("sp", 1:4)))
folds <- fold_map(y, v = 5)
truth <- matrix(runif(200), nrow = 50, dimnames = dimnames(y))
oof <- list(good = 0.8 * y + 0.2 * truth, noise = truth)
ensemble_fit(oof, y, scorable_cells(y, folds), folds)

Description

biomod2's EMcv and EMci, read on the members' predictions under the weights the stack carries. For each target and response: the weighted mean m of the members' predictions; their weighted standard deviation s, the square root of sum(w * (p - m)^2) / (1 - sum(w^2)), which is the sample standard deviation when the weights are equal; the coefficient of variation s / m; and the interval m -+ qt(1 - alpha / 2, n - 1) * s * sqrt(sum(w^2)), n the number of members carrying weight, which is the t interval of a mean of n members when the weights are equal. Under a response head whose predictions are probabilities, the interval is held inside zero and one. An uncertainty map is this on one target per map cell.

Usage

ensemble_spread(stack, preds, alpha = 0.05)

Arguments

Details

A committee's and a median's members are read at equal weight, and a committee's spread is that of the members' predictions rather than of their votes.

Value

A ⁠[target, response, statistic]⁠ array, the statistics being mean, sd, cv, lower and upper. sd, cv and the interval are NA where fewer than two members carry weight.

Examples

set.seed(1)
y <- matrix(rbinom(200, 1, 0.4), nrow = 50,
            dimnames = list(sprintf("p%02d", 1:50), paste0("sp", 1:4)))
folds <- fold_map(y, v = 5)
truth <- matrix(runif(200), nrow = 50, dimnames = dimnames(y))
oof <- list(good = 0.8 * y + 0.2 * truth, fair = 0.6 * y + 0.4 * truth, noise = truth)
st <- ensemble_fit(oof, y, scorable_cells(y, folds), folds, ensemble("mean"))
ensemble_spread(st, oof)[1:3, "sp1", ]

Description

The weights the combiner fitted

Usage

ensemble_weights(fit)

Arguments

Value

A named numeric vector, or NULL where the run fitted no combiner.

Examples

set.seed(1)
y <- matrix(rbinom(200, 1, 0.4), nrow = 50,
            dimnames = list(sprintf("p%02d", 1:50), paste0("sp", 1:4)))
folds <- fold_map(y, v = 5)
truth <- matrix(runif(200), nrow = 50, dimnames = dimnames(y))
oof <- list(good = 0.8 * y + 0.2 * truth, noise = truth)
ensemble_weights(ensemble_fit(oof, y, scorable_cells(y, folds), folds))

Description

One envelope per response, over every bin-by-channel column of the representation: for each column, the quantile and 1 - quantile quantiles of its readings over the units present, and a unit predicted present where every column lies between its two, the ends included. It is biomod2's SRE, and with the same quantile it draws the envelope bm_SRE() draws; the quantile is R's default, type 7.

Usage

envelope(data = NULL, quantile = 0.025)

Arguments

Details

The prediction is zero or one. It enters an ensemble as that, and the combiner weighs it by its held-out loss like any other candidate; a threshold read off it has two values to choose from, so its TSS is the one its own zeros and ones give.

An envelope reads the presences and nothing else. The absences do not move it, and neither do the response head's case weights, which it takes no notice of.

Every column has to agree for a unit to be inside, so the more columns a representation has the fewer units any envelope holds: at the default quantile each column shuts out about one presence in twenty, and a weekly representation of three years has 157 columns per channel. A coarse grain is what an envelope is meant for, and pinning it there with data = grain("season") keeps it there while the rest of the run reads finer ones.

A response with no presence, or with nothing else, is predicted its share among the fitting units, and the fit names it in unfitted. The learner needs a presence-absence response, and a head whose loss is not the binary cross-entropy is refused.

Value

A learner().

Examples

envelope()
envelope(data = grain("season"), quantile = 0.05)

Description

A representation the package did not build, such as a published set of hand-aggregated climate summaries, enters here. It becomes a one-channel ⁠[unit, feature, 1]⁠ array, which is what a learner reads, so a feature table and a temporal grain can be arms of the same grain_ladder() and be scored on the same cells by the same rule.

Usage

feature_matrix(m, label = "features")

Arguments

Details

It carries no time axis, because it has none: the reduction already happened, elsewhere, and what reaches the model is a list of numbers per unit. That is the whole point of comparing against it.

Value

A timesift_matrix of shape ⁠[unit, feature, 1]⁠.

Examples

m <- matrix(rnorm(30), nrow = 10,
            dimnames = list(sprintf("p%02d", 1:10), paste0("bio", 1:3)))
feature_matrix(m)

Description

Fit one learner at one grain

Usage

fit_learner(
  learner,
  x,
  y,
  response = "presence_absence",
  control = NULL,
  group = NULL,
  ...
)

## S3 method for class 'timesift_fit'
predict(object, newdata, ...)

Arguments

Value

A timesift_fit, which stats::predict() takes a new representation.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 120)
units <- sprintf("p%02d", 1:40)
d <- data.frame(plot = rep(units, each = length(t)), t = rep(t, length(units)),
                temp = as.numeric(replicate(length(units), rnorm(length(t)))))
x <- grain_matrix(d, plot, t, temp, grain = "month")
y <- matrix(rbinom(80, 1, 0.4), nrow = 40, dimnames = list(units, c("sp1", "sp2")))
fit <- fit_learner(elasticnet(), x, y)
dim(stats::predict(fit, x))

Description

One fold map, built once and read by everything that scores. Every learner in a ladder is then fitted and scored on identical splits, which is what makes the comparison between them paired rather than a comparison of two clouds of numbers.

Usage

fold_map(y, v = 10L, seed = 1L, strata = 5L, by = NULL, group = NULL)

Arguments

Details

Units are held out singly. Where the input a model reads is measured at the unit itself, as a logger in each plot is, a held-out unit brings its own measured input with it and nothing of its neighbours' reaches the model.

Folds are balanced within strata: units are grouped into strata equal-count groups of the stratifying value, shuffled inside each group, and dealt round-robin by one counter that runs on from each stratum into the next, so each fold carries the same mix and the folds are equal in size to within one unit whatever v is, up to one unit per fold. With a multi-variable response the default stratifies on richness, the number of variables present at a unit, because one fold map has to serve every variable at once and cannot be stratified on any single one of them.

Value

An integer vector of fold numbers named by unit, of class timesift_folds. Any named integer vector of the same shape is accepted wherever this is.

Examples

y <- matrix(rbinom(300, 1, 0.3), nrow = 60,
            dimnames = list(sprintf("p%02d", 1:60), paste0("sp", 1:5)))
f <- fold_map(y, v = 5)
table(f)

Description

One forest per response, over every bin-by-channel column of the representation: a probability forest under a presence-absence head and a regression forest under a head with a squared-error loss or a count head, a count being cut on its variance and a leaf reporting its mean. Trees split on one column at a time and pay nothing for columns that carry nothing, so a forest reads a wide tabular representation without a penalty path and without a selection step, and it finds an interaction between two bins that a linear model would need the product term for.

Usage

forest(
  data = NULL,
  trees = NULL,
  mtry = NULL,
  min_node = NULL,
  balance = FALSE,
  preset = c("default", "bigboss"),
  seed = 1L,
  threads = 1L
)

Arguments

Details

Each tree is grown on a bootstrap draw of the units, as many draws as there are units, and each node is split on the best of mtry columns drawn for it, by the Gini index or the sum of squares tree() splits by. A tree is grown out: a node is split while it holds at least twice min_node units and its responses differ. A leaf reports the share of presences among the draws it holds, or their mean, and the forest the mean over its trees. The forest is grown by the core the Python package calls, and every draw comes from one generator seeded per tree, so the two languages grow the same forest on any number of threads.

balance = TRUE is the down-sampled forest biomod2 fits as RFd: each tree draws as many units from each class as the smaller class holds, so a rare response's presences are half of every tree's draw.

preset says whose defaults the settings left NULL take. "default" is randomForest's own, which is what biomod2's default option set fits: 500 trees, mtry the square root of the column count under presence-absence and a third of it under a squared-error loss, and min_node 1 and 5 under the two. "bigboss" is biomod2's tuned option set: 500 trees, mtry = 2 and min_node = 5. A setting given explicitly beats either, and an mtry above the column count is the column count.

The case weights are the response head's, positive_weights() under presence-absence, and weight the bootstrap draw: a unit is drawn in proportion to its weight, and within its class under balance. Under the shipped presence-absence head those weights are on, so a default forest() is randomForest's specification fitted under them; a head registered without weights fits it unweighted.

Value

A learner().

Examples

forest(trees = 200L)
forest(preset = "bigboss", balance = TRUE)

Description

A mixed model on the per-cell scores of one learner, score ~ grain + (1 | variable) + (1 | fold), fitted by restricted maximum likelihood. Every grain is then compared against the reference by Dunnett's many-to-one procedure, which corrects for the comparisons made without correcting for pairs nobody asked about.

Usage

grain_contrasts(ladder, learner = NULL, reference = NULL, adjust = "mvt")

Arguments

Details

The design is balanced across variables, folds and grains, so each grain is compared within a variable and within a fold and the variation between variables cancels from the comparison. That is what lets a difference of 0.015 hold up where absolute skill ranges across variables by ten times as much.

Taking the best-observed grain as the reference is a choice that favours the reference, so read this beside the paired differences paired_contrast() gives, which single out no grain.

Value

A data frame of one row per grain: the estimated marginal mean difference from the reference, its interval, and the adjusted p-value. Needs lme4, lmerTest and emmeans.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 200)
units <- sprintf("p%03d", 1:60)
warmth <- rnorm(60)
d <- data.frame(
  plot = rep(units, each = length(t)), t = rep(t, length(units)),
  temp = as.numeric(vapply(warmth, function(w) 1.5 * w + rnorm(length(t), sd = 20),
                           numeric(length(t)))))
y <- matrix(rbinom(60 * 6, 1, plogis(3 * as.numeric(outer(warmth, rep(c(1, -1), 3))))),
            ncol = 6, dimnames = list(units, paste0("sp", 1:6)))
x <- grain_matrix(d, plot, t, temp, grain = c("day", "week", "month"))
lad <- grain_ladder(x, y, elasticnet(), folds = fold_map(y, v = 5), verbose = FALSE)
grain_contrasts(lad)

Description

Cross-validates every learner at every grain of a representation set, on one fold map and one mask of scorable cells, and returns the score of each ⁠(grain, learner, variable, fold)⁠ cell. It is the measurement the package exists for: how much of a record a model needs, read off the point where making the record finer stops paying.

Usage

grain_ladder(
  x,
  y,
  learners,
  folds = NULL,
  response = "presence_absence",
  metric = NULL,
  control = train_control(),
  keep_fits = FALSE,
  interval = c("variables", "nested_cv"),
  repeats = 1L,
  seed = 1L,
  verbose = TRUE
)

## S3 method for class 'timesift_ladder'
summary(object, ...)

Arguments

Details

Every arm sees identical splits and is restricted to identical cells, so the arms' means share a denominator and any two of them can be compared cell by cell with paired_contrast().

Value

A data frame of one row per scored cell, of class timesift_ladder, carrying the grain, the learner, the variable, the fold and the score. The held-out prediction of every unit is kept in the predictions attribute, and the scorable-cell mask in cells.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 200)
units <- sprintf("p%02d", 1:60)
warmth <- rnorm(60)
d <- data.frame(
  plot = rep(units, each = length(t)), t = rep(t, length(units)),
  temp = as.numeric(vapply(warmth, function(w) w + sin(seq_along(t) / 300) + rnorm(length(t)),
                           numeric(length(t)))))
y <- matrix(rbinom(120, 1, plogis(c(warmth, -warmth))), nrow = 60,
            dimnames = list(units, c("sp1", "sp2")))
x <- grain_matrix(d, plot, t, temp, grain = c("week", "month"))
lad <- grain_ladder(x, y, elasticnet(), folds = fold_map(y, v = 3), verbose = FALSE)
summary(lad)

Description

Bins each unit's readings by the calendar and summarises every bin by one or more statistics, returning the array a model is fitted on. The reduction is the choice the package exists to make explicit: grain sets how coarse the record becomes, stats sets what survives the reduction, and the two are not interchangeable.

Usage

grain_matrix(
  data,
  id,
  time,
  value,
  grain = "day",
  stats = "mean",
  year_start = "01-01",
  partial = c("keep", "drop")
)

Arguments

Details

Seven statistics are available, and the distinction between an extreme reading, an extreme day and a typical day is deliberate rather than pedantic:

• mean: arithmetic mean of the readings in the bin.

• min, max: coldest and warmest single reading in the bin.

• cold_day, warm_day: coldest and warmest day, each day first reduced to its own mean. Defined for "day" and coarser.

• mean_daily_min, mean_daily_max: the bin's average daily minimum and average daily maximum, each day first reduced to its own extreme. Defined for "day" and coarser.

An extreme day is a state the unit was in; an extreme reading can be one hour; an average daily extreme is the exposure a typical day of the bin brought. On alpine soil temperature the day-level pair carries more predictive signal than the bin mean, and by more the coarser the bin.

Whether a grain is a day or coarser is decided from the bins rather than from the grain's name, so a supplied calendar that cuts inside a day is refused for these four as well, naming the day it splits.

Nothing is standardised here. Scaling belongs to the fold it is computed on, never to the representation, because computing it over all units would leak held-out units into the input.

Value

A numeric array of shape ⁠[unit, bin, channel]⁠, with dimnames giving the sorted unit identifiers, the ISO-8601 start of each bin, and the statistic names. Attributes:

• grain: the grain name.

• stats: the statistic names in channel order.

• year_start: the boundary used.

• bin_start: the instant each bin begins on the calendar, which is earlier than the bin's first reading wherever the record does not reach the boundary.

• bin_end: the last reading instant assigned to each bin.

• bin_n: a ⁠[unit, bin]⁠ matrix of how many readings fell in each bin.

• bin_partial: a logical vector marking the bins the record does not cover for their whole calendar span.

Naming more than one grain returns a timesift_set() of those arrays.

Time zone

Bins follow the calendar the series is carried in, which is the tzone attribute of time; a column with none is read as UTC, and a name the zone database does not know is an error. The zone is resolved once, at the edge: below it the binning works in local time, where a day is 86400 seconds whatever the night did, so a zone that moves its clock at midnight has no midnight to lose. A bin start is a local time, so reporting it back as an instant needs a rule: one the clock skipped resolves to the instant the clock jumped to, one the clock repeated to the first of the two. Instants are read at whole seconds.

In a zone that keeps daylight saving time a day is therefore a wall-clock day: the day the clock goes forward spans 23 hours and the day it goes back 25, and the week and the month holding them one hour less or more. Hourly readings carried in "Europe/Vienna" give a 2021-03-28 day of 23 readings and a 2021-10-31 day of 25. A record kept on a fixed offset, as the Schrankogel deposit is (26,304 hourly readings, 1,096 days of 24), bins into days of 24 hours; a series carried in a daylight-saving zone bins that way too once its tzone is set to the fixed offset of standard time, such as "Etc/GMT-1" for Central European Time, whose days then run from 01:00 to 01:00 on the summer clock.

The "native" grain is the one grain not read on that clock: its bin is the reading itself, so the two readings of an hour a zone repeats are two bins, and the record read at "native" is the same array whichever zone it is carried in.

Bins that do not tile the record

Every unit must reach every bin, and consecutive bins must be one bin apart on the grain's own calendar. A bin no unit reaches is never built, so a month missing from the whole record would otherwise pass as four adjacent monthly bins with one simply gone. Neither the "native" grain, whose bin is the reading itself, nor a supplied calendar, which declares its own bin lengths, is held to the second rule. coverage() lays the same binning out as a count of readings per unit and bin, which is where a refused record's gaps are read off.

Partial bins

A bin is partial when the record does not cover its whole calendar span. Which bins those are follows from where the record starts and stops against the calendar, not from the grain alone: three years of hourly readings from 1 September carry no partial month and no partial season on a "09-01" boundary, but the same record carries a partial week at each end, because 1 September is a Wednesday. A record from an arbitrary deployment date carries one at each end of almost every grain.

A bin is partial if its start precedes the first reading of the record, or if its calendar span runs past the last reading plus the record's own sampling interval, taken as the smallest gap between consecutive distinct reading instants. Only a bin at an end of the record can satisfy either, because every unit is required to span every bin in between. The verdict is returned as the bin_partial attribute whichever way partial is set, so a kept partial bin is labelled rather than silent.

Keeping partial bins is the default because dropping them discards the record's ends: on a seasonal grain that is up to three months of readings at each end. The cost of keeping them is that such a bin's mean is taken over fewer readings and its extremes over fewer days, so cold_day and warm_day there are drawn from a shorter draw and sit closer to the bin mean than a full bin's would. bin_n gives the count the bin was actually reduced from.

A caller-supplied binning declares its own bins, so the package cannot know where the last one was meant to end and takes the record's end as its end. Such a final bin is never reported partial; its leading bin is judged as any other.

Examples

t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 40)
d <- data.frame(plot = rep(c("a", "b"), each = length(t)),
                t = rep(t, 2),
                temp = c(sin(seq_along(t) / 24), cos(seq_along(t) / 24)))
x <- grain_matrix(d, plot, t, temp, grain = "week",
                   stats = c("cold_day", "mean", "warm_day"))
dim(x)
dimnames(x)[[3]]

ladder <- grain_matrix(d, plot, t, temp, grain = c("day", "week"))
names(ladder)

Description

The set a learner without a ⁠data =⁠ of its own is fitted at every member of. grains() names calendar grains, lookbacks() names lookback spans, and either takes its arguments as separate names or as one vector.

Usage

grains(..., stats = "mean", year_start = "01-01")

lookbacks(..., lag = "0 days", bins = 1L, stats = "mean")

as_sift(x)

Arguments

Details

grains("auto") is the set of named grains the record gives at least two bins, in the order native, halfday, day, week, month, season, year, leaving out any grain the requested statistics are not defined at. There is no cap on it: over three years of hourly readings native is 26304 bins, and a caller who does not want that names the grains instead.

Value

A timesift_sift: a named list of representations, labelled by each one's own label where the list was not named.

Examples

grains("day", "week", "month")
grains(c("month", "year"), stats = c("cold_day", "mean", "warm_day"))
lookbacks("30 days", "90 days", bins = 3L)

Description

One model per response: a logistic regression on every bin-by-channel column of the representation, standardised so the prior on a coefficient means the same thing for each, optionally with an intercept for each unit and a Gaussian-process field over the targets' coordinates. A unit then carries what its own record says and what its neighbours' presences say, under the folds every other learner is scored on.

Usage

hierarchical(
  data = NULL,
  spatial = c("none", "nngp", "hsgp"),
  random = FALSE,
  cov = c("exponential", "matern32", "matern52", "gaussian"),
  neighbours = 15L,
  m = 6L,
  boundary = 1.5,
  nodes = 5L,
  threads = 1L
)

Arguments

Details

The coefficients, the intercept included, have a N(0, 2.5^2) prior. A unit's intercept is N(0, sd^2) and the field has a marginal standard deviation and a range, each standard deviation under a penalised-complexity prior (Simpson et al. 2017) with P(sd > 3) = 0.01, and the range under one anchored at a fifth of the coordinates' extent with P(range < anchor) = 0.5 (Fuglstad et al. 2019). Coordinates are centred and divided by one factor, so distances keep their proportions.

The field is a Gaussian process on the two columns named by coords in timesift(). "hsgp" is the Hilbert-space approximation (Solin and Sarkka 2020) with m Laplacian eigenfunctions per axis, and "nngp" the nearest-neighbour process (Datta et al. 2016) over the distinct locations, each conditioned on its neighbours nearest among those before it in lexicographic order of the coordinates, with the covariance cov, whose range is the field's. Its sparse precision is factored once per set of hyperparameters, so the fit scales with the number of locations rather than their square.

Inference is Laplace's method over the coefficients, the intercepts and the field together, and the hyperparameters are integrated over on a grid of nodes points per hyperparameter, centred on the mode of their posterior and weighted by it (Rue, Martino and Chopin 2009). With an intercept alone its standard deviation is set at the mode of its posterior; with neither the fit is the posterior mode of the coefficients, so at spatial = "none" and random = FALSE the learner is a penalised logistic model, every coefficient shrunk towards zero by its prior. A prediction at new units interpolates the field to their coordinates, so predict.timesift() is given targets that carry the same coordinate columns. The head's case weights enter the likelihood in every configuration. A response holding one value is predicted its mean and named in unfitted. The learner fits a presence-absence head.

With random = TRUE each unit, as named by id in timesift(), gets an intercept, which absorbs what its several targets share beyond the record: it is identified where a unit carries more than one target, as an anchored fit's units do. A prediction adds a unit's intercept where the unit was in the fit and leaves it at zero, the population level, where it was not, so a unit held out whole is predicted from its record and its place alone.

Both languages call one C++ core, src/ts_hierarchical.cpp.

Value

A learner().

Examples

hierarchical()
hierarchical(data = grain("month"), spatial = "hsgp")

Description

tss_inflation() maps a population skill to the level a design reports for it. This inverts that map: given a level actually read off a ladder, it solves for the population skill whose expected reported level equals it.

Usage

implied_skill(
  y,
  folds,
  observed,
  grid = seq(0, 0.95, by = 0.05),
  replicates = 200L,
  seed = 1L
)

Arguments

Details

It answers the question a level raises once the inflation is known, under the distribution of predictions tss_inflation() plants; a model whose predictions are distributed otherwise is inflated by a different amount. It does not correct a difference between two arms, whose inflations need not be equal.

Value

A data frame of one row per observed level: the level, the population skill it is consistent with, and whether that sits inside the grid the map was measured on.

Examples

set.seed(1)
y <- matrix(rbinom(1200, 1, 0.15), nrow = 200,
            dimnames = list(sprintf("p%03d", 1:200), paste0("sp", 1:6)))
implied_skill(y, fold_map(y, v = 5), observed = 0.71, replicates = 40)

Description

A chance-corrected agreement rate on a two-by-two table. Read against the observed response it is a skill score beside tss(); read between two models' decisions on the same units it says where the two part company.

Usage

kappa_score(y, p, rule = c("prevalence", "youden", "kappa", "mpa"))

decision_threshold(y, ...)

## Default S3 method:
decision_threshold(
  y,
  p,
  rule = c("youden", "kappa", "prevalence", "mpa"),
  perc = 0.9,
  ...
)

## S3 method for class 'timesift'
decision_threshold(
  y,
  candidate = "ensemble",
  rule = c("youden", "kappa", "prevalence", "mpa"),
  perc = 0.9,
  ...
)

model_agreement(y, p_a, p_b, rule = c("youden", "kappa", "prevalence", "mpa"))

Arguments

Details

Kappa is read at a threshold rather than maximised over one, so the rule that picks the threshold is part of the statistic. "youden" is the operating point tss() is defined at and inherits its selection bias; "kappa" maximises kappa itself and inherits the analogous bias; "prevalence" cuts at the observed presence rate, which selects nothing from the labels and is the rule to read an absolute level at. "mpa" is the minimum predicted area rule: the highest cut that still predicts presence at a share perc of the observed presences, so the labels enter through the presences alone.

Given a timesift() fit, decision_threshold() learns one cut per response from a candidate's out-of-fold predictions of the fit's own targets, which is the cut predict(type = "binary") applies to new targets.

Value

For kappa_score(), one number. For decision_threshold(), the cut itself, applied as p >= threshold, and on a fit one cut per response, named by it, NA where the response's held-out predictions give none. For model_agreement(), a one-row data frame carrying the agreement kappa between two models cut by the same rule, the share of units they decide differently, and how often each is the one that is right there.

Examples

y <- c(0, 0, 0, 1, 1, 1)
kappa_score(y, c(0.1, 0.2, 0.6, 0.4, 0.8, 0.9))
model_agreement(y, c(0.1, 0.2, 0.6, 0.4, 0.8, 0.9), c(0.2, 0.1, 0.3, 0.7, 0.9, 0.8))

Description

A learner is a pair of functions: one that fits a model to a representation and a response, and one that predicts from it. Everything the package fits goes through this pair, so a learner of your own sits beside the ones that ship and needs no change to the ladder, the folds or the scoring.

Usage

learner(
  name,
  fit,
  predict,
  data = NULL,
  reads = c("tabular", "sequence"),
  multi = c("separate", "joint"),
  control = NULL,
  needs = character(),
  params = list()
)

Arguments

Details

A learner declares what it can be handed and how it covers several responses. reads is "tabular" where the bins reach it as a block of predictors and "sequence" where their order in time is what it reads, and multi is "joint" where one fitted model covers every response and "separate" where the learner fits one per response. Either way it is handed the whole response matrix and returns one column per response, so nothing above the learner layer has to know which it was, and the block of predictors is built once for the fit rather than once for every response of it.

data pins a learner to one representation. Left NULL the learner runs across every representation of the run.

Value

A timesift_learner.

Examples

# Every bin of every channel as a predictor, one weighted logistic regression per variable.
flat_glm <- learner(
  "flat_glm",
  fit = function(x, y, weights, ...) {
    f <- as.data.frame(as.matrix(x))
    lapply(seq_len(ncol(y)), function(j) {
      stats::glm(y[, j] ~ ., data = f, weights = weights[, j],
                 family = stats::quasibinomial())
    })
  },
  predict = function(model, x) {
    f <- as.data.frame(as.matrix(x))
    vapply(model, function(m) stats::predict(m, f, type = "response"), numeric(nrow(f)))
  }
)
flat_glm

Description

One generalised linear model per variable over every bin-by-channel column, its terms chosen by Akaike's criterion. The family is the response head's: logistic under a binary cross-entropy loss, Gaussian under a squared-error one and Poisson under a Poisson-deviance one. So are the case weights, so a rare response weighs here what it weighs in every other learner.

Usage

linear(
  data = NULL,
  select = c("both", "forward", "backward", "none"),
  terms = c("power", "column"),
  max_terms = Inf,
  degree = 2L,
  threads = 1L
)

Arguments

Details

The defaults are biomod2's GLM: every column enters as x + I(x^2), and the terms are searched in both directions by AIC as MASS::stepAIC() searches them, with no bound on how many are kept. Under the shipped presence-absence head those weights are on, so a default linear() is biomod2's GLM specification fitted under them; a head registered without weights fits it unweighted.

terms says what one term is. Under "power" each power of a column, x, x^2 and so on, is a term of its own, which is how biomod2 writes a quadratic formula and how stepAIC() walks it. Under "column" a term is a column's orthogonal polynomial of degree degree, so a column enters with its curvature at once and can be non-monotone in the reading the way a niche optimum is. A column holding one value over the fitting units is not a term.

select is the search. "both" starts from the intercept and at every step takes the move that lowers the criterion most, adding a term or dropping one it holds, while one does and the model holds fewer than max_terms. "forward" only adds, and "backward" starts from every term and only drops. "none" fits every term and selects nothing. The two-way and backward searches are stepAIC(), step for step: the model as it stands wins a tie, a term whose removal leaves the rank unchanged is dropped first, and an addition that does not raise the rank is not offered. max_terms bounds what a forward or two-way search adds, and the backward and unselected fits start from every term whatever it is.

Each fit is R's glm.fit: iteratively reweighted least squares, the rank read off the same pivoted decomposition, and the same stopping rule. A move whose fit does not settle within its 25 iterations is refused rather than taken, and the fit names every response whose final model did not settle in stopped.

A model with nothing but the intercept predicts the response's share among the fitting units. Selection happens inside whichever units the learner is handed, so under grain_ladder() it is redone in every fold. Reported beside a penalised fit it also prices discrete selection: choosing a handful of columns out of hundreds is high variance, and that variance is a cost of the selector rather than of the features. An unselected fit over hundreds of columns separates any response it is given and is meant for a coarse grain, data = grain("season").

The search and the fits run on the core the Python package calls, so the two select the same terms and return the same coefficients.

Value

A learner().

Examples

linear()
linear(select = "forward", terms = "column", max_terms = 3)
linear(data = grain("season"), select = "none")

Description

Reads, for every target, a fixed length of record ending a fixed lag before that target's own instant, and summarises it by one or more statistics. It is the reduction a calendar cannot express: two targets on the same unit a fortnight apart read two different stretches of the same series, so the bins are relative to the target rather than to a month or a week.

Usage

lookback_matrix(
  data,
  id,
  time,
  value,
  at,
  span,
  lag = "0 days",
  bins = 1L,
  stats = "mean"
)

Arguments

Details

Bin b of a target anchored at a covers ⁠[a - lag - span + b * step, a - lag - span + (b + 1) * step)⁠, with step the lookback's length divided by bins and b counted from zero. The interval is closed at the left and open at the right, so a reading on a boundary belongs to the later bin, and only the readings of the target's own unit are read.

Every ⁠(target, bin)⁠ cell must hold at least one reading. A bin lying wholly past either end of the record is an error naming the target and the interval, never a padded row: an invented value in front of a model is worse than a target the record cannot answer for. A bin the record reaches only in part is summarised over the readings it holds, and bin_n records how many.

The four day-level statistics reduce each calendar day first, so they are defined only where every day lies whole inside one bin. For a lookback that is two conditions rather than one: step must be a whole number of days, and a - lag - span must fall on a day boundary. Either failing is an error naming the target.

Value

A numeric array of shape ⁠[target, bin, channel]⁠, of class timesift_matrix. Its rows are the rows of at, in at's own order, named by at's row names. Its bins are named by where each one opens relative to the anchor, oldest first, and its channels by the statistic. Attributes:

• grain: "lookback".

• span, lag: the durations, resolved to seconds.

• bins: how many bins the lookback was cut into.

• stats: the statistic names in channel order.

• bin_n: a ⁠[target, bin]⁠ matrix of how many readings fell in each cell.

Durations

span and lag are read from a count and a unit – "30 days", "12 hours", "1 year" – or from a bare number of seconds. A year is 365 days and a month is 30 days here. A lookback of a fixed length is a fixed length, not a calendar step: the point of anchoring on the target is that every target reads the same amount of record, which a February and a leap year would take away. Where the calendar is what matters, grain_matrix() is the call that follows it.

The units are seconds, minutes, hours, days, weeks, months and years, singular or plural.

Time zone

The calendar is the series', taken from the tzone attribute of time as grain_matrix() takes it; a column with none is read as UTC. The anchors are instants and are read as a clock in that same calendar, whatever zone at carries, so one record is binned by one calendar. The span is measured on that clock: a lookback of one day ending at a local midnight holds the whole local day before it, which is 25 hours of record on the night a zone sets its clock back and 23 on the night it sets it forward. That is what keeps a calendar day whole inside a bin for the four day-level statistics; a length fixed in instants could not.

See Also

grain_matrix(), the reduction that follows the calendar instead.

Examples

t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 60)
d <- data.frame(plot = rep(c("a", "b"), each = length(t)),
                t = rep(t, 2),
                temp = c(sin(seq_along(t) / 24), cos(seq_along(t) / 24)))
at <- data.frame(id = c("a", "b"),
                 at = as.POSIXct(c("2021-10-20", "2021-10-25"), tz = "UTC"))
x <- lookback_matrix(d, plot, t, temp, at = at, span = "30 days", bins = 3L,
                   stats = c("cold_day", "mean", "warm_day"))
dim(x)
dimnames(x)[[2]]

Description

One MARS model per response, over every bin-by-channel column of the representation, fitted as the earth package fits it and as biomod2 fits MARS. The forward pass starts from the intercept and at each step multiplies a term already in the model by a pair of hinges on one column, max(0, x - t) and max(0, t - x), taking the parent, the column and the knot t that most reduce the residual sum of squares of a least-squares fit to the response. A knot at a column's least value enters the column linearly. degree bounds how many hinges a term multiplies, so degree = 1 is an additive model and 2 admits pairwise interactions. The pass stops at max_terms terms, when a step raises the R-squared by less than min_gain, or when no term reduces the residuals.

Usage

mars(
  data = NULL,
  degree = 1L,
  penalty = NULL,
  max_terms = NULL,
  min_gain = 0.001,
  minspan = 0L,
  endspan = 0L,
  fast_k = 20L,
  fast_beta = 1,
  prune = TRUE,
  nprune = NULL,
  threads = 1L
)

Arguments

Details

The pruning pass then removes terms one at a time, each time the one whose loss raises the residuals least, and keeps the subset of least generalised cross-validation, which charges penalty for each knot. Under a presence-absence head the kept terms are refitted as a logistic model, as earth's glm = list(family = binomial) refits them and biomod2 asks it to, and the prediction is that model's probability; under a count head they are refitted as a Poisson model with a log link, as earth's glm = list(family = poisson) refits them, and the prediction is its mean count; under a squared-error head they are refitted by least squares.

The defaults are earth's own, which biomod2 uses under its default option set and under "bigboss" alike: degree one, penalty = 2, min_gain = 0.001, ⁠max_terms = min(200, max(20, 2 p)) + 1⁠ for p columns, Friedman's rules for the spans between knots, and Fast MARS over the 20 best parents. The forward pass, the pruning pass and the refit live in the core the Python package calls, pinned against earth in the fixtures, so the two languages keep the same terms and return the same coefficients.

The case weights are the response head's, positive_weights() under presence-absence, and weigh both passes and the refit as earth weighs them. earth refits the whole basis by QR at every candidate knot of a weighted fit; the core reaches the same residual sums with Friedman's running updates, which is what makes a weighted fit over hundreds of columns affordable. Under the shipped presence-absence head those weights are on, so a default mars() is earth's specification fitted under them; a head registered without weights fits it unweighted.

A response holding one value is predicted its mean.

Value

A learner().

Examples

mars()
mars(degree = 2L)
mars(data = grain("season"), nprune = 10L)

Description

One maximum-entropy model per response, over every bin-by-channel column of the representation, in the formulation of the maxnet package (Phillips et al. 2017) and biomod2's MAXNET: maxnet's feature classes, its regularisation of each feature, and a lasso over them, fitted by the penalised core elasticnet() runs on, which the Python package calls too. With the maxnet package's own settings the features and the penalty factors are maxnet's to rounding, and the fit settles at the objective glmnet reaches for maxnet.

Usage

maxent(
  data = NULL,
  classes = NULL,
  regmult = 1,
  formulation = c("background", "absence"),
  type = NULL,
  knots = 50L,
  add_samples = TRUE,
  clamp = TRUE,
  n_inner = 5L,
  s = c("lambda.min", "lambda.1se"),
  tol = 1e-08,
  max_design = 2,
  threads = 1L,
  seed = 1L
)

Arguments

Details

The feature classes are the letters of classes: l the column itself, q its square, p the product of each pair of columns, h forward and reverse hinges at the interior of knots equally spaced points of each column's range, and t thresholds at 49 interior points of it. Left NULL, they follow the response's presence count as maxnet.formula() has them: "l" under 10 presences, "lq" under 15, "lqh" under 80, and "lqph" from 80 on. A column holding one value over the units fitted takes no feature.

formulation says what the absences are. "background" is maxnet's own and what biomod2 fits as MAXNET: every unit is background, each presence joins the background again unless an absence carries the same readings (add_samples), the background is weighted 100 against a presence's 1, and the model is read at the last of maxnet's 200 penalties, which scale with regmult. Its output is maxnet's type, "cloglog" by default as biomod2 predicts it. "absence" reads the absences as absences: a logistic lasso over the same features and penalty factors under the response head's case weights, positive_weights() under presence-absence, with the penalty chosen by an inner cross-validation dealt as the elastic net's is, and a probability as output.

The background formulation takes no case weights, as maxnet takes none and biomod2 passes none: the background weight is what sets a presence's weight there. Either formulation holds each column inside the range it was fitted on, and each feature inside its own, before predicting, as maxnet's predict(clamp = TRUE) does; clamp = FALSE reads them as they are.

A hinge per column per knot makes the design large: a weekly three-channel representation, 471 columns, is 47,100 features under "lqh", and its products under "lqph" 110,685 more. The design is held in memory with a centred copy beside it, and a fit whose design would take more than max_design gigabytes is refused with the size it would have taken. threads fit that many responses at once, each the fit it gets alone, and the designs held at once share max_design. A coarser representation (data = grain("month")), fewer classes or fewer threads is then the way to fit it.

A response with fewer than two presences, or one whose inner training sets cannot each hold two of each outcome under the absence formulation, is predicted its share among the fitting units, and the fit names it in unfitted. A path that does not settle at a penalty ends there, as the elastic net's does, and is read at its last settled point; the fit names every such response in stopped.

The learner needs a presence-absence response: maxnet has no model for a continuous one, and a head whose loss is not the binary cross-entropy is refused.

Value

A learner().

Examples

maxent()
maxent(classes = "lqh", regmult = 2)
maxent(formulation = "absence")

Description

A representation is the way a record is reduced before a learner reads it. It is one of four things: the record as recorded, the record at a calendar grain, several grains side by side in one block of features, or a lookback of fixed length that ends at each target's own instant. A representation holds only its settings. It can therefore be described before any record is seen, name the arm it produced in a fitted object, and be rebuilt for new targets in predict.timesift().

Usage

native(stats = "mean", year_start = "01-01")

grain(grain, stats = "mean", year_start = "01-01")

multigrain(grains = NULL, stats = "mean", year_start = "01-01")

lookback(span, lag = "0 days", bins = 1L, stats = "mean")

Arguments

Details

Every representation carries label, the name it is reported under; kind, one of "grain", "multigrain" and "lookback"; the settings its kind uses; and sequence, which says whether its bins are ordered in time and so mean something to a convolution. native(), grain() and a lookback() of more than one bin are sequences; multigrain() and a one-bin lookback() are blocks of features.

multigrain() flattens each of its grains to one row per target and puts them side by side, so a column of the block names the grain, the statistic and the bin it came from. It is the tabular representation a penalised regression or a random forest reads.

Value

A timesift_representation.

Examples

native()
grain("week", stats = c("cold_day", "mean", "warm_day"))
multigrain(c("month", "season"))
lookback("30 days", bins = 3L)

Description

Holds one bin of the record back at a time, rescores the held-out units, and records the fall in score as that bin's weight. Nothing is refitted: the models kept by grain_ladder(keep_fits = TRUE) are the ones read, so the profile describes the models that produced the reported scores rather than a fresh set of them.

Usage

occlusion(x, ...)

## Default S3 method:
occlusion(x, ...)

## S3 method for class 'timesift_ladder'
occlusion(
  x,
  data,
  y,
  arm,
  over = c("bin", "channel"),
  substitute = c("permute", "fold_mean", "unit_mean"),
  metric = NULL,
  permutations = 20L,
  seed = 1L,
  ...
)

## S3 method for class 'timesift'
occlusion(x, candidate, over = c("bin", "channel"), ...)

Arguments

Details

A model has to be shown something in place of a held-back bin, and what it is shown decides what the weight means. Permuting the bin's values across units keeps the observed readings exactly and cuts only the link between a reading and its unit. Replacing every unit by the fitting-fold mean removes all between-unit variation while keeping the shape of the year. Replacing the bin by each unit's own mean over the record keeps how warm a unit is and removes only that bin's departure from it.

A bin is held back whole: a permutation moves every channel of the bin together, so a unit is shown another unit's week rather than a coldest day from one unit beside a warmest day from another. A channel that is the same for every unit, as calendar_channels() are, says where the bin sits rather than what a unit read there, and is left in place.

Read with over = "channel" the same machinery asks what each statistic of a grain carries, holding one channel back across the whole record instead of one bin across all channels.

A column named in static is a channel of a representation built from a record, constant across its bins, so over = "channel" holds it back; where the run has no record the columns are the bins of one block and over = "bin" does.

Reached through a timesift() run rather than through a ladder, the profile reads the per-fold models the run was told to keep, so every bin is held back from a model that never saw the units it is rescored on. The candidate is named as summary() reports it.

Value

A data frame of one row per held-back part and variable, carrying the mean weight over folds, the score with and without the part, and importance: one minus the correlation between the predictions with and without the part, which is how biomod2 scores a predictor and needs no response. Both are read on the cells the fit was scored on, and importance is NA where either set of predictions is constant.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 120)
units <- sprintf("p%02d", 1:40)
warmth <- rnorm(40)
d <- data.frame(
  plot = rep(units, each = length(t)), t = rep(t, length(units)),
  temp = as.numeric(vapply(warmth, function(w) w + sin(seq_along(t) / 300) + rnorm(length(t)),
                           numeric(length(t)))))
y <- matrix(rbinom(80, 1, plogis(c(warmth, -warmth))), nrow = 40,
            dimnames = list(units, c("sp1", "sp2")))
x <- grain_matrix(d, plot, t, temp, grain = "month")
lad <- grain_ladder(x, y, elasticnet(), folds = fold_map(y, v = 3),
                     keep_fits = TRUE, verbose = FALSE)
head(occlusion(lad, x, y, "month|elasticnet", permutations = 3))

Description

The statistics biomod2 reads a model of an ordinal response by. The response is a column of whole-number classes, the model predicts a number on the same scale, and each prediction is read as the observed class nearest to it, the lower class on a tie. With m[i, j] the units of observed class j predicted as class i, over the k classes observed in the cell:

Usage

ordinal_metric(y, p, metric = names(.ordinal_metrics))

Arguments

Details

A class with no unit, or in which nothing is predicted, adds zero to its mean. The four are registered as ordinal_accuracy, ordinal_recall, ordinal_precision and ordinal_f1.

Value

One number, or NA where the cell defines none.

Examples

ordinal_metric(c(1, 1, 2, 3, 3), c(1.2, 2.4, 2, 2.8, 3.4), "accuracy")

Description

Two arms scored on the same held-out units do not necessarily have the same set of defined cells, so a difference of two marginal means is not a difference between the arms. This takes the difference inside each ⁠(variable, fold)⁠ cell both arms scored, averages it within a variable over its folds, and summarises those per-variable means, the variables being the independent replicates.

Usage

paired_contrast(ladder, a, b, interval = c("variables", "nested_cv"))

Arguments

Details

Pairing removes the variation between variables and the part of a threshold-selected metric's bias that the design sets, but not the part that belongs to each arm. TSS read at the threshold that maximises it is biased upward where presences are thin, and how far depends on how an arm's predictions are distributed as well as on how many presences the cell holds, so two arms of equal skill on the same cell can carry different biases. A difference in TSS can therefore favour one arm with no difference in skill behind it. A threshold-free metric such as roc_auc() has no cut to choose, and a TSS contrast is best read beside the same contrast under it.

An arm is named whole, by its grain and its learner. A learner named alone would have to take its best grain, and that grain is chosen on the held-out scores the contrast is then read off: the difference becomes one between two maxima, favouring whichever learner ran across more grains, and neither the interval nor the p-value accounts for the choice. It is the mechanism tss_inflation() measures one level down, and here pairing does not cancel it. select_grain() chooses a grain on inner folds instead, and its compare argument contrasts the selection with the arms of a ladder on matched cells.

Value

A one-row data frame: the mean per-variable difference; the centre of the interval and its bounds, on Student's t across variables or on nested cross-validation as interval names, which the row carries; the number of variables the difference favours; the paired cells and variables it rests on; a Wilcoxon signed-rank p-value; and p_method, "exact" where the p-value is read off the exact distribution and "normal" where it is the normal approximation with continuity and tie corrections, which it is when the per-variable differences hold a zero or a tie or number fifty or more.

References

Bates, S., Hastie, T. and Tibshirani, R. (2024). Cross-validation: what does it estimate and how well does it do it? Journal of the American Statistical Association 119(546), 1434-1445. doi:10.1080/01621459.2023.2197686

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 200)
units <- sprintf("p%02d", 1:60)
warmth <- rnorm(60)
d <- data.frame(
  plot = rep(units, each = length(t)), t = rep(t, length(units)),
  temp = as.numeric(vapply(warmth, function(w) w + sin(seq_along(t) / 300) + rnorm(length(t)),
                           numeric(length(t)))))
y <- matrix(rbinom(120, 1, plogis(c(warmth, -warmth))), nrow = 60,
            dimnames = list(units, c("sp1", "sp2")))
x <- grain_matrix(d, plot, t, temp, grain = c("week", "month"))
lad <- grain_ladder(x, y, elasticnet(), folds = fold_map(y, v = 3), verbose = FALSE)
paired_contrast(lad, "week|elasticnet", "month|elasticnet")

Description

One network per response, over every bin-by-channel column of the representation, fitted as the nnet package fits it and as biomod2 fits ANN. Each of hidden logistic units takes a bias and every column, and the output takes a bias, every hidden unit and, with skip, every column again. The weights start uniform on ⁠[-range, range]⁠ and are fitted by the variable metric (BFGS) method of Nash (1990), the minimiser nnet uses, on the response head's loss plus decay times the sum of the squared weights, biases included. The fit stops after max_iter iterations, when the objective falls below abs_tol, or when an iteration lowers it by no more than rel_tol of itself.

Usage

perceptron(
  data = NULL,
  hidden = NULL,
  decay = NULL,
  range = NULL,
  max_iter = NULL,
  skip = FALSE,
  standardise = FALSE,
  abs_tol = 1e-04,
  rel_tol = 1e-08,
  preset = c("default", "bigboss"),
  max_hessian = 2,
  threads = 1L,
  seed = 1L
)

Arguments

Details

Under a presence-absence head the output is the logistic function of its sum and the loss the cross-entropy, nnet's entropy = TRUE; under a head with a squared-error loss the output is the sum itself and the loss the sum of squares, nnet's linout = TRUE; under a count head the output is the exponential of the sum and the loss the Poisson deviance, which nnet does not offer. biomod2 leaves nnet's own entropy = FALSE, fitting a presence-absence response by least squares on the logistic output; the learner fits the head's loss, as every learner does.

preset says whose defaults the settings left NULL take. "default" is what biomod2's default option set fits: two hidden units, as biomod2 sets them, and nnet's own decay = 0, range = 0.7 and max_iter = 100. "bigboss" is biomod2's tuned option set: five hidden units, decay = 0.1, range = 0.1 and max_iter = 200. A setting given explicitly beats either.

nnet reads the columns as given, and so does the default: a record in its own units saturates the hidden units sooner the wider its range. With standardise = TRUE each column is centred on its mean and divided by its sample standard deviation over the fitting units, and a prediction centres and scales by the fit's own.

The network, its objective and the minimiser live in the core the Python package calls, pinned against nnet in the fixtures from the same starting weights, so the two languages fit the same network. The starting weights are drawn from the core's own generator, so a fit does not repeat nnet's from the same R seed. threads fit that many responses at once, each network the same as when fitted alone. The minimiser holds an approximate inverse Hessian of one number per pair of weights for every network fitted at once; a fit that would need more than max_hessian gigabytes for them is refused with the size, and a coarser grain, fewer hidden units or fewer threads shrinks it.

The case weights are the response head's, positive_weights() under presence-absence, and weigh each unit's term of the loss as nnet's weights do. Under the shipped presence-absence head those weights are on, so a default perceptron() is biomod2's ANN specification fitted under them; a head registered without weights fits it unweighted.

A response holding one value is predicted its mean.

Value

A learner().

Examples

perceptron()
perceptron(preset = "bigboss")
perceptron(hidden = 4L, decay = 0.01)

Description

One line per learner across the representations it ran on, read the way a ladder is read, and the stack's held-out score drawn across them, its weights fitted inside each outer training fold. The curves are scored on the folds a choice among them would be judged on, so the best of them sits a little high; the ensemble line does not. Where the ensemble line sits above every curve the candidates are carrying different parts of the signal, and where it sits on or below the best curve they are not.

Usage

## S3 method for class 'timesift'
plot(x, col = NULL, interval = TRUE, ...)

Arguments

Value

The table the plot is drawn from, invisibly.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 90)
units <- sprintf("p%02d", 1:30)
warmth <- rnorm(30)
logger <- data.frame(
  plot = rep(units, each = length(t)), datetime = rep(t, 30),
  temp = as.numeric(vapply(warmth, function(w) w + sin(seq_along(t) / 300), numeric(length(t)))))
plots <- data.frame(plot = units,
                    sp_a = rbinom(30, 1, plogis(2 * warmth)),
                    sp_b = rbinom(30, 1, plogis(-2 * warmth)))
fit <- timesift(plots, logger, y = starts_with("sp_"), id = plot, time = datetime,
                sift = grains("week", "month"), resampling = cv(v = 3L), verbose = FALSE)
plot(fit)

Description

One line per learner across the grains, at the across-variable mean of the per-variable score, with a 95 percent interval from its standard error across variables, on Student's t with one degree of freedom fewer than there are variables. An open circle marks each learner's best grain, which is where the curve says the record stops paying for being read more finely.

Usage

## S3 method for class 'timesift_ladder'
plot(x, col = NULL, interval = TRUE, ...)

Arguments

Value

The summary table the plot is drawn from, invisibly.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 200)
units <- sprintf("p%02d", 1:60)
warmth <- rnorm(60)
d <- data.frame(
  plot = rep(units, each = length(t)), t = rep(t, length(units)),
  temp = as.numeric(vapply(warmth, function(w) w + sin(seq_along(t) / 300) + rnorm(length(t)),
                           numeric(length(t)))))
y <- matrix(rbinom(120, 1, plogis(c(warmth, -warmth))), nrow = 60,
            dimnames = list(units, c("sp1", "sp2")))
x <- grain_matrix(d, plot, t, temp, grain = c("day", "week", "month"))
lad <- grain_ladder(x, y, elasticnet(), folds = fold_map(y, v = 3), verbose = FALSE)
plot(lad)

Description

The map of cells lost, kept, absent and gained, for one response.

Usage

## S3 method for class 'timesift_range_change'
plot(x, variable = NULL, ...)

Arguments

Value

x, invisibly.

Examples

now <- terra::rast(nrows = 10, ncols = 10, vals = rep(c(1, 0), each = 50))
names(now) <- "sp1"
later <- terra::rast(now)
terra::values(later) <- rep(c(0, 1, 1, 0), each = 25)
plot(range_change(now, later))

Description

One line per response against the value of the predictor. For two predictors, the prediction of one response as a filled image over the grid of the two.

Usage

## S3 method for class 'timesift_response_curve'
plot(x, variable = NULL, col = NULL, legend = "topright", ...)

Arguments

Value

The data frame drawn, invisibly.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 90)
units <- sprintf("p%02d", 1:30)
warmth <- rnorm(30)
logger <- data.frame(
  plot = rep(units, each = length(t)), datetime = rep(t, 30),
  temp = as.numeric(vapply(warmth, function(w) w + sin(seq_along(t) / 300), numeric(length(t)))))
plots <- data.frame(plot = units,
                    sp_a = rbinom(30, 1, plogis(2 * warmth)),
                    sp_b = rbinom(30, 1, plogis(-2 * warmth)))
fit <- timesift(plots, logger, y = starts_with("sp_"), id = plot, time = datetime,
                sift = grains("month"), resampling = cv(v = 3L), ensemble = FALSE,
                n_inner = NULL, verbose = FALSE)
plot(response_curve(fit, "elasticnet / month", "mean", n = 20))

Description

Every candidate's inner score in every outer fold, one line per outer fold, with an open circle on the candidate that fold selected. A selection that lands on the same candidate each time draws its circles in one column; one that wanders says the grid is flat enough that the choice is arbitrary, which is worth seeing beside the estimate rather than after it.

Usage

## S3 method for class 'timesift_selection'
plot(x, col = NULL, ...)

Arguments

Value

The table of inner scores the plot is drawn from, invisibly, with at, the position of each candidate on the axis, and selected, whether its fold chose it.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 120)
units <- sprintf("p%02d", 1:40)
warmth <- rnorm(40)
d <- data.frame(
  plot = rep(units, each = length(t)), t = rep(t, length(units)),
  temp = as.numeric(vapply(warmth, function(w) w + sin(seq_along(t) / 300) + rnorm(length(t)),
                           numeric(length(t)))))
y <- matrix(rbinom(80, 1, plogis(c(warmth, -warmth))), nrow = 40,
            dimnames = list(units, c("sp1", "sp2")))
x <- grain_matrix(d, plot, t, temp, grain = c("week", "month"))
sel <- select_grain(x, y, elasticnet(), folds = fold_map(y, v = 3), n_inner = 3,
                    verbose = FALSE)
plot(sel)

Description

The weight every learner that ships fits a presence-absence response under: each presence of a response weighs the ratio of absences to presences among the fitting units, capped, and each absence weighs one. A response with a presence in one target of a hundred is otherwise fitted away by any learner that minimises a mean loss, and the encoders, the penalised fit, the forest and the forward search would each have to decide that for themselves.

Usage

positive_weights(y, cap = 50, fitting = NULL)

Arguments

Details

The ratio is read off the units the model is fitted on and the weight applies to every unit handed in. An encoder holds part of its units back as an inner validation set and reads the loss it stops on under the same weights, so the loss that stops the fit is the loss the fit minimises; a unit held back that way is not in the count the ratio is made from.

The weights are the response head's: the shipped presence-absence head carries this function as its weights, and a head registered with weights = function(y, fitting) positive_weights(y, cap = 20, fitting = fitting) weights every learner by that cap instead. A head without weights is fitted unweighted.

Value

A numeric ⁠[unit, variable]⁠ matrix of case weights, one per cell of y.

Examples

y <- cbind(rare = c(1, 0, 0, 0, 0, 0), common = c(1, 1, 1, 0, 0, 0))
positive_weights(y)
positive_weights(y, fitting = c(TRUE, TRUE, TRUE, TRUE, FALSE, FALSE))

Description

Rebuilds every member's representation for the new targets from the settings its own arm was built with, predicts with the model refitted on all targets, and combines them where the ensemble is asked for.

Usage

## S3 method for class 'timesift'
predict(
  object,
  targets,
  series = NULL,
  candidate = "ensemble",
  type = c("response", "binary", "spread"),
  rule = c("youden", "kappa", "prevalence", "mpa"),
  alpha = 0.05,
  ...
)

Arguments

Details

With type = "binary" each response is cut into presence and absence at the threshold decision_threshold() learns from the same candidate's out-of-fold predictions of the fit's own targets: never from predictions of units the model was fitted on, and never from the new targets, which carry no response. rule is the rule that picks it. A binary map of a species is this prediction on one target per map cell.

With type = "spread" the ensemble's members are read side by side rather than combined: ensemble_spread() gives their weighted mean, standard deviation, coefficient of variation and interval at alpha, which is biomod2's EMcv and EMci and, on one target per map cell, an uncertainty map.

Value

A ⁠[target, response]⁠ matrix of predictions, named by target and in the order the fit carries its own targets: sorted by identifier, or the targets' own order where target_time anchors them. Under type = "binary" an integer matrix of 0 and 1, NA for a response whose held-out predictions give no cut. Under type = "spread" the ⁠[target, response, statistic]⁠ array of ensemble_spread().

See Also

decision_threshold() for the cuts themselves.

Examples

set.seed(1)
t <- seq(as.POSIXct("2021-09-01", tz = "UTC"), by = "hour", length.out = 24 * 90)
units <- sprintf("p%02d", 1:30)
warmth <- rnorm(30)
logger <- data.frame(
  plot = rep(units, each = length(t)), datetime = rep(t, 30),
  temp = as.numeric(vapply(warmth, function(w) w + sin(seq_along(t) / 300), numeric(length(t)))))
plots <- data.frame(plot = units,
                    sp_a = rbinom(30, 1, plogis(2 * warmth)),
                    sp_b = rbinom(30, 1, plogis(-2 * warmth)))
fit <- timesift(plots, logger, y = starts_with("sp_"), id = plot, time = datetime,
                sift = grains("week", "month"), resampling = cv(v = 3L), verbose = FALSE)
head(predict(fit, plots, logger), 3)
predict(fit, plots, logger, type = "binary", rule = "prevalence")[1:3, ]