gpboost.GPModel

class gpboost.GPModel(likelihood='gaussian', group_data=None, group_rand_coef_data=None, ind_effect_group_rand_coef=None, drop_intercept_group_rand_effect=None, gp_coords=None, gp_rand_coef_data=None, cov_function='matern', cov_fct_shape=1.5, gp_approx='none', num_parallel_threads=None, GPU_use=False, matrix_inversion_method='default', weights=None, likelihood_learning_rate=1.0, cov_fct_taper_range=1.0, cov_fct_taper_shape=1.0, num_neighbors=None, vecchia_ordering='random', ind_points_selection='kmeans++', num_ind_points=None, cover_tree_radius=1.0, seed=0, cluster_ids=None, num_data=None, likelihood_additional_param=None, free_raw_data=False, model_file=None, model_dict=None, vecchia_approx=None, vecchia_pred_type=None, num_neighbors_pred=None)[source]

Bases: object

Class for random effects model (Gaussian process, grouped random effects, mixed effects models, etc.)

Authors:: Fabio Sigrist

__init__(likelihood='gaussian', group_data=None, group_rand_coef_data=None, ind_effect_group_rand_coef=None, drop_intercept_group_rand_effect=None, gp_coords=None, gp_rand_coef_data=None, cov_function='matern', cov_fct_shape=1.5, gp_approx='none', num_parallel_threads=None, GPU_use=False, matrix_inversion_method='default', weights=None, likelihood_learning_rate=1.0, cov_fct_taper_range=1.0, cov_fct_taper_shape=1.0, num_neighbors=None, vecchia_ordering='random', ind_points_selection='kmeans++', num_ind_points=None, cover_tree_radius=1.0, seed=0, cluster_ids=None, num_data=None, likelihood_additional_param=None, free_raw_data=False, model_file=None, model_dict=None, vecchia_approx=None, vecchia_pred_type=None, num_neighbors_pred=None)[source]

Initialize a GPModel.

Parameters:

likelihood (string, optional (default="gaussian")) –
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
- ”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)
- ”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 can be set via the ‘likelihood_additional_param’ parameter. The default is quantile = 0.5
- ”zero_inflated_gamma”:
  
  Zero-inflated gamma likelihood. The log-transformed mean of the response variable equals the sum of fixed and random effects, E(y) = mu = exp(F(X) + Zb), and the rate parameter equals (1-p0) * gamma / mu, where p0 is the zero-inflation probability and gamma the shape parameter. I.e., the rate parameter depends on F(X) + Zb, and p0 and gamma are (univariate auxiliary) parameters that are estimated. Note that E(y) = mu above refers the the mean of the entire distribution and not just the positive part
- ”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”
- ”gaussian_heteroscedastic”:
  
  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
- 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
group_data (numpy array or pandas DataFrame with numeric or string data or None, optional (default=None)) – Either a vector or a matrix whose columns are categorical grouping variables. The elements are group levels defining grouped random effects. The number of columns corresponds to the number of grouped (intercept) random effects
group_rand_coef_data (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Covariate data for grouped random coefficients
ind_effect_group_rand_coef (list, numpy 1-D array, pandas Series / one-column DataFrame with integer data or None, optional (default=None)) – 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, [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’
drop_intercept_group_rand_effect (list, numpy 1-D array, pandas Series / one-column DataFrame with bool data or None, optional (default=None)) – 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
gp_coords (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Coordinates (= inputs / features) for defining Gaussian processes
gp_rand_coef_data (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Covariate data for Gaussian process random coefficients
cov_function (string, optional (default="exponential")) –
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 dimension / column 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 dimension / column of ‘gp_coords’
- ”powered_exponential”:
  
  Powered exponential covariance function with the exponent specified by ‘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’) = (sigma2 / 2) * ( (s_1^2 + sum_{k=2}^d (s_k / l_k)^2)^H + (s’_1^2 + sum_{k=2}^d (s’_k / l_k)^2)^H - ((s_1 - s’_1)^2 + sum_{k=2}^d ((s_k - s’_k) / l_k)^2)^H )
cov_fct_shape (float, optional (default=0.)) – 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
gp_approx (string, optional (default="none")) –
Specifies the use of a large data approximation for Gaussian processes. Available options:
- ”none”:
  
  No approximation
- ”vecchia”:
  
  Vecchia approximation; see Sigrist (2022, JMLR) for more details
- ”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
num_parallel_threads (integer, optional (default=None)) – The number of parallel threads for OMP. If num_parallel_threads=None, all available threads are used
GPU_use (bool, optional (default=False).) – If TRUE, GPU acceleration will be used if supported.
matrix_inversion_method (string, optional (default="default")) –
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, 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)
weights (list, numpy 1-D array, pandas Series / one-column DataFrame or None, optional (default=None)) – Sample weights. Note that this affects both the random and fixed effects components.
likelihood_learning_rate (float, optional (default=1.)) – A learning rate for the likelihood for generalized Bayesian inference (only non-Gaussian likelihoods)
cov_fct_taper_range (float, optional (default=1.)) – Range parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)
cov_fct_taper_shape (float, optional (default=0.)) – Shape (=smoothness) parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)
num_neighbors (integer, optional) –
Number of neighbors for the Vecchia approximation. Internal default values if None:
- 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
vecchia_ordering (string, optional (default="random")) –
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)
ind_points_selection (string, optional (default="kmeans++")) –
Specifies 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
num_ind_points (integer, optional) –
Number of inducing points / knots for FITC, full_scale_tapering, and VIF approximations. Internal default values if None:
- 500 for gp_approx = “FITC” and gp_approx = “full_scale_tapering”
- 200 for gp_approx = “full_scale_vecchia”
cover_tree_radius (float, optional (default=1.)) – The radius (= “spatial resolution”) for the cover tree algorithm
seed (integer, optional (default=0)) – The seed used for model creation (e.g., random ordering in Vecchia approximation)
cluster_ids (list, numpy 1-D array, pandas Series / one-column DataFrame with numeric or string data or None, optional (default=None)) – The elements indicate independent realizations of random effects / Gaussian processes (same values = same process realization)
num_data (integer, optional (default=None)) – The number of samples. This is only used for iid models.
likelihood_additional_param (float, optional (default=1.)) –
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 = None’, the following internal default values are used:
- df = 2 for ‘likelihood = “t_fix_df”’
- quantile = 0.5 for likelihood = “asymmetric_laplace”
free_raw_data (bool, optional (default=False)) – If True, the data (groups, coordinates, covariate data for random coefficients) is freed in Python after initialization
model_file (string or None, optional (default=None)) – Path to the model file.
model_dict (dict or None, optional (default=None)) – Dict with model file
vecchia_approx (bool or None, discontinued (default=None)) – This is discontinued. Use gp_approx = “none” instead
vecchia_pred_type (string, optional (default=None)) – The type of Vecchia approximation used for making predictions. This is discontinued here. Use the function ‘set_prediction_data’ to specify this
num_neighbors_pred (integer or None, optional (default=None)) – The number of neighbors for making predictions. This is discontinued here. Use the function ‘set_prediction_data’ to specify this

Example

>>> # Grouped random effects model
>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> # Gaussian process model
>>> gp_model = gpb.GPModel(gp_coords=coords, cov_function="matern", cov_fct_shape=1.5, likelihood="gaussian")

Methods

`__init__`([likelihood, group_data, ...])	Initialize a GPModel.
`fit`(y[, X, params, offset, fixed_effects])	Fit / estimate a GPModel by maximizing the marginal likelihood
`get_aux_pars`([format_pandas])	Get (estimated) auxiliary (additional) parameters of the likelihood such as the shape parameter of a gamma or a negative binomial distribution.
`get_coef`([std_err, format_pandas])	Get (estimated) linear regression coefficients
`get_cov_pars`([std_err, format_pandas])	Get (estimated) covariance parameters
`get_current_neg_log_likelihood`()	Get the current value of the negative log-likelihood
`model_to_dict`([include_response_data])	Convert a GPModel to a dict for saving.
`neg_log_likelihood`(cov_pars, y[, ...])	Evaluate the negative log-likelihood.
`predict`([predict_response, predict_var, ...])	Make predictions for a GPModel.
`predict_training_data_random_effects`([...])	Predict ("estimate") training data random effects.
`save_model`(filename)	Save a GPModel to file.
`set_optim_params`(params)	Set parameters for estimation of the covariance parameters.
`set_prediction_data`([vecchia_pred_type, ...])	Set the data required for making predictions with a GPModel.
`summary`([std_err])	Print summary of fitted model parameters.

fit(y, X=None, params=None, offset=None, fixed_effects=None)[source]

Fit / estimate a GPModel by maximizing the marginal likelihood

Parameters:

y (list, numpy 1-D array, pandas Series / one-column DataFrame or None, optional (default=None)) – Response variable data
X (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Covariate data for the fixed effects linear regression term (if there is one)
params (dict or None, optional (default=None)) –
Parameters for the estimation / optimization
- tracebool, optional (default = False)
  If True, information on the progress of the parameter optimization is printed.
- init_cov_parsnumpy array or pandas DataFrame, optional (default = None)
  Initial values for covariance parameters of Gaussian process and random effects (can be None). The order is the 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 range of the Gaussian process. If there are multiple Gaussian processes, then the variances and ranges follow alternatingly. If ‘init_cov_pars = None’, an internatl 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_coefnumpy array or pandas DataFrame, optional (default = None)
  Initial values for the regression coefficients (if there are any, can be None)
- init_aux_parsnumpy array or pandas DataFrame, optional (default = None)
  Initial values for additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative binomial likelihood) (can be None).
- estimate_cov_par_indexlist, numpy 1-D array, pandas Series / one-column DataFrame with integer data or None, optional (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”: [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’).
- estimate_aux_parsbool, (default = True)
  If True, any additional parameters for non-Gaussian likelihoods are also estimated (e.g., shape parameter of a gamma or negative binomial likelihood).
- optimizer_covstring, optional (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_coefstring, optional (default = “wls” for Gaussian data 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.
- maxitinteger, optional (default = 1000)
  Maximal number of iterations for optimization algorithm.
- delta_rel_convdouble, optional (default = 1e-6 except for “nelder_mead” for which the default is 1e-8)
  Convergence tolerance. The algorithm stops if the relative change in eiher the (approximate) log-likelihood or the parameters is below this value. If < 0, internal default values are used. Default = 1e-6 except for “nelder_mead” for which the default is 1e-8.
- cg_max_num_it: integer, optional (default = 1000)
  Maximal number of iterations for conjugate gradient algorithms.
- cg_max_num_it_tridiag: integer, optional (default = 1000)
  Maximal number of iterations for conjugate gradient algorithm when being run as Lanczos algorithm for tridiagonalization.
- cg_delta_conv: double, optional (default = 1e-2)
  Tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithm when being used for parameter estimation.
- num_rand_vec_trace: integer, optional (default = 50)
  Number of random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix.
- reuse_rand_vec_trace: boolean, optional (default = True)
  If true, random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix are sampled only once at the beginning of Newton’s method for finding the mode in the Laplace approximation and are then reused in later trace approximations. Otherwise they are sampled every time a trace is calculated.
- seed_rand_vec_trace: integer, optional (default = 1)
  Seed number to generate random vectors (e.g., Rademacher).
- cg_preconditioner_type: string, optional
  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”: modified predictive process preconditioner for inverting (B^-1 * D * B^-T + W^-1)
  
  ”pivoted_cholesky” (= default): (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, optional
  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). Internal default values if None or < 0:
  
  200 for the FITC preconditioner
  
  50 for the pivoted Cholesky decomposition preconditioner
- convergence_criterionstring, optional (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”.
- lr_covdouble, optional (default = 0.1 for “gradient_descent” and 1. otherwise, only relevant for “gradient_descent”, “fisher_scoring”, and “newton”)
  Initial learning rate for covariance parameters if a gradient-based optimization method is used.
  
  If lr_cov < 0, 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 “lbfgs”, this is divided by the norm of the gradient in the first iteration.
- lr_coefdouble, optional (default = 0.1, only relevant for “gradient_descent”, “fisher_scoring”, and “newton”)
  Learning rate for fixed effect regression coefficients.
- use_nesterov_accbool, optional (default = True, only relevant for “gradient_descent”)
  If True, Nesterov acceleration is used for gradient descent.
- acc_rate_covdouble, optional (default = 0.5, only relevant for “gradient_descent”)
  Acceleration rate for covariance parameters for Nesterov acceleration.
- acc_rate_coefdouble, optional (default = 0.5, only relevant for “gradient_descent”)
  Acceleration rate for regression coefficients (if there are any) for Nesterov acceleration.
- momentum_offsetinteger, optional (default = 2, only relevant for “gradient_descent”)
  Number of iterations for which no momentum is applied in the beginning.
- m_lbfgsinteger, optional (default = 6)
  Number of corrections to approximate the inverse Hessian matrix for the “lbfgs” optimizer
- delta_conv_mode_findingdouble, optional (default = 1e-8)
  Convergence tolerance in mode finding algorithm for Laplace approximation for non-Gaussian likelihoods
offset (numpy 1-D array or None, optional (default=None)) – 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.
fixed_effects (numpy 1-D array or None, optional (default=None)) – This is discontinued. Use the renamed equivalent argument ‘offset’ instead.

Example

>>> # Grouped random effects model
>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> gp_model.fit(y=y, X=X)
>>> # Gaussian process model
>>> gp_model = gpb.GPModel(gp_coords=X, cov_function="matern", cov_fct_shape=1.5, likelihood="gaussian")
>>> gp_model.fit(y=y)

get_aux_pars(format_pandas=True)[source]

Get (estimated) auxiliary (additional) parameters of the likelihood such as the shape parameter of a gamma or a negative binomial distribution. Some likelihoods (e.g., bernoulli_logit or poisson) have no auxiliary parameters

Parameters:: format_pandas (bool (default=True)) – If True, a pandas DataFrame is returned, otherwise a numpy array is returned
Returns:: result – auxiliary (additional) parameters of the likelihood
Return type:: numpy array or pandas DataFrame

Example

>>> gp_model = gpb.GPModel(group_data=group, likelihood="gamma")
>>> gp_model.fit(y=y, X=X)
>>> gp_model.get_aux_pars()

get_coef(std_err=False, format_pandas=True)[source]

Get (estimated) linear regression coefficients

Parameters:

std_err (bool (default=False)) – If True, (approximate) standard errors are calculated
format_pandas (bool (default=True)) – If True, a pandas DataFrame is returned, otherwise a numpy array is returned

Returns:

result – (estimated) linear regression coefficients and standard errors (if std_err=True)

Return type:

numpy array or pandas DataFrame

Example

>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> gp_model.fit(y=y, X=X)
>>> gp_model.get_cov_pars()

get_cov_pars(std_err=False, format_pandas=True)[source]

Get (estimated) covariance parameters

Parameters:

std_err (bool (default=False)) – If True, (approximate) standard errors are calculated
format_pandas (bool (default=True)) – If True, a pandas DataFrame is returned, otherwise a numpy array is returned

Returns:

result – (estimated) covariance parameters and standard errors (if std_err=True)

Return type:

pandas DataFrame

Example

>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> gp_model.fit(y=y, X=X)
>>> gp_model.get_cov_pars()

get_current_neg_log_likelihood()[source]

Get the current value of the negative log-likelihood

Returns:: result
Return type:: the current value of the negative log-likelihood

Example

>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> gp_model.fit(y=y)
>>> gp_model.get_current_neg_log_likelihood()

model_to_dict(include_response_data=True)[source]

Convert a GPModel to a dict for saving.

Parameters:: include_response_data (bool (default=False)) – If true, the response variable data is also included in the dict
Returns:: model_dict – GPModel in dict format.
Return type:: dict

neg_log_likelihood(cov_pars, y, fixed_effects=None, aux_pars=None)[source]

Evaluate the negative log-likelihood. If there is a linear fixed effects predictor term, this needs to be calculated “manually” prior to calling this function (see example below)

Parameters:

cov_pars (list, numpy 1-D array, pandas Series / one-column DataFrame or None, optional (default=None)) – Covariance parameters of Gaussian process and random effects
y (list, numpy 1-D array, pandas Series / one-column DataFrame or None, optional (default=None)) – Response variable data
fixed_effects (numpy 1-D array or None, optional (default=None)) – A vector with fixed effects, e.g., containing a linear predictor. The length of this vector needs to equal the number of training data points.
aux_pars (numpy array or pandas DataFrame, optional (default = None)) – Additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative binomial likelihood) (can be None)

Returns:

result

Return type:

the value of the negative log-likelihood

Example

>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> coef = [0, 0.1]
>>> fixed_effects = X.dot(coef)
>>> gp_model.neg_log_likelihood(y=y, cov_pars=[1.,1.], fixed_effects=fixed_effects)

predict(predict_response=True, predict_var=False, predict_cov_mat=False, sample_posterior=False, sample_prior=False, num_post_samples=100, num_prior_samples=100, y=None, cov_pars=None, group_data_pred=None, group_rand_coef_data_pred=None, gp_coords_pred=None, gp_rand_coef_data_pred=None, cluster_ids_pred=None, X_pred=None, use_saved_data=False, offset=None, offset_pred=None, fixed_effects=None, fixed_effects_pred=None, vecchia_pred_type=None, num_neighbors_pred=None)[source]

Make predictions for a GPModel.

Parameters:

predict_response (bool (default=True)) – If True, the response variable (label) is predicted, otherwise the latent random effects
predict_var (bool (default=False)) – If True, the (posterior) predictive variances are calculated
predict_cov_mat (bool (default=False)) – If True, the (posterior) predictive covariance is calculated in addition to the (posterior) predictive mean
sample_posterior (bool (default=False)) – If True, samples from the posterior are drawn
sample_prior (bool (default=False)) – If True, samples from the prior are drawn
num_post_samples (integer (default=100)) – Number of posterior samples to draw if ‘sample_posterior=True’
num_prior_samples (integer (default=100)) – Number of prior samples to draw if ‘sample_prior=True’
y (list, numpy 1-D array, pandas Series / one-column DataFrame or None, optional (default=None)) – Observed response variable data (can be None, e.g. when the model has been estimated already and the same data is used for making predictions)
cov_pars (numpy array or None, optional (default = None)) – 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 estimated ones
group_data_pred (numpy array or pandas DataFrame with numeric or string data or None, optional (default=None)) – The elements are group levels for which predictions are made (if there are any grouped random effects in the model)
group_rand_coef_data_pred (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Covariate data for grouped random coefficients (if there are some in the model)
gp_coords_pred (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Prediction coordinates (=features) for Gaussian process (if there is a GP in the model)
gp_rand_coef_data_pred (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Covariate data for Gaussian process random coefficients (if there are some in the model)
cluster_ids_pred (list, numpy 1-D array, pandas Series / one-column DataFrame with numeric or string data or None, optional (default=None)) – The elements indicating independent realizations of random effects / Gaussian processes for which predictions are made (set to None if you have not specified this when creating the model)
X_pred (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Prediction covariate data for the fixed effects linear regression term (if there is one)
use_saved_data (bool (default=False)) – If True, predictions are done using a priory set data via the function ‘set_prediction_data’ (this option is not used by users directly)
offset (numpy 1-D array or None, optional (default=None)) – 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.
offset_pred (numpy 1-D array or None, optional (default=None)) – 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.
fixed_effects (numpy 1-D array or None, optional (default=None)) – This is discontinued. Use the renamed equivalent argument ‘offset’ instead
fixed_effects_pred (numpy 1-D array or None, optional (default=None)) – This is discontinued. Use the renamed equivalent argument ‘offset_pred’ instead
vecchia_pred_type (string, optional (default=None)) – The type of Vecchia approximation used for making predictions. This is discontinued here. Use the function ‘set_prediction_data’ to specify this
num_neighbors_pred (integer or None, optional (default=None)) – The number of neighbors for making predictions. This is discontinued here. Use the function ‘set_prediction_data’ to specify this

Returns:

result –

‘mu’ (first entry):

Predictive (=posterior) mean. For (generalized) linear mixed effects models, i.e., models with a linear regression term, this consists of the sum of fixed effects and random effects predictions
’cov’ (second entry):

Predictive (=posterior) covariance matrix. This is None if ‘predict_cov_mat=False’
’var’ (third entry):

Predictive (=posterior) variances. This is None if ‘predict_var=False’

Return type:

a dict with three entries having numpy arrays as values

Example

>>> # Grouped random effects model
>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> gp_model.fit(y=y, X=X)
>>> pred = gp_model.predict(X_pred=X_test, group_data_pred=group_test,
                            predict_var=True, predict_response=False)
>>> print(pred['mu']) # Predicted latent mean
>>> print(pred['var']) # Predicted latent variance
>>> # Gaussian process model
>>> gp_model = gpb.GPModel(gp_coords=X, cov_function="matern", cov_fct_shape=1.5, likelihood="gaussian")
>>> gp_model.fit(y=y)
>>> pred = gp_model.predict(X_pred=X_test, gp_coords_pred=coords_test,
>>>                         predict_var=True, predict_response=False)

predict_training_data_random_effects(predict_var=False)[source]

Predict (“estimate”) training data random effects.

Parameters:: predict_var (bool (default=False)) – If True, the (posterior) predictive variances are calculated
Returns:: result
Return type:: a matrix with predicted (“estimated”) training data random effects

Example

>>> # Grouped random effects model
>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> gp_model.fit(y=y, X=X)
>>> gp_model.predict_training_data_random_effects()
>>> # The function 'predict_training_data_random_effects' returns predicted random effects for all data points.
>>> # Unique random effects for every group can be obtained as follows
>>> first_occurences = [np.where(group==i)[0][0] for i in np.unique(group)]
>>> training_data_random_effects = all_training_data_random_effects.iloc[first_occurences]

save_model(filename)[source]

Save a GPModel to file.

Parameters:: filename (string) – Filename to save a GPModel.
Returns:: self – Returns self.
Return type:: GPModel

Example

>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> gp_model.fit(y=y, X=X)
>>> gp_model.save_model('gp_model.json')

set_optim_params(params)[source]

Set parameters for estimation of the covariance parameters.

Parameters:

params (dict or None, optional (default=None)) –

Parameters for the estimation / optimization

tracebool, optional (default = False)
If True, information on the progress of the parameter optimization is printed.

init_cov_parsnumpy array or pandas DataFrame, optional (default = None)
Initial values for covariance parameters of Gaussian process and random effects (can be None). The order is the 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 range of the Gaussian process. If there are multiple Gaussian processes, then the variances and ranges follow alternatingly. If ‘init_cov_pars = None’, an internatl 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_coefnumpy array or pandas DataFrame, optional (default = None)
Initial values for the regression coefficients (if there are any, can be None)

init_aux_parsnumpy array or pandas DataFrame, optional (default = None)
Initial values for additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative binomial likelihood) (can be None).

estimate_cov_par_indexlist, numpy 1-D array, pandas Series / one-column DataFrame with integer data or None, optional (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”: [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’).

estimate_aux_parsbool, (default = True)
If True, any additional parameters for non-Gaussian likelihoods are also estimated (e.g., shape parameter of a gamma or negative binomial likelihood).

optimizer_covstring, optional (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_coefstring, optional (default = “wls” for Gaussian data 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.

maxitinteger, optional (default = 1000)
Maximal number of iterations for optimization algorithm.

delta_rel_convdouble, optional (default = 1e-6 except for “nelder_mead” for which the default is 1e-8)
Convergence tolerance. The algorithm stops if the relative change in eiher the (approximate) log-likelihood or the parameters is below this value. If < 0, internal default values are used. Default = 1e-6 except for “nelder_mead” for which the default is 1e-8.

cg_max_num_it: integer, optional (default = 1000)
Maximal number of iterations for conjugate gradient algorithms.

cg_max_num_it_tridiag: integer, optional (default = 1000)
Maximal number of iterations for conjugate gradient algorithm when being run as Lanczos algorithm for tridiagonalization.

cg_delta_conv: double, optional (default = 1e-2)
Tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithm when being used for parameter estimation.

num_rand_vec_trace: integer, optional (default = 50)
Number of random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix.

reuse_rand_vec_trace: boolean, optional (default = True)
If true, random vectors (e.g., Rademacher) for stochastic approximation of the trace of a matrix are sampled only once at the beginning of Newton’s method for finding the mode in the Laplace approximation and are then reused in later trace approximations. Otherwise they are sampled every time a trace is calculated.

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

cg_preconditioner_type: string, optional
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”: modified predictive process preconditioner for inverting (B^-1 * D * B^-T + W^-1)

”pivoted_cholesky” (= default): (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, optional
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). Internal default values if None or < 0:

200 for the FITC preconditioner

50 for the pivoted Cholesky decomposition preconditioner

convergence_criterionstring, optional (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”.

lr_covdouble, optional (default = 0.1 for “gradient_descent” and 1. otherwise, only relevant for “gradient_descent”, “fisher_scoring”, and “newton”)
Initial learning rate for covariance parameters if a gradient-based optimization method is used.

If lr_cov < 0, 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 “lbfgs”, this is divided by the norm of the gradient in the first iteration.

lr_coefdouble, optional (default = 0.1, only relevant for “gradient_descent”, “fisher_scoring”, and “newton”)
Learning rate for fixed effect regression coefficients.

use_nesterov_accbool, optional (default = True, only relevant for “gradient_descent”)
If True, Nesterov acceleration is used for gradient descent.

acc_rate_covdouble, optional (default = 0.5, only relevant for “gradient_descent”)
Acceleration rate for covariance parameters for Nesterov acceleration.

acc_rate_coefdouble, optional (default = 0.5, only relevant for “gradient_descent”)
Acceleration rate for regression coefficients (if there are any) for Nesterov acceleration.

momentum_offsetinteger, optional (default = 2, only relevant for “gradient_descent”)
Number of iterations for which no momentum is applied in the beginning.

m_lbfgsinteger, optional (default = 6)
Number of corrections to approximate the inverse Hessian matrix for the “lbfgs” optimizer

delta_conv_mode_findingdouble, optional (default = 1e-8)
Convergence tolerance in mode finding algorithm for Laplace approximation for non-Gaussian likelihoods

Example

>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> gp_model.set_optim_params(params={"optimizer_cov": "nelder_mead", "trace": True})

set_prediction_data(vecchia_pred_type=None, num_neighbors_pred=None, cg_delta_conv_pred=None, nsim_var_pred=None, rank_pred_approx_matrix_lanczos=None, group_data_pred=None, group_rand_coef_data_pred=None, gp_coords_pred=None, gp_rand_coef_data_pred=None, cluster_ids_pred=None, X_pred=None)[source]

Set the data required for making predictions with a GPModel.

Parameters:

vecchia_pred_type (string, optional (default=None)) –
Type of Vecchia approximation used for making predictions

Default value if vecchia_pred_type = None: “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 or 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
num_neighbors_pred (integer or None, optional (default=None)) –
Number of neighbors for the Vecchia approximation for making predictions

Default value if None: num_neighbors_pred = 2 * num_neighbors
cg_delta_conv_pred (double or None, optional (default=None)) –
Tolerance level for L2 norm of residuals for checking convergence in conjugate gradient algorithm when being used for prediction

Default value if None: 1e-3
nsim_var_pred (integer or None, optional (default=None)) –
The number of samples when simulation is used for calculating predictive variances

Internal default values if None:
- 500 for grouped random effects
- 1000 for gp_approx = “vecchia” and gp_approx = “full_scale_tapering”
- 100 for gp_approx = “full_scale_vecchia”
rank_pred_approx_matrix_lanczos (integer or None, optional (default=None)) –
The rank of the matrix for approximating predictive covariances obtained using the Lanczos algorithm

Default value if None: 1000
group_data_pred (numpy array or pandas DataFrame with numeric or string data or None, optional (default=None)) – The elements are group levels for which predictions are made (if there are any grouped random effects in the model)
group_rand_coef_data_pred (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Covariate data for grouped random coefficients (if there are some in the model)
gp_coords_pred (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Prediction coordinates (=features) for Gaussian process (if there is a GP in the model)
gp_rand_coef_data_pred (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Covariate data for Gaussian process random coefficients (if there are some in the model)
cluster_ids_pred (list, numpy 1-D array, pandas Series / one-column DataFrame with numeric or string data)
None (or) – The elements indicating independent realizations of random effects / Gaussian processes for which predictions are made (set to None if you have not specified this when creating the model)
(default=None) (optional) – The elements indicating independent realizations of random effects / Gaussian processes for which predictions are made (set to None if you have not specified this when creating the model)
X_pred (numpy array or pandas DataFrame with numeric data or None, optional (default=None)) – Prediction covariate data for the fixed effects linear regression term (if there is on

Example

>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> pred = gp_model.set_prediction_data(group_data_pred=group_valid)

summary(std_err=True)[source]

Print summary of fitted model parameters.

Parameters:: std_err (bool (default=False)) – If True, (approximate) standard errors are calculated

Example

>>> gp_model = gpb.GPModel(group_data=group, likelihood="gaussian")
>>> gp_model.fit(y=y, X=X)
>>> gp_model.summary()