Finding your way around geosmooth

Prediction and reusable fits

summary(fit) shows the selected settings and the diagnostics stored in the fit, and says so when no convergence information is available. fitted(fit) returns the fitted values as a vector or a matrix, as stored by the fit. residuals(fit) returns the response minus the fitted values, with NA where the response is missing. Residuals exist only at the training points, so a fit that was evaluated at other coordinates has none. They are on the scale of the response, also for binomial fits; they are not Pearson or deviance residuals. Probability mass is a different quantity and is read from fit$rho. See help("model-values", package = "geosmooth") and help("model-summary", package = "geosmooth") for details.

Three further terms: a refit fits a new response while reusing the neighborhoods or the operator of an existing fit. It does not repeat the model selection of the original fit; the help page of each refit method says whether the penalties stay fixed or can be selected again. An operator is the matrix of differences or derivatives that a method penalizes; it is not a fitted response. A spectral basis is a stored set of graph eigenvectors, which allows many responses to be smoothed without recomputing them.

Refit and extract smoother matrices

refit(fit, y = new.y, ...) works for the six kinds of model in the table below. R chooses the method from the class of the fitted object. The arguments named in the last column belong to that method only.

A refit of a metric graph low-pass fit or of a quadratic Hessian fit returns a summary that cannot itself be refitted, so for further responses keep calling refit() on the original fit. The other four kinds return fitted objects that can be refitted again.

smoother.matrix(fit, ...) returns the matrix that maps responses to predictions, for LPS and MALPS fits. For LPS, the fit must have a single candidate configuration; the matrix maps the training responses to the predictions at the evaluation points, and check.tol sets the tolerance of an internal check. For MALPS, the matrix maps training responses to predictions at the training points, given the stored supports and weights; max.n limits the size of the matrix. For a robust MALPS fit, allow.robust = TRUE returns the matrix for the final robust weights held fixed. It does not reflect how those weights depend on the response.

malps.fit <- fit.malps(
  X, y, degree = 1L, support.type = "knn", support.size = 12L
)
new.y <- cos(2 * pi * X[, 1])
malps.refit <- refit(malps.fit, y = new.y)
S <- smoother.matrix(malps.fit)
max(abs(as.vector(S %*% new.y) - predict(malps.refit)))
#> [1] 4.440892e-16

If you are updating code written for version 0.1.0: the six refit.*() functions and the two *.smoother.matrix() functions are no longer exported. Call refit() and smoother.matrix() instead, and in refit calls replace fitted.model = with object = and y.new = with y =. The other arguments and the classes of the results are unchanged. help("refit", package = "geosmooth") and help("smoother.matrix", package = "geosmooth") link to the arguments of each method, and the old help names still work.

Shared model operations

Local regression and prediction synchronization

Model assessment and uncertainty

Generalized cross-validation (GCV) estimates prediction error from the effective degrees of freedom of a smoother, without fitting a separate model for each left-out observation. It is a different criterion from cross-validation with folds, and the two need not select the same settings.

The uncertainty functions for LPS and for MALPS rest on different assumptions. lps.pointwise.band() requires a Gaussian fit made with the R backend and the orthogonal polynomial basis, with a single candidate configuration and a fixed chart dimension, and evaluated at the training points. malps.gcv() and bootstrap.malps() treat the selected supports and weights as fixed and do not repeat the model selection. Neither gives simultaneous coverage, and neither accounts for the uncertainty of choosing the settings from the data.

Lifting trend filters

Both fitting functions accept an operator that was built beforehand, or build one themselves. For a single fit, call the fitting function directly. Build the operator yourself when you want to fit several responses, examine the neighborhoods, or avoid recomputing an expensive operator. Both fitting functions require the optional package genlasso, and both predict only at the training points.

Graph spectral smoothing and graph trend filtering

geosmooth does not construct graphs or compute shortest paths; the dgraphs package does. A graph is passed as an adjacency list, in which element i holds the indices (starting from 1) of the neighbors of vertex i, and a weight list of the same shape that holds a value for each edge.

The meaning of the edge values depends on the method. A length is a distance between two vertices. A conductance is the strength of their connection, and larger values mean a stronger connection. Metric low-pass smoothing expects lengths by default and converts them to conductances; pass weight.type = "conductance" if your values are already conductances. Graph trend filtering with the default recursive operator treats the values as conductances, whereas its path divided-difference operator treats them as lengths.

The heat time eta controls the amount of smoothing in graph low-pass filtering. A basis that contains only part of the spectrum may be unable to represent the filter for small positive heat times. If you build or extend a grid of heat times yourself, consult the spectral diagnostics returned with the basis.

With a complete spectrum, heat smoothing also accepts missing responses, as long as every connected component of the graph has at least one observed response. The observed-data heat guide describes the criterion that is minimized in this case and the optional numerical cross-checks requested with numerical.policy = "verified.fallback". Under that policy, a heat time is not selected automatically when the GCV scores cannot distinguish the candidates, either because the scores are constant by construction or because their differences are smaller than the numerical accuracy. The resolution bands reported there describe numerical accuracy; they are not confidence intervals.

# Each row, vertex, and response refers to the same location.
locations <- matrix(seq(0, 1, length.out = 20), ncol = 1,
                    dimnames = list(sprintf("site-%02d", 1:20), "position"))
response <- sin(2 * pi * locations[, 1])
graph <- dgraphs::create.rknn.graph(locations, radius = .06,
                                    graph.detail = "minimal")
# dgraphs' new dgraph representation has public accessors. The second branch
# supports the pinned transitional dependency used by this development release.
if (inherits(graph, "dgraph")) {
  adjacency <- dgraphs::graph.adjacency(graph)
  edge.lengths <- dgraphs::graph.lengths(graph)
} else {
  adjacency <- graph$adj_list
  edge.lengths <- graph$weight_list
}
graph.fit <- fit.metric.graph.lowpass(
  adjacency, edge.lengths, response,
  n.eigenpairs = nrow(locations), eigen.solver = "dense",
  conductance.rule = "inverse.length.power", conductance.alpha = 1,
  conductance.epsilon = 1e-8, eta.grid = c(.001, .01, .1, 1)
)
head(data.frame(site = rownames(locations), response,
                fitted = fitted(graph.fit)))
#>            site  response      fitted
#> site-01 site-01 0.0000000 0.006105107
#> site-02 site-02 0.3246995 0.324089090
#> site-03 site-03 0.6142127 0.612949746
#> site-04 site-04 0.8371665 0.835444574
#> site-05 site-05 0.9694003 0.967406378
#> site-06 site-06 0.9965845 0.994534692
graph.fit$gcv[c("eta.optimal", "gcv.optimal", "effective.df")]
#> $eta.optimal
#> [1] 0.001
#> 
#> $gcv.optimal
#> [1] 0.004647749
#> 
#> $effective.df
#> [1] 19.2978
# Reuse the fitted object's spectral basis, with the selected smoothing fixed.
second.fit <- refit(graph.fit, y = cos(2 * pi * locations[, 1]))
head(fitted(second.fit))
#> [1]  0.99896194  0.94388145  0.78751744  0.54582318  0.24498057 -0.08240949

One-dimensional Gaussian and logistic trend filtering

For a continuous outcome, fit the model over a list of penalties that you supply, and then assess both the association with the predictor and the held-out prediction:

x.continuous <- seq(0, 1, length.out = 40)
y.continuous <- sin(2*pi*x.continuous) + 0.2*cos(seq_along(x.continuous))
gaussian.fit <- fit.gaussian.trend.filtering(
  x.continuous, y.continuous, lambda = c(Inf, 0.01), backend = "banded")
# Nine permutations keep this API example short; use more for inference.
gaussian.assessment <- assess.gaussian.trend.filtering(
  gaussian.fit, permutations = 9, inner.folds = 3, outer.folds = 3, seed = 42,
  additional.tests = c("absolute.deviation", "total.variation"))
gaussian.assessment$tests
#>                 test  statistic   status requested completed failed exceedances
#> 1      squared.error 19.3392875 complete         9         9      0           0
#> 2 absolute.deviation  0.6093144 complete         9         9      0           0
#> 3    total.variation  3.6858303 complete         9         9      0           0
#>   comparison.tolerance resolution p.value p.value.lower p.value.upper
#> 1         1.933929e-09        0.1     0.1           0.1           0.1
#> 2         1.000000e-10        0.1     0.1           0.1           0.1
#> 3         3.685830e-10        0.1     0.1           0.1           0.1
#>   p.value.holm
#> 1          0.3
#> 2          0.3
#> 3          0.3
gaussian.assessment$prediction$scores
#>             model        mse       mae
#> 1   Constant mean 0.55796292 0.6414749
#> 2 Trend filtering 0.03138962 0.1546267

The test statistic is the reduction in the residual sum of squares relative to the sample mean. The tuning is repeated for every permutation. The test assumes that the outcomes are exchangeable under the null hypothesis; equal means with variances that depend on the predictor are not enough. The nested cross-validation is a separate assessment. It compares the smoother with the mean of the outer training data, using MSE and MAE on the same held-out rows. The curve summaries are in the units of the response. Note that lambda = Inf gives a quadratic curve, whereas the baseline of the assessment is a constant. The example shows how the functions are used; with nine permutations it says nothing about the calibration or the power of the test.

The logistic trend-filtering vignette walks through synthetic probabilities, trial weights and held-out Pima predictions.

fit.logistic.trend.filtering() accepts binary observations, or proportions together with the numbers of trials as weights. You must supply the penalties in lambda. They refer to the predictor after it has been rescaled to the range of the training data, not to its original units. Observations with exactly equal predictor values share one value of the curve but keep their own covariates and offsets. To select a penalty by cross-validation, supply foldid; give repeated observations of one subject the same fold label.

Always inspect fit$diagnostics. A fit that reached the iteration limit, detected separation, or failed numerically is not a successful fit. predict() applies the transformations learned from the training data and holds the curve constant outside the training range. No ridge penalty is added unless you request one.

Two solvers are available. backend = "dense" is the reference solver for small problems. The banded solver does not form dense matrices over all knots, yet its stopping checks still examine every knot coefficient. backend = "auto" chooses the banded solver when there are more than 256 distinct predictor values or when the dense solver would need too much memory. Both solvers check an estimate of their memory use before allocating.

Cross-validation for Gaussian, logistic, local polynomial lifting, and graph trend filtering uses cv.penalty.scale = "training_mass" by default. A fold is fitted to less data than the full fit, so its penalty is multiplied by the ratio of the training weight to the full weight. This keeps the penalty in the same proportion to the mean loss. Set cv.penalty.scale = "none" to use the same numerical penalty in every fold, as earlier versions did. The one-dimensional fits scale explicit ridge penalties in the same way and record the settings used in each fold. The assessment functions apply the same rule to their outer and inner training sets. Fits saved by earlier versions keep the earlier convention, and fits with a fixed penalty are not affected.

Hessian-energy regression

The Hessian of a function is its matrix of second derivatives. These methods estimate it locally from the points near each observation and penalize its size. A quadratic penalty reduces curvature everywhere by a moderate amount, whereas an absolute-value penalty can set some second derivatives exactly to zero. The two penalties therefore give different estimates.

Cross-validation works when some points have no observed response (semi-supervised regression). GCV requires a single response vector with no missing values and strictly positive observation weights. The trace needed for GCV can be computed exactly or estimated with random probes; the exact computation costs more, and the estimate depends on a random seed. For the absolute-value penalty, solver = "admm" uses the built-in alternating direction method of multipliers, and the genlasso solver requires that optional package.

To select the penalties of the quadratic method, call the same fitting function with lambda.selection = "cv" or "gcv", and put the settings of the selection in cv.control or gcv.control. A single penalty value and a grid of candidates cannot be given together.

hessian.fit <- fit.ssrhe.hessian.regression(
  X, y, k = 12L, tangent.dim = 1L,
  lambda.selection = "cv", lambda1.grid = c(0.01, 0.1),
  cv.control = list(foldid = rep(1:5, length.out = nrow(X)))
)
hessian.fit$lambda
#> $lambda1
#> [1] 0.01
#> 
#> $lambda2
#> [1] 0
#> 
#> $ridge
#> [1] 0

The separate .cv() and .gcv() fitting functions of version 0.1.0 have been removed. Their old help names now open the help page of fit.ssrhe.hessian.regression(), which explains how to update a call. Calls with fixed lambda1 and optional lambda2 work as before.

Harmonic interpolation

Harmonic smoothing replaces the values inside a region of the graph by values that are, at each vertex, the weighted average of the neighbors’ values, while the values on the boundary of the region stay fixed. A boundary vertex is a vertex of the region that has a neighbor outside it. A vertex of degree one is therefore not a boundary vertex when the region is the whole graph.

The two functions behave differently when the region has no boundary, so they are not interchangeable. perform.harmonic.smoothing() returns its input unchanged with status = "no_boundary", whereas harmonic.smoother() still runs the iteration. Both report numerical convergence only when the change between iterations and the harmonic residual are both small.

harmonic.smoother() returns a harmonic_smoother object with print(), summary(), and plot() methods. It also reports whether the local extrema of the values stopped changing: stability_detected says whether a stable stretch of iterations was found, and stable_iteration identifies it, or is NA if none was found. This is a separate diagnostic from convergence. For example, on a path whose values alternate, recording every second iteration can look stable while the values are still oscillating. Check converged and max_residual, and do not take stable extrema as a sign that the amount of smoothing is right.

Probability mass and subject occupation

support <- matrix(seq(0, 1, length.out = 6), ncol = 1)
occupation <- fit.subject.od(
  support, subject.index = c(1L, 3L, 3L, 6L), method = "empirical"
)
occupation$rho
#> [1] 0.25 0.00 0.50 0.00 0.00 0.25
sum(occupation$rho)
#> [1] 1

In this example the third support point is visited twice and the first and sixth once each. rho is a vector of probability masses, one for each support point, that sums to one. It is not a density per unit volume. normalize.density() converts a fitted field to such a vector; it does not fit anything. If the field has negative values, they are set to zero before normalizing, and the result reports how much was changed.

fit.subject.od() offers more methods than fit.density(), because it can also fit one of the regression smoothers and normalize the result. Two kinds of cross-validation occur here and should not be confused. Visit-level cross-validation holds out visits and takes its folds from visit.foldid. The cross-validation of a regression smoother holds out rows of the support and takes its folds from foldid, or from the argument of that name in the chosen method.

Experimental transport methods

Derivatives estimated at neighboring points are expressed in different local coordinate frames. Parallel transport aligns the frames so that the derivatives can be compared. The functions in this table are experimental, and they return detailed diagnostics about the frames, the operator, and the solver. transported.graph.hessian.operator() serves a different purpose from pttf.operator(): it builds a diagnostic operator and does not lead to a fitted response.

Synthetic datasets: specifications, responses, and reproducibility

The functions whose names begin with synthetic. only describe a data set; no random numbers are drawn until you call materialize.synthetic(). The result is a synthetic_dataset object with print(), plot(), and as.data.frame() methods. vignette("synthetic-datasets", package = "geosmooth") gives complete examples. For new work, use the random-number policy "named.stream.v1". The other seed and algorithm settings exist so that data sets generated by earlier versions can be reproduced exactly; leave them at their defaults unless that is your aim.

Geometry constructors and helpers supplied by dgraphs

The following functions were exported by geosmooth 0.1.0 and are now provided by dgraphs. In existing code, replace geosmooth:: with dgraphs::. Their help pages are in dgraphs; see, for example, help("synthetic.circle", package = "dgraphs").

Sampling constructors supplied by dgraphs

Two sampling descriptions of similar regions do not necessarily draw the same random numbers, because each description also fixes the order of the draws and the algorithm.