---
title: "Mixed effects machine learning with GPBoost for grouped and areal spatial econometric data"
subtitle: "A demo using European GDP data"
author: "Fabio Sigrist"
date: "2023-06-16"
output:
  html_document:
    toc: yes
    number_sections: yes
    theme: united
  # html_notebook: #figures are not shown when using this option
  #   toc: yes
  #   number_sections: yes
  #   theme: united
---

```{r setup, include=FALSE}
cache <- FALSE # caching does not work properly...
knitr::opts_chunk$set(echo = TRUE, cache = cache, dev = "ragg_png")
run_python <- FALSE
library(reticulate)
reticulate::use_python("C:/ProgramData/Anaconda3", require=TRUE)
```

# Introduction

The [GPBoost algorithm](https://www.jmlr.org/papers/v23/20-322.html) extends linear mixed effects and Gaussian process models by replacing the linear fixed effects function with a non-parametric non-linear function modeled using tree-boosting. This article shows how the GPBoost algorithm implemented in the [``GPBoost`` library](https://github.com/fabsig/GPBoost) can be used for modeling data with a spatial and grouped structure. We demonstrate the functionality of the ``GPBoost`` library using European GDP data which is an example of areal spatial econometric data. All results are obtained using ``GPBoost`` version 1.2.1. This demo uses the R package, but the corresponding Python package provides the same functionality.

Applying a GPBoost model (= combined tree-boosting and random effects / GP models) involves the following main steps:

1. Define a ``GPModel`` in which you specify the following:

    - A random effects model (e.g., spatial random effects, grouped random effects, combined spatial and grouped, etc.)
    - The likelihood (= distribution of the response variable conditional on fixed and random effects)

2. Create a ``gpb.Dataset`` with the response variable and fixed effects predictor variables
3. Choose tuning parameters, e.g., using the function ``gpb.grid.search.tune.parameters``
4. Train the model with the ``gpboost`` / ``gpb.train`` function
5. Interpret the trained model an/or make predictions

## Data description

The data used in this demo is European gross domestic product (GDP) data. It can downloaded from https://raw.githubusercontent.com/fabsig/GPBoost/master/data/gdp_european_regions.csv. The data was collected by Massimo Giannini, University of Rome Tor Vergata, from Eurostat and kindly provided to Fabio Sigrist for a talk at the University of Rome Tor Vergata on June 16, 2023. 

Data was collected for 242 European regions for the two years 2000 and 2021. I.e., the total number of data points is 484. The response variable is

- (log) GDP / capita.

There are four predictor variables:

- L: (log) share of employment (empl/pop)
- K: (log) fixed capital/population
- Pop: log(population)
- Edu: share of tertiary education

Further, there are centroid spatial coordinates or the regions:

- Long: region longitude
- Lat: region latitude

A spatial region ID:

- Group: the region ID (from 1 to 242)

A spatial cluster ID:

- cl: identifies the cluster the region belongs (there are two clusters)

## Data loading and short visualization

We first load the data and create a map illustrating the (log) GDP / capita over space. We create two maps: one with all data and another one when excluding some remote islands. In the commented code below, we also show how to create a spatial plot when no shape file for spatial polygons is available.

<!-- ***R*** -->
```{r load_data_r, echo=TRUE, fig.height=4, fig.width=10, message=FALSE, warning=FALSE, results=FALSE}
library(gpboost)
library(ggplot2)
library(gridExtra)
library(viridis)
library(sf)

# Load data
data <- read.csv("https://raw.githubusercontent.com/fabsig/GPBoost/master/data/gdp_european_regions.csv")
FID <- data$FID
data <- as.matrix(data[,names(data)!="FID"]) # convert to matrix since the boosting part currently does not support data.frames
covars <- c("L", "K", "pop", "edu")
# Load shape file for spatial plots
cur_tempfile <- tempfile()
download.file(url = "https://raw.githubusercontent.com/fabsig/GPBoost/master/data/shape_European_regions.zip", destfile = cur_tempfile)
out_directory <- tempfile()
unzip(cur_tempfile, exdir = out_directory)
shape <- st_read(dsn = out_directory)

# Create spatial plot of GDP
data_plot <- merge(shape, data.frame(FID = FID, y = data[,"y"]), by="FID")
p1 <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="GDP / capita (log)", option = "B") + 
   ggtitle("GDP / capita (log)")
# Sample plot excluding islands
p2 <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="GDP / capita (log)", option = "B") +
  xlim(2700000,6400000) + ylim(1500000,5200000) + 
  ggtitle("GDP / capita (log) -- excluding islands")
grid.arrange(p1, p2, ncol=2)

# # Spatial plot without a shape file
# p1 <- ggplot(data = data.frame(Lat=data[,"Lat"], Long=data[,"Long"],
#                                GDP=data[,"y"]), aes(x=Long,y=Lat,color=GDP)) +
#   geom_point(size=2, alpha=0.5) + scale_color_viridis(option = "B") + 
#   ggtitle("GDP / capita (log)")
# p2 <- ggplot(data = data.frame(Lat=data[,"Lat"], Long=data[,"Long"], 
#                                GDP=data[,"y"]), aes(x=Long,y=Lat,color=GDP)) +
#   geom_point(size=3, alpha=0.5) + scale_color_viridis(option = "B") + 
#   ggtitle("GDP / capita (log) -- Europe excluding islands") + xlim(-10,28) + ylim(35,67)
# grid.arrange(p1, p2, ncol=2)
```

# Training a GPBoost model
In the following, we use a Gaussian process model with an exponential covariance function to model spatial random effects. Additionally, we include grouped random effects for the cluster variable cl. In the ``GPBoost`` library, Gaussian process random effects are defined by the ``gp_coords`` argument and grouped random effects via the ``group_data`` argument of the ``GPModel`` constructor. The above-mentioned predictor variables are used in the fixed effects tree-ensemble function. We fit the model using the ``gpboost``, or equivalently the ``gpb.train``, function. Note that we use tuning parameters that are selected below using cross-validation.

```{r train_r, results=FALSE, cache=cache}
gp_model <- GPModel(group_data = data[, c("cl")], 
                    gp_coords = data[, c("Long", "Lat")],
                    likelihood = "gaussian", cov_function = "exponential")
params <- list(learning_rate = 0.01, max_depth = 2, num_leaves = 2^10,
               min_data_in_leaf = 10, lambda_l2 = 0)
nrounds <- 37
# gp_model$set_optim_params(params = list(trace=TRUE)) # To monitor hyperparameter estimation
boost_data <- gpb.Dataset(data = data[, covars], label = data[, "y"])
gpboost_model <- gpboost(data = boost_data, gp_model = gp_model, 
                         nrounds = nrounds, params = params, 
                         verbose = 1) # same as gpb.train
```

# Choosing tuning parameters
It is important that tuning parameters are appropriately chosen for boosting. There are no universal default values and every data set will likely need different tuning parameters. Below we show how tuning parameters can be chosen using the `` gpb.grid.search.tune.parameters`` function. We use the mean square error (``mse``) as prediction accuracy measure on the validation data. Alternatively, one can also use, e.g.,  the test negative log-likelihood (``test_neg_log_likelihood`` = default value if nothing is specified) which also takes prediction uncertainty into account. *Depending on the data set and the grid size, this can take some time. Instead of a deterministic grid search as below, one can also do a random grid search to speed things up (see ``num_try_random``).*

```{r tune_pars_r, eval=FALSE, include=TRUE}
gp_model <- GPModel(group_data = data[, "cl"], 
                    gp_coords = data[, c("Long", "Lat")],
                    likelihood = "gaussian", cov_function = "exponential")
boost_data <- gpb.Dataset(data = data[, covars], label = data[, "y"])
param_grid = list("learning_rate" = c(1,0.1,0.01), 
                  "min_data_in_leaf" = c(10,100,1000),
                  "max_depth" = c(1,2,3,5,10),
                  "lambda_l2" = c(0,1,10))
other_params <- list(num_leaves = 2^10)
set.seed(1)
opt_params <- gpb.grid.search.tune.parameters(param_grid = param_grid, params = other_params,
                                              num_try_random = NULL, nfold = 4,
                                              data = boost_data, gp_model = gp_model, 
                                              nrounds = 1000, early_stopping_rounds = 10,
                                              verbose_eval = 1, metric = "mse") # metric = "test_neg_log_likelihood"
opt_params
# ***** New best test score (l2 = 0.0255393919591794) found for the following parameter combination: learning_rate: 0.01, min_data_in_leaf: 10, max_depth: 2, lambda_l2: 0, nrounds: 37
```

# Model interpretation

## Estimated random effects model
Information on the estimated random effects model can be obtained by calling the ``summary`` function on the ``gp_model``. For Gaussian processes, ``GP_var`` is the marginal variance, i.e., the amount of spatial correlation or structure spatial variation, and ``GP_range`` is the range, or scale, parameter that measures how fast correlation decays over space. For an exponential covariance function, three times this value (approx. `r signif(3 * gp_model$get_cov_pars()[4], digits=2)` here) is the distance where the (residual) spatial correlation is essentially zero (below 5%). As the results below show, the amount of spatial correlation is relatively small since the marginal variance of `r signif(gp_model$get_cov_pars()[3], digits=2)` is small compared to the total variance of the response variable which is approx. `r signif(var(data[,"y"]), digits=2)`. Additionally the error term and the cl grouped random effects also have small variances (`r signif(gp_model$get_cov_pars()[1], digits=2)` and `r signif(gp_model$get_cov_pars()[2], digits=2)`). We thus conclude that most of the variance in the response variable is explained by the fixed effects predictor variables. 
```{r summary_r}
summary(gp_model)
```

## Spatial effect map
We can plot the estimated ("predicted") spatial random effects at the training data locations by calling the ``predict`` function on the training data; see the code below. Such a plot shows the spatial effect when factoring out the effect of the fixed effects predictor variables. Note that these spatial effects take into account both the spatial Gaussian process and the grouped region cluster random effects. If one wants to obtain only spatial random effects from the Gaussian process part, one can use the ``predict_training_data_random_effects`` function (see the commented code below). Alternatively, one can also do spatial extrapolation (="Krigging"), but this makes not a lot of sense for areal data.
```{r spatial_map_r, fig.width=10, fig.height=4, echo=TRUE, message=FALSE, warning=FALSE}
pred <- predict(gpboost_model, group_data_pred = data[1:242, c("cl")], 
                gp_coords_pred = data[1:242, c("Long", "Lat")],
                data = data[1:242, covars], predict_var = TRUE, pred_latent = TRUE)
data_plot <- merge(shape, data.frame(FID = FID, y = pred$random_effect_mean), by="FID")
plot_mu <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="Spatial effect", option = "B") +
  ggtitle("Spatial effect (mean)") + xlim(2700000,6400000) + ylim(1500000,5200000)
data_plot <- merge(shape, data.frame(FID = FID, y = pred$random_effect_cov), by="FID")
plot_sd <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="Std. dev.", option = "B") +
  ggtitle("Uncertainty (std. dev.)") + xlim(2700000,6400000) + ylim(1500000,5200000)
grid.arrange(plot_mu, plot_sd, ncol=2)

# Only spatial effecst from the Gaussian process
# rand_effs <- predict_training_data_random_effects(gp_model, predict_var = TRUE)[1:242, c("GP", "GP_var")]
# data_plot <- merge(shape, data.frame(FID = FID, y = rand_effs[,1]), by="FID")
# plot_mu <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
#   scale_fill_viridis(name="Spatial effect (mean)", option = "B") +
#   ggtitle("Spatial effect from Gausian process (mean)") + xlim(2700000,6400000) + ylim(1500000,5200000)
# data_plot <- merge(shape, data.frame(FID = FID, y = rand_effs[,2]), by="FID")
# plot_sd <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
#   scale_fill_viridis(name="Uncertainty (std. dev.)", option = "B") +
#   ggtitle("Uncertainty (std. dev.)") + xlim(2700000,6400000) + ylim(1500000,5200000)
# grid.arrange(plot_mu, plot_sd, ncol=2)
```

<!-- #   # Spatial extrapolation ("Krigging") to produce a continuos spatial map -->
<!--  #   # Define coordinates for making spatial map -->
<!--  # nx <- 50 # number locations on each axis for prediction -->
<!--  # x2 <- x1 <- rep((1:nx)/nx,nx) -->
<!--  # for(i in 1:nx) x2[((i-1)*nx+1):(i*nx)]=i/nx -->
<!--  # coord_pred <- cbind(x1,x2) -->
<!--  # # min_c <- apply(data[, c("Long", "Lat")], 2, min) # all data -->
<!--  # # max_c <- apply(data[, c("Long", "Lat")], 2, max) -->
<!--  # min_c <- c(-12, 33) # only "core" Europe, no islands -->
<!--  # max_c <- c(30, 67) -->
<!--  # coord_pred[,1] <- min_c[1] + (max_c[1] - min_c[1]) * coord_pred[,1] -->
<!--  # coord_pred[,2] <- min_c[2] + (max_c[2] - min_c[2]) * coord_pred[,2] -->
<!--  # # Other data required for making predictions -->
<!--  # X_med <- matrix(rep(apply(data[, covars],2,median), dim(coord_pred)[1]),  -->
<!--  #                 ncol=length(covars), byrow = TRUE) -->
<!--  # group_data_med <- matrix(rep(median(data[, "cl"]), dim(coord_pred)[1]),  -->
<!--  #                          ncol=1, byrow = TRUE) # Note: we ignore the region identifier for making the spatial map -->
<!--  # # Create spatial predictive map -->
<!--  # preds <- predict(gpboost_model, data = X_med, gp_coords_pred = coord_pred, -->
<!--  #                  group_data_pred = group_data_med, predict_var = TRUE, pred_latent = TRUE) -->
<!--  # plot_mu <- ggplot(data = data.frame(Long=coord_pred[,1], Lat=coord_pred[,2],  -->
<!--  #                                     spat_eff=preds$random_effect_mean), aes(x=Long,y=Lat,color=spat_eff)) + -->
<!--  #   geom_point(size=2, shape=15) + scale_color_viridis(option = "B") +  -->
<!--  #   ggtitle("Spatial effect map (mean)") -->
<!--  # plot_sd <- ggplot(data = data.frame(Long=coord_pred[,1] ,Lat=coord_pred[,2],  -->
<!--  #                                     std_dev=sqrt(preds$random_effect_cov)), aes(x=Long,y=Lat,color=std_dev)) + -->
<!--  #   geom_point(size=2, shape=15) + scale_color_viridis(option = "B") +  -->
<!--  #   ggtitle("Uncertainty (std. dev.)") -->
<!--  # grid.arrange(plot_mu, plot_sd, ncol=2) -->

## Understanding the fixed effects function
There are several tools for understanding the form of the fixed effects function. Below we consider variable importance measures, interaction measures, and dependence plots. Specifically, we look at 

- SHAP values
- SHAP dependence plots
- SHAP force plot
- Split-based variable importance
- Friedman's H-statistic
- Partial dependence plots (one and two-dimensional)

As the results below show, the information obtained from SHAP values and SHAP dependence plots is the same as when looking at traditional variable importance measures and partial dependence plots. The most important variables are 'K' and 'edu'. From the dependence plots, we conclude that there are non-linearities. For instance, the effect of K is almost flat for large and small values of K and increasing in-between. Further, the effect of edu is steeper for small values of edu and flattens for larger values of edu. Friedman's H-statistic indicates that there are some interactions. For the two variables with the largest amount of interaction, L and pop, we create a two-dimensional partial dependence plot below. Further, the SHAP force plot shows that the effect of the predictor variables is different for the year 2000 compared to the year 2021.

```{r interpret_r, fig.width=7, fig.height=5, cache=cache, message=FALSE, warning=FALSE, results=FALSE}
# SHAP values
library(SHAPforxgboost)
shap.plot.summary.wrap1(gpboost_model, X = data[,covars]) + ggtitle("SHAP values")

# SHAP dependence plots
shap_long <- shap.prep(gpboost_model, X_train = data[,covars])
shap.plot.dependence(data_long = shap_long, x = covars[2], 
                     color_feature = covars[4], smooth = FALSE, size = 2) + 
  ggtitle("SHAP dependence plot for K")
shap.plot.dependence(data_long = shap_long, x = covars[4], 
                     color_feature = covars[2], smooth = FALSE, size = 2) + 
  ggtitle("SHAP dependence plot for edu")

# SHAP force plot
shap_contrib <- shap.values(gpboost_model, X_train = data[,covars])
plot_data <- cbind(shap_contrib$shap_score, ID = 1:dim(data)[1])
shap.plot.force_plot(plot_data, zoom_in=FALSE, id = "ID") + 
  scale_fill_discrete(name="Variable") + xlab("Regions") + 
  geom_vline(xintercept=242.5, linewidth=1) + ggtitle("SHAP force plot") + 
  geom_text(x=100,y=-1.2,label="2000") + geom_text(x=400,y=-1.2,label="2021")

# Split-based feature importances
feature_importances <- gpb.importance(gpboost_model, percentage = TRUE)
gpb.plot.importance(feature_importances, top_n = 25, measure = "Gain", 
                    main = "Split-based variable importances")

# H-statistic for interactions
library(flashlight)
fl <- flashlight(model = gpboost_model, data = data.frame(y = data[,"y"], data[,covars]), 
                 y = "y", label = "gpb",
                 predict_fun = function(m, X) predict(m, data.matrix(X[,covars]), 
                                                      gp_coords_pred = matrix(-100, ncol = 2, nrow = dim(X)[1]),
                                                      group_data_pred = matrix(-1, ncol = 1, nrow = dim(X)[1]),
                                                      pred_latent = TRUE)$fixed_effect)
plot(imp <- light_interaction(fl, v = covars, pairwise = TRUE)) + 
  ggtitle("H interaction statistic") # takes a few seconds

# Partial dependence plots
gpb.plot.partial.dependence(gpboost_model, data[,covars], variable = 2, xlab = covars[2], 
                            ylab = "gdp", main = "Partial dependence plot for K" )
gpb.plot.partial.dependence(gpboost_model, data[,covars], variable = 4, xlab = covars[4], 
                            ylab = "gdp", main = "Partial dependence plot edu" )

# Two-dimensional partial dependence plot (to visualize interactions)
i = 1; j = 3;# i vs j
gpb.plot.part.dep.interact(gpboost_model, data[,covars], variables = c(i,j), xlab = covars[i], 
                           ylab = covars[j], main = "Pairwise partial dependence plot for pop and L")

```

# Extensions

## Separate random effects for different time periods

In the above model, we have used the same random effects for both years 2000 and 2021. Alternatively, one can also use independent spatial and grouped random effects for different time periods (*independent under the prior, conditional on the data, there is dependence ...*). In the ``GPBoost`` library, this can be done via the ``cluster_ids`` argument. ``cluster_ids`` needs to be a vector of length the sample size, and every entry indicates the "cluster" to which an observation belongs to. Different clusters then have independent spatial and grouped random effects, but the hyperparameters (e.g., marginal variance, variance components, etc.) and the fixed effects function are the same across clusters. 

Below we show how we can fit such a model and create two separate spatial maps. As the results below show, the spatial effects are quite different for the two years 2000 and 2021. In particular, there is less (residual) spatial variation (or heterogeneity) for the year 2021 compared to 2000. This is confirmed by looking at standard deviations of the predicted spatial random effects which is almost twice as large for 2000 compared to the year 2021 (see below). A further conclusion is that in 2000, there were more regions with "low" GDP numbers (spatial effects below 0), and this is no longer the case for 2021.

```{r cluster_ids_r, echo=TRUE, fig.height=4, fig.width=10, message=FALSE, warning=FALSE}
gp_model <- GPModel(group_data = data[, c("cl")], gp_coords = data[, c("Long", "Lat")],
                    likelihood = "gaussian", cov_function = "exponential",
                    cluster_ids  = c(rep(1,242), rep(2,242)))
boost_data <- gpb.Dataset(data = data[, covars], label = data[, "y"])
params <- list(learning_rate = 0.01, max_depth = 1, num_leaves = 2^10,
               min_data_in_leaf = 10, lambda_l2 = 1) 
# Note: we use the same tuning parameters as above. Ideally, the would have to be chosen again
gpboost_model <- gpboost(data = boost_data, gp_model = gp_model, nrounds = nrounds,
                         params = params, verbose = 0)
# Separate spatial maps for the years 2000 and 2021
pred <- predict(gpboost_model, group_data_pred = data[, c("cl")], 
                gp_coords_pred = data[, c("Long", "Lat")],
                data = data[, covars], predict_var = TRUE, pred_latent = TRUE,
                cluster_ids_pred = c(rep(1,242), rep(2,242)))
data_plot <- merge(shape, data.frame(FID = FID, y = pred$random_effect_mean[1:242]), by="FID")
plot_mu_2000 <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="Spatial effect", option = "B") +
  ggtitle("Spatial effect for 2000 (mean)") + xlim(2700000,6400000) + ylim(1500000,5200000)
data_plot <- merge(shape, data.frame(FID = FID, y = pred$random_effect_mean[243:484]), by="FID")
plot_mu_2021 <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="Spatial effect", option = "B") +
  ggtitle("Spatial effect for 2021 (mean)") + xlim(2700000,6400000) + ylim(1500000,5200000)
grid.arrange(plot_mu_2000, plot_mu_2021, ncol=2)
sd(pred$random_effect_mean[1:242]) # 0.2321267
sd(pred$random_effect_mean[243:484]) # 0.1286398
```

## Interaction between space and fixed effects predictor variables

In the above model, there is no interaction between the random effects and the fixed effects predictor variables. I.e., there is no interaction between the spatial coordinates and the fixed effects predictor variables. Such interaction can be modeled by additionally including the random effects input variables (= the spatial coordinates or the categorical grouping variable) in the fixed effects function. The code below shows how this can be done. As the variable importance plot below shows, the coordinates are not used in the tree-ensemble, and there is thus no such interaction detectable for this data set.

```{r train_interact_r, results=FALSE}
gp_model <- GPModel(group_data = data[, c("cl")], gp_coords = data[, c("Long", "Lat")],
                    likelihood = "gaussian", cov_function = "exponential")
covars_interact <- c(c("Long", "Lat"), covars) ## add spatial coordinates to fixed effects predictor variables
boost_data <- gpb.Dataset(data = data[, covars_interact], label = data[, "y"])
params <- list(learning_rate = 0.01, max_depth = 1, num_leaves = 2^10,
               min_data_in_leaf = 10, lambda_l2 = 1) 
# Note: we use the same tuning parameters as above. Ideally, the would have to be chosen again
gpboost_model <- gpboost(data = boost_data, gp_model = gp_model, nrounds = nrounds,
                         params = params, verbose = 0)
feature_importances <- gpb.importance(gpboost_model, percentage = TRUE)
gpb.plot.importance(feature_importances, top_n = 25, measure = "Gain", 
                    main = "Var. imp. when including coordinates in the fixed effects")
```

## Large data

For large data sets, calculations with Gaussian processes are slow, and one has to use an approximation. The ``GPBoost`` library implements several ones. For instance, setting ``gp_approx="vecchia"`` in the ``GPModel`` constructor will use a Vecchia approximation. The data set used in this article is relatively small, and we can do all calculations exactly. 

## Other spatial random effects models

Above, we have used a Gaussian process to model spatial random effects. Since the data is areal data, another option is to use models that rely on neighborhood information such as CAR and SAR models. These models are currently not yet implemented in the ``GPBoost`` library (*might be added in the future -> contact the author*).

Another option is to use a grouped random effects model defined by the categorical region ID variable for modeling spatial effects. The code below shows how this can be done in ``GPBoost``. However, this model essentially ignores the spatial structure.

```{r grouped_spat_r, eval=FALSE, include=TRUE}
gp_model <- GPModel(group_data = data[, c("group", "cl")], likelihood = "gaussian")
```

## (Generalized) linear mixed effects and Gaussian process models

(Generalized) linear mixed effects and Gaussian process models can also be run in the ``GPBoost`` library. The code below shows how this can be done with the ``fitGPModel`` function. Note that one needs to manually add a column of 1's to include an intercept. With the option ``params = list(std_dev=TRUE)``, one obtains standard deviations and p-values in the ``summary`` function.

```{r linear_r}
X_incl_1 = cbind(Intercept = rep(1,dim(data)[1]), data[, covars])
gp_model <- fitGPModel(group_data = data[, c("cl")], gp_coords = data[, c("Long", "Lat")],
                       likelihood = "gaussian", cov_function = "exponential",
                       y = data[, "y"], X = X_incl_1, params = list(std_dev=TRUE))
summary(gp_model)
```
# References

- Sigrist Fabio. "[Gaussian Process Boosting](https://www.jmlr.org/papers/v23/20-322.html)". *Journal of Machine Learning Research* (2022).
- Sigrist Fabio. "[Latent Gaussian Model Boosting](https://ieeexplore.ieee.org/document/9759834)". *IEEE Transactions on Pattern Analysis and Machine Intelligence* (2023).
- https://github.com/fabsig/GPBoost