Package {gpboost}

A string specifying the likelihood function (distribution) of the response variable. The currently supported likelihoods, grouped by response type, are:

The table lists the main likelihood names. Aliases, heteroscedastic and varying-shape variants, and the "_regression" variants of the two-part likelihoods are documented in the detailed descriptions below.

Continuous response: y in (-inf, inf)

• "gaussian": Gaussian likelihood

• "t": t-distribution (e.g., for robust regression). The default approximation is Fisher-Laplace, i.e., the Fisher information is used instead of the observed Hessian in the Laplace approximation of the marginal likelihood.

• "t_fix_df": t-distribution with the degrees-of-freedom (df) held fixed and not estimated

  • The degrees-of-freedom (df) can be set via the likelihood_additional_param parameter. The default is df = 2

• "quantile_regression" / "asymmetric_laplace" : an asymmetric Laplace likelihood for quantile regression. Both names are accepted, as is the alias "quantile"

  • The quantile must be supplied through the likelihood_additional_param parameter and must be strictly between 0 and 1

  The default approximation is Fisher-Laplace, i.e., the Fisher information is used instead of the observed Hessian in the Laplace approximation of the marginal likelihood. The triangular-kernel-curvature (TKC) approximation can be enabled by appending "_triangular_kernel_curvature" or the shorthand "_tkc", for example, "quantile_regression_tkc" or "asymmetric_laplace_triangular_kernel_curvature".

  • experimental: the log-likelihood is not differentiable at the kinks where 'y_i' equals the location parameter, and the quasi-Newton mode finding can therefore stall at a point that is not the exact posterior mode. Appending "_ssn_alm", for example "quantile_regression_ssn_alm", enables an exact check of the (non-smooth) optimality conditions after the mode finding and, if the check fails, a refinement of the mode with a semismooth Newton method applied to the subproblems of an augmented Lagrangian method. This makes the approximate marginal likelihood a well-defined function of the parameters, i.e., independent of the path taken by the optimizer and of the number of threads, at the price of additional linear solves during mode finding. This matters for standard errors and for likelihood-based model comparison, but it does not necessarily improve predictive accuracy; increasing max_num_restarts_lbfgs can be a cheaper way of avoiding poor optima. Appending "_admm_ssn_alm" instead warm starts the refinement with an ADMM phase, which typically halves its cost when matrix_inversion_method = "cholesky" and is not recommended for iterative methods

• "gaussian_heteroscedastic": Gaussian likelihood where the mean is related to fixed and random effects and the log-error variance is related to fixed effects only (linear predictor or GPBoost algorithm). The estimated coefficients of the log-variance model are returned alongside the mean-model coefficients (with the suffix "_scale"). Fisher-Laplace is the default and currently the only implemented approximation.

• "gaussian_heteroscedastic_fixed_and_random": Gaussian likelihood where both the mean and the variance are related to fixed and random effects. This is currently only implemented for GPs with a 'vecchia' approximation. Fisher-Laplace is the default and currently the only implemented approximation. By default, Cholesky factorizations are used. If init_cov_pars is not given, the initial values are obtained from a homoscedastic Gaussian process model (unless some weights are zero).

Positive continuous response: y in (0, inf)

• "gamma": Gamma likelihood with a log link function

• "gamma_varying_shape": As "gamma", but the gamma shape varies across observations and is modeled by an additional fixed-effects-only predictor: log(shape) = F_s(X) (F_s(X) = linear predictor or the GPBoost algorithm), while the log mean log(mu) = F(X) + Zb is related to both fixed and random effects. Note that the shape also governs the dispersion, var(y) = mu^2 / shape. The estimated coefficients of the log-shape model are returned alongside the mean-model coefficients (with the suffix "_shape"). See also the corresponding two-part likelihoods "hurdle_gamma_varying_shape" and "hurdle_regression_gamma_varying_shape" below

• "lognormal": Log-normal likelihood with a log link function

• "gpd": Generalized Pareto likelihood for finite positive responses. The log scale parameter equals the latent predictor eta (sum of fixed and random effects), sigma = exp(eta), and the estimated auxiliary parameter is shape > -0.5

• "egpd_power": Naveau power-carrier extended generalized Pareto likelihood with auxiliary parameters shape and kappa

• "egpd_power_mixture": Naveau ordered power-mixture carrier with auxiliary parameters shape, kappa1, delta_kappa, and p

• "egpd_beta": Naveau beta-carrier extended generalized Pareto likelihood with auxiliary parameters shape and delta

• "egpd_power_beta": Naveau power-beta carrier with auxiliary parameters shape, delta, and kappa

Non-negative continuous / semicontinuous response: y in [0, inf) with point mass at 0

• "tweedie": Compound Poisson–Gamma Tweedie likelihood with a log link, variance phi * mu^p, and 1.01 < p < 1.99. The dispersion phi and power p are estimated

• "tweedie_fixed_p": The same Tweedie likelihood with p fixed through likelihood_additional_param; only phi is estimated. The fixed power must satisfy 1.01 < p < 1.99

• "tweedie_joint", "tweedie_joint_fixed_p": The same Tweedie model, but the observed number of events (e.g., claims) N, given in the first column of additional_likelihood_data, is used jointly with the aggregate response y. The likelihood is the joint density of (y, N) of the compound Poisson–Gamma representation, N ~ Poisson(mu^(2-p) / (phi * (2-p))) and y | N = n ~ Gamma(n * (2-p) / (p-1), scale = phi * (p-1) * mu^(p-1)), see Jorgensen and de Souza (1994, Scandinavian Actuarial Journal). Given phi and p, the mean model is the same as for "tweedie"; N adds information for estimating phi and p. N must be a non-negative integer that is 0 if and only if y is 0. Predictions are the same as for "tweedie", and the "test_neg_log_likelihood" metric uses the marginal density of y since N is not known for new data

• "tweedie_varying_dispersion", "tweedie_varying_dispersion_fixed_p": As "tweedie" and "tweedie_fixed_p", but the dispersion phi varies across observations: log(phi) = F_d(X) is related to fixed effects only (linear predictor or GPBoost algorithm), while log(mu) = F(X) + Zb is related to both fixed and random effects. The power p is then the only (auxiliary) parameter. The estimated coefficients of the log-dispersion model are returned alongside the mean-model coefficients (with the suffix "_dispersion")

• "tweedie_joint_varying_dispersion", "tweedie_joint_varying_dispersion_fixed_p": The joint (y, N) likelihood of "tweedie_joint" with the varying dispersion of "tweedie_varying_dispersion"

• "hurdle_<base>": Two-part likelihoods for non-negative response variables with an excess probability 'p0' of exact zeros. They combine a point mass 'p0' at zero with a base distribution with support 'y > 0' for the remaining mass '1 - p0'. The fixed effects 'F(X)' and random effects 'Zb' enter only through the base component: 'exp(F(X) + Zb)' is its mean or scale parameter, not the unconditional response mean, which is E(y) = (1 - p0) * base_mean. The structural-zero probability 'p0' is estimated jointly with the base auxiliary parameters. Currently supported bases: "hurdle_gamma", "hurdle_lognormal", and the extreme-value bases "hurdle_gpd", "hurdle_egpd_power", "hurdle_egpd_power_mixture", "hurdle_egpd_beta", "hurdle_egpd_power_beta". For all these bases, the "zero_inflated" prefix is accepted as an alias, e.g. "zero_inflated_gamma" maps to "hurdle_gamma". See "zero_inflated_<base>" under count responses for the corresponding count-data family

• "hurdle_regression_<base>": As "hurdle_<base>", but the structural-zero probability is modeled as a logistic regression on the covariates instead of as a constant (e.g. "hurdle_regression_gamma", "hurdle_regression_lognormal", "hurdle_regression_gpd", "hurdle_regression_egpd_power", "hurdle_regression_egpd_power_mixture", "hurdle_regression_egpd_beta", "hurdle_regression_egpd_power_beta"). The zero probability is then pi_i = 1 / (1 + exp(-x_i'alpha)), modeled through a second fixed-effects-only predictor that reuses the same design matrix X as the response model (the response predictor carries the random effects, the zero predictor does not). The estimated zero-model coefficients alpha are returned by get_coef() alongside the response-model coefficients, with the suffix "_zero"

• "hurdle_gamma_varying_shape", "hurdle_regression_gamma_varying_shape": As "hurdle_gamma" and "hurdle_regression_gamma", but with the observation-specific gamma shape of "gamma_varying_shape". For "hurdle_regression_gamma_varying_shape" there are three predictors: the response mean, the structural-zero logit (coefficients with the suffix "_zero"), and log(shape)

• "zero_censored_power_transformed_normal": Likelihood of a censored and power-transformed normal variable for modeling data with a point mass at 0 and a continuous distribution for y > 0. The model used is Y = max(0,X)^lambda, X ~ N(mu, sigma^2), where mu = F(X) + Zb, and sigma and lambda are (auxiliary) parameters that are estimated. For more details on this model, see Sigrist et al. (2012, AOAS) "A dynamic nonstationary spatio-temporal model for short term prediction of precipitation"

• "zero_censored_power_transformed_normal_heteroscedastic": As "zero_censored_power_transformed_normal", but the standard deviation sigma of the latent normal variable varies across observations: log(sigma) = F_2(X) is related to fixed effects only (linear predictor or GPBoost algorithm), while mu = F(X) + Zb is related to both fixed and random effects. lambda is then the only (auxiliary) parameter that is estimated

• "zero_censored_shifted_gamma": Zero-censored shifted gamma likelihood for modeling data with a point mass at 0 and a continuous distribution for y > 0. The model used is Y = max(Z - xi, 0), where Z follows a gamma distribution with mean mu = exp(F(X) + Zb) and shape k. The shape k and shift xi are (auxiliary) parameters that are estimated. This is the version of "zero_one_censored_shifted_gamma" (Sigrist and Stahel, 2011) without the upper censoring at 1

• "zero_censored_shifted_gamma_varying_shape": As "zero_censored_shifted_gamma", but the shape k varies across observations: log(k) = F_s(X) is related to fixed effects only (covariates and / or the GPBoost tree-boosting algorithm; no random effects / GPs for the shape), while log(mu) = F(X) + Zb is related to both fixed and random effects. The shift xi is then the only (auxiliary) parameter that is estimated. The estimated coefficients of the log-shape model are returned alongside the mean-model coefficients (with the suffix "_shape")

Count response: y in {0, 1, 2, ...}

• "poisson": Poisson likelihood with a log link function

• "negative_binomial": negative binomial likelihood with a log link function (aka "nbinom2", "negative_binomial_2"). The variance is mu * (mu + r) / r, mu = mean, r = shape, with this parametrization

• "negative_binomial_1": Negative binomial 1 (aka "nbinom1") likelihood with a log link function. The variance is mu * (1 + phi), mu = mean, phi = dispersion, with this parametrization

• "zero_inflated_<base>": Two-part count likelihoods that combine a point mass 'p0' at zero with a count base distribution for the remaining mass '1 - p0' (the base can itself generate additional zeros). As for the hurdle likelihoods, 'exp(F(X) + Zb)' is the mean of the (non-structural) base component - not the unconditional response mean - so E(y) = (1 - p0) * base_mean, and 'p0' is estimated jointly with the base auxiliary parameters. Currently supported bases: "zero_inflated_poisson", "zero_inflated_negative_binomial" (aux: shape; aka "zero_inflated_nbinom2"), "zero_inflated_negative_binomial_1" (aux: dispersion; aka "zero_inflated_nbinom1")

• "zero_inflated_regression_<base>": As "zero_inflated_<base>", but with the logistic structural-zero model described under "hurdle_regression_<base>" above: "zero_inflated_regression_poisson", "zero_inflated_regression_negative_binomial", and "zero_inflated_regression_negative_binomial_1"

Binary response: y in {0, 1}

• "bernoulli_logit": Bernoulli likelihood with a logit link function for binary classification. Aliases: "binary", "binary_logit"

• "bernoulli_probit": Bernoulli likelihood with a probit link function for binary classification. Aliases: "binary_probit"

Proportion / fractional / bounded response: y in [0, 1]

• "quasi_bernoulli_logit": quasi-Bernoulli likelihood with a logit link function for y in [0,1]. Aliases: "quasi_binary", "quasi_binary_logit"

• "quasi_bernoulli_probit": quasi-Bernoulli likelihood with a probit link function for y in [0,1]. Aliases: "quasi_binary_probit"

• "binomial_logit": Binomial likelihood with a logit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials. Aliases: "binomial"

• "binomial_probit": Binomial likelihood with a probit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials

• "beta_binomial": Beta-binomial likelihood with a logit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials. Aliases: "betabinomial", "beta-binomial"

• "beta" : Beta likelihood with a logit link function (parametrization of Ferrari and Cribari-Neto, 2004). The response must lie strictly in (0,1)

• "zoctn": Zero-one censored transformed normal likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). The model used is T ~ N(mu, sigma^2), W = max(min(T,1),0), and Y = g(W), where g(x) = expit(a + b * logit(x)) for x in (0,1), mu = F(X) + Z_RE u, u denotes the random effects, Z_RE is their design matrix, and sigma, a, and b are (auxiliary) parameters that are estimated. For more details on this model, see Qiang and Sigrist (2026)

• "zero_one_censored_transformed_beta": Zero-one censored transformed beta likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). If T follows a beta distribution with mean mu = expit(F(X) + Zb) and precision phi, the observed response is obtained by applying the linear transformation Y = (1 + 2u) * T - u and censoring the result to [0,1]. The precision phi and shift u are (auxiliary) parameters that are estimated. For more details on this model, see Kosmidis and Zeileis (2025)

• "zero_one_censored_shifted_gamma": Zero-one censored shifted gamma likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). The model used is Y = min(max(Z - xi, 0), 1), where Z follows a gamma distribution with mean mu = exp(F(X) + Zb) and shape k. The shape k and shift xi are (auxiliary) parameters that are estimated. For more details on this model, see Sigrist and Stahel (2011)

Notes

• The first lines in the likelihoods source file contain additional comments on the specific parametrizations used

• Other likelihoods can be implemented upon request

A vector or matrix whose columns are categorical grouping variables. The elements being group levels defining grouped random effects. The elements of 'group_data' can be integer, double, or character. The number of columns corresponds to the number of grouped (intercept) random effects

A vector or matrix with numeric covariate data for grouped random coefficients

A vector with integer indices that indicate the corresponding categorical grouping variable (=columns) in 'group_data' for every covariate in 'group_rand_coef_data'. Counting starts at 1. The length of this index vector must equal the number of covariates in 'group_rand_coef_data'. For instance, c(1,1,2) means that the first two covariates (=first two columns) in 'group_rand_coef_data' have random coefficients corresponding to the first categorical grouping variable (=first column) in 'group_data', and the third covariate (=third column) in 'group_rand_coef_data' has a random coefficient corresponding to the second grouping variable (=second column) in 'group_data'

A vector of type logical (boolean). Indicates whether intercept random effects are dropped (only for random coefficients). If drop_intercept_group_rand_effect[k] is TRUE, the intercept random effect number k is dropped / not included. Only random effects with random slopes can be dropped.

A matrix with numeric coordinates (= inputs / features) for defining Gaussian processes

A vector or matrix with numeric covariate data for Gaussian process random coefficients

A string specifying the covariance function for the Gaussian process. Available options:

• "matern": Matern covariance function with the smoothness specified by the cov_fct_shape parameter (using the parametrization of Rasmussen and Williams, 2006)

• "matern_estimate_shape": same as "matern" but the smoothness parameter is also estimated

• "matern_space_time": Spatio-temporal Matern covariance function with different range parameters for space and time. Note that the first column in gp_coords must correspond to the time dimension

• "space_time_gneiting": Spatio-temporal covariance function given in Eq. (16) of Gneiting (2002). Note that the first column in gp_coords must correspond to the time dimension. This covariance has seven parameters (in the following order: sigma2, a, c, alpha, nu, beta, delta) which are all estimated by default. You can disable the estimation of some of these parameter using the 'estimate_cov_par_index' argument of the params argument in either the fit function of a gp_model object or the set_optim_params function prior to estimation.

• "matern_ard": anisotropic Matern covariance function with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of gp_coords

• "matern_ard_estimate_shape": same as "matern_ard" but the smoothness parameter is also estimated

• "exponential": Exponential covariance function (using the parametrization of Diggle and Ribeiro, 2007)

• "gaussian": Gaussian, aka squared exponential, covariance function (using the parametrization of Diggle and Ribeiro, 2007)

• "gaussian_ard": anisotropic Gaussian, aka squared exponential, covariance function with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of gp_coords

• "powered_exponential": powered exponential covariance function with the exponent specified by the cov_fct_shape parameter (using the parametrization of Diggle and Ribeiro, 2007)

• "wendland": Compactly supported Wendland covariance function (using the parametrization of Bevilacqua et al., 2019, AOS)

• "linear": linear covariance function. This corresponds to a Bayesian linear regression model with a Gaussian prior on the coefficients with a constant variance diagonal prior covariance, and the prior variance is estimated using empirical Bayes.

• "hurst": Hurst covariance function cov(s, s') = (sigma2 / 2) * ( ||s||^(2H) + ||s'||^(2H) - ||s - s'||^(2H) ), 0 < H < 1. For H = 0.5, this corresponds to Brownian motion in one dimension (-> see the 'estimate_cov_par_index' argument). This is the covariance for cov_fct_order = 1 (default). For cov_fct_order = 2, the second-order Hurst covariance cov(s, s') = sigma2 / (2 (2H - 1)) * ( ||s - s'||^(2H) - ||s||^(2H) - ||s'||^(2H) + 2H (s^T s') (||s||^(2H-2) + ||s'||^(2H-2)) ), 1 < H < 2, is used. In one dimension and for s, s' >= 0, H = 1.5 corresponds to an integrated Wiener process. See cov_fct_order for more details

• "hurst_ard": Hurst covariance function with with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of “gp_coords“ except for the first coordinate which has a range parameter of 1 due to identifiability with the marginal variance: cov(s, s') = (\sigma^2/2)\left[ \left(s_1^2 + \sum_{k=2}^d (s_k/l_k)^2\right)^H + \left({s'}_1^2 + \sum_{k=2}^d ({s'}_k/l_k)^2\right)^H - \left((s_1-{s'}_1)^2 + \sum_{k=2}^d ((s_k-{s'}_k)/l_k)^2\right)^H \right] . For cov_fct_order = 2, the second-order Hurst covariance (see "hurst") is applied to the scaled coordinates (s_1, s_2/l_2, ..., s_d/l_d)

• "ar1_mf_<base>": two-level autoregressive multifidelity covariance defined by f_H(x)=\rho f_L(x)+\delta(x), where f_L and \delta are independent Gaussian processes using the same base covariance type but separate parameter vectors.

  • For example, use "ar1_mf_matern", "ar1_mf_matern_ard", or "ar1_mf_matern_estimate_shape".

  • The last column of gp_coords must be 0 for low fidelity or 1 for high fidelity. All preceding columns are input coordinates for the GPs.

  • Covariance parameters are ordered as [low-fidelity base parameters, discrepancy base parameters, rho]. The two base blocks use the ordinary parameter ordering of the base covariance, and rho is unrestricted and can be negative.

  • Any supported base covariance except "wendland" can be used. Correlation tapering and Gaussian-process random coefficients are currently not supported for this model.

A numeric specifying the shape parameter of the covariance function (e.g., smoothness parameter for Matern and Wendland covariance) This parameter is irrelevant for some covariance functions such as the exponential or Gaussian

A string specifying the large data approximation for Gaussian processes. Available options:

• "none": No approximation

• "vecchia": Vecchia approximation; see Sigrist (2022, JMLR) for more details. For "matern_space_time" and the anisotropic ARD covariance functions ("matern_ard", "gaussian_ard", "matern_ard_estimate_shape"), neighbors are selected by absolute correlation, i.e., using distances in coordinates scaled by the range parameters, and they are redetermined during parameter estimation. For "space_time_gneiting" and "ar1_mf_<base>", neighbors are also selected by absolute correlation by default. Use gp_approx = "vecchia_euclidean" for Euclidean-distance selection.

• "full_scale_vecchia": Vecchia-inducing points full-scale (VIF) approximation; see Gyger, Furrer, and Sigrist (2026, JMLR) for more details

• "tapering": The covariance function is multiplied by a compactly supported Wendland correlation function

• "fitc": Fully Independent Training Conditional approximation aka modified predictive process approximation

• "full_scale_tapering": Full-scale approximation combining an inducing point / predictive process approximation with tapering on the residual process; see Gyger, Furrer, and Sigrist (2026, JUQ) for more details

• "vecchia_latent": similar as "vecchia" but a Vecchia approximation is applied to the latent Gaussian process for likelihood == "gaussian". For likelihood != "gaussian", "vecchia" and "vecchia_latent" are equivalent

An integer specifying the number of parallel threads for OMP. If num_parallel_threads = NULL, the default number of threads of the session is used, see gpb.get.default.num.threads. Initially, this is the number of physical performance cores, which is determined automatically as follows. On CPUs whose cores have different speeds (e.g., the performance and efficiency cores of recent Intel CPUs or Apple silicon), the fastest cores are counted, and hyperthreads are counted only once, since slow threads can slow down an entire model. The next fastest cores are added if the fastest ones alone would leave only a single thread. On Linux, a CPU bandwidth limit of a control group (e.g., of a container) is respected as well. If the environment variable OMP_NUM_THREADS is set, it takes precedence over both the topology of the CPU and the limit of the control group, and the number of threads that OpenMP is configured to use is taken as it is if the cores of the CPU cannot be determined. In either case, the limits of the OpenMP runtime itself still apply, in particular the limit of the contention group (OMP_THREAD_LIMIT). For ordinary use, leave num_parallel_threads unspecified to use this default. The default of the session can be changed with gpb.tune.num.threads, which benchmarks different numbers of threads, or with gpb.set.default.num.threads. Setting num_parallel_threads = 1 disables the OpenMP parallelization of the model and can substantially increase the runtime. A single thread should be chosen deliberately, e.g., to distribute the resources among concurrent model fits or for a specific test; it is not required for correctness or reproducibility

A boolean. If TRUE, GPU acceleration will be used if supported

A string specifying the method used for inverting covariance matrices. Available options:

• "default": iterative methods where possible, otherwise Cholesky factorization. Cholesky factorization is also used for likelihood = "gaussian_heteroscedastic_fixed_and_random"

• "cholesky": Cholesky factorization

• "iterative": iterative methods. A combination of the conjugate gradient, the Lanczos algorithm, and other methods.

  This is currently only supported for the following cases:

  • grouped random effects with more than one level

  • likelihood != "gaussian" and gp_approx == "vecchia" (non-Gaussian likelihoods with a Vecchia-Laplace approximation)

  • likelihood != "gaussian" and gp_approx == "full_scale_vecchia" (non-Gaussian likelihoods with a VIF approximation)

  • likelihood == "gaussian" and gp_approx == "full_scale_tapering" (Gaussian likelihood with a full-scale tapering approximation)

A vector with sample weights. For a Gaussian likelihood, the error variance ("nugget") for observation i is divided by weights[i]. For non-Gaussian likelihoods, the conditional log-likelihood contribution of observation i is multiplied by weights[i]. Consequently, weights affect the estimation of both random and fixed effects. Note that a Gaussian likelihood is calculated via the Laplace approximation when gp_approx = "vecchia_latent", when likelihood = "gaussian_latent", and when grouped random effects are combined with a Vecchia-approximated Gaussian process. In these cases, the weights act as for a non-Gaussian likelihood, i.e., the Gaussian log-likelihood contribution of observation i is multiplied by weights[i] instead of the error variance being divided by it.

A vector or matrix with observation-level data that some likelihoods require in addition to the response variable y (one row per data point). Currently, this is only used by the joint Tweedie likelihoods ("tweedie_joint" and its variants), for which its first column contains the observed number of events (e.g., claims) N. N must be a non-negative integer that is 0 if and only if y is 0

A numeric with a learning rate for the likelihood for generalized Bayesian inference (only non-Gaussian likelihoods)

A numeric specifying the range parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)

A numeric specifying the shape (=smoothness) parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)

An integer specifying the number of neighbors for the Vecchia and VIF approximations. Internal default values if NULL:

• 20 for gp_approx = "vecchia"

• 30 for gp_approx = "full_scale_vecchia"

Note: for prediction, the number of neighbors can be set through the 'num_neighbors_pred' parameter in the 'set_prediction_data' function. By default, num_neighbors_pred = 2 * num_neighbors. Further, the type of Vecchia approximation used for making predictions is set through the 'vecchia_pred_type' parameter in the 'set_prediction_data' function

A string specifying the ordering used in the Vecchia approximation. Available options:

• "none": the default ordering in the data is used

• "random": a random ordering

• "time": ordering accorrding to time (only for space-time models)

• "time_random_space": ordering according to time and randomly for all spatial points with the same time points (only for space-time models)

A string specifying the method for choosing inducing points Available options:

• "kmeans++: the k-means++ algorithm

• "cover_tree": the cover tree algorithm

• "random": random selection from data points

An integer specifying the number of inducing points / knots for FITC, full_scale_tapering, and VIF approximations. Internal default values if NULL:

• 500 for gp_approx = "FITC" and gp_approx = "full_scale_tapering"

• 200 for gp_approx = "full_scale_vecchia"

A numeric specifying the radius (= "spatial resolution") for the cover tree algorithm

An integer specifying the seed used for model creation (e.g., random ordering in Vecchia approximation)

A vector with elements indicating independent realizations of random effects / Gaussian processes (same values = same process realization). The elements of 'cluster_ids' can be integer, double, or character.

A numeric specifying an additional parameter for the likelihood which cannot be estimated for this likelihood (e.g., degrees of freedom for likelihood = "t_fix_df"). This is not to be confused with any auxiliary parameters that can be estimated and accessed through the function get_aux_pars after estimation. Note that this likelihood_additional_param parameter is irrelevant for many likelihoods. If likelihood_additional_param = NULL, the following default and mandatory-value rules apply:

• df = 2 for likelihood = "t_fix_df"

• No default is used for likelihood = "asymmetric_laplace"; a quantile strictly between 0 and 1 is required

• No default is used for the Tweedie likelihoods with a fixed power ("tweedie_fixed_p" and the other "_fixed_p" variants); a power strictly between 1.01 and 1.99 is required

A logical. For an ar1_mf_<base> covariance, whether marginal mean models are fitted independently for low and high fidelity. If TRUE (the default), linear covariates supplied to fit() are internally expanded into separate low- and high-fidelity coefficient blocks, and the fidelity indicator is automatically added as a feature when the model is used in the GPBoost algorithm. Ignored for other covariance functions.

An integer specifying the order m of the "hurst" and "hurst_ard" covariance functions (also when used as base covariance in "ar1_mf_hurst" and "ar1_mf_hurst_ard"). Currently, the orders 1 (default) and 2 are supported. The Hurst exponent H of the order m satisfies m - 1 < H < m, and H = m - 0.5 is used as initial value. Order 1 is a fractional Brownian motion / field anchored at the origin (b(0) = 0), and order 2 is a second-order Hurst process / field anchored at the origin such that both the process and its gradient are zero there (b(0) = 0, grad b(0) = 0). The origin of the coordinates thus is a part of the model, and the coordinates are not centered internally. The order is not estimated. It can be chosen, e.g., by comparing the marginal likelihoods or by cross-validation.

• Continuous-time RW1 prior: use cov_function = "hurst", cov_fct_order = 1, and fix H = 0.5. This is a Brownian motion with b(0) = 0 and b'(t) = sqrt(q) W'(t), where W' denotes white noise and q = sigma2

• Continuous-time RW2 prior: use cov_function = "hurst", cov_fct_order = 2, and fix H = 1.5. This is an integrated Wiener process with b(0) = b'(0) = 0 and b”(t) = sqrt(q) W'(t), where q = 3 * sigma2

For these priors, use a one-dimensional time coordinate t >= 0 whose origin t = 0 is the desired anchor (shift the time points if necessary). H can be fixed using the 'init_cov_pars' and 'estimate_cov_par_index' arguments of the params argument, see the examples of fitGPModel. If the polynomial null space of the intrinsic model should not be penalized, include an intercept (RW1) or an intercept and a linear time trend (RW2) as fixed effects. Note that the order 2 model with H = 1.5 is the continuous-time RW2 prior evaluated at the time points and not the usual discrete intrinsic RW2 prior on a lattice: for equally spaced time points, adjacent second differences are correlated (correlation 1/4) whereas they are independent for the discrete RW2

A numeric with the number of samples. This is only used for iid models

A boolean. If TRUE, the data (groups, coordinates, covariate data for random coefficients, cluster IDs, weights, and additional likelihood data) is freed in R after initialization

Discontinued. Use the argument gp_approx instead

A string specifying the type of Vecchia approximation used for making predictions. This is discontinued here. Use the function 'set_prediction_data' to specify this

an integer specifying the number of neighbors for making predictions. This is discontinued here. Use the function 'set_prediction_data' to specify this

A string specifying the likelihood function (distribution) of the response variable. The currently supported likelihoods, grouped by response type, are:

The table lists the main likelihood names. Aliases, heteroscedastic and varying-shape variants, and the "_regression" variants of the two-part likelihoods are documented in the detailed descriptions below.

Continuous response: y in (-inf, inf)

• "gaussian": Gaussian likelihood

• "t": t-distribution (e.g., for robust regression). The default approximation is Fisher-Laplace, i.e., the Fisher information is used instead of the observed Hessian in the Laplace approximation of the marginal likelihood.

• "t_fix_df": t-distribution with the degrees-of-freedom (df) held fixed and not estimated

  • The degrees-of-freedom (df) can be set via the likelihood_additional_param parameter. The default is df = 2

• "quantile_regression" / "asymmetric_laplace" : an asymmetric Laplace likelihood for quantile regression. Both names are accepted, as is the alias "quantile"

  • The quantile must be supplied through the likelihood_additional_param parameter and must be strictly between 0 and 1

  The default approximation is Fisher-Laplace, i.e., the Fisher information is used instead of the observed Hessian in the Laplace approximation of the marginal likelihood. The triangular-kernel-curvature (TKC) approximation can be enabled by appending "_triangular_kernel_curvature" or the shorthand "_tkc", for example, "quantile_regression_tkc" or "asymmetric_laplace_triangular_kernel_curvature".

  • experimental: the log-likelihood is not differentiable at the kinks where 'y_i' equals the location parameter, and the quasi-Newton mode finding can therefore stall at a point that is not the exact posterior mode. Appending "_ssn_alm", for example "quantile_regression_ssn_alm", enables an exact check of the (non-smooth) optimality conditions after the mode finding and, if the check fails, a refinement of the mode with a semismooth Newton method applied to the subproblems of an augmented Lagrangian method. This makes the approximate marginal likelihood a well-defined function of the parameters, i.e., independent of the path taken by the optimizer and of the number of threads, at the price of additional linear solves during mode finding. This matters for standard errors and for likelihood-based model comparison, but it does not necessarily improve predictive accuracy; increasing max_num_restarts_lbfgs can be a cheaper way of avoiding poor optima. Appending "_admm_ssn_alm" instead warm starts the refinement with an ADMM phase, which typically halves its cost when matrix_inversion_method = "cholesky" and is not recommended for iterative methods

• "gaussian_heteroscedastic": Gaussian likelihood where the mean is related to fixed and random effects and the log-error variance is related to fixed effects only (linear predictor or GPBoost algorithm). The estimated coefficients of the log-variance model are returned alongside the mean-model coefficients (with the suffix "_scale"). Fisher-Laplace is the default and currently the only implemented approximation.

• "gaussian_heteroscedastic_fixed_and_random": Gaussian likelihood where both the mean and the variance are related to fixed and random effects. This is currently only implemented for GPs with a 'vecchia' approximation. Fisher-Laplace is the default and currently the only implemented approximation. By default, Cholesky factorizations are used. If init_cov_pars is not given, the initial values are obtained from a homoscedastic Gaussian process model (unless some weights are zero).

Positive continuous response: y in (0, inf)

• "gamma": Gamma likelihood with a log link function

• "gamma_varying_shape": As "gamma", but the gamma shape varies across observations and is modeled by an additional fixed-effects-only predictor: log(shape) = F_s(X) (F_s(X) = linear predictor or the GPBoost algorithm), while the log mean log(mu) = F(X) + Zb is related to both fixed and random effects. Note that the shape also governs the dispersion, var(y) = mu^2 / shape. The estimated coefficients of the log-shape model are returned alongside the mean-model coefficients (with the suffix "_shape"). See also the corresponding two-part likelihoods "hurdle_gamma_varying_shape" and "hurdle_regression_gamma_varying_shape" below

• "lognormal": Log-normal likelihood with a log link function

• "gpd": Generalized Pareto likelihood for finite positive responses. The log scale parameter equals the latent predictor eta (sum of fixed and random effects), sigma = exp(eta), and the estimated auxiliary parameter is shape > -0.5

• "egpd_power": Naveau power-carrier extended generalized Pareto likelihood with auxiliary parameters shape and kappa

• "egpd_power_mixture": Naveau ordered power-mixture carrier with auxiliary parameters shape, kappa1, delta_kappa, and p

• "egpd_beta": Naveau beta-carrier extended generalized Pareto likelihood with auxiliary parameters shape and delta

• "egpd_power_beta": Naveau power-beta carrier with auxiliary parameters shape, delta, and kappa

Non-negative continuous / semicontinuous response: y in [0, inf) with point mass at 0

• "tweedie": Compound Poisson–Gamma Tweedie likelihood with a log link, variance phi * mu^p, and 1.01 < p < 1.99. The dispersion phi and power p are estimated

• "tweedie_fixed_p": The same Tweedie likelihood with p fixed through likelihood_additional_param; only phi is estimated. The fixed power must satisfy 1.01 < p < 1.99

• "tweedie_joint", "tweedie_joint_fixed_p": The same Tweedie model, but the observed number of events (e.g., claims) N, given in the first column of additional_likelihood_data, is used jointly with the aggregate response y. The likelihood is the joint density of (y, N) of the compound Poisson–Gamma representation, N ~ Poisson(mu^(2-p) / (phi * (2-p))) and y | N = n ~ Gamma(n * (2-p) / (p-1), scale = phi * (p-1) * mu^(p-1)), see Jorgensen and de Souza (1994, Scandinavian Actuarial Journal). Given phi and p, the mean model is the same as for "tweedie"; N adds information for estimating phi and p. N must be a non-negative integer that is 0 if and only if y is 0. Predictions are the same as for "tweedie", and the "test_neg_log_likelihood" metric uses the marginal density of y since N is not known for new data

• "tweedie_varying_dispersion", "tweedie_varying_dispersion_fixed_p": As "tweedie" and "tweedie_fixed_p", but the dispersion phi varies across observations: log(phi) = F_d(X) is related to fixed effects only (linear predictor or GPBoost algorithm), while log(mu) = F(X) + Zb is related to both fixed and random effects. The power p is then the only (auxiliary) parameter. The estimated coefficients of the log-dispersion model are returned alongside the mean-model coefficients (with the suffix "_dispersion")

• "tweedie_joint_varying_dispersion", "tweedie_joint_varying_dispersion_fixed_p": The joint (y, N) likelihood of "tweedie_joint" with the varying dispersion of "tweedie_varying_dispersion"

• "hurdle_<base>": Two-part likelihoods for non-negative response variables with an excess probability 'p0' of exact zeros. They combine a point mass 'p0' at zero with a base distribution with support 'y > 0' for the remaining mass '1 - p0'. The fixed effects 'F(X)' and random effects 'Zb' enter only through the base component: 'exp(F(X) + Zb)' is its mean or scale parameter, not the unconditional response mean, which is E(y) = (1 - p0) * base_mean. The structural-zero probability 'p0' is estimated jointly with the base auxiliary parameters. Currently supported bases: "hurdle_gamma", "hurdle_lognormal", and the extreme-value bases "hurdle_gpd", "hurdle_egpd_power", "hurdle_egpd_power_mixture", "hurdle_egpd_beta", "hurdle_egpd_power_beta". For all these bases, the "zero_inflated" prefix is accepted as an alias, e.g. "zero_inflated_gamma" maps to "hurdle_gamma". See "zero_inflated_<base>" under count responses for the corresponding count-data family

• "hurdle_regression_<base>": As "hurdle_<base>", but the structural-zero probability is modeled as a logistic regression on the covariates instead of as a constant (e.g. "hurdle_regression_gamma", "hurdle_regression_lognormal", "hurdle_regression_gpd", "hurdle_regression_egpd_power", "hurdle_regression_egpd_power_mixture", "hurdle_regression_egpd_beta", "hurdle_regression_egpd_power_beta"). The zero probability is then pi_i = 1 / (1 + exp(-x_i'alpha)), modeled through a second fixed-effects-only predictor that reuses the same design matrix X as the response model (the response predictor carries the random effects, the zero predictor does not). The estimated zero-model coefficients alpha are returned by get_coef() alongside the response-model coefficients, with the suffix "_zero"

• "hurdle_gamma_varying_shape", "hurdle_regression_gamma_varying_shape": As "hurdle_gamma" and "hurdle_regression_gamma", but with the observation-specific gamma shape of "gamma_varying_shape". For "hurdle_regression_gamma_varying_shape" there are three predictors: the response mean, the structural-zero logit (coefficients with the suffix "_zero"), and log(shape)

• "zero_censored_power_transformed_normal": Likelihood of a censored and power-transformed normal variable for modeling data with a point mass at 0 and a continuous distribution for y > 0. The model used is Y = max(0,X)^lambda, X ~ N(mu, sigma^2), where mu = F(X) + Zb, and sigma and lambda are (auxiliary) parameters that are estimated. For more details on this model, see Sigrist et al. (2012, AOAS) "A dynamic nonstationary spatio-temporal model for short term prediction of precipitation"

• "zero_censored_power_transformed_normal_heteroscedastic": As "zero_censored_power_transformed_normal", but the standard deviation sigma of the latent normal variable varies across observations: log(sigma) = F_2(X) is related to fixed effects only (linear predictor or GPBoost algorithm), while mu = F(X) + Zb is related to both fixed and random effects. lambda is then the only (auxiliary) parameter that is estimated

• "zero_censored_shifted_gamma": Zero-censored shifted gamma likelihood for modeling data with a point mass at 0 and a continuous distribution for y > 0. The model used is Y = max(Z - xi, 0), where Z follows a gamma distribution with mean mu = exp(F(X) + Zb) and shape k. The shape k and shift xi are (auxiliary) parameters that are estimated. This is the version of "zero_one_censored_shifted_gamma" (Sigrist and Stahel, 2011) without the upper censoring at 1

• "zero_censored_shifted_gamma_varying_shape": As "zero_censored_shifted_gamma", but the shape k varies across observations: log(k) = F_s(X) is related to fixed effects only (covariates and / or the GPBoost tree-boosting algorithm; no random effects / GPs for the shape), while log(mu) = F(X) + Zb is related to both fixed and random effects. The shift xi is then the only (auxiliary) parameter that is estimated. The estimated coefficients of the log-shape model are returned alongside the mean-model coefficients (with the suffix "_shape")

Count response: y in {0, 1, 2, ...}

• "poisson": Poisson likelihood with a log link function

• "negative_binomial": negative binomial likelihood with a log link function (aka "nbinom2", "negative_binomial_2"). The variance is mu * (mu + r) / r, mu = mean, r = shape, with this parametrization

• "negative_binomial_1": Negative binomial 1 (aka "nbinom1") likelihood with a log link function. The variance is mu * (1 + phi), mu = mean, phi = dispersion, with this parametrization

• "zero_inflated_<base>": Two-part count likelihoods that combine a point mass 'p0' at zero with a count base distribution for the remaining mass '1 - p0' (the base can itself generate additional zeros). As for the hurdle likelihoods, 'exp(F(X) + Zb)' is the mean of the (non-structural) base component - not the unconditional response mean - so E(y) = (1 - p0) * base_mean, and 'p0' is estimated jointly with the base auxiliary parameters. Currently supported bases: "zero_inflated_poisson", "zero_inflated_negative_binomial" (aux: shape; aka "zero_inflated_nbinom2"), "zero_inflated_negative_binomial_1" (aux: dispersion; aka "zero_inflated_nbinom1")

• "zero_inflated_regression_<base>": As "zero_inflated_<base>", but with the logistic structural-zero model described under "hurdle_regression_<base>" above: "zero_inflated_regression_poisson", "zero_inflated_regression_negative_binomial", and "zero_inflated_regression_negative_binomial_1"

Binary response: y in {0, 1}

• "bernoulli_logit": Bernoulli likelihood with a logit link function for binary classification. Aliases: "binary", "binary_logit"

• "bernoulli_probit": Bernoulli likelihood with a probit link function for binary classification. Aliases: "binary_probit"

Proportion / fractional / bounded response: y in [0, 1]

• "quasi_bernoulli_logit": quasi-Bernoulli likelihood with a logit link function for y in [0,1]. Aliases: "quasi_binary", "quasi_binary_logit"

• "quasi_bernoulli_probit": quasi-Bernoulli likelihood with a probit link function for y in [0,1]. Aliases: "quasi_binary_probit"

• "binomial_logit": Binomial likelihood with a logit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials. Aliases: "binomial"

• "binomial_probit": Binomial likelihood with a probit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials

• "beta_binomial": Beta-binomial likelihood with a logit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials. Aliases: "betabinomial", "beta-binomial"

• "beta" : Beta likelihood with a logit link function (parametrization of Ferrari and Cribari-Neto, 2004). The response must lie strictly in (0,1)

• "zoctn": Zero-one censored transformed normal likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). The model used is T ~ N(mu, sigma^2), W = max(min(T,1),0), and Y = g(W), where g(x) = expit(a + b * logit(x)) for x in (0,1), mu = F(X) + Z_RE u, u denotes the random effects, Z_RE is their design matrix, and sigma, a, and b are (auxiliary) parameters that are estimated. For more details on this model, see Qiang and Sigrist (2026)

• "zero_one_censored_transformed_beta": Zero-one censored transformed beta likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). If T follows a beta distribution with mean mu = expit(F(X) + Zb) and precision phi, the observed response is obtained by applying the linear transformation Y = (1 + 2u) * T - u and censoring the result to [0,1]. The precision phi and shift u are (auxiliary) parameters that are estimated. For more details on this model, see Kosmidis and Zeileis (2025)

• "zero_one_censored_shifted_gamma": Zero-one censored shifted gamma likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). The model used is Y = min(max(Z - xi, 0), 1), where Z follows a gamma distribution with mean mu = exp(F(X) + Zb) and shape k. The shape k and shift xi are (auxiliary) parameters that are estimated. For more details on this model, see Sigrist and Stahel (2011)

Notes

• The first lines in the likelihoods source file contain additional comments on the specific parametrizations used

• Other likelihoods can be implemented upon request

A numeric specifying an additional parameter for the likelihood which cannot be estimated for this likelihood (e.g., degrees of freedom for likelihood = "t_fix_df"). This is not to be confused with any auxiliary parameters that can be estimated and accessed through the function get_aux_pars after estimation. Note that this likelihood_additional_param parameter is irrelevant for many likelihoods. If likelihood_additional_param = NULL, the following default and mandatory-value rules apply:

• df = 2 for likelihood = "t_fix_df"

• No default is used for likelihood = "asymmetric_laplace"; a quantile strictly between 0 and 1 is required

• No default is used for the Tweedie likelihoods with a fixed power ("tweedie_fixed_p" and the other "_fixed_p" variants); a power strictly between 1.01 and 1.99 is required

A vector or matrix whose columns are categorical grouping variables. The elements being group levels defining grouped random effects. The elements of 'group_data' can be integer, double, or character. The number of columns corresponds to the number of grouped (intercept) random effects

A vector or matrix with numeric covariate data for grouped random coefficients

A vector with integer indices that indicate the corresponding categorical grouping variable (=columns) in 'group_data' for every covariate in 'group_rand_coef_data'. Counting starts at 1. The length of this index vector must equal the number of covariates in 'group_rand_coef_data'. For instance, c(1,1,2) means that the first two covariates (=first two columns) in 'group_rand_coef_data' have random coefficients corresponding to the first categorical grouping variable (=first column) in 'group_data', and the third covariate (=third column) in 'group_rand_coef_data' has a random coefficient corresponding to the second grouping variable (=second column) in 'group_data'

A vector of type logical (boolean). Indicates whether intercept random effects are dropped (only for random coefficients). If drop_intercept_group_rand_effect[k] is TRUE, the intercept random effect number k is dropped / not included. Only random effects with random slopes can be dropped.

A matrix with numeric coordinates (= inputs / features) for defining Gaussian processes

A vector or matrix with numeric covariate data for Gaussian process random coefficients

A string specifying the covariance function for the Gaussian process. Available options:

• "matern": Matern covariance function with the smoothness specified by the cov_fct_shape parameter (using the parametrization of Rasmussen and Williams, 2006)

• "matern_estimate_shape": same as "matern" but the smoothness parameter is also estimated

• "matern_space_time": Spatio-temporal Matern covariance function with different range parameters for space and time. Note that the first column in gp_coords must correspond to the time dimension

• "space_time_gneiting": Spatio-temporal covariance function given in Eq. (16) of Gneiting (2002). Note that the first column in gp_coords must correspond to the time dimension. This covariance has seven parameters (in the following order: sigma2, a, c, alpha, nu, beta, delta) which are all estimated by default. You can disable the estimation of some of these parameter using the 'estimate_cov_par_index' argument of the params argument in either the fit function of a gp_model object or the set_optim_params function prior to estimation.

• "matern_ard": anisotropic Matern covariance function with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of gp_coords

• "matern_ard_estimate_shape": same as "matern_ard" but the smoothness parameter is also estimated

• "exponential": Exponential covariance function (using the parametrization of Diggle and Ribeiro, 2007)

• "gaussian": Gaussian, aka squared exponential, covariance function (using the parametrization of Diggle and Ribeiro, 2007)

• "gaussian_ard": anisotropic Gaussian, aka squared exponential, covariance function with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of gp_coords

• "powered_exponential": powered exponential covariance function with the exponent specified by the cov_fct_shape parameter (using the parametrization of Diggle and Ribeiro, 2007)

• "wendland": Compactly supported Wendland covariance function (using the parametrization of Bevilacqua et al., 2019, AOS)

• "linear": linear covariance function. This corresponds to a Bayesian linear regression model with a Gaussian prior on the coefficients with a constant variance diagonal prior covariance, and the prior variance is estimated using empirical Bayes.

• "hurst": Hurst covariance function cov(s, s') = (sigma2 / 2) * ( ||s||^(2H) + ||s'||^(2H) - ||s - s'||^(2H) ), 0 < H < 1. For H = 0.5, this corresponds to Brownian motion in one dimension (-> see the 'estimate_cov_par_index' argument). This is the covariance for cov_fct_order = 1 (default). For cov_fct_order = 2, the second-order Hurst covariance cov(s, s') = sigma2 / (2 (2H - 1)) * ( ||s - s'||^(2H) - ||s||^(2H) - ||s'||^(2H) + 2H (s^T s') (||s||^(2H-2) + ||s'||^(2H-2)) ), 1 < H < 2, is used. In one dimension and for s, s' >= 0, H = 1.5 corresponds to an integrated Wiener process. See cov_fct_order for more details

• "hurst_ard": Hurst covariance function with with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of “gp_coords“ except for the first coordinate which has a range parameter of 1 due to identifiability with the marginal variance: cov(s, s') = (\sigma^2/2)\left[ \left(s_1^2 + \sum_{k=2}^d (s_k/l_k)^2\right)^H + \left({s'}_1^2 + \sum_{k=2}^d ({s'}_k/l_k)^2\right)^H - \left((s_1-{s'}_1)^2 + \sum_{k=2}^d ((s_k-{s'}_k)/l_k)^2\right)^H \right] . For cov_fct_order = 2, the second-order Hurst covariance (see "hurst") is applied to the scaled coordinates (s_1, s_2/l_2, ..., s_d/l_d)

• "ar1_mf_<base>": two-level autoregressive multifidelity covariance defined by f_H(x)=\rho f_L(x)+\delta(x), where f_L and \delta are independent Gaussian processes using the same base covariance type but separate parameter vectors.

  • For example, use "ar1_mf_matern", "ar1_mf_matern_ard", or "ar1_mf_matern_estimate_shape".

  • The last column of gp_coords must be 0 for low fidelity or 1 for high fidelity. All preceding columns are input coordinates for the GPs.

  • Covariance parameters are ordered as [low-fidelity base parameters, discrepancy base parameters, rho]. The two base blocks use the ordinary parameter ordering of the base covariance, and rho is unrestricted and can be negative.

  • Any supported base covariance except "wendland" can be used. Correlation tapering and Gaussian-process random coefficients are currently not supported for this model.

A logical. For an ar1_mf_<base> covariance, whether marginal mean models are fitted independently for low and high fidelity. If TRUE (the default), linear covariates supplied to fit() are internally expanded into separate low- and high-fidelity coefficient blocks, and the fidelity indicator is automatically added as a feature when the model is used in the GPBoost algorithm. Ignored for other covariance functions.

A numeric specifying the shape parameter of the covariance function (e.g., smoothness parameter for Matern and Wendland covariance) This parameter is irrelevant for some covariance functions such as the exponential or Gaussian

An integer specifying the order m of the "hurst" and "hurst_ard" covariance functions (also when used as base covariance in "ar1_mf_hurst" and "ar1_mf_hurst_ard"). Currently, the orders 1 (default) and 2 are supported. The Hurst exponent H of the order m satisfies m - 1 < H < m, and H = m - 0.5 is used as initial value. Order 1 is a fractional Brownian motion / field anchored at the origin (b(0) = 0), and order 2 is a second-order Hurst process / field anchored at the origin such that both the process and its gradient are zero there (b(0) = 0, grad b(0) = 0). The origin of the coordinates thus is a part of the model, and the coordinates are not centered internally. The order is not estimated. It can be chosen, e.g., by comparing the marginal likelihoods or by cross-validation.

• Continuous-time RW1 prior: use cov_function = "hurst", cov_fct_order = 1, and fix H = 0.5. This is a Brownian motion with b(0) = 0 and b'(t) = sqrt(q) W'(t), where W' denotes white noise and q = sigma2

• Continuous-time RW2 prior: use cov_function = "hurst", cov_fct_order = 2, and fix H = 1.5. This is an integrated Wiener process with b(0) = b'(0) = 0 and b”(t) = sqrt(q) W'(t), where q = 3 * sigma2

For these priors, use a one-dimensional time coordinate t >= 0 whose origin t = 0 is the desired anchor (shift the time points if necessary). H can be fixed using the 'init_cov_pars' and 'estimate_cov_par_index' arguments of the params argument, see the examples of fitGPModel. If the polynomial null space of the intrinsic model should not be penalized, include an intercept (RW1) or an intercept and a linear time trend (RW2) as fixed effects. Note that the order 2 model with H = 1.5 is the continuous-time RW2 prior evaluated at the time points and not the usual discrete intrinsic RW2 prior on a lattice: for equally spaced time points, adjacent second differences are correlated (correlation 1/4) whereas they are independent for the discrete RW2

A string specifying the large data approximation for Gaussian processes. Available options:

• "none": No approximation

• "vecchia": Vecchia approximation; see Sigrist (2022, JMLR) for more details. For "matern_space_time" and the anisotropic ARD covariance functions ("matern_ard", "gaussian_ard", "matern_ard_estimate_shape"), neighbors are selected by absolute correlation, i.e., using distances in coordinates scaled by the range parameters, and they are redetermined during parameter estimation. For "space_time_gneiting" and "ar1_mf_<base>", neighbors are also selected by absolute correlation by default. Use gp_approx = "vecchia_euclidean" for Euclidean-distance selection.

• "full_scale_vecchia": Vecchia-inducing points full-scale (VIF) approximation; see Gyger, Furrer, and Sigrist (2026, JMLR) for more details

• "tapering": The covariance function is multiplied by a compactly supported Wendland correlation function

• "fitc": Fully Independent Training Conditional approximation aka modified predictive process approximation

• "full_scale_tapering": Full-scale approximation combining an inducing point / predictive process approximation with tapering on the residual process; see Gyger, Furrer, and Sigrist (2026, JUQ) for more details

• "vecchia_latent": similar as "vecchia" but a Vecchia approximation is applied to the latent Gaussian process for likelihood == "gaussian". For likelihood != "gaussian", "vecchia" and "vecchia_latent" are equivalent

An integer specifying the number of parallel threads for OMP. If num_parallel_threads = NULL, the default number of threads of the session is used, see gpb.get.default.num.threads. Initially, this is the number of physical performance cores, which is determined automatically as follows. On CPUs whose cores have different speeds (e.g., the performance and efficiency cores of recent Intel CPUs or Apple silicon), the fastest cores are counted, and hyperthreads are counted only once, since slow threads can slow down an entire model. The next fastest cores are added if the fastest ones alone would leave only a single thread. On Linux, a CPU bandwidth limit of a control group (e.g., of a container) is respected as well. If the environment variable OMP_NUM_THREADS is set, it takes precedence over both the topology of the CPU and the limit of the control group, and the number of threads that OpenMP is configured to use is taken as it is if the cores of the CPU cannot be determined. In either case, the limits of the OpenMP runtime itself still apply, in particular the limit of the contention group (OMP_THREAD_LIMIT). For ordinary use, leave num_parallel_threads unspecified to use this default. The default of the session can be changed with gpb.tune.num.threads, which benchmarks different numbers of threads, or with gpb.set.default.num.threads. Setting num_parallel_threads = 1 disables the OpenMP parallelization of the model and can substantially increase the runtime. A single thread should be chosen deliberately, e.g., to distribute the resources among concurrent model fits or for a specific test; it is not required for correctness or reproducibility

A boolean. If TRUE, GPU acceleration will be used if supported

A numeric specifying the range parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)

A numeric specifying the shape (=smoothness) parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)

An integer specifying the number of neighbors for the Vecchia and VIF approximations. Internal default values if NULL:

• 20 for gp_approx = "vecchia"

• 30 for gp_approx = "full_scale_vecchia"

Note: for prediction, the number of neighbors can be set through the 'num_neighbors_pred' parameter in the 'set_prediction_data' function. By default, num_neighbors_pred = 2 * num_neighbors. Further, the type of Vecchia approximation used for making predictions is set through the 'vecchia_pred_type' parameter in the 'set_prediction_data' function

A string specifying the ordering used in the Vecchia approximation. Available options:

• "none": the default ordering in the data is used

• "random": a random ordering

• "time": ordering accorrding to time (only for space-time models)

• "time_random_space": ordering according to time and randomly for all spatial points with the same time points (only for space-time models)

A string specifying the method for choosing inducing points Available options:

• "kmeans++: the k-means++ algorithm

• "cover_tree": the cover tree algorithm

• "random": random selection from data points

An integer specifying the number of inducing points / knots for FITC, full_scale_tapering, and VIF approximations. Internal default values if NULL:

• 500 for gp_approx = "FITC" and gp_approx = "full_scale_tapering"

• 200 for gp_approx = "full_scale_vecchia"

A numeric specifying the radius (= "spatial resolution") for the cover tree algorithm

A string specifying the method used for inverting covariance matrices. Available options:

• "default": iterative methods where possible, otherwise Cholesky factorization. Cholesky factorization is also used for likelihood = "gaussian_heteroscedastic_fixed_and_random"

• "cholesky": Cholesky factorization

• "iterative": iterative methods. A combination of the conjugate gradient, the Lanczos algorithm, and other methods.

  This is currently only supported for the following cases:

  • grouped random effects with more than one level

  • likelihood != "gaussian" and gp_approx == "vecchia" (non-Gaussian likelihoods with a Vecchia-Laplace approximation)

  • likelihood != "gaussian" and gp_approx == "full_scale_vecchia" (non-Gaussian likelihoods with a VIF approximation)

  • likelihood == "gaussian" and gp_approx == "full_scale_tapering" (Gaussian likelihood with a full-scale tapering approximation)

An integer specifying the seed used for model creation (e.g., random ordering in Vecchia approximation)

A string specifying the type of Vecchia approximation used for making predictions. Default value if vecchia_pred_type = NULL: "order_obs_first_cond_obs_only". Available options:

• "order_obs_first_cond_obs_only": Vecchia approximation for the observable process and observed training data is ordered first and the neighbors are only observed training data points

• "order_obs_first_cond_all": Vecchia approximation for the observable process and observed training data is ordered first and the neighbors are selected among all points (training + prediction)

• "latent_order_obs_first_cond_obs_only": Vecchia approximation for the latent process and observed data is ordered first and neighbors are only observed points

• "latent_order_obs_first_cond_all": Vecchia approximation for the latent process and observed data is ordered first and neighbors are selected among all points

• "order_pred_first": Vecchia approximation for the observable process and prediction data is ordered first for making predictions. This option is only available for Gaussian likelihoods

an integer specifying the number of neighbors for the Vecchia approximation for making predictions. Default value if NULL: num_neighbors_pred = 2 * num_neighbors

a numeric specifying the tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithms when being used for prediction. This is the tolerance of the "absolute" criterion only, the "relative" criterion is controlled by cg_rel_tol_pred and cg_abs_tol_pred. Default value if NULL: 1e-3

a string specifying the stopping rule of conjugate gradient algorithms when being used for prediction: "absolute" (stop when ||r||_2 < cg_delta_conv_pred) or "relative" (stop when ||r||_2 <= max(cg_abs_tol_pred, cg_rel_tol_pred * ||b||_2)). Default value if NULL: the value of cg_convergence_criterion used for the parameter estimation

a numeric specifying the relative tolerance of the "relative" stopping rule when being used for prediction. Default value if NULL: the value of cg_rel_tol used for the parameter estimation. This is inherited independently of cg_convergence_criterion_pred

a numeric specifying the absolute tolerance floor of the "relative" stopping rule when being used for prediction. This makes the rule robust for right-hand sides that are zero or very small. Default value if NULL: the value of cg_abs_tol used for the parameter estimation. This is inherited independently of the other two prediction options

an integer specifying the number of samples when simulation is used for calculating predictive variances Internal default values if NULL:

• 500 for grouped random effects

• 1000 for gp_approx = "vecchia" and gp_approx = "full_scale_tapering"

• 100 for gp_approx = "full_scale_vecchia"

an integer specifying the rank of the matrix for approximating predictive covariances obtained using the Lanczos algorithm Default value if NULL: 1000

A vector with elements indicating independent realizations of random effects / Gaussian processes (same values = same process realization). The elements of 'cluster_ids' can be integer, double, or character.

A vector with sample weights. For a Gaussian likelihood, the error variance ("nugget") for observation i is divided by weights[i]. For non-Gaussian likelihoods, the conditional log-likelihood contribution of observation i is multiplied by weights[i]. Consequently, weights affect the estimation of both random and fixed effects. Note that a Gaussian likelihood is calculated via the Laplace approximation when gp_approx = "vecchia_latent", when likelihood = "gaussian_latent", and when grouped random effects are combined with a Vecchia-approximated Gaussian process. In these cases, the weights act as for a non-Gaussian likelihood, i.e., the Gaussian log-likelihood contribution of observation i is multiplied by weights[i] instead of the error variance being divided by it.

A vector or matrix with observation-level data that some likelihoods require in addition to the response variable y (one row per data point). Currently, this is only used by the joint Tweedie likelihoods ("tweedie_joint" and its variants), for which its first column contains the observed number of events (e.g., claims) N. N must be a non-negative integer that is 0 if and only if y is 0

A numeric with a learning rate for the likelihood for generalized Bayesian inference (only non-Gaussian likelihoods)

A boolean. If TRUE, the data (groups, coordinates, covariate data for random coefficients, cluster IDs, weights, and additional likelihood data) is freed in R after initialization

A vector with response variable data

A matrix with numeric covariate data for the fixed effects linear regression term (if there is one)

A list with parameters for the estimation / optimization

• trace: boolean (default = FALSE). If TRUE, information on the progress of the parameter optimization is printed

• init_cov_pars: vector with numeric elements (default = NULL). Initial values for covariance parameters of Gaussian process and random effects (can be NULL). The order is same as the order of the parameters in the summary function: first is the error variance (only for "gaussian" likelihood), next follow the variances of the grouped random effects (if there are any, in the order provided in 'group_data'), and then follow the marginal variance and the ranges of the Gaussian process. If there are multiple Gaussian processes, then the variances and ranges follow alternatingly. If 'init_cov_pars = NULL', an internal choice is used that depends on the likelihood and the random effects type and covariance function. If you select the option 'trace = TRUE' in the 'params' argument, you will see the first initial covariance parameters in iteration 0.

• init_coef: vector with numeric elements (default = NULL). Initial values for the regression coefficients (if there are any, can be NULL)

• init_aux_pars: vector with numeric elements (default = NULL). Initial values for additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative_binomial likelihood)

• init_coef_aux_pars_from_iid_model: boolean (default = TRUE). If TRUE, regression coefficients and auxiliary parameters are initialized from an iid model (only for models with a linear regression term). This option is ignored if init_coef is provided. If init_aux_pars is provided but init_coef is not, only regression coefficients are initialized from an iid model. For likelihood = "gaussian_heteroscedastic_fixed_and_random", if init_cov_pars is not provided, the initial covariance parameters are obtained from a homoscedastic Gaussian process model, and the regression coefficients are then initialized from this model instead: the coefficients of the mean from its coefficients and the intercept of the log-error variance from the logarithm of its error variance.

• estimate_cov_par_index: vector with integer (default = -1). This allows for disabling the estimation of some (or all) covariance parameters. If 'estimate_cov_par_index' = -1, all covariance parameters are estimated. If estimate_cov_par_index != -1, this should be a vector with length equal to the number of covariance parameters, and estimate_cov_par_index[i] should be of bool type indicating whether parameter number i is estimated or not. For instance, estimate_cov_par_index = c(1,1,0) means that the first two covariance parameters are estimated and the last one not. Parameters that are not estimated are kept at their initial values (see 'init_cov_pars'). They are treated as known constants when calculating the standard errors of the estimated parameters, and their own standard errors are NaN.

• estimate_aux_pars: boolean (default = TRUE). If TRUE, additional parameters for non-Gaussian likelihoods are also estimated (e.g., shape parameter of a gamma or negative_binomial likelihood)

• optimizer_cov: string (default = "lbfgs"). Optimizer used for estimating covariance parameters. Options: "lbfgs", "gradient_descent", "fisher_scoring", "newton", "nelder_mead". If there are additional auxiliary parameters for non-Gaussian likelihoods, 'optimizer_cov' is also used for those

• optimizer_coef: string (default = "wls" for Gaussian likelihoods and "lbfgs" for other likelihoods). Optimizer used for estimating linear regression coefficients, if there are any (for the GPBoost algorithm there are usually none). Options: "gradient_descent", "lbfgs", "wls", "nelder_mead". Gradient descent steps are done simultaneously with gradient descent steps for the covariance parameters. "wls" refers to doing coordinate descent for the regression coefficients using weighted least squares. If 'optimizer_cov' is set to "nelder_mead" or "lbfgs", 'optimizer_coef' is automatically also set to the same value.

• maxit: integer (default = 1000). Maximal number of iterations for optimization algorithm. If maxit = -999, internal default values are used

• delta_rel_conv: numeric (default = 1E-6 except for "nelder_mead" for which the default is 1E-8). Convergence tolerance. The algorithm stops if the relative change in either the (approximate) log-likelihood or the parameters is below this value. If delta_rel_conv = -999, internal default values are used

• cg_max_num_it: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithms. If cg_max_num_it = -999, internal default values are used

• cg_max_num_it_tridiag: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithm when being run as Lanczos algorithm for tridiagonalization. If cg_max_num_it_tridiag = -999, internal default values are used

• cg_delta_conv: numeric (default = 1E-2). Tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithm when being used for parameter estimation. This is the tolerance of the "absolute" criterion only, the "relative" criterion is controlled by cg_rel_tol and cg_abs_tol. If cg_delta_conv = -999, internal default values are used

• cg_convergence_criterion: string (default = "absolute"). Stopping rule of conjugate gradient algorithms.

  • "absolute": stop when ||r||_2 < cg_delta_conv

  • "relative": stop when ||r||_2 <= max(cg_abs_tol, cg_rel_tol * ||b||_2), where b is the right-hand side of the linear system. The cg_abs_tol floor makes this robust for right-hand sides that are zero or very small

• cg_rel_tol: numeric (default = 1E-6). Relative tolerance of the "relative" stopping rule. If cg_rel_tol = -999, internal default values are used

• cg_abs_tol: numeric (default = 1E-8). Absolute tolerance floor of the "relative" stopping rule, which makes it robust for right-hand sides that are zero or very small. This is independent of cg_delta_conv, so that the two criteria can be varied separately. If cg_abs_tol = -999, internal default values are used

• cg_multi_rhs_convergence: string (default = "average"). How the stopping rule is aggregated over the columns of a linear system with several right-hand sides, such as the stochastic Lanczos quadrature. For cg_convergence_criterion = "relative", every right-hand side is normalized by its own norm, i.e. q_j = ||r_j||_2 / max(cg_abs_tol, cg_rel_tol * ||b_j||_2):

  • "average": stop when mean(q_j) <= 1

  • "max": stop when max(q_j) <= 1

  • "per_rhs": stop iterating on column j individually as soon as q_j <= 1

  For cg_convergence_criterion = "absolute" there is no per-column normalization:

  • "average": stop when mean(||r_j||_2) < cg_delta_conv. This is the historic behavior and the default

  • "max": stop when ||r_j||_2 < cg_delta_conv for all j

  • "per_rhs": stop iterating on column j as soon as ||r_j||_2 < cg_delta_conv

• num_rand_vec_trace: integer (default = 50). Number of random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix. If num_rand_vec_trace = -999, internal default values are used

• reuse_rand_vec_trace: boolean (default = TRUE). If true, random vectors (e.g., Rademacher) for stochastic approximations of the trace of a matrix are sampled only once at the beginning of the parameter estimation and reused in later trace approximations. Otherwise they are sampled every time a trace is calculated

• seed_rand_vec_trace: integer (default = 1). Seed number to generate random vectors (e.g., Rademacher)

• cg_preconditioner_type (string): Type of preconditioner used for conjugate gradient algorithms.

  • Options for grouped random effects:

    • "ssor" (= default): SSOR preconditioner

    • "incomplete_cholesky": zero fill-in incomplete Cholesky factorization

  • Options for likelihood != "gaussian" and gp_approx == "vecchia" or likelihood == "gaussian" and gp_approx == "vecchia_latent":

    • "vadu" (= default): (B^T * (D^-1 + W) * B) as preconditioner for inverting (B^T * D^-1 * B + W), where B^T * D^-1 * B approx= Sigma^-1

    • "fitc": FITC / modified predictive process preconditioner for inverting (B^-1 * D * B^-T + W^-1)

    • "pivoted_cholesky": (Lk * Lk^T + W^-1) as preconditioner for inverting (B^-1 * D * B^-T + W^-1), where Lk is a low-rank pivoted Cholesky approximation for Sigma and B^-1 * D * B^-T approx= Sigma

    • "incomplete_cholesky": zero fill-in incomplete (reverse) Cholesky factorization of (B^T * D^-1 * B + W) using the sparsity pattern of B^T * D^-1 * B approx= Sigma^-1

  • Options for likelihood != "gaussian" and gp_approx == "full_scale_vecchia":

    • "fitc" ( = default): FITC / modified predictive process preconditioner

    • "vifdu": VIF with diagonal update preconditioner

  • Options for likelihood == "gaussian" and gp_approx == "full_scale_tapering":

    • "fitc" (= default): modified predictive process preconditioner

    • "none": no preconditioner

• fitc_piv_chol_preconditioner_rank (integer ): Rank of the FITC and pivoted Cholesky decomposition preconditioners for iterative methods for Vecchia and VIF approximations (for full_scale_tapering, the same inducing points as in the approximation as used). If fitc_piv_chol_preconditioner_rank = -999, internal default values are used:

  • 200 for the FITC preconditioner

  • 50 for the pivoted Cholesky decomposition preconditioner

• convergence_criterion: string (default = "relative_change_in_log_likelihood", only relevant for "gradient_descent", "fisher_scoring", and "newton"). The convergence criterion used for terminating the optimization algorithm. Options: "relative_change_in_log_likelihood" or "relative_change_in_parameters". If convergence_criterion = "default", internal default values are used

• lr_cov: numeric (default = 0.1 for "gradient_descent" and 1. otherwise, relevant for "gradient_descent", "fisher_scoring", "newton", and "lbfgs"). Initial learning rate for covariance parameters if a gradient-based optimization method is used

  • If lr_cov = -999, internal default values are used (0.1 for "gradient_descent" and 1. otherwise)

  • If there are additional auxiliary parameters for non-Gaussian likelihoods, 'lr_cov' is also used for those

  • For "gradient_descent", "fisher_scoring", and "newton", this is the learning rate by which the search direction is multiplied. It can be decreased, and increased again, during the optimization

  • For "lbfgs" and "lbfgs_linesearch_nocedal_wright", this is not a learning rate by which every update is multiplied. It only determines the first trial step, which is the negative gradient divided by its norm and multiplied by 'lr_cov', i.e., a step of length 'lr_cov'. The subsequent steps are determined by the approximate Hessian of the L-BFGS algorithm and the line search. 'lr_cov' has no effect when the approximate Hessian of a previous call is reused (see 'reuse_learning_rates_gp_model' of the GPBoost algorithm)

• lr_coef: numeric (default = 0.1, only relevant for "gradient_descent", "fisher_scoring", and "newton"). Learning rate for fixed effect regression coefficients if gradient descent is used. If lr_coef = -999, internal default values are used

• use_nesterov_acc: boolean (default = TRUE, only relevant for "gradient_descent"). If TRUE Nesterov acceleration is used. This is used only for gradient descent

• nesterov_schedule_version: integer (Default = 0, only relevant for "gradient_descent"). Which version of Nesterov schedule should be used. If nesterov_schedule_version = -999, internal default values are used.

• acc_rate_coef: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for regression coefficients (if there are any) for Nesterov acceleration. If acc_rate_coef = -999, internal default values are used

• acc_rate_cov: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for covariance parameters for Nesterov acceleration. If acc_rate_cov = -999, internal default values are used

• momentum_offset: integer (Default = 2, only relevant for "gradient_descent"). Number of iterations for which no momentum is applied in the beginning. If momentum_offset = -999, internal default values are used.

• m_lbfgs: integer (Default = 6). Number of corrections to approximate the inverse Hessian matrix for the "lbfgs" optimizer. If m_lbfgs = -999, internal default values are used

• cold_restart_lbfgs: boolean (Default = TRUE). If TRUE, restarts of the "lbfgs" optimizers are "cold" restarts: the regression coefficients, the auxiliary parameters, and the modes of the Laplace approximations are reset to their initial values and only the covariance parameters are kept. Such restarts can escape a local optimum in which the optimizer is stuck. If FALSE, "warm" restarts are done: the optimization simply continues from the current parameters with a re-initialized approximate Hessian. Only relevant if max_num_restarts_lbfgs > 0. Note that every cold restart is a full optimization with up to maxit iterations, and the number of iterations reported is the sum over all restarts

• max_num_restarts_lbfgs: integer (Default = 0). Maximal number of restarts of the "lbfgs" optimizers after they have terminated. A restart continues from the current parameters with a re-initialized approximate Hessian. What else is reset depends on cold_restart_lbfgs. Restarts are only done as long as the negative log-likelihood still decreases by more than delta_rel_conv and as long as the total number of iterations is below maxit. This can help when the line search of "lbfgs" fails and the optimizer thus terminates prematurely without having converged. This can happen, e.g., for the "asymmetric_laplace" (= quantile regression) likelihood, whose approximate marginal likelihood is not smooth

• reset_modes_line_search_lbfgs: boolean (Default = FALSE). If TRUE, the mode finding of the Laplace approximation starts at the modes of the current parameters for every trial point of the line search of the "lbfgs" optimizer, and not at the modes of the previous trial point. Only relevant for non-Gaussian likelihoods and optimizer_cov = "lbfgs" or "lbfgs_linesearch_nocedal_wright"

• delta_conv_mode_finding: numeric (Default = 1E-8, and 1E-10 for the likelihood "gaussian_heteroscedastic_fixed_and_random"). Convergence tolerance in mode finding algorithm for Laplace approximation for non-Gaussian likelihoods. If delta_conv_mode_finding = -999, internal default values are used

A numeric vector with additional fixed effects contributions that are added to the linear predictor (= offset). The length of this vector needs to equal the number of training data points times the number of fixed-effect sets.

This is discontinued. Use the renamed equivalent argument offset instead

A vector or matrix with elements being group levels for which predictions are made (if there are grouped random effects in the GPModel)

A vector or matrix with covariate data for grouped random coefficients (if there are some in the GPModel)

A matrix with prediction coordinates (=features) for Gaussian process (if there is a GP in the GPModel)

A vector or matrix with covariate data for Gaussian process random coefficients (if there are some in the GPModel)

A vector with elements indicating the realizations of random effects / Gaussian processes for which predictions are made (set to NULL if you have not specified this when creating the GPModel)

A matrix with prediction covariate data for the fixed effects linear regression term (if there is one in the GPModel)

A boolean. If TRUE, the (posterior) predictive covariance is calculated in addition to the (posterior) predictive mean

A boolean. If TRUE, the (posterior) predictive variances are calculated

A boolean. If TRUE, samples from the posterior are drawn

A boolean. If TRUE, samples from the prior are drawn

A numeric with the number of posterior samples to draw if 'sample_posterior=TRUE'

A numeric with the number of prior samples to draw if 'sample_prior=TRUE'

A boolean. If TRUE, (approximate) standard errors are calculated (= square root of diagonal of the inverse Fisher information for Gaussian likelihoods and square root of diagonal of a numerically approximated inverse Hessian for non-Gaussian likelihoods). For non-Gaussian likelihoods, no standard errors of the covariance and auxiliary parameters are calculated if the model is used in the GPBoost algorithm, since they require the fixed effects given by the tree ensemble

Discontinued. Use the argument gp_approx instead

A numeric with the number of samples. This is only used for iid models

Object of class gpb.Dataset

other parameters

object of class gpb.Dataset

a list of two elements: the first one is ignored and the second one is column names

a GPModel

A vector with response variable data

A matrix with numeric covariate data for the fixed effects linear regression term (if there is one)

A list with parameters for the estimation / optimization

• trace: boolean (default = FALSE). If TRUE, information on the progress of the parameter optimization is printed

• init_cov_pars: vector with numeric elements (default = NULL). Initial values for covariance parameters of Gaussian process and random effects (can be NULL). The order is same as the order of the parameters in the summary function: first is the error variance (only for "gaussian" likelihood), next follow the variances of the grouped random effects (if there are any, in the order provided in 'group_data'), and then follow the marginal variance and the ranges of the Gaussian process. If there are multiple Gaussian processes, then the variances and ranges follow alternatingly. If 'init_cov_pars = NULL', an internal choice is used that depends on the likelihood and the random effects type and covariance function. If you select the option 'trace = TRUE' in the 'params' argument, you will see the first initial covariance parameters in iteration 0.

• init_coef: vector with numeric elements (default = NULL). Initial values for the regression coefficients (if there are any, can be NULL)

• init_aux_pars: vector with numeric elements (default = NULL). Initial values for additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative_binomial likelihood)

• init_coef_aux_pars_from_iid_model: boolean (default = TRUE). If TRUE, regression coefficients and auxiliary parameters are initialized from an iid model (only for models with a linear regression term). This option is ignored if init_coef is provided. If init_aux_pars is provided but init_coef is not, only regression coefficients are initialized from an iid model. For likelihood = "gaussian_heteroscedastic_fixed_and_random", if init_cov_pars is not provided, the initial covariance parameters are obtained from a homoscedastic Gaussian process model, and the regression coefficients are then initialized from this model instead: the coefficients of the mean from its coefficients and the intercept of the log-error variance from the logarithm of its error variance.

• estimate_cov_par_index: vector with integer (default = -1). This allows for disabling the estimation of some (or all) covariance parameters. If 'estimate_cov_par_index' = -1, all covariance parameters are estimated. If estimate_cov_par_index != -1, this should be a vector with length equal to the number of covariance parameters, and estimate_cov_par_index[i] should be of bool type indicating whether parameter number i is estimated or not. For instance, estimate_cov_par_index = c(1,1,0) means that the first two covariance parameters are estimated and the last one not. Parameters that are not estimated are kept at their initial values (see 'init_cov_pars'). They are treated as known constants when calculating the standard errors of the estimated parameters, and their own standard errors are NaN.

• estimate_aux_pars: boolean (default = TRUE). If TRUE, additional parameters for non-Gaussian likelihoods are also estimated (e.g., shape parameter of a gamma or negative_binomial likelihood)

• optimizer_cov: string (default = "lbfgs"). Optimizer used for estimating covariance parameters. Options: "lbfgs", "gradient_descent", "fisher_scoring", "newton", "nelder_mead". If there are additional auxiliary parameters for non-Gaussian likelihoods, 'optimizer_cov' is also used for those

• optimizer_coef: string (default = "wls" for Gaussian likelihoods and "lbfgs" for other likelihoods). Optimizer used for estimating linear regression coefficients, if there are any (for the GPBoost algorithm there are usually none). Options: "gradient_descent", "lbfgs", "wls", "nelder_mead". Gradient descent steps are done simultaneously with gradient descent steps for the covariance parameters. "wls" refers to doing coordinate descent for the regression coefficients using weighted least squares. If 'optimizer_cov' is set to "nelder_mead" or "lbfgs", 'optimizer_coef' is automatically also set to the same value.

• maxit: integer (default = 1000). Maximal number of iterations for optimization algorithm. If maxit = -999, internal default values are used

• delta_rel_conv: numeric (default = 1E-6 except for "nelder_mead" for which the default is 1E-8). Convergence tolerance. The algorithm stops if the relative change in either the (approximate) log-likelihood or the parameters is below this value. If delta_rel_conv = -999, internal default values are used

• cg_max_num_it: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithms. If cg_max_num_it = -999, internal default values are used

• cg_max_num_it_tridiag: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithm when being run as Lanczos algorithm for tridiagonalization. If cg_max_num_it_tridiag = -999, internal default values are used

• cg_delta_conv: numeric (default = 1E-2). Tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithm when being used for parameter estimation. This is the tolerance of the "absolute" criterion only, the "relative" criterion is controlled by cg_rel_tol and cg_abs_tol. If cg_delta_conv = -999, internal default values are used

• cg_convergence_criterion: string (default = "absolute"). Stopping rule of conjugate gradient algorithms.

  • "absolute": stop when ||r||_2 < cg_delta_conv

  • "relative": stop when ||r||_2 <= max(cg_abs_tol, cg_rel_tol * ||b||_2), where b is the right-hand side of the linear system. The cg_abs_tol floor makes this robust for right-hand sides that are zero or very small

• cg_rel_tol: numeric (default = 1E-6). Relative tolerance of the "relative" stopping rule. If cg_rel_tol = -999, internal default values are used

• cg_abs_tol: numeric (default = 1E-8). Absolute tolerance floor of the "relative" stopping rule, which makes it robust for right-hand sides that are zero or very small. This is independent of cg_delta_conv, so that the two criteria can be varied separately. If cg_abs_tol = -999, internal default values are used

• cg_multi_rhs_convergence: string (default = "average"). How the stopping rule is aggregated over the columns of a linear system with several right-hand sides, such as the stochastic Lanczos quadrature. For cg_convergence_criterion = "relative", every right-hand side is normalized by its own norm, i.e. q_j = ||r_j||_2 / max(cg_abs_tol, cg_rel_tol * ||b_j||_2):

  • "average": stop when mean(q_j) <= 1

  • "max": stop when max(q_j) <= 1

  • "per_rhs": stop iterating on column j individually as soon as q_j <= 1

  For cg_convergence_criterion = "absolute" there is no per-column normalization:

  • "average": stop when mean(||r_j||_2) < cg_delta_conv. This is the historic behavior and the default

  • "max": stop when ||r_j||_2 < cg_delta_conv for all j

  • "per_rhs": stop iterating on column j as soon as ||r_j||_2 < cg_delta_conv

• num_rand_vec_trace: integer (default = 50). Number of random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix. If num_rand_vec_trace = -999, internal default values are used

• reuse_rand_vec_trace: boolean (default = TRUE). If true, random vectors (e.g., Rademacher) for stochastic approximations of the trace of a matrix are sampled only once at the beginning of the parameter estimation and reused in later trace approximations. Otherwise they are sampled every time a trace is calculated

• seed_rand_vec_trace: integer (default = 1). Seed number to generate random vectors (e.g., Rademacher)

• cg_preconditioner_type (string): Type of preconditioner used for conjugate gradient algorithms.

  • Options for grouped random effects:

    • "ssor" (= default): SSOR preconditioner

    • "incomplete_cholesky": zero fill-in incomplete Cholesky factorization

  • Options for likelihood != "gaussian" and gp_approx == "vecchia" or likelihood == "gaussian" and gp_approx == "vecchia_latent":

    • "vadu" (= default): (B^T * (D^-1 + W) * B) as preconditioner for inverting (B^T * D^-1 * B + W), where B^T * D^-1 * B approx= Sigma^-1

    • "fitc": FITC / modified predictive process preconditioner for inverting (B^-1 * D * B^-T + W^-1)

    • "pivoted_cholesky": (Lk * Lk^T + W^-1) as preconditioner for inverting (B^-1 * D * B^-T + W^-1), where Lk is a low-rank pivoted Cholesky approximation for Sigma and B^-1 * D * B^-T approx= Sigma

    • "incomplete_cholesky": zero fill-in incomplete (reverse) Cholesky factorization of (B^T * D^-1 * B + W) using the sparsity pattern of B^T * D^-1 * B approx= Sigma^-1

  • Options for likelihood != "gaussian" and gp_approx == "full_scale_vecchia":

    • "fitc" ( = default): FITC / modified predictive process preconditioner

    • "vifdu": VIF with diagonal update preconditioner

  • Options for likelihood == "gaussian" and gp_approx == "full_scale_tapering":

    • "fitc" (= default): modified predictive process preconditioner

    • "none": no preconditioner

• fitc_piv_chol_preconditioner_rank (integer ): Rank of the FITC and pivoted Cholesky decomposition preconditioners for iterative methods for Vecchia and VIF approximations (for full_scale_tapering, the same inducing points as in the approximation as used). If fitc_piv_chol_preconditioner_rank = -999, internal default values are used:

  • 200 for the FITC preconditioner

  • 50 for the pivoted Cholesky decomposition preconditioner

• convergence_criterion: string (default = "relative_change_in_log_likelihood", only relevant for "gradient_descent", "fisher_scoring", and "newton"). The convergence criterion used for terminating the optimization algorithm. Options: "relative_change_in_log_likelihood" or "relative_change_in_parameters". If convergence_criterion = "default", internal default values are used

• lr_cov: numeric (default = 0.1 for "gradient_descent" and 1. otherwise, relevant for "gradient_descent", "fisher_scoring", "newton", and "lbfgs"). Initial learning rate for covariance parameters if a gradient-based optimization method is used

  • If lr_cov = -999, internal default values are used (0.1 for "gradient_descent" and 1. otherwise)

  • If there are additional auxiliary parameters for non-Gaussian likelihoods, 'lr_cov' is also used for those

  • For "gradient_descent", "fisher_scoring", and "newton", this is the learning rate by which the search direction is multiplied. It can be decreased, and increased again, during the optimization

  • For "lbfgs" and "lbfgs_linesearch_nocedal_wright", this is not a learning rate by which every update is multiplied. It only determines the first trial step, which is the negative gradient divided by its norm and multiplied by 'lr_cov', i.e., a step of length 'lr_cov'. The subsequent steps are determined by the approximate Hessian of the L-BFGS algorithm and the line search. 'lr_cov' has no effect when the approximate Hessian of a previous call is reused (see 'reuse_learning_rates_gp_model' of the GPBoost algorithm)

• lr_coef: numeric (default = 0.1, only relevant for "gradient_descent", "fisher_scoring", and "newton"). Learning rate for fixed effect regression coefficients if gradient descent is used. If lr_coef = -999, internal default values are used

• use_nesterov_acc: boolean (default = TRUE, only relevant for "gradient_descent"). If TRUE Nesterov acceleration is used. This is used only for gradient descent

• nesterov_schedule_version: integer (Default = 0, only relevant for "gradient_descent"). Which version of Nesterov schedule should be used. If nesterov_schedule_version = -999, internal default values are used.

• acc_rate_coef: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for regression coefficients (if there are any) for Nesterov acceleration. If acc_rate_coef = -999, internal default values are used

• acc_rate_cov: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for covariance parameters for Nesterov acceleration. If acc_rate_cov = -999, internal default values are used

• momentum_offset: integer (Default = 2, only relevant for "gradient_descent"). Number of iterations for which no momentum is applied in the beginning. If momentum_offset = -999, internal default values are used.

• m_lbfgs: integer (Default = 6). Number of corrections to approximate the inverse Hessian matrix for the "lbfgs" optimizer. If m_lbfgs = -999, internal default values are used

• cold_restart_lbfgs: boolean (Default = TRUE). If TRUE, restarts of the "lbfgs" optimizers are "cold" restarts: the regression coefficients, the auxiliary parameters, and the modes of the Laplace approximations are reset to their initial values and only the covariance parameters are kept. Such restarts can escape a local optimum in which the optimizer is stuck. If FALSE, "warm" restarts are done: the optimization simply continues from the current parameters with a re-initialized approximate Hessian. Only relevant if max_num_restarts_lbfgs > 0. Note that every cold restart is a full optimization with up to maxit iterations, and the number of iterations reported is the sum over all restarts

• max_num_restarts_lbfgs: integer (Default = 0). Maximal number of restarts of the "lbfgs" optimizers after they have terminated. A restart continues from the current parameters with a re-initialized approximate Hessian. What else is reset depends on cold_restart_lbfgs. Restarts are only done as long as the negative log-likelihood still decreases by more than delta_rel_conv and as long as the total number of iterations is below maxit. This can help when the line search of "lbfgs" fails and the optimizer thus terminates prematurely without having converged. This can happen, e.g., for the "asymmetric_laplace" (= quantile regression) likelihood, whose approximate marginal likelihood is not smooth

• reset_modes_line_search_lbfgs: boolean (Default = FALSE). If TRUE, the mode finding of the Laplace approximation starts at the modes of the current parameters for every trial point of the line search of the "lbfgs" optimizer, and not at the modes of the previous trial point. Only relevant for non-Gaussian likelihoods and optimizer_cov = "lbfgs" or "lbfgs_linesearch_nocedal_wright"

• delta_conv_mode_finding: numeric (Default = 1E-8, and 1E-10 for the likelihood "gaussian_heteroscedastic_fixed_and_random"). Convergence tolerance in mode finding algorithm for Laplace approximation for non-Gaussian likelihoods. If delta_conv_mode_finding = -999, internal default values are used

A numeric vector with additional fixed effects contributions that are added to the linear predictor (= offset). The length of this vector needs to equal the number of training data points times the number of fixed-effect sets.

This is discontinued. Use the renamed equivalent argument offset instead

a GPModel

A vector with response variable data

A matrix with numeric covariate data for the fixed effects linear regression term (if there is one)

A list with parameters for the estimation / optimization

• trace: boolean (default = FALSE). If TRUE, information on the progress of the parameter optimization is printed

• init_cov_pars: vector with numeric elements (default = NULL). Initial values for covariance parameters of Gaussian process and random effects (can be NULL). The order is same as the order of the parameters in the summary function: first is the error variance (only for "gaussian" likelihood), next follow the variances of the grouped random effects (if there are any, in the order provided in 'group_data'), and then follow the marginal variance and the ranges of the Gaussian process. If there are multiple Gaussian processes, then the variances and ranges follow alternatingly. If 'init_cov_pars = NULL', an internal choice is used that depends on the likelihood and the random effects type and covariance function. If you select the option 'trace = TRUE' in the 'params' argument, you will see the first initial covariance parameters in iteration 0.

• init_coef: vector with numeric elements (default = NULL). Initial values for the regression coefficients (if there are any, can be NULL)

• init_aux_pars: vector with numeric elements (default = NULL). Initial values for additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative_binomial likelihood)

• init_coef_aux_pars_from_iid_model: boolean (default = TRUE). If TRUE, regression coefficients and auxiliary parameters are initialized from an iid model (only for models with a linear regression term). This option is ignored if init_coef is provided. If init_aux_pars is provided but init_coef is not, only regression coefficients are initialized from an iid model. For likelihood = "gaussian_heteroscedastic_fixed_and_random", if init_cov_pars is not provided, the initial covariance parameters are obtained from a homoscedastic Gaussian process model, and the regression coefficients are then initialized from this model instead: the coefficients of the mean from its coefficients and the intercept of the log-error variance from the logarithm of its error variance.

• estimate_cov_par_index: vector with integer (default = -1). This allows for disabling the estimation of some (or all) covariance parameters. If 'estimate_cov_par_index' = -1, all covariance parameters are estimated. If estimate_cov_par_index != -1, this should be a vector with length equal to the number of covariance parameters, and estimate_cov_par_index[i] should be of bool type indicating whether parameter number i is estimated or not. For instance, estimate_cov_par_index = c(1,1,0) means that the first two covariance parameters are estimated and the last one not. Parameters that are not estimated are kept at their initial values (see 'init_cov_pars'). They are treated as known constants when calculating the standard errors of the estimated parameters, and their own standard errors are NaN.

• estimate_aux_pars: boolean (default = TRUE). If TRUE, additional parameters for non-Gaussian likelihoods are also estimated (e.g., shape parameter of a gamma or negative_binomial likelihood)

• optimizer_cov: string (default = "lbfgs"). Optimizer used for estimating covariance parameters. Options: "lbfgs", "gradient_descent", "fisher_scoring", "newton", "nelder_mead". If there are additional auxiliary parameters for non-Gaussian likelihoods, 'optimizer_cov' is also used for those

• optimizer_coef: string (default = "wls" for Gaussian likelihoods and "lbfgs" for other likelihoods). Optimizer used for estimating linear regression coefficients, if there are any (for the GPBoost algorithm there are usually none). Options: "gradient_descent", "lbfgs", "wls", "nelder_mead". Gradient descent steps are done simultaneously with gradient descent steps for the covariance parameters. "wls" refers to doing coordinate descent for the regression coefficients using weighted least squares. If 'optimizer_cov' is set to "nelder_mead" or "lbfgs", 'optimizer_coef' is automatically also set to the same value.

• maxit: integer (default = 1000). Maximal number of iterations for optimization algorithm. If maxit = -999, internal default values are used

• delta_rel_conv: numeric (default = 1E-6 except for "nelder_mead" for which the default is 1E-8). Convergence tolerance. The algorithm stops if the relative change in either the (approximate) log-likelihood or the parameters is below this value. If delta_rel_conv = -999, internal default values are used

• cg_max_num_it: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithms. If cg_max_num_it = -999, internal default values are used

• cg_max_num_it_tridiag: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithm when being run as Lanczos algorithm for tridiagonalization. If cg_max_num_it_tridiag = -999, internal default values are used

• cg_delta_conv: numeric (default = 1E-2). Tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithm when being used for parameter estimation. This is the tolerance of the "absolute" criterion only, the "relative" criterion is controlled by cg_rel_tol and cg_abs_tol. If cg_delta_conv = -999, internal default values are used

• cg_convergence_criterion: string (default = "absolute"). Stopping rule of conjugate gradient algorithms.

  • "absolute": stop when ||r||_2 < cg_delta_conv

  • "relative": stop when ||r||_2 <= max(cg_abs_tol, cg_rel_tol * ||b||_2), where b is the right-hand side of the linear system. The cg_abs_tol floor makes this robust for right-hand sides that are zero or very small

• cg_rel_tol: numeric (default = 1E-6). Relative tolerance of the "relative" stopping rule. If cg_rel_tol = -999, internal default values are used

• cg_abs_tol: numeric (default = 1E-8). Absolute tolerance floor of the "relative" stopping rule, which makes it robust for right-hand sides that are zero or very small. This is independent of cg_delta_conv, so that the two criteria can be varied separately. If cg_abs_tol = -999, internal default values are used

• cg_multi_rhs_convergence: string (default = "average"). How the stopping rule is aggregated over the columns of a linear system with several right-hand sides, such as the stochastic Lanczos quadrature. For cg_convergence_criterion = "relative", every right-hand side is normalized by its own norm, i.e. q_j = ||r_j||_2 / max(cg_abs_tol, cg_rel_tol * ||b_j||_2):

  • "average": stop when mean(q_j) <= 1

  • "max": stop when max(q_j) <= 1

  • "per_rhs": stop iterating on column j individually as soon as q_j <= 1

  For cg_convergence_criterion = "absolute" there is no per-column normalization:

  • "average": stop when mean(||r_j||_2) < cg_delta_conv. This is the historic behavior and the default

  • "max": stop when ||r_j||_2 < cg_delta_conv for all j

  • "per_rhs": stop iterating on column j as soon as ||r_j||_2 < cg_delta_conv

• num_rand_vec_trace: integer (default = 50). Number of random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix. If num_rand_vec_trace = -999, internal default values are used

• reuse_rand_vec_trace: boolean (default = TRUE). If true, random vectors (e.g., Rademacher) for stochastic approximations of the trace of a matrix are sampled only once at the beginning of the parameter estimation and reused in later trace approximations. Otherwise they are sampled every time a trace is calculated

• seed_rand_vec_trace: integer (default = 1). Seed number to generate random vectors (e.g., Rademacher)

• cg_preconditioner_type (string): Type of preconditioner used for conjugate gradient algorithms.

  • Options for grouped random effects:

    • "ssor" (= default): SSOR preconditioner

    • "incomplete_cholesky": zero fill-in incomplete Cholesky factorization

  • Options for likelihood != "gaussian" and gp_approx == "vecchia" or likelihood == "gaussian" and gp_approx == "vecchia_latent":

    • "vadu" (= default): (B^T * (D^-1 + W) * B) as preconditioner for inverting (B^T * D^-1 * B + W), where B^T * D^-1 * B approx= Sigma^-1

    • "fitc": FITC / modified predictive process preconditioner for inverting (B^-1 * D * B^-T + W^-1)

    • "pivoted_cholesky": (Lk * Lk^T + W^-1) as preconditioner for inverting (B^-1 * D * B^-T + W^-1), where Lk is a low-rank pivoted Cholesky approximation for Sigma and B^-1 * D * B^-T approx= Sigma

    • "incomplete_cholesky": zero fill-in incomplete (reverse) Cholesky factorization of (B^T * D^-1 * B + W) using the sparsity pattern of B^T * D^-1 * B approx= Sigma^-1

  • Options for likelihood != "gaussian" and gp_approx == "full_scale_vecchia":

    • "fitc" ( = default): FITC / modified predictive process preconditioner

    • "vifdu": VIF with diagonal update preconditioner

  • Options for likelihood == "gaussian" and gp_approx == "full_scale_tapering":

    • "fitc" (= default): modified predictive process preconditioner

    • "none": no preconditioner

• fitc_piv_chol_preconditioner_rank (integer ): Rank of the FITC and pivoted Cholesky decomposition preconditioners for iterative methods for Vecchia and VIF approximations (for full_scale_tapering, the same inducing points as in the approximation as used). If fitc_piv_chol_preconditioner_rank = -999, internal default values are used:

  • 200 for the FITC preconditioner

  • 50 for the pivoted Cholesky decomposition preconditioner

• convergence_criterion: string (default = "relative_change_in_log_likelihood", only relevant for "gradient_descent", "fisher_scoring", and "newton"). The convergence criterion used for terminating the optimization algorithm. Options: "relative_change_in_log_likelihood" or "relative_change_in_parameters". If convergence_criterion = "default", internal default values are used

• lr_cov: numeric (default = 0.1 for "gradient_descent" and 1. otherwise, relevant for "gradient_descent", "fisher_scoring", "newton", and "lbfgs"). Initial learning rate for covariance parameters if a gradient-based optimization method is used

  • If lr_cov = -999, internal default values are used (0.1 for "gradient_descent" and 1. otherwise)

  • If there are additional auxiliary parameters for non-Gaussian likelihoods, 'lr_cov' is also used for those

  • For "gradient_descent", "fisher_scoring", and "newton", this is the learning rate by which the search direction is multiplied. It can be decreased, and increased again, during the optimization

  • For "lbfgs" and "lbfgs_linesearch_nocedal_wright", this is not a learning rate by which every update is multiplied. It only determines the first trial step, which is the negative gradient divided by its norm and multiplied by 'lr_cov', i.e., a step of length 'lr_cov'. The subsequent steps are determined by the approximate Hessian of the L-BFGS algorithm and the line search. 'lr_cov' has no effect when the approximate Hessian of a previous call is reused (see 'reuse_learning_rates_gp_model' of the GPBoost algorithm)

• lr_coef: numeric (default = 0.1, only relevant for "gradient_descent", "fisher_scoring", and "newton"). Learning rate for fixed effect regression coefficients if gradient descent is used. If lr_coef = -999, internal default values are used

• use_nesterov_acc: boolean (default = TRUE, only relevant for "gradient_descent"). If TRUE Nesterov acceleration is used. This is used only for gradient descent

• nesterov_schedule_version: integer (Default = 0, only relevant for "gradient_descent"). Which version of Nesterov schedule should be used. If nesterov_schedule_version = -999, internal default values are used.

• acc_rate_coef: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for regression coefficients (if there are any) for Nesterov acceleration. If acc_rate_coef = -999, internal default values are used

• acc_rate_cov: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for covariance parameters for Nesterov acceleration. If acc_rate_cov = -999, internal default values are used

• momentum_offset: integer (Default = 2, only relevant for "gradient_descent"). Number of iterations for which no momentum is applied in the beginning. If momentum_offset = -999, internal default values are used.

• m_lbfgs: integer (Default = 6). Number of corrections to approximate the inverse Hessian matrix for the "lbfgs" optimizer. If m_lbfgs = -999, internal default values are used

• cold_restart_lbfgs: boolean (Default = TRUE). If TRUE, restarts of the "lbfgs" optimizers are "cold" restarts: the regression coefficients, the auxiliary parameters, and the modes of the Laplace approximations are reset to their initial values and only the covariance parameters are kept. Such restarts can escape a local optimum in which the optimizer is stuck. If FALSE, "warm" restarts are done: the optimization simply continues from the current parameters with a re-initialized approximate Hessian. Only relevant if max_num_restarts_lbfgs > 0. Note that every cold restart is a full optimization with up to maxit iterations, and the number of iterations reported is the sum over all restarts

• max_num_restarts_lbfgs: integer (Default = 0). Maximal number of restarts of the "lbfgs" optimizers after they have terminated. A restart continues from the current parameters with a re-initialized approximate Hessian. What else is reset depends on cold_restart_lbfgs. Restarts are only done as long as the negative log-likelihood still decreases by more than delta_rel_conv and as long as the total number of iterations is below maxit. This can help when the line search of "lbfgs" fails and the optimizer thus terminates prematurely without having converged. This can happen, e.g., for the "asymmetric_laplace" (= quantile regression) likelihood, whose approximate marginal likelihood is not smooth

• reset_modes_line_search_lbfgs: boolean (Default = FALSE). If TRUE, the mode finding of the Laplace approximation starts at the modes of the current parameters for every trial point of the line search of the "lbfgs" optimizer, and not at the modes of the previous trial point. Only relevant for non-Gaussian likelihoods and optimizer_cov = "lbfgs" or "lbfgs_linesearch_nocedal_wright"

• delta_conv_mode_finding: numeric (Default = 1E-8, and 1E-10 for the likelihood "gaussian_heteroscedastic_fixed_and_random"). Convergence tolerance in mode finding algorithm for Laplace approximation for non-Gaussian likelihoods. If delta_conv_mode_finding = -999, internal default values are used

A numeric vector with additional fixed effects contributions that are added to the linear predictor (= offset). The length of this vector needs to equal the number of training data points times the number of fixed-effect sets.

This is discontinued. Use the renamed equivalent argument offset instead

A string specifying the likelihood function (distribution) of the response variable. The currently supported likelihoods, grouped by response type, are:

The table lists the main likelihood names. Aliases, heteroscedastic and varying-shape variants, and the "_regression" variants of the two-part likelihoods are documented in the detailed descriptions below.

Continuous response: y in (-inf, inf)

• "gaussian": Gaussian likelihood

• "t": t-distribution (e.g., for robust regression). The default approximation is Fisher-Laplace, i.e., the Fisher information is used instead of the observed Hessian in the Laplace approximation of the marginal likelihood.

• "t_fix_df": t-distribution with the degrees-of-freedom (df) held fixed and not estimated

  • The degrees-of-freedom (df) can be set via the likelihood_additional_param parameter. The default is df = 2

• "quantile_regression" / "asymmetric_laplace" : an asymmetric Laplace likelihood for quantile regression. Both names are accepted, as is the alias "quantile"

  • The quantile must be supplied through the likelihood_additional_param parameter and must be strictly between 0 and 1

  The default approximation is Fisher-Laplace, i.e., the Fisher information is used instead of the observed Hessian in the Laplace approximation of the marginal likelihood. The triangular-kernel-curvature (TKC) approximation can be enabled by appending "_triangular_kernel_curvature" or the shorthand "_tkc", for example, "quantile_regression_tkc" or "asymmetric_laplace_triangular_kernel_curvature".

  • experimental: the log-likelihood is not differentiable at the kinks where 'y_i' equals the location parameter, and the quasi-Newton mode finding can therefore stall at a point that is not the exact posterior mode. Appending "_ssn_alm", for example "quantile_regression_ssn_alm", enables an exact check of the (non-smooth) optimality conditions after the mode finding and, if the check fails, a refinement of the mode with a semismooth Newton method applied to the subproblems of an augmented Lagrangian method. This makes the approximate marginal likelihood a well-defined function of the parameters, i.e., independent of the path taken by the optimizer and of the number of threads, at the price of additional linear solves during mode finding. This matters for standard errors and for likelihood-based model comparison, but it does not necessarily improve predictive accuracy; increasing max_num_restarts_lbfgs can be a cheaper way of avoiding poor optima. Appending "_admm_ssn_alm" instead warm starts the refinement with an ADMM phase, which typically halves its cost when matrix_inversion_method = "cholesky" and is not recommended for iterative methods

• "gaussian_heteroscedastic": Gaussian likelihood where the mean is related to fixed and random effects and the log-error variance is related to fixed effects only (linear predictor or GPBoost algorithm). The estimated coefficients of the log-variance model are returned alongside the mean-model coefficients (with the suffix "_scale"). Fisher-Laplace is the default and currently the only implemented approximation.

• "gaussian_heteroscedastic_fixed_and_random": Gaussian likelihood where both the mean and the variance are related to fixed and random effects. This is currently only implemented for GPs with a 'vecchia' approximation. Fisher-Laplace is the default and currently the only implemented approximation. By default, Cholesky factorizations are used. If init_cov_pars is not given, the initial values are obtained from a homoscedastic Gaussian process model (unless some weights are zero).

Positive continuous response: y in (0, inf)

• "gamma": Gamma likelihood with a log link function

• "gamma_varying_shape": As "gamma", but the gamma shape varies across observations and is modeled by an additional fixed-effects-only predictor: log(shape) = F_s(X) (F_s(X) = linear predictor or the GPBoost algorithm), while the log mean log(mu) = F(X) + Zb is related to both fixed and random effects. Note that the shape also governs the dispersion, var(y) = mu^2 / shape. The estimated coefficients of the log-shape model are returned alongside the mean-model coefficients (with the suffix "_shape"). See also the corresponding two-part likelihoods "hurdle_gamma_varying_shape" and "hurdle_regression_gamma_varying_shape" below

• "lognormal": Log-normal likelihood with a log link function

• "gpd": Generalized Pareto likelihood for finite positive responses. The log scale parameter equals the latent predictor eta (sum of fixed and random effects), sigma = exp(eta), and the estimated auxiliary parameter is shape > -0.5

• "egpd_power": Naveau power-carrier extended generalized Pareto likelihood with auxiliary parameters shape and kappa

• "egpd_power_mixture": Naveau ordered power-mixture carrier with auxiliary parameters shape, kappa1, delta_kappa, and p

• "egpd_beta": Naveau beta-carrier extended generalized Pareto likelihood with auxiliary parameters shape and delta

• "egpd_power_beta": Naveau power-beta carrier with auxiliary parameters shape, delta, and kappa

Non-negative continuous / semicontinuous response: y in [0, inf) with point mass at 0

• "tweedie": Compound Poisson–Gamma Tweedie likelihood with a log link, variance phi * mu^p, and 1.01 < p < 1.99. The dispersion phi and power p are estimated

• "tweedie_fixed_p": The same Tweedie likelihood with p fixed through likelihood_additional_param; only phi is estimated. The fixed power must satisfy 1.01 < p < 1.99

• "tweedie_joint", "tweedie_joint_fixed_p": The same Tweedie model, but the observed number of events (e.g., claims) N, given in the first column of additional_likelihood_data, is used jointly with the aggregate response y. The likelihood is the joint density of (y, N) of the compound Poisson–Gamma representation, N ~ Poisson(mu^(2-p) / (phi * (2-p))) and y | N = n ~ Gamma(n * (2-p) / (p-1), scale = phi * (p-1) * mu^(p-1)), see Jorgensen and de Souza (1994, Scandinavian Actuarial Journal). Given phi and p, the mean model is the same as for "tweedie"; N adds information for estimating phi and p. N must be a non-negative integer that is 0 if and only if y is 0. Predictions are the same as for "tweedie", and the "test_neg_log_likelihood" metric uses the marginal density of y since N is not known for new data

• "tweedie_varying_dispersion", "tweedie_varying_dispersion_fixed_p": As "tweedie" and "tweedie_fixed_p", but the dispersion phi varies across observations: log(phi) = F_d(X) is related to fixed effects only (linear predictor or GPBoost algorithm), while log(mu) = F(X) + Zb is related to both fixed and random effects. The power p is then the only (auxiliary) parameter. The estimated coefficients of the log-dispersion model are returned alongside the mean-model coefficients (with the suffix "_dispersion")

• "tweedie_joint_varying_dispersion", "tweedie_joint_varying_dispersion_fixed_p": The joint (y, N) likelihood of "tweedie_joint" with the varying dispersion of "tweedie_varying_dispersion"

• "hurdle_<base>": Two-part likelihoods for non-negative response variables with an excess probability 'p0' of exact zeros. They combine a point mass 'p0' at zero with a base distribution with support 'y > 0' for the remaining mass '1 - p0'. The fixed effects 'F(X)' and random effects 'Zb' enter only through the base component: 'exp(F(X) + Zb)' is its mean or scale parameter, not the unconditional response mean, which is E(y) = (1 - p0) * base_mean. The structural-zero probability 'p0' is estimated jointly with the base auxiliary parameters. Currently supported bases: "hurdle_gamma", "hurdle_lognormal", and the extreme-value bases "hurdle_gpd", "hurdle_egpd_power", "hurdle_egpd_power_mixture", "hurdle_egpd_beta", "hurdle_egpd_power_beta". For all these bases, the "zero_inflated" prefix is accepted as an alias, e.g. "zero_inflated_gamma" maps to "hurdle_gamma". See "zero_inflated_<base>" under count responses for the corresponding count-data family

• "hurdle_regression_<base>": As "hurdle_<base>", but the structural-zero probability is modeled as a logistic regression on the covariates instead of as a constant (e.g. "hurdle_regression_gamma", "hurdle_regression_lognormal", "hurdle_regression_gpd", "hurdle_regression_egpd_power", "hurdle_regression_egpd_power_mixture", "hurdle_regression_egpd_beta", "hurdle_regression_egpd_power_beta"). The zero probability is then pi_i = 1 / (1 + exp(-x_i'alpha)), modeled through a second fixed-effects-only predictor that reuses the same design matrix X as the response model (the response predictor carries the random effects, the zero predictor does not). The estimated zero-model coefficients alpha are returned by get_coef() alongside the response-model coefficients, with the suffix "_zero"

• "hurdle_gamma_varying_shape", "hurdle_regression_gamma_varying_shape": As "hurdle_gamma" and "hurdle_regression_gamma", but with the observation-specific gamma shape of "gamma_varying_shape". For "hurdle_regression_gamma_varying_shape" there are three predictors: the response mean, the structural-zero logit (coefficients with the suffix "_zero"), and log(shape)

• "zero_censored_power_transformed_normal": Likelihood of a censored and power-transformed normal variable for modeling data with a point mass at 0 and a continuous distribution for y > 0. The model used is Y = max(0,X)^lambda, X ~ N(mu, sigma^2), where mu = F(X) + Zb, and sigma and lambda are (auxiliary) parameters that are estimated. For more details on this model, see Sigrist et al. (2012, AOAS) "A dynamic nonstationary spatio-temporal model for short term prediction of precipitation"

• "zero_censored_power_transformed_normal_heteroscedastic": As "zero_censored_power_transformed_normal", but the standard deviation sigma of the latent normal variable varies across observations: log(sigma) = F_2(X) is related to fixed effects only (linear predictor or GPBoost algorithm), while mu = F(X) + Zb is related to both fixed and random effects. lambda is then the only (auxiliary) parameter that is estimated

• "zero_censored_shifted_gamma": Zero-censored shifted gamma likelihood for modeling data with a point mass at 0 and a continuous distribution for y > 0. The model used is Y = max(Z - xi, 0), where Z follows a gamma distribution with mean mu = exp(F(X) + Zb) and shape k. The shape k and shift xi are (auxiliary) parameters that are estimated. This is the version of "zero_one_censored_shifted_gamma" (Sigrist and Stahel, 2011) without the upper censoring at 1

• "zero_censored_shifted_gamma_varying_shape": As "zero_censored_shifted_gamma", but the shape k varies across observations: log(k) = F_s(X) is related to fixed effects only (covariates and / or the GPBoost tree-boosting algorithm; no random effects / GPs for the shape), while log(mu) = F(X) + Zb is related to both fixed and random effects. The shift xi is then the only (auxiliary) parameter that is estimated. The estimated coefficients of the log-shape model are returned alongside the mean-model coefficients (with the suffix "_shape")

Count response: y in {0, 1, 2, ...}

• "poisson": Poisson likelihood with a log link function

• "negative_binomial": negative binomial likelihood with a log link function (aka "nbinom2", "negative_binomial_2"). The variance is mu * (mu + r) / r, mu = mean, r = shape, with this parametrization

• "negative_binomial_1": Negative binomial 1 (aka "nbinom1") likelihood with a log link function. The variance is mu * (1 + phi), mu = mean, phi = dispersion, with this parametrization

• "zero_inflated_<base>": Two-part count likelihoods that combine a point mass 'p0' at zero with a count base distribution for the remaining mass '1 - p0' (the base can itself generate additional zeros). As for the hurdle likelihoods, 'exp(F(X) + Zb)' is the mean of the (non-structural) base component - not the unconditional response mean - so E(y) = (1 - p0) * base_mean, and 'p0' is estimated jointly with the base auxiliary parameters. Currently supported bases: "zero_inflated_poisson", "zero_inflated_negative_binomial" (aux: shape; aka "zero_inflated_nbinom2"), "zero_inflated_negative_binomial_1" (aux: dispersion; aka "zero_inflated_nbinom1")

• "zero_inflated_regression_<base>": As "zero_inflated_<base>", but with the logistic structural-zero model described under "hurdle_regression_<base>" above: "zero_inflated_regression_poisson", "zero_inflated_regression_negative_binomial", and "zero_inflated_regression_negative_binomial_1"

Binary response: y in {0, 1}

• "bernoulli_logit": Bernoulli likelihood with a logit link function for binary classification. Aliases: "binary", "binary_logit"

• "bernoulli_probit": Bernoulli likelihood with a probit link function for binary classification. Aliases: "binary_probit"

Proportion / fractional / bounded response: y in [0, 1]

• "quasi_bernoulli_logit": quasi-Bernoulli likelihood with a logit link function for y in [0,1]. Aliases: "quasi_binary", "quasi_binary_logit"

• "quasi_bernoulli_probit": quasi-Bernoulli likelihood with a probit link function for y in [0,1]. Aliases: "quasi_binary_probit"

• "binomial_logit": Binomial likelihood with a logit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials. Aliases: "binomial"

• "binomial_probit": Binomial likelihood with a probit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials

• "beta_binomial": Beta-binomial likelihood with a logit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials. Aliases: "betabinomial", "beta-binomial"

• "beta" : Beta likelihood with a logit link function (parametrization of Ferrari and Cribari-Neto, 2004). The response must lie strictly in (0,1)

• "zoctn": Zero-one censored transformed normal likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). The model used is T ~ N(mu, sigma^2), W = max(min(T,1),0), and Y = g(W), where g(x) = expit(a + b * logit(x)) for x in (0,1), mu = F(X) + Z_RE u, u denotes the random effects, Z_RE is their design matrix, and sigma, a, and b are (auxiliary) parameters that are estimated. For more details on this model, see Qiang and Sigrist (2026)

• "zero_one_censored_transformed_beta": Zero-one censored transformed beta likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). If T follows a beta distribution with mean mu = expit(F(X) + Zb) and precision phi, the observed response is obtained by applying the linear transformation Y = (1 + 2u) * T - u and censoring the result to [0,1]. The precision phi and shift u are (auxiliary) parameters that are estimated. For more details on this model, see Kosmidis and Zeileis (2025)

• "zero_one_censored_shifted_gamma": Zero-one censored shifted gamma likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). The model used is Y = min(max(Z - xi, 0), 1), where Z follows a gamma distribution with mean mu = exp(F(X) + Zb) and shape k. The shape k and shift xi are (auxiliary) parameters that are estimated. For more details on this model, see Sigrist and Stahel (2011)

Notes

• The first lines in the likelihoods source file contain additional comments on the specific parametrizations used

• Other likelihoods can be implemented upon request

A vector or matrix whose columns are categorical grouping variables. The elements being group levels defining grouped random effects. The elements of 'group_data' can be integer, double, or character. The number of columns corresponds to the number of grouped (intercept) random effects

A vector or matrix with numeric covariate data for grouped random coefficients

A vector with integer indices that indicate the corresponding categorical grouping variable (=columns) in 'group_data' for every covariate in 'group_rand_coef_data'. Counting starts at 1. The length of this index vector must equal the number of covariates in 'group_rand_coef_data'. For instance, c(1,1,2) means that the first two covariates (=first two columns) in 'group_rand_coef_data' have random coefficients corresponding to the first categorical grouping variable (=first column) in 'group_data', and the third covariate (=third column) in 'group_rand_coef_data' has a random coefficient corresponding to the second grouping variable (=second column) in 'group_data'

A vector of type logical (boolean). Indicates whether intercept random effects are dropped (only for random coefficients). If drop_intercept_group_rand_effect[k] is TRUE, the intercept random effect number k is dropped / not included. Only random effects with random slopes can be dropped.

A matrix with numeric coordinates (= inputs / features) for defining Gaussian processes

A vector or matrix with numeric covariate data for Gaussian process random coefficients

A string specifying the covariance function for the Gaussian process. Available options:

• "matern": Matern covariance function with the smoothness specified by the cov_fct_shape parameter (using the parametrization of Rasmussen and Williams, 2006)

• "matern_estimate_shape": same as "matern" but the smoothness parameter is also estimated

• "matern_space_time": Spatio-temporal Matern covariance function with different range parameters for space and time. Note that the first column in gp_coords must correspond to the time dimension

• "space_time_gneiting": Spatio-temporal covariance function given in Eq. (16) of Gneiting (2002). Note that the first column in gp_coords must correspond to the time dimension. This covariance has seven parameters (in the following order: sigma2, a, c, alpha, nu, beta, delta) which are all estimated by default. You can disable the estimation of some of these parameter using the 'estimate_cov_par_index' argument of the params argument in either the fit function of a gp_model object or the set_optim_params function prior to estimation.

• "matern_ard": anisotropic Matern covariance function with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of gp_coords

• "matern_ard_estimate_shape": same as "matern_ard" but the smoothness parameter is also estimated

• "exponential": Exponential covariance function (using the parametrization of Diggle and Ribeiro, 2007)

• "gaussian": Gaussian, aka squared exponential, covariance function (using the parametrization of Diggle and Ribeiro, 2007)

• "gaussian_ard": anisotropic Gaussian, aka squared exponential, covariance function with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of gp_coords

• "powered_exponential": powered exponential covariance function with the exponent specified by the cov_fct_shape parameter (using the parametrization of Diggle and Ribeiro, 2007)

• "wendland": Compactly supported Wendland covariance function (using the parametrization of Bevilacqua et al., 2019, AOS)

• "linear": linear covariance function. This corresponds to a Bayesian linear regression model with a Gaussian prior on the coefficients with a constant variance diagonal prior covariance, and the prior variance is estimated using empirical Bayes.

• "hurst": Hurst covariance function cov(s, s') = (sigma2 / 2) * ( ||s||^(2H) + ||s'||^(2H) - ||s - s'||^(2H) ), 0 < H < 1. For H = 0.5, this corresponds to Brownian motion in one dimension (-> see the 'estimate_cov_par_index' argument). This is the covariance for cov_fct_order = 1 (default). For cov_fct_order = 2, the second-order Hurst covariance cov(s, s') = sigma2 / (2 (2H - 1)) * ( ||s - s'||^(2H) - ||s||^(2H) - ||s'||^(2H) + 2H (s^T s') (||s||^(2H-2) + ||s'||^(2H-2)) ), 1 < H < 2, is used. In one dimension and for s, s' >= 0, H = 1.5 corresponds to an integrated Wiener process. See cov_fct_order for more details

• "hurst_ard": Hurst covariance function with with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of “gp_coords“ except for the first coordinate which has a range parameter of 1 due to identifiability with the marginal variance: cov(s, s') = (\sigma^2/2)\left[ \left(s_1^2 + \sum_{k=2}^d (s_k/l_k)^2\right)^H + \left({s'}_1^2 + \sum_{k=2}^d ({s'}_k/l_k)^2\right)^H - \left((s_1-{s'}_1)^2 + \sum_{k=2}^d ((s_k-{s'}_k)/l_k)^2\right)^H \right] . For cov_fct_order = 2, the second-order Hurst covariance (see "hurst") is applied to the scaled coordinates (s_1, s_2/l_2, ..., s_d/l_d)

• "ar1_mf_<base>": two-level autoregressive multifidelity covariance defined by f_H(x)=\rho f_L(x)+\delta(x), where f_L and \delta are independent Gaussian processes using the same base covariance type but separate parameter vectors.

  • For example, use "ar1_mf_matern", "ar1_mf_matern_ard", or "ar1_mf_matern_estimate_shape".

  • The last column of gp_coords must be 0 for low fidelity or 1 for high fidelity. All preceding columns are input coordinates for the GPs.

  • Covariance parameters are ordered as [low-fidelity base parameters, discrepancy base parameters, rho]. The two base blocks use the ordinary parameter ordering of the base covariance, and rho is unrestricted and can be negative.

  • Any supported base covariance except "wendland" can be used. Correlation tapering and Gaussian-process random coefficients are currently not supported for this model.

A numeric specifying the shape parameter of the covariance function (e.g., smoothness parameter for Matern and Wendland covariance) This parameter is irrelevant for some covariance functions such as the exponential or Gaussian

A string specifying the large data approximation for Gaussian processes. Available options:

• "none": No approximation

• "vecchia": Vecchia approximation; see Sigrist (2022, JMLR) for more details. For "matern_space_time" and the anisotropic ARD covariance functions ("matern_ard", "gaussian_ard", "matern_ard_estimate_shape"), neighbors are selected by absolute correlation, i.e., using distances in coordinates scaled by the range parameters, and they are redetermined during parameter estimation. For "space_time_gneiting" and "ar1_mf_<base>", neighbors are also selected by absolute correlation by default. Use gp_approx = "vecchia_euclidean" for Euclidean-distance selection.

• "full_scale_vecchia": Vecchia-inducing points full-scale (VIF) approximation; see Gyger, Furrer, and Sigrist (2026, JMLR) for more details

• "tapering": The covariance function is multiplied by a compactly supported Wendland correlation function

• "fitc": Fully Independent Training Conditional approximation aka modified predictive process approximation

• "full_scale_tapering": Full-scale approximation combining an inducing point / predictive process approximation with tapering on the residual process; see Gyger, Furrer, and Sigrist (2026, JUQ) for more details

• "vecchia_latent": similar as "vecchia" but a Vecchia approximation is applied to the latent Gaussian process for likelihood == "gaussian". For likelihood != "gaussian", "vecchia" and "vecchia_latent" are equivalent

An integer specifying the number of parallel threads for OMP. If num_parallel_threads = NULL, the default number of threads of the session is used, see gpb.get.default.num.threads. Initially, this is the number of physical performance cores, which is determined automatically as follows. On CPUs whose cores have different speeds (e.g., the performance and efficiency cores of recent Intel CPUs or Apple silicon), the fastest cores are counted, and hyperthreads are counted only once, since slow threads can slow down an entire model. The next fastest cores are added if the fastest ones alone would leave only a single thread. On Linux, a CPU bandwidth limit of a control group (e.g., of a container) is respected as well. If the environment variable OMP_NUM_THREADS is set, it takes precedence over both the topology of the CPU and the limit of the control group, and the number of threads that OpenMP is configured to use is taken as it is if the cores of the CPU cannot be determined. In either case, the limits of the OpenMP runtime itself still apply, in particular the limit of the contention group (OMP_THREAD_LIMIT). For ordinary use, leave num_parallel_threads unspecified to use this default. The default of the session can be changed with gpb.tune.num.threads, which benchmarks different numbers of threads, or with gpb.set.default.num.threads. Setting num_parallel_threads = 1 disables the OpenMP parallelization of the model and can substantially increase the runtime. A single thread should be chosen deliberately, e.g., to distribute the resources among concurrent model fits or for a specific test; it is not required for correctness or reproducibility

A boolean. If TRUE, GPU acceleration will be used if supported

A string specifying the method used for inverting covariance matrices. Available options:

• "default": iterative methods where possible, otherwise Cholesky factorization. Cholesky factorization is also used for likelihood = "gaussian_heteroscedastic_fixed_and_random"

• "cholesky": Cholesky factorization

• "iterative": iterative methods. A combination of the conjugate gradient, the Lanczos algorithm, and other methods.

  This is currently only supported for the following cases:

  • grouped random effects with more than one level

  • likelihood != "gaussian" and gp_approx == "vecchia" (non-Gaussian likelihoods with a Vecchia-Laplace approximation)

  • likelihood != "gaussian" and gp_approx == "full_scale_vecchia" (non-Gaussian likelihoods with a VIF approximation)

  • likelihood == "gaussian" and gp_approx == "full_scale_tapering" (Gaussian likelihood with a full-scale tapering approximation)

A vector with sample weights. For a Gaussian likelihood, the error variance ("nugget") for observation i is divided by weights[i]. For non-Gaussian likelihoods, the conditional log-likelihood contribution of observation i is multiplied by weights[i]. Consequently, weights affect the estimation of both random and fixed effects. Note that a Gaussian likelihood is calculated via the Laplace approximation when gp_approx = "vecchia_latent", when likelihood = "gaussian_latent", and when grouped random effects are combined with a Vecchia-approximated Gaussian process. In these cases, the weights act as for a non-Gaussian likelihood, i.e., the Gaussian log-likelihood contribution of observation i is multiplied by weights[i] instead of the error variance being divided by it.

A vector or matrix with observation-level data that some likelihoods require in addition to the response variable y (one row per data point). Currently, this is only used by the joint Tweedie likelihoods ("tweedie_joint" and its variants), for which its first column contains the observed number of events (e.g., claims) N. N must be a non-negative integer that is 0 if and only if y is 0

A numeric with a learning rate for the likelihood for generalized Bayesian inference (only non-Gaussian likelihoods)

A numeric specifying the range parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)

A numeric specifying the shape (=smoothness) parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)

An integer specifying the number of neighbors for the Vecchia and VIF approximations. Internal default values if NULL:

• 20 for gp_approx = "vecchia"

• 30 for gp_approx = "full_scale_vecchia"

Note: for prediction, the number of neighbors can be set through the 'num_neighbors_pred' parameter in the 'set_prediction_data' function. By default, num_neighbors_pred = 2 * num_neighbors. Further, the type of Vecchia approximation used for making predictions is set through the 'vecchia_pred_type' parameter in the 'set_prediction_data' function

A string specifying the ordering used in the Vecchia approximation. Available options:

• "none": the default ordering in the data is used

• "random": a random ordering

• "time": ordering accorrding to time (only for space-time models)

• "time_random_space": ordering according to time and randomly for all spatial points with the same time points (only for space-time models)

A string specifying the method for choosing inducing points Available options:

• "kmeans++: the k-means++ algorithm

• "cover_tree": the cover tree algorithm

• "random": random selection from data points

An integer specifying the number of inducing points / knots for FITC, full_scale_tapering, and VIF approximations. Internal default values if NULL:

• 500 for gp_approx = "FITC" and gp_approx = "full_scale_tapering"

• 200 for gp_approx = "full_scale_vecchia"

A numeric specifying the radius (= "spatial resolution") for the cover tree algorithm

An integer specifying the seed used for model creation (e.g., random ordering in Vecchia approximation)

A vector with elements indicating independent realizations of random effects / Gaussian processes (same values = same process realization). The elements of 'cluster_ids' can be integer, double, or character.

A logical. For an ar1_mf_<base> covariance, whether marginal mean models are fitted independently for low and high fidelity. If TRUE (the default), linear covariates supplied to fit() are internally expanded into separate low- and high-fidelity coefficient blocks, and the fidelity indicator is automatically added as a feature when the model is used in the GPBoost algorithm. Ignored for other covariance functions.

An integer specifying the order m of the "hurst" and "hurst_ard" covariance functions (also when used as base covariance in "ar1_mf_hurst" and "ar1_mf_hurst_ard"). Currently, the orders 1 (default) and 2 are supported. The Hurst exponent H of the order m satisfies m - 1 < H < m, and H = m - 0.5 is used as initial value. Order 1 is a fractional Brownian motion / field anchored at the origin (b(0) = 0), and order 2 is a second-order Hurst process / field anchored at the origin such that both the process and its gradient are zero there (b(0) = 0, grad b(0) = 0). The origin of the coordinates thus is a part of the model, and the coordinates are not centered internally. The order is not estimated. It can be chosen, e.g., by comparing the marginal likelihoods or by cross-validation.

• Continuous-time RW1 prior: use cov_function = "hurst", cov_fct_order = 1, and fix H = 0.5. This is a Brownian motion with b(0) = 0 and b'(t) = sqrt(q) W'(t), where W' denotes white noise and q = sigma2

• Continuous-time RW2 prior: use cov_function = "hurst", cov_fct_order = 2, and fix H = 1.5. This is an integrated Wiener process with b(0) = b'(0) = 0 and b”(t) = sqrt(q) W'(t), where q = 3 * sigma2

For these priors, use a one-dimensional time coordinate t >= 0 whose origin t = 0 is the desired anchor (shift the time points if necessary). H can be fixed using the 'init_cov_pars' and 'estimate_cov_par_index' arguments of the params argument, see the examples of fitGPModel. If the polynomial null space of the intrinsic model should not be penalized, include an intercept (RW1) or an intercept and a linear time trend (RW2) as fixed effects. Note that the order 2 model with H = 1.5 is the continuous-time RW2 prior evaluated at the time points and not the usual discrete intrinsic RW2 prior on a lattice: for equally spaced time points, adjacent second differences are correlated (correlation 1/4) whereas they are independent for the discrete RW2

A boolean. If TRUE, the data (groups, coordinates, covariate data for random coefficients, cluster IDs, weights, and additional likelihood data) is freed in R after initialization

A vector with response variable data

A matrix with numeric covariate data for the fixed effects linear regression term (if there is one)

A list with parameters for the estimation / optimization

• trace: boolean (default = FALSE). If TRUE, information on the progress of the parameter optimization is printed

• init_cov_pars: vector with numeric elements (default = NULL). Initial values for covariance parameters of Gaussian process and random effects (can be NULL). The order is same as the order of the parameters in the summary function: first is the error variance (only for "gaussian" likelihood), next follow the variances of the grouped random effects (if there are any, in the order provided in 'group_data'), and then follow the marginal variance and the ranges of the Gaussian process. If there are multiple Gaussian processes, then the variances and ranges follow alternatingly. If 'init_cov_pars = NULL', an internal choice is used that depends on the likelihood and the random effects type and covariance function. If you select the option 'trace = TRUE' in the 'params' argument, you will see the first initial covariance parameters in iteration 0.

• init_coef: vector with numeric elements (default = NULL). Initial values for the regression coefficients (if there are any, can be NULL)

• init_aux_pars: vector with numeric elements (default = NULL). Initial values for additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative_binomial likelihood)

• init_coef_aux_pars_from_iid_model: boolean (default = TRUE). If TRUE, regression coefficients and auxiliary parameters are initialized from an iid model (only for models with a linear regression term). This option is ignored if init_coef is provided. If init_aux_pars is provided but init_coef is not, only regression coefficients are initialized from an iid model. For likelihood = "gaussian_heteroscedastic_fixed_and_random", if init_cov_pars is not provided, the initial covariance parameters are obtained from a homoscedastic Gaussian process model, and the regression coefficients are then initialized from this model instead: the coefficients of the mean from its coefficients and the intercept of the log-error variance from the logarithm of its error variance.

• estimate_cov_par_index: vector with integer (default = -1). This allows for disabling the estimation of some (or all) covariance parameters. If 'estimate_cov_par_index' = -1, all covariance parameters are estimated. If estimate_cov_par_index != -1, this should be a vector with length equal to the number of covariance parameters, and estimate_cov_par_index[i] should be of bool type indicating whether parameter number i is estimated or not. For instance, estimate_cov_par_index = c(1,1,0) means that the first two covariance parameters are estimated and the last one not. Parameters that are not estimated are kept at their initial values (see 'init_cov_pars'). They are treated as known constants when calculating the standard errors of the estimated parameters, and their own standard errors are NaN.

• estimate_aux_pars: boolean (default = TRUE). If TRUE, additional parameters for non-Gaussian likelihoods are also estimated (e.g., shape parameter of a gamma or negative_binomial likelihood)

• optimizer_cov: string (default = "lbfgs"). Optimizer used for estimating covariance parameters. Options: "lbfgs", "gradient_descent", "fisher_scoring", "newton", "nelder_mead". If there are additional auxiliary parameters for non-Gaussian likelihoods, 'optimizer_cov' is also used for those

• optimizer_coef: string (default = "wls" for Gaussian likelihoods and "lbfgs" for other likelihoods). Optimizer used for estimating linear regression coefficients, if there are any (for the GPBoost algorithm there are usually none). Options: "gradient_descent", "lbfgs", "wls", "nelder_mead". Gradient descent steps are done simultaneously with gradient descent steps for the covariance parameters. "wls" refers to doing coordinate descent for the regression coefficients using weighted least squares. If 'optimizer_cov' is set to "nelder_mead" or "lbfgs", 'optimizer_coef' is automatically also set to the same value.

• maxit: integer (default = 1000). Maximal number of iterations for optimization algorithm. If maxit = -999, internal default values are used

• delta_rel_conv: numeric (default = 1E-6 except for "nelder_mead" for which the default is 1E-8). Convergence tolerance. The algorithm stops if the relative change in either the (approximate) log-likelihood or the parameters is below this value. If delta_rel_conv = -999, internal default values are used

• cg_max_num_it: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithms. If cg_max_num_it = -999, internal default values are used

• cg_max_num_it_tridiag: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithm when being run as Lanczos algorithm for tridiagonalization. If cg_max_num_it_tridiag = -999, internal default values are used

• cg_delta_conv: numeric (default = 1E-2). Tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithm when being used for parameter estimation. This is the tolerance of the "absolute" criterion only, the "relative" criterion is controlled by cg_rel_tol and cg_abs_tol. If cg_delta_conv = -999, internal default values are used

• cg_convergence_criterion: string (default = "absolute"). Stopping rule of conjugate gradient algorithms.

  • "absolute": stop when ||r||_2 < cg_delta_conv

  • "relative": stop when ||r||_2 <= max(cg_abs_tol, cg_rel_tol * ||b||_2), where b is the right-hand side of the linear system. The cg_abs_tol floor makes this robust for right-hand sides that are zero or very small

• cg_rel_tol: numeric (default = 1E-6). Relative tolerance of the "relative" stopping rule. If cg_rel_tol = -999, internal default values are used

• cg_abs_tol: numeric (default = 1E-8). Absolute tolerance floor of the "relative" stopping rule, which makes it robust for right-hand sides that are zero or very small. This is independent of cg_delta_conv, so that the two criteria can be varied separately. If cg_abs_tol = -999, internal default values are used

• cg_multi_rhs_convergence: string (default = "average"). How the stopping rule is aggregated over the columns of a linear system with several right-hand sides, such as the stochastic Lanczos quadrature. For cg_convergence_criterion = "relative", every right-hand side is normalized by its own norm, i.e. q_j = ||r_j||_2 / max(cg_abs_tol, cg_rel_tol * ||b_j||_2):

  • "average": stop when mean(q_j) <= 1

  • "max": stop when max(q_j) <= 1

  • "per_rhs": stop iterating on column j individually as soon as q_j <= 1

  For cg_convergence_criterion = "absolute" there is no per-column normalization:

  • "average": stop when mean(||r_j||_2) < cg_delta_conv. This is the historic behavior and the default

  • "max": stop when ||r_j||_2 < cg_delta_conv for all j

  • "per_rhs": stop iterating on column j as soon as ||r_j||_2 < cg_delta_conv

• num_rand_vec_trace: integer (default = 50). Number of random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix. If num_rand_vec_trace = -999, internal default values are used

• reuse_rand_vec_trace: boolean (default = TRUE). If true, random vectors (e.g., Rademacher) for stochastic approximations of the trace of a matrix are sampled only once at the beginning of the parameter estimation and reused in later trace approximations. Otherwise they are sampled every time a trace is calculated

• seed_rand_vec_trace: integer (default = 1). Seed number to generate random vectors (e.g., Rademacher)

• cg_preconditioner_type (string): Type of preconditioner used for conjugate gradient algorithms.

  • Options for grouped random effects:

    • "ssor" (= default): SSOR preconditioner

    • "incomplete_cholesky": zero fill-in incomplete Cholesky factorization

  • Options for likelihood != "gaussian" and gp_approx == "vecchia" or likelihood == "gaussian" and gp_approx == "vecchia_latent":

    • "vadu" (= default): (B^T * (D^-1 + W) * B) as preconditioner for inverting (B^T * D^-1 * B + W), where B^T * D^-1 * B approx= Sigma^-1

    • "fitc": FITC / modified predictive process preconditioner for inverting (B^-1 * D * B^-T + W^-1)

    • "pivoted_cholesky": (Lk * Lk^T + W^-1) as preconditioner for inverting (B^-1 * D * B^-T + W^-1), where Lk is a low-rank pivoted Cholesky approximation for Sigma and B^-1 * D * B^-T approx= Sigma

    • "incomplete_cholesky": zero fill-in incomplete (reverse) Cholesky factorization of (B^T * D^-1 * B + W) using the sparsity pattern of B^T * D^-1 * B approx= Sigma^-1

  • Options for likelihood != "gaussian" and gp_approx == "full_scale_vecchia":

    • "fitc" ( = default): FITC / modified predictive process preconditioner

    • "vifdu": VIF with diagonal update preconditioner

  • Options for likelihood == "gaussian" and gp_approx == "full_scale_tapering":

    • "fitc" (= default): modified predictive process preconditioner

    • "none": no preconditioner

• fitc_piv_chol_preconditioner_rank (integer ): Rank of the FITC and pivoted Cholesky decomposition preconditioners for iterative methods for Vecchia and VIF approximations (for full_scale_tapering, the same inducing points as in the approximation as used). If fitc_piv_chol_preconditioner_rank = -999, internal default values are used:

  • 200 for the FITC preconditioner

  • 50 for the pivoted Cholesky decomposition preconditioner

• convergence_criterion: string (default = "relative_change_in_log_likelihood", only relevant for "gradient_descent", "fisher_scoring", and "newton"). The convergence criterion used for terminating the optimization algorithm. Options: "relative_change_in_log_likelihood" or "relative_change_in_parameters". If convergence_criterion = "default", internal default values are used

• lr_cov: numeric (default = 0.1 for "gradient_descent" and 1. otherwise, relevant for "gradient_descent", "fisher_scoring", "newton", and "lbfgs"). Initial learning rate for covariance parameters if a gradient-based optimization method is used

  • If lr_cov = -999, internal default values are used (0.1 for "gradient_descent" and 1. otherwise)

  • If there are additional auxiliary parameters for non-Gaussian likelihoods, 'lr_cov' is also used for those

  • For "gradient_descent", "fisher_scoring", and "newton", this is the learning rate by which the search direction is multiplied. It can be decreased, and increased again, during the optimization

  • For "lbfgs" and "lbfgs_linesearch_nocedal_wright", this is not a learning rate by which every update is multiplied. It only determines the first trial step, which is the negative gradient divided by its norm and multiplied by 'lr_cov', i.e., a step of length 'lr_cov'. The subsequent steps are determined by the approximate Hessian of the L-BFGS algorithm and the line search. 'lr_cov' has no effect when the approximate Hessian of a previous call is reused (see 'reuse_learning_rates_gp_model' of the GPBoost algorithm)

• lr_coef: numeric (default = 0.1, only relevant for "gradient_descent", "fisher_scoring", and "newton"). Learning rate for fixed effect regression coefficients if gradient descent is used. If lr_coef = -999, internal default values are used

• use_nesterov_acc: boolean (default = TRUE, only relevant for "gradient_descent"). If TRUE Nesterov acceleration is used. This is used only for gradient descent

• nesterov_schedule_version: integer (Default = 0, only relevant for "gradient_descent"). Which version of Nesterov schedule should be used. If nesterov_schedule_version = -999, internal default values are used.

• acc_rate_coef: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for regression coefficients (if there are any) for Nesterov acceleration. If acc_rate_coef = -999, internal default values are used

• acc_rate_cov: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for covariance parameters for Nesterov acceleration. If acc_rate_cov = -999, internal default values are used

• momentum_offset: integer (Default = 2, only relevant for "gradient_descent"). Number of iterations for which no momentum is applied in the beginning. If momentum_offset = -999, internal default values are used.

• m_lbfgs: integer (Default = 6). Number of corrections to approximate the inverse Hessian matrix for the "lbfgs" optimizer. If m_lbfgs = -999, internal default values are used

• cold_restart_lbfgs: boolean (Default = TRUE). If TRUE, restarts of the "lbfgs" optimizers are "cold" restarts: the regression coefficients, the auxiliary parameters, and the modes of the Laplace approximations are reset to their initial values and only the covariance parameters are kept. Such restarts can escape a local optimum in which the optimizer is stuck. If FALSE, "warm" restarts are done: the optimization simply continues from the current parameters with a re-initialized approximate Hessian. Only relevant if max_num_restarts_lbfgs > 0. Note that every cold restart is a full optimization with up to maxit iterations, and the number of iterations reported is the sum over all restarts

• max_num_restarts_lbfgs: integer (Default = 0). Maximal number of restarts of the "lbfgs" optimizers after they have terminated. A restart continues from the current parameters with a re-initialized approximate Hessian. What else is reset depends on cold_restart_lbfgs. Restarts are only done as long as the negative log-likelihood still decreases by more than delta_rel_conv and as long as the total number of iterations is below maxit. This can help when the line search of "lbfgs" fails and the optimizer thus terminates prematurely without having converged. This can happen, e.g., for the "asymmetric_laplace" (= quantile regression) likelihood, whose approximate marginal likelihood is not smooth

• reset_modes_line_search_lbfgs: boolean (Default = FALSE). If TRUE, the mode finding of the Laplace approximation starts at the modes of the current parameters for every trial point of the line search of the "lbfgs" optimizer, and not at the modes of the previous trial point. Only relevant for non-Gaussian likelihoods and optimizer_cov = "lbfgs" or "lbfgs_linesearch_nocedal_wright"

• delta_conv_mode_finding: numeric (Default = 1E-8, and 1E-10 for the likelihood "gaussian_heteroscedastic_fixed_and_random"). Convergence tolerance in mode finding algorithm for Laplace approximation for non-Gaussian likelihoods. If delta_conv_mode_finding = -999, internal default values are used

Discontinued. Use the argument gp_approx instead

A string specifying the type of Vecchia approximation used for making predictions. This is discontinued here. Use the function 'set_prediction_data' to specify this

an integer specifying the number of neighbors for making predictions. This is discontinued here. Use the function 'set_prediction_data' to specify this

A numeric vector with additional fixed effects contributions that are added to the linear predictor (= offset). The length of this vector needs to equal the number of training data points times the number of fixed-effect sets.

This is discontinued. Use the renamed equivalent argument offset instead

A numeric specifying an additional parameter for the likelihood which cannot be estimated for this likelihood (e.g., degrees of freedom for likelihood = "t_fix_df"). This is not to be confused with any auxiliary parameters that can be estimated and accessed through the function get_aux_pars after estimation. Note that this likelihood_additional_param parameter is irrelevant for many likelihoods. If likelihood_additional_param = NULL, the following default and mandatory-value rules apply:

• df = 2 for likelihood = "t_fix_df"

• No default is used for likelihood = "asymmetric_laplace"; a quantile strictly between 0 and 1 is required

• No default is used for the Tweedie likelihoods with a fixed power ("tweedie_fixed_p" and the other "_fixed_p" variants); a power strictly between 1.01 and 1.99 is required

A GPModel

A boolean. If TRUE, (approximate) standard errors are calculated (= square root of diagonal of the inverse Fisher information for Gaussian likelihoods and square root of diagonal of a numerically approximated inverse Hessian for non-Gaussian likelihoods). For non-Gaussian likelihoods, no standard errors of the covariance and auxiliary parameters are calculated if the model is used in the GPBoost algorithm, since they require the fixed effects given by the tree ensemble

A GPModel

A boolean. If TRUE, (approximate) standard errors are calculated (= square root of diagonal of the inverse Fisher information for Gaussian likelihoods and square root of diagonal of a numerically approximated inverse Hessian for non-Gaussian likelihoods). For non-Gaussian likelihoods, no standard errors of the covariance and auxiliary parameters are calculated if the model is used in the GPBoost algorithm, since they require the fixed effects given by the tree ensemble

A GPModel

A boolean. If TRUE, (approximate) standard errors are calculated (= square root of diagonal of the inverse Fisher information for Gaussian likelihoods and square root of diagonal of a numerically approximated inverse Hessian for non-Gaussian likelihoods). For non-Gaussian likelihoods, no standard errors of the covariance and auxiliary parameters are calculated if the model is used in the GPBoost algorithm, since they require the fixed effects given by the tree ensemble

A GPModel

A boolean. If TRUE, (approximate) standard errors are calculated (= square root of diagonal of the inverse Fisher information for Gaussian likelihoods and square root of diagonal of a numerically approximated inverse Hessian for non-Gaussian likelihoods). For non-Gaussian likelihoods, no standard errors of the covariance and auxiliary parameters are calculated if the model is used in the GPBoost algorithm, since they require the fixed effects given by the tree ensemble

A GPModel

A boolean. If TRUE, (approximate) standard errors are calculated (= square root of diagonal of the inverse Fisher information for Gaussian likelihoods and square root of diagonal of a numerically approximated inverse Hessian for non-Gaussian likelihoods). For non-Gaussian likelihoods, no standard errors of the covariance and auxiliary parameters are calculated if the model is used in the GPBoost algorithm, since they require the fixed effects given by the tree ensemble

A GPModel

A boolean. If TRUE, (approximate) standard errors are calculated (= square root of diagonal of the inverse Fisher information for Gaussian likelihoods and square root of diagonal of a numerically approximated inverse Hessian for non-Gaussian likelihoods). For non-Gaussian likelihoods, no standard errors of the covariance and auxiliary parameters are calculated if the model is used in the GPBoost algorithm, since they require the fixed effects given by the tree ensemble

A vector containing the outer categorical grouping variable within which the inner_var is nested in. Can be of type integer, double, or character.

A vector containing the inner nested categorical grouping variable

Object of class gpb.Dataset

other parameters

the name of the information field to get (see details)

a matrix object, a dgCMatrix object or a character representing a filename

a list of parameters. See the "Dataset Parameters" section of the parameter documentation for a list of parameters and valid values.

reference dataset. When GPBoost creates a Dataset, it does some preprocessing like binning continuous features into histograms. If you want to apply the same bin boundaries from an existing dataset to new data, pass that existing Dataset to this argument.

names of columns

categorical features. This can either be a character vector of feature names or an integer vector with the indices of the features (e.g. c(1L, 10L) to say "the first and tenth columns").

GPBoost constructs its data format, called a "Dataset", from tabular data. By default, this Dataset object on the R side does keep a copy of the raw data. If you set free_raw_data = TRUE, no copy of the raw data is kept (this reduces memory usage)

a list of information of the gpb.Dataset object

other information to pass to info or parameters pass to params

Object of class gpb.Dataset

gpb.Dataset object, training data

a matrix object, a dgCMatrix object or a character representing a filename

a list of information of the gpb.Dataset object

other information to pass to info.

object of class gpb.Dataset

object filename of output file

object of class gpb.Dataset

categorical features. This can either be a character vector of feature names or an integer vector with the indices of the features (e.g. c(1L, 10L) to say "the first and tenth columns").

object of class gpb.Dataset

object of class gpb.Dataset

A data.frame or data.table to prepare.

A set of rules from the data preparator, if already used. This should be an R list, where names are column names in data and values are named character vectors whose names are column values and whose values are new values to replace them with.

list of "tuning" parameters. See the parameter documentation for more information. A few key parameters:

• learning_rate: The learning rate, also called shrinkage or damping parameter (default = 0.1). An important tuning parameter for boosting. Lower values usually lead to higher predictive accuracy but more boosting iterations are needed

• num_leaves: Number of leaves in a tree. Tuning parameter for tree-boosting (default = 31)

• max_depth: Maximal depth of a tree. Tuning parameter for tree-boosting (default = no limit)

• min_data_in_leaf: Minimal number of samples per leaf. Tuning parameter for tree-boosting (default = 20)

• lambda_l2: L2 regularization (default = 0)

• lambda_l1: L1 regularization (default = 0)

• max_bin: Maximal number of bins that feature values will be bucketed in (default = 255)

• line_search_step_length (default = FALSE): If TRUE, a line search is done to find the optimal step length for every boosting update (see, e.g., Friedman 2001). This is then multiplied by the learning rate

• train_gp_model_cov_pars (default = TRUE): If TRUE, the covariance parameters of the Gaussian process are estimated in every boosting iterations, otherwise the gp_model parameters are not estimated. In the latter case, you need to either estimate them beforehand or provide values via the 'init_cov_pars' parameter when creating the gp_model

• use_gp_model_for_validation (default = TRUE): If TRUE, the Gaussian process is also used (in addition to the tree model) for calculating predictions on the validation data

• leaves_newton_update (default = FALSE): Set this to TRUE to do a Newton update step for the tree leaves after the gradient step. Applies only to Gaussian process boosting (GPBoost algorithm)

• num_threads: Number of threads. For the best speed, set this to the number of real CPU cores(parallel::detectCores(logical = FALSE)), not the number of threads (most CPU using hyper-threading to generate 2 threads per CPU core). For ordinary use, leave this unspecified to use the default number of threads. Setting num_threads = 1 disables the parallelization and can substantially increase the runtime. A single thread should be chosen deliberately, e.g., to distribute the resources among concurrent model fits or for a specific test; it is not required for correctness or reproducibility.

a gpb.Dataset object, used for training. Some functions, such as gpb.cv, may allow you to pass other types of data like matrix and then separately supply label as a keyword argument.

A GPModel object that contains the random effects (Gaussian process and / or grouped random effects) model

number of boosting iterations (= number of trees). This is the most important tuning parameter for boosting

int. Activates early stopping. Requires at least one validation data and one metric. When this parameter is non-null, training will stop if the evaluation of any metric on any validation set fails to improve for early_stopping_rounds consecutive boosting rounds. If training stops early, the returned model will have attribute best_iter set to the iteration number of the best iteration.

(list) This provides the possibility to use a list of pre-defined CV folds (each element must be a vector of test fold's indices). When folds are supplied, the nfold and stratified parameters are ignored.

The original dataset is randomly partitioned into nfold equal-sized subsamples. This allows for doing nfold-CV

Evaluation metric to be monitored when doing CV and parameter tuning. Can be a character string or vector of character strings. If not NULL, the metric in params will be overridden.

• Default ="test_neg_log_likelihood" if there is a GPModel

• Non-exhaustive list of supported metrics: "test_neg_log_likelihood", "mse", "rmse", "mae", "crps_gaussian", "auc", "average_precision", "binary_logloss", "binary_error"

• See the "metric" section of the parameter documentation for a complete list of valid metrics

verbosity for output, if <= 0, also will disable the print of evaluation during training

Boolean. If TRUE, the gp_model (Gaussian process and/or random effects) is also used (in addition to the tree model) for calculating predictions on the validation data. If FALSE, the gp_model (random effects part) is ignored for making predictions and only the tree ensemble is used for making predictions for calculating the validation / test error.

Boolean (default = FALSE). If TRUE, the covariance parameters of the gp_model model are estimated using the out-of-sample (OOS) predictions on the validation data using the optimal number of iterations (after performing the CV). This corresponds to the GPBoostOOS algorithm.

Boolean. If TRUE, the covariance parameters of the gp_model (Gaussian process and/or random effects) are estimated in every boosting iterations, otherwise the gp_model parameters are not estimated. In the latter case, you need to either estimate them beforehand or provide the values via the init_cov_pars parameter when creating the gp_model

Vector of labels, used if data is not an gpb.Dataset

vector of response values. If not NULL, will set to dataset

(character) The distribution of the response variable (=label) conditional on fixed and random effects. This only needs to be set when doing independent boosting without random effects / Gaussian processes.

Evaluation metric to be monitored when doing CV and parameter tuning. This can be a string, function, or list with a mixture of strings and functions.

• a. character vector: Non-exhaustive list of supported metrics: "test_neg_log_likelihood", "mse", "rmse", "mae", "auc", "average_precision", "binary_logloss", "binary_error" See the "metric" section of the parameter documentation for a complete list of valid metrics.

• b. function: You can provide a custom evaluation function. This should accept the keyword arguments preds and dtrain and should return a named list with three elements:

  • name: A string with the name of the metric, used for printing and storing results.

  • value: A single number indicating the value of the metric for the given predictions and true values

  • higher_better: A boolean indicating whether higher values indicate a better fit. For example, this would be FALSE for metrics like MAE or RMSE.

• c. list: If a list is given, it should only contain character vectors and functions. These should follow the requirements from the descriptions above.

Boolean, TRUE will record iteration message to booster$record_evals

evaluation output frequency, only effect when verbose > 0

boolean, whether to show standard deviation of cross validation. This parameter defaults to TRUE.

a boolean indicating whether sampling of folds should be stratified by the values of outcome labels.

path of model file of gpb.Booster object, will continue training from this model

feature names, if not null, will use this to overwrite the names in dataset

categorical features. This can either be a character vector of feature names or an integer vector with the indices of the features (e.g. c(1L, 10L) to say "the first and tenth columns").

List of callback functions that are applied at each iteration.

Boolean, setting it to TRUE (not the default value) will transform the booster model into a predictor model which frees up memory and the original datasets

Boolean, setting it to TRUE (not the default value) will delete the boosters of the individual folds

other parameters, see Parameters.rst for more information.

Object of class gpb.Booster

number of iteration want to predict with, NULL or <= 0 means use best iteration

Object of class gpb.Booster

Name of the dataset to return evaluation results for.

Name of the evaluation metric to return results for.

An integer vector of iterations you want to get evaluation results for. If NULL (the default), evaluation results for all iterations will be returned.

TRUE will return evaluation error instead

list with candidate parameters defining the grid over which a search is done

integer with number of random trial on parameter grid. If NULL, a deterministic search is done

a gpb.Dataset object, used for training. Some functions, such as gpb.cv, may allow you to pass other types of data like matrix and then separately supply label as a keyword argument.

A GPModel object that contains the random effects (Gaussian process and / or grouped random effects) model

list with other parameters not included in param_grid

number of boosting iterations (= number of trees). This is the most important tuning parameter for boosting

int. Activates early stopping. Requires at least one validation data and one metric. When this parameter is non-null, training will stop if the evaluation of any metric on any validation set fails to improve for early_stopping_rounds consecutive boosting rounds. If training stops early, the returned model will have attribute best_iter set to the iteration number of the best iteration.

(list) This provides the possibility to use a list of pre-defined CV folds (each element must be a vector of test fold's indices). When folds are supplied, the nfold and stratified parameters are ignored.

The original dataset is randomly partitioned into nfold equal-sized subsamples. This allows for doing nfold-CV

Evaluation metric to be monitored when doing CV and parameter tuning. Can be a character string or vector of character strings. If not NULL, the metric in params will be overridden.

• Default ="test_neg_log_likelihood" if there is a GPModel

• Non-exhaustive list of supported metrics: "test_neg_log_likelihood", "mse", "rmse", "mae", "crps_gaussian", "auc", "average_precision", "binary_logloss", "binary_error"

• See the "metric" section of the parameter documentation for a complete list of valid metrics

integer. Whether to display information on the progress of tuning parameter choice. If None or 0, verbose is of. If = 1, summary progress information is displayed for every parameter combination. If >= 2, detailed progress is displayed at every boosting stage for every parameter combination.

Seed for generating folds when doing nfold CV

Boolean. If TRUE, the gp_model (Gaussian process and/or random effects) is also used (in addition to the tree model) for calculating predictions on the validation data. If FALSE, the gp_model (random effects part) is ignored for making predictions and only the tree ensemble is used for making predictions for calculating the validation / test error.

Boolean. If TRUE, the covariance parameters of the gp_model (Gaussian process and/or random effects) are estimated in every boosting iterations, otherwise the gp_model parameters are not estimated. In the latter case, you need to either estimate them beforehand or provide the values via the init_cov_pars parameter when creating the gp_model

Vector of labels, used if data is not an gpb.Dataset

vector of response values. If not NULL, will set to dataset

(character) The distribution of the response variable (=label) conditional on fixed and random effects. This only needs to be set when doing independent boosting without random effects / Gaussian processes.

Evaluation metric to be monitored when doing CV and parameter tuning. This can be a string, function, or list with a mixture of strings and functions.

• a. character vector: Non-exhaustive list of supported metrics: "test_neg_log_likelihood", "mse", "rmse", "mae", "auc", "average_precision", "binary_logloss", "binary_error" See the "metric" section of the parameter documentation for a complete list of valid metrics.

• b. function: You can provide a custom evaluation function. This should accept the keyword arguments preds and dtrain and should return a named list with three elements:

  • name: A string with the name of the metric, used for printing and storing results.

  • value: A single number indicating the value of the metric for the given predictions and true values

  • higher_better: A boolean indicating whether higher values indicate a better fit. For example, this would be FALSE for metrics like MAE or RMSE.

• c. list: If a list is given, it should only contain character vectors and functions. These should follow the requirements from the descriptions above.

a boolean indicating whether sampling of folds should be stratified by the values of outcome labels.

path of model file of gpb.Booster object, will continue training from this model

feature names, if not null, will use this to overwrite the names in dataset

categorical features. This can either be a character vector of feature names or an integer vector with the indices of the features (e.g. c(1L, 10L) to say "the first and tenth columns").

List of callback functions that are applied at each iteration.

a boolean indicating whether all tried parameter combinations are returned

other parameters, see Parameters.rst for more information.

object of class gpb.Booster.

whether to show importance in relative percentage.

object of class gpb.Booster.

a matrix object or a dgCMatrix object.

an integer vector of indices of rows needed.

number of iteration want to predict with, NULL or <= 0 means use best iteration.

path of model file

a str containing the model

object of class gpb.Booster

number of iterations you want to predict with. NULL or <= 0 means use best iteration

a data.table returned by gpb.importance.

maximal number of top features to include into the plot.

the name of importance measure to plot, can be "Gain", "Cover" or "Frequency".

(base R barplot) allows to adjust the left margin size to fit feature names.

(base R barplot) passed as cex.names parameter to barplot. Set a number smaller than 1.0 to make the bar labels smaller than R's default and values greater than 1.0 to make them larger.

other parameters passed to graphics::barplot

a data.table returned by gpb.interprete.

maximal number of top features to include into the plot.

the column numbers of layout, will be used only for multiclass classification feature contribution.

(base R barplot) allows to adjust the left margin size to fit feature names.

(base R barplot) passed as cex.names parameter to barplot.

A gpb.Booster model object

A matrix with data for creating partial dependence plots

A vector of length two of type string with names of the columns or integer with indices of the columns in data for which an interaction dependence plot is created

If TRUE, the plot is done on the scale of the tree-ensemble. If FALSE, the plot is done on the (potentially transformed) response variable scale (e.g. probabilities for classification)

Number of grid points per variable (used only if a variable is not discrete) For continuous variables, the two-dimensional grid for the interaction plot has dimension c(n.pt.per.var, n.pt.per.var)

Fraction of random samples in data to be used for calculating the partial dependence plot

A vector of length two of type boolean. If an entry is TRUE, the evaluation grid of the corresponding variable is set to the unique values of the variable

An integer indicating the class in multi-class classification (value from 0 to num_class - 1)

A character string indicating the type of the plot. Supported values: "filled.contour" and "contour"

Parameter passed to the filled.contour or contour function

Parameter passed to the filled.contour or contour function

Parameter passed to the filled.contour or contour function

Parameter passed to the filled.contour or contour function

Parameter passed to the filled.contour or contour function

A boolean. If TRUE, the data for creating the partial dependence plot is returned

Additional parameters passed to the filled.contour or contour function

A gpb.Booster model object

A matrix with data for creating partial dependence plots

A string with a name of the column or an integer with an index of the column in data for which a dependence plot is created

If TRUE, the plot is done on the scale of the tree-ensemble. If FALSE, the plot is done on the (potentially transformed) response variable scale (e.g. probabilities for classification)

Evaluation grid size (used only if x is not discrete)

Fraction of random samples in data to be used for calculating the partial dependence plot

A boolean. If TRUE, the evaluation grid is set to the unique values of x

An integer indicating the class in multi-class classification (value from 0 to num_class - 1)

Parameter passed to plot

Parameter passed to plot

Parameter passed to plot

Parameter passed to plot

A boolean. If TRUE, the data for creating the partial dependence plot is returned

Additional parameters passed to plot

Object of class gpb.Booster

saved filename

int or NULL, optional (default=NULL) Start index of the iteration to predict. If NULL or <= 0, starts from the first iteration.

int or NULL, optional (default=NULL) Limit number of iterations in the prediction. If NULL, if the best iteration exists and start_iteration is NULL or <= 0, the best iteration is used; otherwise, all iterations from start_iteration are used. If <= 0, all iterations from start_iteration are used (no limits).

If TRUE, the raw data (predictor / covariate data) for the Booster is also saved. Enable this option if you want to change start_iteration or num_iteration at prediction time after loading.

Additional named arguments passed to the predict() method of the gpb.Booster object passed to object. This is only used when there is a gp_model and when save_raw_data=FALSE

An integer specifying the number of threads. It is limited by the largest number of threads that GPBoost uses on its own: the number of threads that OMP uses when GPBoost determines its default (usually the number of logical processors, or the value of the environment variable OMP_NUM_THREADS if it is set), the limit of the contention group of OpenMP (OMP_THREAD_LIMIT), and, unless OMP_NUM_THREADS is set, a CPU bandwidth limit of a control group on Linux. A number of threads that is specified for an individual model via num_parallel_threads is not limited by this. If num_threads is not positive, the automatically selected number of threads is used again

An integer specifying the number of threads. If num_threads is not positive, the default number of threads of the session is used, see gpb.get.default.num.threads. Initially, this is the automatically selected number of threads (normally based on the number of physical performance cores, subject to the settings of OpenMP and of the system, see the num_parallel_threads argument of GPModel), and it is the number of threads set by gpb.tune.num.threads or gpb.set.default.num.threads once one of them has been called

list of "tuning" parameters. See the parameter documentation for more information. A few key parameters:

• learning_rate: The learning rate, also called shrinkage or damping parameter (default = 0.1). An important tuning parameter for boosting. Lower values usually lead to higher predictive accuracy but more boosting iterations are needed

• num_leaves: Number of leaves in a tree. Tuning parameter for tree-boosting (default = 31)

• max_depth: Maximal depth of a tree. Tuning parameter for tree-boosting (default = no limit)

• min_data_in_leaf: Minimal number of samples per leaf. Tuning parameter for tree-boosting (default = 20)

• lambda_l2: L2 regularization (default = 0)

• lambda_l1: L1 regularization (default = 0)

• max_bin: Maximal number of bins that feature values will be bucketed in (default = 255)

• line_search_step_length (default = FALSE): If TRUE, a line search is done to find the optimal step length for every boosting update (see, e.g., Friedman 2001). This is then multiplied by the learning rate

• train_gp_model_cov_pars (default = TRUE): If TRUE, the covariance parameters of the Gaussian process are estimated in every boosting iterations, otherwise the gp_model parameters are not estimated. In the latter case, you need to either estimate them beforehand or provide values via the 'init_cov_pars' parameter when creating the gp_model

• use_gp_model_for_validation (default = TRUE): If TRUE, the Gaussian process is also used (in addition to the tree model) for calculating predictions on the validation data

• leaves_newton_update (default = FALSE): Set this to TRUE to do a Newton update step for the tree leaves after the gradient step. Applies only to Gaussian process boosting (GPBoost algorithm)

• num_threads: Number of threads. For the best speed, set this to the number of real CPU cores(parallel::detectCores(logical = FALSE)), not the number of threads (most CPU using hyper-threading to generate 2 threads per CPU core). For ordinary use, leave this unspecified to use the default number of threads. Setting num_threads = 1 disables the parallelization and can substantially increase the runtime. A single thread should be chosen deliberately, e.g., to distribute the resources among concurrent model fits or for a specific test; it is not required for correctness or reproducibility.

a gpb.Dataset object, used for training. Some functions, such as gpb.cv, may allow you to pass other types of data like matrix and then separately supply label as a keyword argument.

number of boosting iterations (= number of trees). This is the most important tuning parameter for boosting

A GPModel object that contains the random effects (Gaussian process and / or grouped random effects) model

Boolean. If TRUE, the gp_model (Gaussian process and/or random effects) is also used (in addition to the tree model) for calculating predictions on the validation data. If FALSE, the gp_model (random effects part) is ignored for making predictions and only the tree ensemble is used for making predictions for calculating the validation / test error.

Boolean. If TRUE, the covariance parameters of the gp_model (Gaussian process and/or random effects) are estimated in every boosting iterations, otherwise the gp_model parameters are not estimated. In the latter case, you need to either estimate them beforehand or provide the values via the init_cov_pars parameter when creating the gp_model

a list of gpb.Dataset objects, used for validation

(character) The distribution of the response variable (=label) conditional on fixed and random effects. This only needs to be set when doing independent boosting without random effects / Gaussian processes.

Evaluation metric to be monitored when doing CV and parameter tuning. This can be a string, function, or list with a mixture of strings and functions.

• a. character vector: Non-exhaustive list of supported metrics: "test_neg_log_likelihood", "mse", "rmse", "mae", "auc", "average_precision", "binary_logloss", "binary_error" See the "metric" section of the parameter documentation for a complete list of valid metrics.

• b. function: You can provide a custom evaluation function. This should accept the keyword arguments preds and dtrain and should return a named list with three elements:

  • name: A string with the name of the metric, used for printing and storing results.

  • value: A single number indicating the value of the metric for the given predictions and true values

  • higher_better: A boolean indicating whether higher values indicate a better fit. For example, this would be FALSE for metrics like MAE or RMSE.

• c. list: If a list is given, it should only contain character vectors and functions. These should follow the requirements from the descriptions above.

verbosity for output, if <= 0, also will disable the print of evaluation during training

Boolean, TRUE will record iteration message to booster$record_evals

evaluation output frequency, only effect when verbose > 0

path of model file of gpb.Booster object, will continue training from this model

feature names, if not null, will use this to overwrite the names in dataset

categorical features. This can either be a character vector of feature names or an integer vector with the indices of the features (e.g. c(1L, 10L) to say "the first and tenth columns").

int. Activates early stopping. Requires at least one validation data and one metric. When this parameter is non-null, training will stop if the evaluation of any metric on any validation set fails to improve for early_stopping_rounds consecutive boosting rounds. If training stops early, the returned model will have attribute best_iter set to the iteration number of the best iteration.

List of callback functions that are applied at each iteration.

Boolean, setting it to TRUE (not the default value) will transform the booster model into a predictor model which frees up memory and the original datasets

other parameters, see the parameter documentation for more information.

A character vector specifying the workloads that are benchmarked:

• "all" (= default): all workloads below

• "grouped_re": a single-level grouped random effect model. The calculations are cheap per data point, which makes the overhead of the parallelization visible

• "vecchia_non_gaussian": a non-Gaussian Gaussian process model with a Vecchia approximation and iterative methods

• "crossed_re_iterative": crossed grouped random effects with iterative methods

Restrict the workloads to the model class that you mainly use if you know it

An integer vector with the numbers of threads that are benchmarked. If NULL, powers of two, the automatically selected number of threads and the largest number of threads that GPBoost uses on its own are benchmarked

An integer specifying the number of repeated measurements per workload and number of threads. The median of the repetitions is used. Every workload is measured at least twice if n_rep is at least two, since a single measurement has no spread from which the noise could be estimated

A numeric specifying the relative difference in runtime that is considered negligible. The smallest number of threads whose aggregated runtime is within this tolerance of the best aggregated runtime is selected, so that more threads are not used for a negligible gain. If the measurements are noisier than this, the observed noise is used instead

A numeric specifying approximately how many seconds the benchmark may take. Repetitions are dropped if the measurements take longer. This is a soft limit: every workload simulates its data, creates its models and is measured at least once, and the remaining time is divided among the workloads that are still to be measured

A numeric specifying how much slower than the fastest number of threads a number of threads may be for a single workload and still be selected. This prevents that a number of threads which is much slower for one model class is selected because the aggregate over all workloads looks acceptable

A logical. If TRUE, the selected number of threads is set as the default of the session, see gpb.set.default.num.threads. Set this to FALSE to only measure

A logical. If TRUE, the progress and the results are printed

A string, either "default" or "small". The small workloads run in a few seconds, but they are too small for the number of threads to matter, so they are only useful for checking that the benchmark runs and the default of the session is never changed for them

List of callback functions that are applied at each iteration.

a gpb.Dataset object, used for training. Some functions, such as gpb.cv, may allow you to pass other types of data like matrix and then separately supply label as a keyword argument.

(list) This provides the possibility to use a list of pre-defined CV folds (each element must be a vector of test fold's indices). When folds are supplied, the nfold and stratified parameters are ignored.

The original dataset is randomly partitioned into nfold equal-sized subsamples. This allows for doing nfold-CV

Seed for generating folds when doing nfold CV

int. Activates early stopping. Requires at least one validation data and one metric. When this parameter is non-null, training will stop if the evaluation of any metric on any validation set fails to improve for early_stopping_rounds consecutive boosting rounds. If training stops early, the returned model will have attribute best_iter set to the iteration number of the best iteration.

Evaluation metric to be monitored when doing CV and parameter tuning. Can be a character string or vector of character strings. If not NULL, the metric in params will be overridden.

• Default ="test_neg_log_likelihood" if there is a GPModel

• Non-exhaustive list of supported metrics: "test_neg_log_likelihood", "mse", "rmse", "mae", "crps_gaussian", "auc", "average_precision", "binary_logloss", "binary_error"

• See the "metric" section of the parameter documentation for a complete list of valid metrics

integer. Whether to display information on the progress of tuning parameter choice. If None or 0, verbose is of. If = 1, summary progress information is displayed for every parameter combination. If >= 2, detailed progress is displayed at every boosting stage for every parameter combination.

Evaluation metric to be monitored when doing CV and parameter tuning. This can be a string, function, or list with a mixture of strings and functions.

• a. character vector: Non-exhaustive list of supported metrics: "test_neg_log_likelihood", "mse", "rmse", "mae", "auc", "average_precision", "binary_logloss", "binary_error" See the "metric" section of the parameter documentation for a complete list of valid metrics.

• b. function: You can provide a custom evaluation function. This should accept the keyword arguments preds and dtrain and should return a named list with three elements:

  • name: A string with the name of the metric, used for printing and storing results.

  • value: A single number indicating the value of the metric for the given predictions and true values

  • higher_better: A boolean indicating whether higher values indicate a better fit. For example, this would be FALSE for metrics like MAE or RMSE.

• c. list: If a list is given, it should only contain character vectors and functions. These should follow the requirements from the descriptions above.

evaluation output frequency, only effect when verbose > 0

a list of gpb.Dataset objects, used for validation

Boolean, TRUE will record iteration message to booster$record_evals

feature names, if not null, will use this to overwrite the names in dataset

categorical features. This can either be a character vector of feature names or an integer vector with the indices of the features (e.g. c(1L, 10L) to say "the first and tenth columns").

path of model file of gpb.Booster object, will continue training from this model

number of boosting iterations (= number of trees). This is the most important tuning parameter for boosting

(character) The distribution of the response variable (=label) conditional on fixed and random effects. This only needs to be set when doing independent boosting without random effects / Gaussian processes.

list of "tuning" parameters. See the parameter documentation for more information. A few key parameters:

• learning_rate: The learning rate, also called shrinkage or damping parameter (default = 0.1). An important tuning parameter for boosting. Lower values usually lead to higher predictive accuracy but more boosting iterations are needed

• num_leaves: Number of leaves in a tree. Tuning parameter for tree-boosting (default = 31)

• max_depth: Maximal depth of a tree. Tuning parameter for tree-boosting (default = no limit)

• min_data_in_leaf: Minimal number of samples per leaf. Tuning parameter for tree-boosting (default = 20)

• lambda_l2: L2 regularization (default = 0)

• lambda_l1: L1 regularization (default = 0)

• max_bin: Maximal number of bins that feature values will be bucketed in (default = 255)

• line_search_step_length (default = FALSE): If TRUE, a line search is done to find the optimal step length for every boosting update (see, e.g., Friedman 2001). This is then multiplied by the learning rate

• train_gp_model_cov_pars (default = TRUE): If TRUE, the covariance parameters of the Gaussian process are estimated in every boosting iterations, otherwise the gp_model parameters are not estimated. In the latter case, you need to either estimate them beforehand or provide values via the 'init_cov_pars' parameter when creating the gp_model

• use_gp_model_for_validation (default = TRUE): If TRUE, the Gaussian process is also used (in addition to the tree model) for calculating predictions on the validation data

• leaves_newton_update (default = FALSE): Set this to TRUE to do a Newton update step for the tree leaves after the gradient step. Applies only to Gaussian process boosting (GPBoost algorithm)

• num_threads: Number of threads. For the best speed, set this to the number of real CPU cores(parallel::detectCores(logical = FALSE)), not the number of threads (most CPU using hyper-threading to generate 2 threads per CPU core). For ordinary use, leave this unspecified to use the default number of threads. Setting num_threads = 1 disables the parallelization and can substantially increase the runtime. A single thread should be chosen deliberately, e.g., to distribute the resources among concurrent model fits or for a specific test; it is not required for correctness or reproducibility.

verbosity for output, if <= 0, also will disable the print of evaluation during training

A GPModel object that contains the random effects (Gaussian process and / or grouped random effects) model

Boolean. If TRUE, a line search is done to find the optimal step length for every boosting update (see, e.g., Friedman 2001). This is then multiplied by the learning_rate. Applies only to the GPBoost algorithm

Boolean. If TRUE, the gp_model (Gaussian process and/or random effects) is also used (in addition to the tree model) for calculating predictions on the validation data. If FALSE, the gp_model (random effects part) is ignored for making predictions and only the tree ensemble is used for making predictions for calculating the validation / test error.

Boolean. If TRUE, the covariance parameters of the gp_model (Gaussian process and/or random effects) are estimated in every boosting iterations, otherwise the gp_model parameters are not estimated. In the latter case, you need to either estimate them beforehand or provide the values via the init_cov_pars parameter when creating the gp_model

a gpb.Dataset object, used for training. Some functions, such as gpb.cv, may allow you to pass other types of data like matrix and then separately supply label as a keyword argument.

Vector of response values / labels, used if data is not an gpb.Dataset

Vector of weights. The GPBoost algorithm currently does not support weights

list of "tuning" parameters. See the parameter documentation for more information. A few key parameters:

• learning_rate: The learning rate, also called shrinkage or damping parameter (default = 0.1). An important tuning parameter for boosting. Lower values usually lead to higher predictive accuracy but more boosting iterations are needed

• num_leaves: Number of leaves in a tree. Tuning parameter for tree-boosting (default = 31)

• max_depth: Maximal depth of a tree. Tuning parameter for tree-boosting (default = no limit)

• min_data_in_leaf: Minimal number of samples per leaf. Tuning parameter for tree-boosting (default = 20)

• lambda_l2: L2 regularization (default = 0)

• lambda_l1: L1 regularization (default = 0)

• max_bin: Maximal number of bins that feature values will be bucketed in (default = 255)

• line_search_step_length (default = FALSE): If TRUE, a line search is done to find the optimal step length for every boosting update (see, e.g., Friedman 2001). This is then multiplied by the learning rate

• train_gp_model_cov_pars (default = TRUE): If TRUE, the covariance parameters of the Gaussian process are estimated in every boosting iterations, otherwise the gp_model parameters are not estimated. In the latter case, you need to either estimate them beforehand or provide values via the 'init_cov_pars' parameter when creating the gp_model

• use_gp_model_for_validation (default = TRUE): If TRUE, the Gaussian process is also used (in addition to the tree model) for calculating predictions on the validation data

• leaves_newton_update (default = FALSE): Set this to TRUE to do a Newton update step for the tree leaves after the gradient step. Applies only to Gaussian process boosting (GPBoost algorithm)

• num_threads: Number of threads. For the best speed, set this to the number of real CPU cores(parallel::detectCores(logical = FALSE)), not the number of threads (most CPU using hyper-threading to generate 2 threads per CPU core). For ordinary use, leave this unspecified to use the default number of threads. Setting num_threads = 1 disables the parallelization and can substantially increase the runtime. A single thread should be chosen deliberately, e.g., to distribute the resources among concurrent model fits or for a specific test; it is not required for correctness or reproducibility.

number of boosting iterations (= number of trees). This is the most important tuning parameter for boosting

A GPModel object that contains the random effects (Gaussian process and / or grouped random effects) model

Boolean. If TRUE, the gp_model (Gaussian process and/or random effects) is also used (in addition to the tree model) for calculating predictions on the validation data. If FALSE, the gp_model (random effects part) is ignored for making predictions and only the tree ensemble is used for making predictions for calculating the validation / test error.

Boolean. If TRUE, the covariance parameters of the gp_model (Gaussian process and/or random effects) are estimated in every boosting iterations, otherwise the gp_model parameters are not estimated. In the latter case, you need to either estimate them beforehand or provide the values via the init_cov_pars parameter when creating the gp_model

a list of gpb.Dataset objects, used for validation

(character) The distribution of the response variable (=label) conditional on fixed and random effects. This only needs to be set when doing independent boosting without random effects / Gaussian processes.

Evaluation metric to be monitored when doing CV and parameter tuning. This can be a string, function, or list with a mixture of strings and functions.

• a. character vector: Non-exhaustive list of supported metrics: "test_neg_log_likelihood", "mse", "rmse", "mae", "auc", "average_precision", "binary_logloss", "binary_error" See the "metric" section of the parameter documentation for a complete list of valid metrics.

• b. function: You can provide a custom evaluation function. This should accept the keyword arguments preds and dtrain and should return a named list with three elements:

  • name: A string with the name of the metric, used for printing and storing results.

  • value: A single number indicating the value of the metric for the given predictions and true values

  • higher_better: A boolean indicating whether higher values indicate a better fit. For example, this would be FALSE for metrics like MAE or RMSE.

• c. list: If a list is given, it should only contain character vectors and functions. These should follow the requirements from the descriptions above.

verbosity for output, if <= 0, also will disable the print of evaluation during training

Boolean, TRUE will record iteration message to booster$record_evals

evaluation output frequency, only effect when verbose > 0

int. Activates early stopping. Requires at least one validation data and one metric. When this parameter is non-null, training will stop if the evaluation of any metric on any validation set fails to improve for early_stopping_rounds consecutive boosting rounds. If training stops early, the returned model will have attribute best_iter set to the iteration number of the best iteration.

path of model file of gpb.Booster object, will continue training from this model

feature names, if not null, will use this to overwrite the names in dataset

categorical features. This can either be a character vector of feature names or an integer vector with the indices of the features (e.g. c(1L, 10L) to say "the first and tenth columns").

List of callback functions that are applied at each iteration.

Additional arguments passed to gpb.train. For example

• valids: a list of gpb.Dataset objects, used for validation

• eval: evaluation function, can be (a list of) character or custom eval function

• record: Boolean, TRUE will record iteration message to booster$record_evals

• colnames: feature names, if not null, will use this to overwrite the names in dataset

• categorical_feature: categorical features. This can either be a character vector of feature names or an integer vector with the indices of the features (e.g. c(1L, 10L) to say "the first and tenth columns").

• reset_data: Boolean, setting it to TRUE (not the default value) will transform the booster model into a predictor model which frees up memory and the original datasets

filename for loading

A GPModel

A vector with numeric elements. Covariance parameters of Gaussian process and random effects

A vector with response variable data

A numeric vector with fixed effects, e.g., containing a linear predictor. The length of this vector needs to equal the number of training data points times the number of fixed-effect sets.

A vector with numeric elements. Additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative_binomial likelihood)

A GPModel

A vector with numeric elements. Covariance parameters of Gaussian process and random effects

A vector with response variable data

A numeric vector with fixed effects, e.g., containing a linear predictor. The length of this vector needs to equal the number of training data points times the number of fixed-effect sets.

A vector with numeric elements. Additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative_binomial likelihood)

a GPModel

A boolean. If TRUE, the response variable (label) is predicted, otherwise the latent random effects

A boolean. If TRUE, the (posterior) predictive variances are calculated

A boolean. If TRUE, the (posterior) predictive covariance is calculated in addition to the (posterior) predictive mean

A boolean. If TRUE, samples from the posterior are drawn

A boolean. If TRUE, samples from the prior are drawn

A numeric with the number of posterior samples to draw if 'sample_posterior=TRUE'

A numeric with the number of prior samples to draw if 'sample_prior=TRUE'

Observed data (can be NULL, e.g. when the model has been estimated already and the same data is used for making predictions)

A vector containing covariance parameters which are used if the GPModel has not been trained or if predictions should be made for other parameters than the trained ones

A vector or matrix with elements being group levels for which predictions are made (if there are grouped random effects in the GPModel)

A vector or matrix with covariate data for grouped random coefficients (if there are some in the GPModel)

A matrix with prediction coordinates (=features) for Gaussian process (if there is a GP in the GPModel)

A vector or matrix with covariate data for Gaussian process random coefficients (if there are some in the GPModel)

A vector with elements indicating the realizations of random effects / Gaussian processes for which predictions are made (set to NULL if you have not specified this when creating the GPModel)

A matrix with prediction covariate data for the fixed effects linear regression term (if there is one in the GPModel)

A boolean. If TRUE, predictions are done using a priory set data via the function '$set_prediction_data' (this option is not used by users directly)

A numeric vector with additional fixed effects contributions that are added to the linear predictor (= offset). The length of this vector needs to equal the number of training data points times the number of fixed-effect sets.

A numeric vector with additional fixed effects contributions that are added to the linear predictor for the prediction points (= offset). The length of this vector needs to equal the number of prediction points times the number of fixed-effect sets.

This is discontinued. Use the renamed equivalent argument offset instead

This is discontinued. Use the renamed equivalent argument offset_pred instead

A string specifying the type of Vecchia approximation used for making predictions. This is discontinued here. Use the function 'set_prediction_data' to specify this

an integer specifying the number of neighbors for making predictions. This is discontinued here. Use the function 'set_prediction_data' to specify this

(not used, ignore this, simply here that there is no CRAN warning)

Object of class gpb.Booster

a matrix object, a dgCMatrix object or a character representing a filename

int or NULL, optional (default=NULL) Start index of the iteration to predict. If NULL or <= 0, starts from the first iteration.

int or NULL, optional (default=NULL) Limit number of iterations in the prediction. If NULL, if the best iteration exists and start_iteration is NULL or <= 0, the best iteration is used; otherwise, all iterations from start_iteration are used. If <= 0, all iterations from start_iteration are used (no limits).

If TRUE latent variables, both fixed effects (tree-ensemble) and random effects (gp_model) are predicted. Otherwise, the response variable (label) is predicted. Depending on how the argument 'pred_latent' is set, different values are returned from this function; see the 'Value' section for more details. If there is no gp_model, this argument corresponds to 'raw_score' in LightGBM.

whether predict leaf index instead.

return per-feature contributions for each record.

only used for prediction for text file. True if text file has header

whether to reshape the vector of predictions to a matrix form when there are several prediction outputs per case.

A vector or matrix with elements being group levels for which predictions are made (if there are grouped random effects in the GPModel)

A vector or matrix with covariate data for grouped random coefficients (if there are some in the GPModel)

A matrix with prediction coordinates (=features) for Gaussian process (if there is a GP in the GPModel)

A vector or matrix with covariate data for Gaussian process random coefficients (if there are some in the GPModel)

A vector with elements indicating the realizations of random effects / Gaussian processes for which predictions are made (set to NULL if you have not specified this when creating the GPModel)

A boolean. If TRUE, the (posterior) predictive covariance is calculated in addition to the (posterior) predictive mean

A boolean. If TRUE, the (posterior) predictive variances are calculated

A boolean. If TRUE, samples from the posterior are drawn

A numeric with the number of posterior samples to draw if 'sample_posterior=TRUE'

A vector containing covariance parameters which are used if the gp_model has not been trained or if predictions should be made for other parameters than the trained ones

A numeric vector. Offsets for prediction: additional fixed effects contributions that are added to the predictor for the prediction points. The length of this vector needs to equal the number of prediction points times the number of fixed-effect sets.

A boolean. If TRUE, predictions are only made for the tree ensemble part and the gp_model is ignored

This is discontinued. Use the renamed equivalent argument pred_latent instead

A string specifying the type of Vecchia approximation used for making predictions. This is discontinued here. Use the function 'set_prediction_data' to specify this

an integer specifying the number of neighbors for making predictions. This is discontinued here. Use the function 'set_prediction_data' to specify this

Additional named arguments passed to the predict() method of the gpb.Booster object passed to object.

A GPModel or a gpb.Booster with a GPModel

A boolean. If TRUE, the (posterior) predictive variances are calculated

Additional arguments passed to methods.

A GPModel or a gpb.Booster with a GPModel

A boolean. If TRUE, the (posterior) predictive variances are calculated

A numeric vector with additional fixed effects contributions that are added to the linear predictor (= offset). The length of this vector needs to equal the number of training data points times the number of fixed-effect sets.

Additional arguments passed to methods.

a connection or the name of the file where the R object is saved to or read from.

a hook function for handling reference objects.

a GPModel

filename for saving

R object to serialize.

a connection or the name of the file where the R object is saved to or read from.

a logical. If TRUE or NA, an ASCII representation is written; otherwise (default), a binary one is used. See the comments in the help for save.

the workspace format version to use. NULL specifies the current default version (2). Versions prior to 2 are not supported, so this will only be relevant when there are later versions.

a logical specifying whether saving to a named file is to use "gzip" compression, or one of "gzip", "bzip2" or "xz" to indicate the type of compression to be used. Ignored if file is a connection.

a hook function for handling reference objects.

whether to save the model in a raw variable or not, recommended to leave it to TRUE.

A GPModel

A list with parameters for the estimation / optimization

• trace: boolean (default = FALSE). If TRUE, information on the progress of the parameter optimization is printed

• init_cov_pars: vector with numeric elements (default = NULL). Initial values for covariance parameters of Gaussian process and random effects (can be NULL). The order is same as the order of the parameters in the summary function: first is the error variance (only for "gaussian" likelihood), next follow the variances of the grouped random effects (if there are any, in the order provided in 'group_data'), and then follow the marginal variance and the ranges of the Gaussian process. If there are multiple Gaussian processes, then the variances and ranges follow alternatingly. If 'init_cov_pars = NULL', an internal choice is used that depends on the likelihood and the random effects type and covariance function. If you select the option 'trace = TRUE' in the 'params' argument, you will see the first initial covariance parameters in iteration 0.

• init_coef: vector with numeric elements (default = NULL). Initial values for the regression coefficients (if there are any, can be NULL)

• init_aux_pars: vector with numeric elements (default = NULL). Initial values for additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative_binomial likelihood)

• init_coef_aux_pars_from_iid_model: boolean (default = TRUE). If TRUE, regression coefficients and auxiliary parameters are initialized from an iid model (only for models with a linear regression term). This option is ignored if init_coef is provided. If init_aux_pars is provided but init_coef is not, only regression coefficients are initialized from an iid model. For likelihood = "gaussian_heteroscedastic_fixed_and_random", if init_cov_pars is not provided, the initial covariance parameters are obtained from a homoscedastic Gaussian process model, and the regression coefficients are then initialized from this model instead: the coefficients of the mean from its coefficients and the intercept of the log-error variance from the logarithm of its error variance.

• estimate_cov_par_index: vector with integer (default = -1). This allows for disabling the estimation of some (or all) covariance parameters. If 'estimate_cov_par_index' = -1, all covariance parameters are estimated. If estimate_cov_par_index != -1, this should be a vector with length equal to the number of covariance parameters, and estimate_cov_par_index[i] should be of bool type indicating whether parameter number i is estimated or not. For instance, estimate_cov_par_index = c(1,1,0) means that the first two covariance parameters are estimated and the last one not. Parameters that are not estimated are kept at their initial values (see 'init_cov_pars'). They are treated as known constants when calculating the standard errors of the estimated parameters, and their own standard errors are NaN.

• estimate_aux_pars: boolean (default = TRUE). If TRUE, additional parameters for non-Gaussian likelihoods are also estimated (e.g., shape parameter of a gamma or negative_binomial likelihood)

• optimizer_cov: string (default = "lbfgs"). Optimizer used for estimating covariance parameters. Options: "lbfgs", "gradient_descent", "fisher_scoring", "newton", "nelder_mead". If there are additional auxiliary parameters for non-Gaussian likelihoods, 'optimizer_cov' is also used for those

• optimizer_coef: string (default = "wls" for Gaussian likelihoods and "lbfgs" for other likelihoods). Optimizer used for estimating linear regression coefficients, if there are any (for the GPBoost algorithm there are usually none). Options: "gradient_descent", "lbfgs", "wls", "nelder_mead". Gradient descent steps are done simultaneously with gradient descent steps for the covariance parameters. "wls" refers to doing coordinate descent for the regression coefficients using weighted least squares. If 'optimizer_cov' is set to "nelder_mead" or "lbfgs", 'optimizer_coef' is automatically also set to the same value.

• maxit: integer (default = 1000). Maximal number of iterations for optimization algorithm. If maxit = -999, internal default values are used

• delta_rel_conv: numeric (default = 1E-6 except for "nelder_mead" for which the default is 1E-8). Convergence tolerance. The algorithm stops if the relative change in either the (approximate) log-likelihood or the parameters is below this value. If delta_rel_conv = -999, internal default values are used

• cg_max_num_it: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithms. If cg_max_num_it = -999, internal default values are used

• cg_max_num_it_tridiag: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithm when being run as Lanczos algorithm for tridiagonalization. If cg_max_num_it_tridiag = -999, internal default values are used

• cg_delta_conv: numeric (default = 1E-2). Tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithm when being used for parameter estimation. This is the tolerance of the "absolute" criterion only, the "relative" criterion is controlled by cg_rel_tol and cg_abs_tol. If cg_delta_conv = -999, internal default values are used

• cg_convergence_criterion: string (default = "absolute"). Stopping rule of conjugate gradient algorithms.

  • "absolute": stop when ||r||_2 < cg_delta_conv

  • "relative": stop when ||r||_2 <= max(cg_abs_tol, cg_rel_tol * ||b||_2), where b is the right-hand side of the linear system. The cg_abs_tol floor makes this robust for right-hand sides that are zero or very small

• cg_rel_tol: numeric (default = 1E-6). Relative tolerance of the "relative" stopping rule. If cg_rel_tol = -999, internal default values are used

• cg_abs_tol: numeric (default = 1E-8). Absolute tolerance floor of the "relative" stopping rule, which makes it robust for right-hand sides that are zero or very small. This is independent of cg_delta_conv, so that the two criteria can be varied separately. If cg_abs_tol = -999, internal default values are used

• cg_multi_rhs_convergence: string (default = "average"). How the stopping rule is aggregated over the columns of a linear system with several right-hand sides, such as the stochastic Lanczos quadrature. For cg_convergence_criterion = "relative", every right-hand side is normalized by its own norm, i.e. q_j = ||r_j||_2 / max(cg_abs_tol, cg_rel_tol * ||b_j||_2):

  • "average": stop when mean(q_j) <= 1

  • "max": stop when max(q_j) <= 1

  • "per_rhs": stop iterating on column j individually as soon as q_j <= 1

  For cg_convergence_criterion = "absolute" there is no per-column normalization:

  • "average": stop when mean(||r_j||_2) < cg_delta_conv. This is the historic behavior and the default

  • "max": stop when ||r_j||_2 < cg_delta_conv for all j

  • "per_rhs": stop iterating on column j as soon as ||r_j||_2 < cg_delta_conv

• num_rand_vec_trace: integer (default = 50). Number of random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix. If num_rand_vec_trace = -999, internal default values are used

• reuse_rand_vec_trace: boolean (default = TRUE). If true, random vectors (e.g., Rademacher) for stochastic approximations of the trace of a matrix are sampled only once at the beginning of the parameter estimation and reused in later trace approximations. Otherwise they are sampled every time a trace is calculated

• seed_rand_vec_trace: integer (default = 1). Seed number to generate random vectors (e.g., Rademacher)

• cg_preconditioner_type (string): Type of preconditioner used for conjugate gradient algorithms.

  • Options for grouped random effects:

    • "ssor" (= default): SSOR preconditioner

    • "incomplete_cholesky": zero fill-in incomplete Cholesky factorization

  • Options for likelihood != "gaussian" and gp_approx == "vecchia" or likelihood == "gaussian" and gp_approx == "vecchia_latent":

    • "vadu" (= default): (B^T * (D^-1 + W) * B) as preconditioner for inverting (B^T * D^-1 * B + W), where B^T * D^-1 * B approx= Sigma^-1

    • "fitc": FITC / modified predictive process preconditioner for inverting (B^-1 * D * B^-T + W^-1)

    • "pivoted_cholesky": (Lk * Lk^T + W^-1) as preconditioner for inverting (B^-1 * D * B^-T + W^-1), where Lk is a low-rank pivoted Cholesky approximation for Sigma and B^-1 * D * B^-T approx= Sigma

    • "incomplete_cholesky": zero fill-in incomplete (reverse) Cholesky factorization of (B^T * D^-1 * B + W) using the sparsity pattern of B^T * D^-1 * B approx= Sigma^-1

  • Options for likelihood != "gaussian" and gp_approx == "full_scale_vecchia":

    • "fitc" ( = default): FITC / modified predictive process preconditioner

    • "vifdu": VIF with diagonal update preconditioner

  • Options for likelihood == "gaussian" and gp_approx == "full_scale_tapering":

    • "fitc" (= default): modified predictive process preconditioner

    • "none": no preconditioner

• fitc_piv_chol_preconditioner_rank (integer ): Rank of the FITC and pivoted Cholesky decomposition preconditioners for iterative methods for Vecchia and VIF approximations (for full_scale_tapering, the same inducing points as in the approximation as used). If fitc_piv_chol_preconditioner_rank = -999, internal default values are used:

  • 200 for the FITC preconditioner

  • 50 for the pivoted Cholesky decomposition preconditioner

• convergence_criterion: string (default = "relative_change_in_log_likelihood", only relevant for "gradient_descent", "fisher_scoring", and "newton"). The convergence criterion used for terminating the optimization algorithm. Options: "relative_change_in_log_likelihood" or "relative_change_in_parameters". If convergence_criterion = "default", internal default values are used

• lr_cov: numeric (default = 0.1 for "gradient_descent" and 1. otherwise, relevant for "gradient_descent", "fisher_scoring", "newton", and "lbfgs"). Initial learning rate for covariance parameters if a gradient-based optimization method is used

  • If lr_cov = -999, internal default values are used (0.1 for "gradient_descent" and 1. otherwise)

  • If there are additional auxiliary parameters for non-Gaussian likelihoods, 'lr_cov' is also used for those

  • For "gradient_descent", "fisher_scoring", and "newton", this is the learning rate by which the search direction is multiplied. It can be decreased, and increased again, during the optimization

  • For "lbfgs" and "lbfgs_linesearch_nocedal_wright", this is not a learning rate by which every update is multiplied. It only determines the first trial step, which is the negative gradient divided by its norm and multiplied by 'lr_cov', i.e., a step of length 'lr_cov'. The subsequent steps are determined by the approximate Hessian of the L-BFGS algorithm and the line search. 'lr_cov' has no effect when the approximate Hessian of a previous call is reused (see 'reuse_learning_rates_gp_model' of the GPBoost algorithm)

• lr_coef: numeric (default = 0.1, only relevant for "gradient_descent", "fisher_scoring", and "newton"). Learning rate for fixed effect regression coefficients if gradient descent is used. If lr_coef = -999, internal default values are used

• use_nesterov_acc: boolean (default = TRUE, only relevant for "gradient_descent"). If TRUE Nesterov acceleration is used. This is used only for gradient descent

• nesterov_schedule_version: integer (Default = 0, only relevant for "gradient_descent"). Which version of Nesterov schedule should be used. If nesterov_schedule_version = -999, internal default values are used.

• acc_rate_coef: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for regression coefficients (if there are any) for Nesterov acceleration. If acc_rate_coef = -999, internal default values are used

• acc_rate_cov: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for covariance parameters for Nesterov acceleration. If acc_rate_cov = -999, internal default values are used

• momentum_offset: integer (Default = 2, only relevant for "gradient_descent"). Number of iterations for which no momentum is applied in the beginning. If momentum_offset = -999, internal default values are used.

• m_lbfgs: integer (Default = 6). Number of corrections to approximate the inverse Hessian matrix for the "lbfgs" optimizer. If m_lbfgs = -999, internal default values are used

• cold_restart_lbfgs: boolean (Default = TRUE). If TRUE, restarts of the "lbfgs" optimizers are "cold" restarts: the regression coefficients, the auxiliary parameters, and the modes of the Laplace approximations are reset to their initial values and only the covariance parameters are kept. Such restarts can escape a local optimum in which the optimizer is stuck. If FALSE, "warm" restarts are done: the optimization simply continues from the current parameters with a re-initialized approximate Hessian. Only relevant if max_num_restarts_lbfgs > 0. Note that every cold restart is a full optimization with up to maxit iterations, and the number of iterations reported is the sum over all restarts

• max_num_restarts_lbfgs: integer (Default = 0). Maximal number of restarts of the "lbfgs" optimizers after they have terminated. A restart continues from the current parameters with a re-initialized approximate Hessian. What else is reset depends on cold_restart_lbfgs. Restarts are only done as long as the negative log-likelihood still decreases by more than delta_rel_conv and as long as the total number of iterations is below maxit. This can help when the line search of "lbfgs" fails and the optimizer thus terminates prematurely without having converged. This can happen, e.g., for the "asymmetric_laplace" (= quantile regression) likelihood, whose approximate marginal likelihood is not smooth

• reset_modes_line_search_lbfgs: boolean (Default = FALSE). If TRUE, the mode finding of the Laplace approximation starts at the modes of the current parameters for every trial point of the line search of the "lbfgs" optimizer, and not at the modes of the previous trial point. Only relevant for non-Gaussian likelihoods and optimizer_cov = "lbfgs" or "lbfgs_linesearch_nocedal_wright"

• delta_conv_mode_finding: numeric (Default = 1E-8, and 1E-10 for the likelihood "gaussian_heteroscedastic_fixed_and_random"). Convergence tolerance in mode finding algorithm for Laplace approximation for non-Gaussian likelihoods. If delta_conv_mode_finding = -999, internal default values are used

A GPModel

A list with parameters for the estimation / optimization

• trace: boolean (default = FALSE). If TRUE, information on the progress of the parameter optimization is printed

• init_cov_pars: vector with numeric elements (default = NULL). Initial values for covariance parameters of Gaussian process and random effects (can be NULL). The order is same as the order of the parameters in the summary function: first is the error variance (only for "gaussian" likelihood), next follow the variances of the grouped random effects (if there are any, in the order provided in 'group_data'), and then follow the marginal variance and the ranges of the Gaussian process. If there are multiple Gaussian processes, then the variances and ranges follow alternatingly. If 'init_cov_pars = NULL', an internal choice is used that depends on the likelihood and the random effects type and covariance function. If you select the option 'trace = TRUE' in the 'params' argument, you will see the first initial covariance parameters in iteration 0.

• init_coef: vector with numeric elements (default = NULL). Initial values for the regression coefficients (if there are any, can be NULL)

• init_aux_pars: vector with numeric elements (default = NULL). Initial values for additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative_binomial likelihood)

• init_coef_aux_pars_from_iid_model: boolean (default = TRUE). If TRUE, regression coefficients and auxiliary parameters are initialized from an iid model (only for models with a linear regression term). This option is ignored if init_coef is provided. If init_aux_pars is provided but init_coef is not, only regression coefficients are initialized from an iid model. For likelihood = "gaussian_heteroscedastic_fixed_and_random", if init_cov_pars is not provided, the initial covariance parameters are obtained from a homoscedastic Gaussian process model, and the regression coefficients are then initialized from this model instead: the coefficients of the mean from its coefficients and the intercept of the log-error variance from the logarithm of its error variance.

• estimate_cov_par_index: vector with integer (default = -1). This allows for disabling the estimation of some (or all) covariance parameters. If 'estimate_cov_par_index' = -1, all covariance parameters are estimated. If estimate_cov_par_index != -1, this should be a vector with length equal to the number of covariance parameters, and estimate_cov_par_index[i] should be of bool type indicating whether parameter number i is estimated or not. For instance, estimate_cov_par_index = c(1,1,0) means that the first two covariance parameters are estimated and the last one not. Parameters that are not estimated are kept at their initial values (see 'init_cov_pars'). They are treated as known constants when calculating the standard errors of the estimated parameters, and their own standard errors are NaN.

• estimate_aux_pars: boolean (default = TRUE). If TRUE, additional parameters for non-Gaussian likelihoods are also estimated (e.g., shape parameter of a gamma or negative_binomial likelihood)

• optimizer_cov: string (default = "lbfgs"). Optimizer used for estimating covariance parameters. Options: "lbfgs", "gradient_descent", "fisher_scoring", "newton", "nelder_mead". If there are additional auxiliary parameters for non-Gaussian likelihoods, 'optimizer_cov' is also used for those

• optimizer_coef: string (default = "wls" for Gaussian likelihoods and "lbfgs" for other likelihoods). Optimizer used for estimating linear regression coefficients, if there are any (for the GPBoost algorithm there are usually none). Options: "gradient_descent", "lbfgs", "wls", "nelder_mead". Gradient descent steps are done simultaneously with gradient descent steps for the covariance parameters. "wls" refers to doing coordinate descent for the regression coefficients using weighted least squares. If 'optimizer_cov' is set to "nelder_mead" or "lbfgs", 'optimizer_coef' is automatically also set to the same value.

• maxit: integer (default = 1000). Maximal number of iterations for optimization algorithm. If maxit = -999, internal default values are used

• delta_rel_conv: numeric (default = 1E-6 except for "nelder_mead" for which the default is 1E-8). Convergence tolerance. The algorithm stops if the relative change in either the (approximate) log-likelihood or the parameters is below this value. If delta_rel_conv = -999, internal default values are used

• cg_max_num_it: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithms. If cg_max_num_it = -999, internal default values are used

• cg_max_num_it_tridiag: integer (default = 1000). Maximal number of iterations for conjugate gradient algorithm when being run as Lanczos algorithm for tridiagonalization. If cg_max_num_it_tridiag = -999, internal default values are used

• cg_delta_conv: numeric (default = 1E-2). Tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithm when being used for parameter estimation. This is the tolerance of the "absolute" criterion only, the "relative" criterion is controlled by cg_rel_tol and cg_abs_tol. If cg_delta_conv = -999, internal default values are used

• cg_convergence_criterion: string (default = "absolute"). Stopping rule of conjugate gradient algorithms.

  • "absolute": stop when ||r||_2 < cg_delta_conv

  • "relative": stop when ||r||_2 <= max(cg_abs_tol, cg_rel_tol * ||b||_2), where b is the right-hand side of the linear system. The cg_abs_tol floor makes this robust for right-hand sides that are zero or very small

• cg_rel_tol: numeric (default = 1E-6). Relative tolerance of the "relative" stopping rule. If cg_rel_tol = -999, internal default values are used

• cg_abs_tol: numeric (default = 1E-8). Absolute tolerance floor of the "relative" stopping rule, which makes it robust for right-hand sides that are zero or very small. This is independent of cg_delta_conv, so that the two criteria can be varied separately. If cg_abs_tol = -999, internal default values are used

• cg_multi_rhs_convergence: string (default = "average"). How the stopping rule is aggregated over the columns of a linear system with several right-hand sides, such as the stochastic Lanczos quadrature. For cg_convergence_criterion = "relative", every right-hand side is normalized by its own norm, i.e. q_j = ||r_j||_2 / max(cg_abs_tol, cg_rel_tol * ||b_j||_2):

  • "average": stop when mean(q_j) <= 1

  • "max": stop when max(q_j) <= 1

  • "per_rhs": stop iterating on column j individually as soon as q_j <= 1

  For cg_convergence_criterion = "absolute" there is no per-column normalization:

  • "average": stop when mean(||r_j||_2) < cg_delta_conv. This is the historic behavior and the default

  • "max": stop when ||r_j||_2 < cg_delta_conv for all j

  • "per_rhs": stop iterating on column j as soon as ||r_j||_2 < cg_delta_conv

• num_rand_vec_trace: integer (default = 50). Number of random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix. If num_rand_vec_trace = -999, internal default values are used

• reuse_rand_vec_trace: boolean (default = TRUE). If true, random vectors (e.g., Rademacher) for stochastic approximations of the trace of a matrix are sampled only once at the beginning of the parameter estimation and reused in later trace approximations. Otherwise they are sampled every time a trace is calculated

• seed_rand_vec_trace: integer (default = 1). Seed number to generate random vectors (e.g., Rademacher)

• cg_preconditioner_type (string): Type of preconditioner used for conjugate gradient algorithms.

  • Options for grouped random effects:

    • "ssor" (= default): SSOR preconditioner

    • "incomplete_cholesky": zero fill-in incomplete Cholesky factorization

  • Options for likelihood != "gaussian" and gp_approx == "vecchia" or likelihood == "gaussian" and gp_approx == "vecchia_latent":

    • "vadu" (= default): (B^T * (D^-1 + W) * B) as preconditioner for inverting (B^T * D^-1 * B + W), where B^T * D^-1 * B approx= Sigma^-1

    • "fitc": FITC / modified predictive process preconditioner for inverting (B^-1 * D * B^-T + W^-1)

    • "pivoted_cholesky": (Lk * Lk^T + W^-1) as preconditioner for inverting (B^-1 * D * B^-T + W^-1), where Lk is a low-rank pivoted Cholesky approximation for Sigma and B^-1 * D * B^-T approx= Sigma

    • "incomplete_cholesky": zero fill-in incomplete (reverse) Cholesky factorization of (B^T * D^-1 * B + W) using the sparsity pattern of B^T * D^-1 * B approx= Sigma^-1

  • Options for likelihood != "gaussian" and gp_approx == "full_scale_vecchia":

    • "fitc" ( = default): FITC / modified predictive process preconditioner

    • "vifdu": VIF with diagonal update preconditioner

  • Options for likelihood == "gaussian" and gp_approx == "full_scale_tapering":

    • "fitc" (= default): modified predictive process preconditioner

    • "none": no preconditioner

• fitc_piv_chol_preconditioner_rank (integer ): Rank of the FITC and pivoted Cholesky decomposition preconditioners for iterative methods for Vecchia and VIF approximations (for full_scale_tapering, the same inducing points as in the approximation as used). If fitc_piv_chol_preconditioner_rank = -999, internal default values are used:

  • 200 for the FITC preconditioner

  • 50 for the pivoted Cholesky decomposition preconditioner

• convergence_criterion: string (default = "relative_change_in_log_likelihood", only relevant for "gradient_descent", "fisher_scoring", and "newton"). The convergence criterion used for terminating the optimization algorithm. Options: "relative_change_in_log_likelihood" or "relative_change_in_parameters". If convergence_criterion = "default", internal default values are used

• lr_cov: numeric (default = 0.1 for "gradient_descent" and 1. otherwise, relevant for "gradient_descent", "fisher_scoring", "newton", and "lbfgs"). Initial learning rate for covariance parameters if a gradient-based optimization method is used

  • If lr_cov = -999, internal default values are used (0.1 for "gradient_descent" and 1. otherwise)

  • If there are additional auxiliary parameters for non-Gaussian likelihoods, 'lr_cov' is also used for those

  • For "gradient_descent", "fisher_scoring", and "newton", this is the learning rate by which the search direction is multiplied. It can be decreased, and increased again, during the optimization

  • For "lbfgs" and "lbfgs_linesearch_nocedal_wright", this is not a learning rate by which every update is multiplied. It only determines the first trial step, which is the negative gradient divided by its norm and multiplied by 'lr_cov', i.e., a step of length 'lr_cov'. The subsequent steps are determined by the approximate Hessian of the L-BFGS algorithm and the line search. 'lr_cov' has no effect when the approximate Hessian of a previous call is reused (see 'reuse_learning_rates_gp_model' of the GPBoost algorithm)

• lr_coef: numeric (default = 0.1, only relevant for "gradient_descent", "fisher_scoring", and "newton"). Learning rate for fixed effect regression coefficients if gradient descent is used. If lr_coef = -999, internal default values are used

• use_nesterov_acc: boolean (default = TRUE, only relevant for "gradient_descent"). If TRUE Nesterov acceleration is used. This is used only for gradient descent

• nesterov_schedule_version: integer (Default = 0, only relevant for "gradient_descent"). Which version of Nesterov schedule should be used. If nesterov_schedule_version = -999, internal default values are used.

• acc_rate_coef: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for regression coefficients (if there are any) for Nesterov acceleration. If acc_rate_coef = -999, internal default values are used

• acc_rate_cov: numeric (default = 0.5, only relevant for "gradient_descent"). Acceleration rate for covariance parameters for Nesterov acceleration. If acc_rate_cov = -999, internal default values are used

• momentum_offset: integer (Default = 2, only relevant for "gradient_descent"). Number of iterations for which no momentum is applied in the beginning. If momentum_offset = -999, internal default values are used.

• m_lbfgs: integer (Default = 6). Number of corrections to approximate the inverse Hessian matrix for the "lbfgs" optimizer. If m_lbfgs = -999, internal default values are used

• cold_restart_lbfgs: boolean (Default = TRUE). If TRUE, restarts of the "lbfgs" optimizers are "cold" restarts: the regression coefficients, the auxiliary parameters, and the modes of the Laplace approximations are reset to their initial values and only the covariance parameters are kept. Such restarts can escape a local optimum in which the optimizer is stuck. If FALSE, "warm" restarts are done: the optimization simply continues from the current parameters with a re-initialized approximate Hessian. Only relevant if max_num_restarts_lbfgs > 0. Note that every cold restart is a full optimization with up to maxit iterations, and the number of iterations reported is the sum over all restarts

• max_num_restarts_lbfgs: integer (Default = 0). Maximal number of restarts of the "lbfgs" optimizers after they have terminated. A restart continues from the current parameters with a re-initialized approximate Hessian. What else is reset depends on cold_restart_lbfgs. Restarts are only done as long as the negative log-likelihood still decreases by more than delta_rel_conv and as long as the total number of iterations is below maxit. This can help when the line search of "lbfgs" fails and the optimizer thus terminates prematurely without having converged. This can happen, e.g., for the "asymmetric_laplace" (= quantile regression) likelihood, whose approximate marginal likelihood is not smooth

• reset_modes_line_search_lbfgs: boolean (Default = FALSE). If TRUE, the mode finding of the Laplace approximation starts at the modes of the current parameters for every trial point of the line search of the "lbfgs" optimizer, and not at the modes of the previous trial point. Only relevant for non-Gaussian likelihoods and optimizer_cov = "lbfgs" or "lbfgs_linesearch_nocedal_wright"

• delta_conv_mode_finding: numeric (Default = 1E-8, and 1E-10 for the likelihood "gaussian_heteroscedastic_fixed_and_random"). Convergence tolerance in mode finding algorithm for Laplace approximation for non-Gaussian likelihoods. If delta_conv_mode_finding = -999, internal default values are used

A GPModel

A string specifying the type of Vecchia approximation used for making predictions. Default value if vecchia_pred_type = NULL: "order_obs_first_cond_obs_only". Available options:

• "order_obs_first_cond_obs_only": Vecchia approximation for the observable process and observed training data is ordered first and the neighbors are only observed training data points

• "order_obs_first_cond_all": Vecchia approximation for the observable process and observed training data is ordered first and the neighbors are selected among all points (training + prediction)

• "latent_order_obs_first_cond_obs_only": Vecchia approximation for the latent process and observed data is ordered first and neighbors are only observed points

• "latent_order_obs_first_cond_all": Vecchia approximation for the latent process and observed data is ordered first and neighbors are selected among all points

• "order_pred_first": Vecchia approximation for the observable process and prediction data is ordered first for making predictions. This option is only available for Gaussian likelihoods

an integer specifying the number of neighbors for the Vecchia approximation for making predictions. Default value if NULL: num_neighbors_pred = 2 * num_neighbors

a numeric specifying the tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithms when being used for prediction. This is the tolerance of the "absolute" criterion only, the "relative" criterion is controlled by cg_rel_tol_pred and cg_abs_tol_pred. Default value if NULL: 1e-3

a string specifying the stopping rule of conjugate gradient algorithms when being used for prediction: "absolute" (stop when ||r||_2 < cg_delta_conv_pred) or "relative" (stop when ||r||_2 <= max(cg_abs_tol_pred, cg_rel_tol_pred * ||b||_2)). Default value if NULL: the value of cg_convergence_criterion used for the parameter estimation

a numeric specifying the relative tolerance of the "relative" stopping rule when being used for prediction. Default value if NULL: the value of cg_rel_tol used for the parameter estimation. This is inherited independently of cg_convergence_criterion_pred

a numeric specifying the absolute tolerance floor of the "relative" stopping rule when being used for prediction. This makes the rule robust for right-hand sides that are zero or very small. Default value if NULL: the value of cg_abs_tol used for the parameter estimation. This is inherited independently of the other two prediction options

an integer specifying the number of samples when simulation is used for calculating predictive variances Internal default values if NULL:

• 500 for grouped random effects

• 1000 for gp_approx = "vecchia" and gp_approx = "full_scale_tapering"

• 100 for gp_approx = "full_scale_vecchia"

an integer specifying the rank of the matrix for approximating predictive covariances obtained using the Lanczos algorithm Default value if NULL: 1000

A vector or matrix with elements being group levels for which predictions are made (if there are grouped random effects in the GPModel)

A vector or matrix with covariate data for grouped random coefficients (if there are some in the GPModel)

A matrix with prediction coordinates (=features) for Gaussian process (if there is a GP in the GPModel)

A vector or matrix with covariate data for Gaussian process random coefficients (if there are some in the GPModel)

A vector with elements indicating the realizations of random effects / Gaussian processes for which predictions are made (set to NULL if you have not specified this when creating the GPModel)

A matrix with prediction covariate data for the fixed effects linear regression term (if there is one in the GPModel)

A GPModel

A string specifying the type of Vecchia approximation used for making predictions. Default value if vecchia_pred_type = NULL: "order_obs_first_cond_obs_only". Available options:

• "order_obs_first_cond_obs_only": Vecchia approximation for the observable process and observed training data is ordered first and the neighbors are only observed training data points

• "order_obs_first_cond_all": Vecchia approximation for the observable process and observed training data is ordered first and the neighbors are selected among all points (training + prediction)

• "latent_order_obs_first_cond_obs_only": Vecchia approximation for the latent process and observed data is ordered first and neighbors are only observed points

• "latent_order_obs_first_cond_all": Vecchia approximation for the latent process and observed data is ordered first and neighbors are selected among all points

• "order_pred_first": Vecchia approximation for the observable process and prediction data is ordered first for making predictions. This option is only available for Gaussian likelihoods

an integer specifying the number of neighbors for the Vecchia approximation for making predictions. Default value if NULL: num_neighbors_pred = 2 * num_neighbors

a numeric specifying the tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithms when being used for prediction. This is the tolerance of the "absolute" criterion only, the "relative" criterion is controlled by cg_rel_tol_pred and cg_abs_tol_pred. Default value if NULL: 1e-3

a string specifying the stopping rule of conjugate gradient algorithms when being used for prediction: "absolute" (stop when ||r||_2 < cg_delta_conv_pred) or "relative" (stop when ||r||_2 <= max(cg_abs_tol_pred, cg_rel_tol_pred * ||b||_2)). Default value if NULL: the value of cg_convergence_criterion used for the parameter estimation

a numeric specifying the relative tolerance of the "relative" stopping rule when being used for prediction. Default value if NULL: the value of cg_rel_tol used for the parameter estimation. This is inherited independently of cg_convergence_criterion_pred

a numeric specifying the absolute tolerance floor of the "relative" stopping rule when being used for prediction. This makes the rule robust for right-hand sides that are zero or very small. Default value if NULL: the value of cg_abs_tol used for the parameter estimation. This is inherited independently of the other two prediction options

an integer specifying the number of samples when simulation is used for calculating predictive variances Internal default values if NULL:

• 500 for grouped random effects

• 1000 for gp_approx = "vecchia" and gp_approx = "full_scale_tapering"

• 100 for gp_approx = "full_scale_vecchia"

an integer specifying the rank of the matrix for approximating predictive covariances obtained using the Lanczos algorithm Default value if NULL: 1000

A vector or matrix with elements being group levels for which predictions are made (if there are grouped random effects in the GPModel)

A vector or matrix with covariate data for grouped random coefficients (if there are some in the GPModel)

A matrix with prediction coordinates (=features) for Gaussian process (if there is a GP in the GPModel)

A vector or matrix with covariate data for Gaussian process random coefficients (if there are some in the GPModel)

A vector with elements indicating the realizations of random effects / Gaussian processes for which predictions are made (set to NULL if you have not specified this when creating the GPModel)

A matrix with prediction covariate data for the fixed effects linear regression term (if there is one in the GPModel)

Object of class gpb.Dataset

other parameters

the name of the field to get

the specific field of information to set

Object of class gpb.Dataset

other parameters (currently not used)

an integer vector of indices of rows needed

a GPModel

A boolean. If TRUE, (approximate) standard errors are calculated (= square root of diagonal of the inverse Fisher information for Gaussian likelihoods and square root of diagonal of a numerically approximated inverse Hessian for non-Gaussian likelihoods). For non-Gaussian likelihoods, no standard errors of the covariance and auxiliary parameters are calculated if the model is used in the GPBoost algorithm, since they require the fixed effects given by the tree ensemble

(not used, ignore this, simply here that there is no CRAN warning)