Exact Gaussian processes and their nugget

gp(x) in a formula fits an exact Gaussian process, and gp(x, k = ) its Hilbert-space approximation, with the arguments of brms::gp(). This page states where the exact form is not brms's model.

Value

gp() is a formula term, not a free-standing function: bf() reads it at parse time. Its estimates, sdgp and the length scales, are in the Gaussian Process Hyperparameters section of summary(). This page itself documents the term and returns nothing.

The nugget

A squared exponential kernel is near singular, so a nugget on its diagonal keeps it factorable. frm() and frm_sample() read one tape, the field sdgp * L z, with L the Cholesky factor of the kernel plus 1e-12 on the correlation: a nugget of 1e-12 * sdgp^2. brms adds an absolute 1e-12. At sdgp near 1 the two are one number. A relative nugget does not depend on the units of y. brms's does: at y times 1000 its kernel does not factor at the estimates (its run reads Rhat 3.77), and at y divided by 1000 its nugget floors sigma.

The two nuggets give one model until the noise reaches the nugget, near sigma = 1e-6 * sdgp:

With one row per position and an estimated sigma, sigma can read low.
Where positions repeat, or sigma is known through se(), or the family has no sigma, sdgp and the length scale can read low.

The warning

frm() warns with class frmtmb_gp_nugget, on top of frmtmb_warning, when the 90th percentile over the positions of the nugget's variance over the noise at a position reaches sqrt(2 / m) for m positions, which is about one standard error of sigma. The noise at a position is read from the likelihood's curvature, so the rule is the same for every family and route. frm() does not check an importance = or integrate = fit. frm_sample(fit) on a fit with this warning starts at random, because chains started at that mode stay pinned.

What to do about it

frm_sample(laplace = TRUE) integrates the field out, and it mixes fast at low noise, for the outer parameters only (sigma, sdgp, the length scales): it refuses posterior_epred(), posterior_linpred(), posterior_predict(), fitted(), predict(), log_lik(), loo() and conditional_effects(). The drawn route, the default, is the one for the field and for predictions, and at low noise it needs long chains (at sigma 1e-3 and 15 rows, Rhat on sdgp read 1.04 to 2.37 on short chains, and 1.03 on sigma and 1.07 on sdgp after an hour, and not yet converged). gp(x, k = ) has no nugget and reaches brms's numbers when k is large enough for the field. Too small a k is a different failure, with no warning and with convergence 0: sin(8 * x) on 0 to 6 with noise 1e-3 at k = 30 has a basis that stops at frequency 6.5, so it reads sigma 720 times high; k = 60 reads it. Compare several k. A smooth trend written as fixed effects also lowers what the field has to carry.

The order of the field

frm_sample()'s b[i] columns and fit$estimates$b of an exact term hold z, in a farthest-first order of the positions. ranef(), coef() and getME(fit, "b") read the field, at the sorted distinct positions, where brms keeps first appearance. None of these is in the order of the rows of the data. Read a field at the data through posterior_epred() or frm_linpred(), which are right row by row.

Examples

set.seed(1)
d <- data.frame(x = seq(0, 6, length.out = 25))
d$y <- sin(d$x) + rnorm(25, 0, 0.3)
fit <- frm(bf(y ~ gp(x)), data = d)
summary(fit)
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ gp(x) 
#>    Data: d (Number of observations: 25) 
#>  Method: ML   logLik: -12.5176   AIC: 33.0352   BIC: 37.9108 
#> 
#> Gaussian Process Hyperparameters:
#>             Estimate Est.Error l-95% CI u-95% CI
#> sdgp(gpx)       0.73      0.30     0.32     1.65
#> lscale(gpx)     0.23      0.07     0.13     0.41
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept    -0.09      0.47    -1.01     0.82    -0.2     0.84
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     0.29      0.05     0.21      0.4

# the field at the rows of the data, in their order
head(fitted(fit))
#>        Estimate Est.Error        Q2.5     Q97.5
#> [1,] 0.02490738 0.2019663 -0.37093926 0.4207540
#> [2,] 0.23793064 0.1557666 -0.06736626 0.5432275
#> [3,] 0.46481417 0.1347154  0.20077693 0.7288514
#> [4,] 0.68489542 0.1286750  0.43269703 0.9370938
#> [5,] 0.87562173 0.1283887  0.62398456 1.1272589
#> [6,] 1.01552011 0.1303555  0.76002794 1.2710123

# noise far below the field reaches the nugget, and frm() says so
d$y_low <- sin(d$x) + rnorm(25, 0, 1e-7)
fit_low <- withCallingHandlers(
  frm(bf(y_low ~ gp(x)), data = d),
  frmtmb_gp_nugget = function(w) {
    message("nugget warning: ", conditionMessage(w))
    invokeRestart("muffleWarning")
  }
)
#> nugget warning: The noise has reached the nugget of the exact gp() term gp(x) (1.5e+05 of the noise at a position): white noise of sd 1e-6 * sdgp, which the likelihood cannot tell from the data's. Where each position has one row, sigma can read low. Where positions repeat, or sigma is known, or the family has none, sdgp and the length scale can read low. The fitted means at observed positions hold; standard errors and new positions may not.