GPModel objectCreate a GPModel which contains a Gaussian process and / or mixed effects model with grouped random effects
GPModel(likelihood = "gaussian", group_data = NULL,
group_rand_coef_data = NULL, ind_effect_group_rand_coef = NULL,
drop_intercept_group_rand_effect = NULL, gp_coords = NULL,
gp_rand_coef_data = NULL, cov_function = "matern", cov_fct_shape = 1.5,
gp_approx = "none", num_parallel_threads = NULL, GPU_use = FALSE,
matrix_inversion_method = "default", weights = NULL,
likelihood_learning_rate = 1, cov_fct_taper_range = 1,
cov_fct_taper_shape = 1, num_neighbors = NULL,
vecchia_ordering = "random", ind_points_selection = "kmeans++",
num_ind_points = NULL, cover_tree_radius = 1, seed = 0L,
cluster_ids = NULL, likelihood_additional_param = NULL,
fidelity_specific_mean = TRUE, num_data = NULL, free_raw_data = FALSE,
vecchia_approx = NULL, vecchia_pred_type = NULL,
num_neighbors_pred = NULL)A GPModel containing ontains a Gaussian process and / or mixed effects model with grouped random effects
A string specifying the likelihood function (distribution) of the response variable.
Available options:
"gaussian"
"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"
"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"
"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
"gamma": Gamma likelihood with a log link function
"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
"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 The GPD/EGPD likelihoods require finite y > 0. Response means exist for shape < 1 and response variances for shape < 0.5.
"lognormal": Log-normal likelihood with a log link function
"beta" : Beta likelihood with a logit link function (parametrization of Ferrari and Cribari-Neto, 2004)
"t": t-distribution (e.g., for robust regression). The default approximation is Fisher-Laplace: Fisher information is used for both mode finding and determinant evaluation.
"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, aliases: "asymmetric_laplace", "quantile_regression"
The quantile must be supplied through the likelihood_additional_param parameter and must be strictly between 0 and 1
"hurdle_<base>" and "zero_inflated_<base>": Two-part likelihoods for response variables with an excess probability 'p0' of exact zeros. They combine a point mass 'p0' at zero with a base distribution for the remaining mass '1 - p0'. Use "hurdle_<base>" when the base has support 'y > 0' (positive continuous responses) and "zero_inflated_<base>" for counts (where the base can itself generate additional zeros). In both cases 'exp(F(X) + Zb)' is the mean or scale parameter of the (non-structural) base component - not the unconditional response mean - so E(y) = (1 - p0) * base_mean. The structural-zero probability 'p0' is estimated jointly with the base auxiliary parameters. Currently supported variants:
Hurdle (positive continuous base): "hurdle_gamma" (aux: shape), "hurdle_lognormal" (aux: log_variance), and the extreme-value bases "hurdle_gpd", "hurdle_egpd_power", "hurdle_egpd_power_mixture", "hurdle_egpd_beta", "hurdle_egpd_power_beta" (same aux parameters as the corresponding non-hurdle GPD/EGPD likelihoods, plus 'p0'). The alias "zero_inflated_gamma" maps to "hurdle_gamma".
Zero-inflated counts (integer y >= 0): "zero_inflated_poisson", "zero_inflated_negative_binomial" (aux: shape; aka "zero_inflated_nbinom2"), "zero_inflated_negative_binomial_1" (aux: dispersion; aka "zero_inflated_nbinom1"). The unsuffixed names default to combined Fisher-Laplace: the exact log-likelihood score and Fisher information (quasi-Fisher information for NB1) are used for mode finding, while the observed Hessian and its derivatives are used for the Laplace determinant. Append "_laplace" for observed-Hessian Newton mode finding, or "_fisher_laplace" to use Fisher (quasi-Fisher for NB1) information for both mode finding and determinant evaluation.
For both the hurdle (positive continuous) and the zero-inflated count families, the structural-zero probability can alternatively be modeled as a logistic regression on the covariates by inserting "regression" after the family prefix in the likelihood name (e.g. "hurdle_regression_gamma", "hurdle_regression_lognormal", "hurdle_regression_gpd", "zero_inflated_regression_poisson", "zero_inflated_regression_negative_binomial"). 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".
"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 (covariates and / or the GPBoost tree-boosting algorithm; no random effects / GPs for sigma), while mu = F(X) + Zb is related to both fixed and random effects. lambda is then the only (auxiliary) parameter that is estimated
"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)
"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.
"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 (covariates and / or the GPBoost tree-boosting algorithm; no random effects / GPs for the variance). Fisher-Laplace is the default and currently the only implemented approximation.
Note: the first lines in the likelihoods source file contain additional comments on the specific parametrizations used
Note: 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) ). For H = 0.5, this corresponds to Brownian motion (-> see the 'estimate_cov_par_index' argument)
"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] \)
"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 "space_time_gneiting" and "ar1_mf_<base>", neighbors are 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 (2025) 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; see Gyger, Furrer, and Sigrist (2024) for more details
"full_scale_tapering": Full-scale approximation combining an inducing point / predictive process approximation with tapering on the residual process; see Gyger, Furrer, and Sigrist (2024) 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, all available threads are used
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": 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.
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 likelihood = "tweedie_fixed_p"; 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.
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)
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
Fabio Sigrist
# See https://github.com/fabsig/GPBoost/tree/master/R-package for more examples
data(GPBoost_data, package = "gpboost")
#--------------------Grouped random effects model: single-level random effect----------------
gp_model <- GPModel(group_data = group_data[,1], likelihood="gaussian")
#--------------------Gaussian process model----------------
gp_model <- GPModel(gp_coords = coords, cov_function = "matern", cov_fct_shape = 1.5,
likelihood="gaussian")
#--------------------Combine Gaussian process with grouped random effects----------------
gp_model <- GPModel(group_data = group_data,
gp_coords = coords, cov_function = "matern", cov_fct_shape = 1.5,
likelihood="gaussian")
Run the code above in your browser using DataLab