---
title: "catalytic_glm_gaussian"
author: Yitong Wu <ywu039@e.ntu.edu.sg>
bibliography: assets/catalytic_refs.bib
link-citations: TRUE
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{catalytic_glm_gaussian}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  fig.width = 7,
  fig.height = 6,
  cache = TRUE
)
```

# Introduction{#Introduction}

This vignette provides an _overview_ of how to use the functions in the 
__catalytic__ package that focuses on GLM Linear Regression. The other __catalytic__ vignettes go into other model-estimating functions.

The goal of the __catalytic__ package is to build framework for catalytic prior distributions. Stabilizing high-dimensional working models by shrinking them towards simplified models. This is achieved by supplementing observed data with weighted synthetic data generated from a predictive distribution under the simpler model. For more information, see [@catalytic_glm].

The two steps of using catalytic package for GLM Linear Regression are

1. Initialization: The `cat_glm_initialization` function constructs a `cat_init` object based on the formula provided by the user to generate synthetic data. The resulting `cat_init` object is tailored to facilitate further analysis, and is integral for subsequent modeling steps in the `catalytic` package.

2. Choose Method(s): Users have the flexibility to choose from four main functions within the catalytic package: `cat_glm`, `cat_glm_tune`, `cat_glm_bayes`, and `cat_glm_bayes_joint`. Each function serves a specific purpose in modeling with catalytic priors and offers distinct capabilities tailored to different modeling scenarios for GLM Linear Regression. This approach enables users to seamlessly incorporate synthetic data with varying weights from different method into GLM Linear Regression analyses, providing flexibility and control over the modeling process.

# Data Preparation{#DataPreparation}

Creating a high-dimensional dataset with a low data size. This step involves increasing the number of features (dimensions) while keeping the number of observations (data size) relatively small. This is useful for testing the performance of `catalytic` models in high-dimensional settings.

A randomly generated dataset is split into training (`train_data`) and test (`test_data`) datasets.

```{r}
library(catalytic)

set.seed(1)

n <- 20 # Number of observations
p <- 5 # Number of predictors

obs_x <- matrix(rnorm(n * (p - 1)), ncol = (p - 1)) # Observation covariates
true_coefs <- rnorm(p) # True coefficient
noise <- rnorm(n) # Noise for more response variability
obs_y <- true_coefs[1] + obs_x %*% true_coefs[-1] + noise # Observation response

obs_data <- as.data.frame(cbind(obs_x, obs_y))
names(obs_data) <- c(paste0("X", 1:(p - 1)), "Y")

# Seperate observation data into train and test data
train_idx <- sample(n, 10)
train_data <- obs_data[train_idx, ]
test_data <- obs_data[-train_idx, ]

print(dim(train_data))
```

In this section, we explore the foundational steps of fitting a Linear regression model (GLM) using the `stats::glm` function with the gaussian family. 

```{r}
# Fit a Linear regression model (GLM)
glm_model <- stats::glm(
  formula = Y ~ .,
  family = gaussian,
  data = train_data
)

predicted_y <- predict(
  glm_model,
  newdata = test_data
)

cat(
  "MLE GLM gaussian Model - Mean Square Error (Data):",
  mean((predicted_y - test_data$Y)^2)
)

cat(
  "\nMLE GLM gaussian Model - Sum Square Error (Coefficients):",
  sum((coef(glm_model) - true_coefs)^2)
)
```

Let us check the scatter plot of the `predicted_y` from `glm_model` versus the `test_data$Y`, this can be a great way to visually assess the accuracy and performance of the model.

```{r}
plot(test_data$Y,
  predicted_y,
  main = "Scatter Plot of true Y vs Predicted Y",
  xlab = "true Y",
  ylab = "Predicted Y",
  pch = 19,
  col = "blue"
)

# Add a 45-degree line for reference
abline(a = 0, b = 1, col = "red", lwd = 2)
```


# Usage of `catalytic`

## Step 1: Initialization{#Initialization}

To initialize data for GLM Linear Regression using the `catalytic` package, the `cat_glm_initialization` function is employed. This function facilitates the setup by preparing synthetic data tailored for modeling purposes.

Here's a breakdown of the parameters used:

* `formula`: A formula specifying the GLMs. Should include response and predictor variables.
* `family`: The type of GLM family. Defaults to Gaussian.
* `data`: A data frame containing the data for modeling.
* `syn_size`: An integer specifying the size of the synthetic dataset to be generated. Default is four times the number of predictor columns.
* `custom_variance`: A custom variance value to be applied if using a Gaussian model. Defaults to `NULL`.
* `gaussian_known_variance`: A logical value indicating whether the data variance is known. Defaults to `FALSE`. Only applicable to Gaussian family.
* `x_degree`: A numeric vector indicating the degree for polynomial expansion of predictors. Default is 1 for each predictor.
* `resample_only`: A logical indicating whether to perform resampling only. Default is FALSE.
* `na_replace`: A function to handle NA values in the data. Default is `stats::na.omit`.

```{r}
cat_init <- cat_glm_initialization(
  formula = Y ~ 1,
  family = gaussian,
  data = train_data,
  syn_size = 50,
  custom_variance = NULL,
  gaussian_known_variance = FALSE,
  x_degree = NULL,
  resample_only = FALSE,
  na_replace = stats::na.omit
)

cat_init
```

Here shows how users can simplify the input for `cat_glm_initialization`. User do not have to specify `syn_size` and other parameters, as they have default values, which mentioned above. `cat_init` objects contain a list of attributes, which is typically generated from above function [`cat_glm_initialization`]{#Initialization}. These attributes provide comprehensive information for below modeling tasks or for user check. 

Here's a breakdown of all attributes except the input parameters:

* `function_name`: The name of this function, which is `cat_glm_initialization`.

* `obs_size`: The number of observations (rows) in the original dataset (`obs_data`).

* `obs_data`: The original dataset (`data`).

* `obs_x`: The covariates from original dataset (`obs_data`).

* `obs_y`: The response from original dataset (`obs_data`).

* `syn_size`: The size of synthetic data generated.

* `syn_data`: The synthetic data created for modeling purposes, based on the original dataset (`obs_data`) characteristics.

* `syn_x`: The covariates from synthetic dataset (`syn_data`).

* `syn_y`: The response from synthetic dataset (`syn_data`).

* `syn_x_resample_inform`: The information detailing the process of resampling synthetic data.

* `size`: The total size of the combined dataset (`obs_size` and `syn_size`).

* `data`: The combined dataset (`obs_data` and `syn_data`).

* `x`: The combined covariates (`obs_x` and `syn_x`)

* `y`: The combined response (`obs_y` and `syn_y`)

For more details, please check `?cat_glm_initialization`.

```{r}
names(cat_init)
```

And of course, user can extract items mentioned above from `cat_glm_initialization` object.

```{r}
# The number of observations (rows) in the original dataset (`obs_data`)
cat_init$obs_size

# The information detailing the process of resampling synthetic data
cat_init$syn_x_resample_inform
```


## Step 2.1: Choose Method(s) - Estimation with Fixed tau{#cat_glm}

The `cat_glm` function fits a Generalized Linear Model (GLM) with a catalytic prior on the regression coefficients. It utilizes information from the `cat_init` object generated during the initialization step, which includes both observed and synthetic data, plus other relevant information. 

The GLM model is then fitted using the specified formula, family, and a single `tau`(synthetic data down-weight factor). The resulting `cat_glm` object encapsulates the fitted model, including estimated coefficients and family information, facilitating further analysis.

Here's a breakdown of the parameters used:

* `formula`: This parameter specifies the model formula used in the GLM (Generalized Linear Model). It defines the relationship between the response variable and the predictors. Alternatively, besides using in format `RESPONSE ~ COVARIATES`, user can also use `~ COVARIATES` without specifying the response name, since the response name is defined in the initialization step.

* `cat_init`: This parameter is essential and represents the initialization object (`cat_init`) created by using [`cat_glm_initialization`](#Initialization). It contains both observed and synthetic data, plus other relevant information, necessary for model fitting.

* `tau`: This parameter determines the down-weight assigned to synthetic data relative to observed data. It influences the influence of synthetic data in the model fitting process. If not specified (`NULL`), it defaults to a one forth of the number of predictors.

```{r}
cat_glm_model <- cat_glm(
  formula = Y ~ ., # Same as `~ .`
  cat_init = cat_init, # Output object from `cat_glm_initialization`
  tau = 10 # Defaults to the number of predictors / 4
)

cat_glm_model
```

Here shows how users can simplify the input for `cat_glm`. User do not have to specify `tau`, as `tau` has default value , which mentioned above.

Let's check the prediction error. 

```{r}
cat_glm_predicted_y <- predict(
  cat_glm_model,
  newdata = test_data
)

cat(
  "Catalytic `cat_glm` - Mean Square Error (Data):",
  mean((cat_glm_predicted_y - test_data$Y)^2)
)

cat(
  "\nCatalytic `cat_glm` - Sum Square Error (Coefficients):",
  sum((coef(cat_glm_model) - true_coefs)^2)
)
```

Let us check the scatter plot of the `cat_glm_predicted_y` from `cat_glm_model` versus the `test_data$Y`, this can be a great way to visually assess the accuracy and performance of the model.

```{r}
plot(test_data$Y,
  cat_glm_predicted_y,
  main = "Scatter Plot of true Y vs Predicted Y (cat_glm)",
  xlab = "true Y",
  ylab = "Predicted Y (cat_glm)",
  pch = 19,
  col = "blue"
)

# Add a 45-degree line for reference
abline(a = 0, b = 1, col = "red", lwd = 2)
```

Both `cat_glm_model` objects are outputs from the `cat_glm` function, providing a list of attributes for further analysis or user inspection.

Here's a breakdown of all attributes except the input parameters:

* `function_name`: The name of the function (`cat_glm`).

* `model`: The fitted GLM model object obtained from `stats::glm`, with `tau`.

* `coefficients`: The estimated coefficients from the fitted GLM model `model`.

For more details, please check `?cat_glm`.

```{r}
names(cat_glm_model)
```

User can extract items mentioned above from `cat_glm` object.

```{r}
# The formula used for modeling
cat_glm_model$formula

# The fitted GLM model object obtained from `stats::glm` with `tau`
cat_glm_model$coefficients
```

## Step 2.2: Choose Method(s) - Estimation with Selective tau{#cat_glm_tune}

The `cat_glm_tune` function fits a GLM with a catalytic prior on the regression coefficients and provides options for optimizing model performance over a range of tau values(`tau_seq`).

These methods empower users to fit and optimize GLM models with catalytic priors, leveraging both observed and synthetic data to enhance model performance and robustness in various statistical analyses.

### Cross-validation (risk_estimate_method = "cross_validation"){#CrossValidation}

This method computes the partial likelihood across a specified range of tau values (`tau_seq`). It iterates through each tau value, evaluating its performance based on cross-validation folds (`cv_fold_num`) to select the optimal tau that minimizes the discrepancy error.

Here's a breakdown of the parameters used:

* `formula` and `cat_init` are same as [above](#cat_glm).

* `risk_estimate_method`: Method for risk estimation, chosen from "parametric_bootstrap", "cross_validation", "mallowian_estimate", "steinian_estimate". In this example, "cross_validation" is used.

* `discrepancy_method`: Method for discrepancy calculation, chosen from "mean_square_error", "mean_classification_error", "logistic_deviance". In this example, "mean_square_error" is used because the family is `gaussian`.

* `tau_seq`: Vector of positive numeric values for down-weighting synthetic data. Defaults to a sequence around one fourth of the number of predictors.

* `cross_validation_fold_num`: Number of folds for cross-validation. Defaults to 5.

```{r}
cat_glm_tune_cv <- cat_glm_tune(
  formula = Y ~ ., # Same as `~ .`
  cat_init = cat_init,
  risk_estimate_method = "cross_validation", # Default auto-select based on the data size
  discrepancy_method = "mean_square_error", # Default auto-select based on family
  tau_seq = seq(0, 5, 0.5), # Default is a numeric sequence around the number of predictors / 4
  cross_validation_fold_num = 2 # Default: 5
)

cat_glm_tune_cv
```

User can plot the tau_seq (x) against discrepancy error (y) using the `plot()` function. This plot will show the lowest discrepancy error at the optimal tau value.

```{r}
plot(cat_glm_tune_cv, legend_pos = "topright", text_pos = 3)
```

### Bootstrap (risk_estimate_method = "parametric_bootstrap"){#Bootstrap}

This method estimates tau using bootstrap resampling, refining the model through iterative sampling to enhance robustness and accuracy.

Here's a breakdown of other parameters used:

* `tau_0`: Initial tau value used for discrepancy calculation in risk estimation. Defaults to one fourth of the number of predictors for binomial and 1 for gaussian.

* `parametric_bootstrap_iteration_times`: Number of bootstrap iterations for "parametric_bootstrap" risk estimation. Defaults to 100.

For the breakdown of other input parameters, please check section [Cross Validation](#CrossValidation)

```{r}
cat_glm_tune_boots <- cat_glm_tune(
  formula = ~., # Same as `Y ~ .`
  cat_init = cat_init,
  risk_estimate_method = "parametric_bootstrap", # Default auto-select based on the data size
  discrepancy_method = "mean_square_error", # Default auto-select based on family
  tau_0 = 2, # Default: 1
  parametric_bootstrap_iteration_times = 5, # Default: 100
)

cat_glm_tune_boots
```

### Mallowian Estimate (risk_estimate_method = " mallowian_estimate"){#MallowianEstimate}

This method computes the risk estimate using a mallowian estimate approach, optimizing the model based on observed and synthetic data.

For the breakdown of the input parameters, please check section [Cross Validation](#CrossValidation) and [Bootstrap](#Bootstrap))

```{r}
cat_glm_tune_mallowian <- cat_glm_tune(
  formula = ~., # Same as `Y ~ .`
  cat_init = cat_init,
  risk_estimate_method = "mallowian_estimate", # Default auto-select based on the data size
  discrepancy_method = "mean_square_error", # Default auto-select based on family
)

cat_glm_tune_mallowian
```

### Recommendations for Choosing `risk_estimate_method` and `discrepancy_method`{#Recommendations}

Choosing the appropriate risk_estimate_method and discrepancy_method depends on the data size, model complexity, and the specific requirements of user's analysis.

1. Choosing `risk_estimate_method`

  * __Small to Medium Data Size (observation data size <= 200)__
    * Default: "parametric_bootstrap"
    * For smaller datasets, parametric bootstrap is generally chosen. It efficiently estimates the model's performance using resampling techniques.

  * __Large Data Size (observation data size > 200)__:
    * Default: "cross_validation"
    * For larger datasets, cross-validation is preferred. It provides a robust estimate of model performance by splitting the data into multiple folds and averaging the results.

2. Setting `discrepancy_method`

  * __GLM Family is gaussian (family == "gaussian")__
    * For gaussian GLM families, "mean_square_error" is the only option. It calculates the discrepancy based on the linear likelihood function, suitable for linear outcomes.

  * __GLM Family is not gaussian (family != "gaussian")__:
    * Default: "logistic_deviance"
    * For other GLM families (e.g. binomial), "logistic_deviance" is typically used. Besides "logistic_deviance", user can also choose "classification_error" for binary outcomes. Please check `catalytic_glm_binomial` for more details.

### Automatic Parameter Selection{#AutoParamSelect}

Of course, user don't need to worry about specifying these parameters explicitly, and they just need to simply provide the `cat_init` object and the `formula`. then `cat_glm_tune` will automatically select `risk_estimate_method` and `discrepancy_method` based on the dataset size and GLM family type. 

In this example, it is `risk_estimate_method = "parametric_bootstrap"` and `discrepancy_method = "square_error"`.

For the breakdown of the input parameters, please check section [Cross Validation](#CrossValidation) and [Bootstrap](#Bootstrap)

```{r}
cat_glm_tune_auto <- cat_glm_tune(
  formula = ~., # Same as `Y ~ .`
  cat_init = cat_init
)

cat_glm_tune_auto
```

Let's check the prediction error. 

```{r, eval=FALSE}
cat_glm_tune_predicted_y <- predict(
  cat_glm_tune_auto,
  newdata = test_data
)

cat(
  "Catalytic `cat_glm_tune` - Mean Square Error (Data):",
  mean((cat_glm_tune_predicted_y - test_data$Y)^2)
)

cat(
  "\nCatalytic `cat_glm_tune` - Sum Square Error (Coefficients):",
  sum((coef(cat_glm_tune_auto) - true_coefs)^2)
)
```

Let us check the scatter plot of the `cat_glm_tune_predicted_y` from `cat_glm_tune_auto` versus the `test_data$Y`, this can be a great way to visually assess the accuracy and performance of the model.

```{r, eval=FALSE}
plot(test_data$Y, cat_glm_tune_predicted_y,
  main = "Scatter Plot of true Y vs Predicted Y (cat_glm_tune)",
  xlab = "true Y",
  ylab = "Predicted Y (cat_glm_tune)",
  pch = 19,
  col = "blue"
)

# Add a 45-degree line for reference
abline(a = 0, b = 1, col = "red", lwd = 2)
```

All above objects in this section including `cat_glm_tune_auto` objects are outputs from the `cat_glm_tune` function, providing a list of attributes for further analysis or user inspection.

Here's a breakdown of all attributes except the input parameters:

* `function_name`: The name of the function (`cat_glm_tune`).

* `tau`: Selected optimal tau value from `tau_seq` that minimizes discrepancy error.

* `model`: The fitted GLM model object obtained from `stats::glm`, with the selected tau (`tau`).

* `coefficients`: The estimated coefficients from the fitted GLM model (`model`).

* `risk_estimate_list`: Collected risk estimates across different tau values.

For more details, please check `?cat_glm_tune`.

```{r, eval=FALSE}
names(cat_glm_tune_auto)
```

User can extract items mentioned above from `cat_glm_tune` object.

```{r, eval=FALSE}
# The method used for risk estimation
cat_glm_tune_auto$risk_estimate_method

# Selected optimal tau value from `tau_seq` that minimizes discrepancy error
cat_glm_tune_auto$tau
```

## Step 2.3: Choose Method(s) - Bayesian Posterior Sampling with Fixed tau{#cat_glm_bayes}

Now, we will explore advanced Bayesian modeling techniques tailored for GLM gaussian using the `catalytic` package. Bayesian inference offers a powerful framework to estimate model parameters and quantify uncertainty by integrating prior knowledge with observed data. 

Below functions enable Bayesian inference for GLM Linear Regression Model with catalytic priors. This function utilizes Markov chain Monte Carlo (MCMC) methods, implemented using the `rstan` package, to sample from the posterior distribution of model parameters. Users can specify various options such as the number of MCMC chains (`chains`), iterations (`iter`), warmup steps (`warmup`), and the MCMC algorithm (`algorithm`). User could also apply other attributes to `rstan::sampling`, like `refresh` and `control`.

In this section, we explore Bayesian approaches using the `cat_glm_bayes` function from the `catalytic` package. This function can fit a Bayesian Generalized Linear Model (GLM) using a fixed tau value. The MCMC sampling process will generate posterior distributions for the coefficients based on the specified tau.

Here's a breakdown of the parameters used:

* `formula`, `cat_init` and `tau` are same as [above](#cat_glm).

* `chains`: Number of Markov chains to run during MCMC sampling in `rstan`. Defaults to 4.

* `iter`: Total number of iterations per chain in the MCMC sampling process in `rstan`. Defaults to 2000.

* `warmup`: Number of warmup iterations in the MCMC sampling process in `rstan`, discarded as burn-in. Defaults to 1000.

* `algorithm`: Specifies the sampling algorithm used in `rstan`. Defaults to "NUTS" (No-U-Turn Sampler).

* `gaussian_variance_alpha`: The shape parameter for the inverse-gamma prior on variance if the variance is unknown in Gaussian models. Defaults to the number of predictors. 

* `gaussian_variance_beta`: The scale parameter for the inverse-gamma prior on variance if the variance is unknown in Gaussian models. Defaults to the number of predictors times variance of observation response.

* `...` (ellipsis): Denotes additional arguments that can be passed directly to the underlying `rstan::sampling` function used within `cat_glm_bayes` to fit the Bayesian GLM model. These arguments allow for customization of the Bayesian GLM fitting process, such as `control`, `refresh`, or other model-specific settings.

For more details, please refer to `?cat_glm_bayes`.

```{r, eval = FALSE}
cat_glm_bayes_model <- cat_glm_bayes(
  formula = Y ~ ., # Same as `~ .`
  cat_init = cat_init,
  tau = 50, # Default: number of predictors / 4
  chains = 1, # Default: 4
  iter = 100, # Default: 2000
  warmup = 50, # Default: 1000
  algorithm = "NUTS", # Default: NUTS
  gaussian_variance_alpha = 1, # Default: number of predictors
  gaussian_variance_beta = 1 # Default: number of predictors times variance of observation response
)

cat_glm_bayes_model
```

Here shows how users can simplify the input for `cat_glm_bayes`. User do not have to specify `tau` and other attributes, as `tau` and other attributes have default value, which mentioned above. Here we assign lower value to `chains`, `iter` and `warmup` for quicker processing time.

User can also get the traceplot of the `rstan` model by using `traceplot()` directly into the output from `cat_glm_bayes`.

```{r, eval = FALSE}
traceplot(cat_glm_bayes_model)
```

Plus, user can use this `catlaytic::traceplot` just like the `rstan::traceplot`, user can add parameters used in `rstan::traceplot`, like `include` and `inc_warmup`.

```{r, eval = FALSE}
traceplot(cat_glm_bayes_model, inc_warmup = TRUE)
```

Let's check the prediction error. 

```{r, eval = FALSE}
cat_glm_bayes_predicted_y <- predict(
  cat_glm_bayes_model,
  newdata = test_data
)

cat(
  "Catalytic cat_glm_bayes - Mean Square Error (Data):",
  mean((cat_glm_bayes_predicted_y - test_data$Y)^2)
)

cat(
  "\nCatalytic cat_glm_bayes - Sum Square Error (Coefficients):",
  sum((coef(cat_glm_bayes_model) - true_coefs)^2)
)
```

Let us check the scatter plot of the `cat_glm_bayes_predicted_y` from `cat_glm_bayes` versus the `test_data$Y`, this can be a great way to visually assess the accuracy and performance of the model.

```{r, eval = FALSE}
plot(test_data$Y, cat_glm_bayes_predicted_y,
  main = "Scatter Plot of true Y vs Predicted Y (cat_glm_bayes)",
  xlab = "true Y",
  ylab = "Predicted Y (cat_glm_bayes)",
  pch = 19,
  col = "blue"
)

# Add a 45-degree line for reference
abline(a = 0, b = 1, col = "red", lwd = 2)
```

Both `cat_glm_bayes_model` and `cat_glm_bayes_model` objects are outputs from the `cat_glm_bayes` function, providing a list of attributes for further analysis or user inspection.

Here's a breakdown of all attributes except the input parameters:

* `function_name`: The name of the function (`cat_glm_bayes`).

* `stan_data`: The data list used for `rstan::sampling`.

* `stan_model`: The `rstan::stan_model` object used for Bayesian modeling.

* `stan_sample_model`: The result object obtained from `rstan::sampling`, encapsulating the MCMC sampling results.

* `coefficients`: The mean estimated coefficients from the Bayesian GLM model, extracted from `rstan::summary(stan_sample_model)$summary`.

For more details, please refer to `?cat_glm_bayes`.

```{r, eval = FALSE}
names(cat_glm_bayes_model)
```

User can extract items mentioned above from `cat_glm_bayes` object.

```{r, eval = FALSE}
# The number of Markov chains used during MCMC sampling in `rstan`.
cat_glm_bayes_model$chain

# The mean estimated coefficients from the Bayesian GLM model
cat_glm_bayes_model$coefficients
```

## Step 2.4: Choose Method(s) - Bayesian Posterior Sampling with Adaptive Tau{#cat_glm_bayes_joint}

In this section, we delve into Bayesian methodologies employing the `cat_glm_bayes_joint` function within the `catalytic` package. Unlike its non-adaptive counterpart (`cat_glm_bayes`), this method employs a joint tau prior approach where tau is treated as a parameter within the MCMC sampling process, improving the robustness and accuracy of parameter estimation in Bayesian gaussian modeling.

In this section, we explore Bayesian approaches using the `cat_glm_bayes_joint` function from the `catalytic` package. These functions are similar to their non-adaptive (non-joint) version [`cat_glm_bayes`](#cat_glm_bayes), but corporate `tau` into the MCMC sampling process.

Here's a breakdown of the parameters used:

* `formula`, and `cat_init` are same as [above](#cat_glm).

* `chains`, `iter`, `warmup`, `algorithm`, `gaussian_variance_alpha` and
 `gaussian_variance_beta` are same in section [Bayesian Posterior Sampling with Fixed Tau](#cat_glm_bayes).

* `tau_alpha`: Alpha parameter controlling degrees of freedom for distribution in the joint tau approach. Default is 2.

* `tau_gamma`: Gamma parameter in the joint tau approach. Default is 1.

* `binomial_joint_theta`: Logical. If TRUE, use theta parameter in the binomial model. Default is FALSE. More explanation in `catalytic_glm_binomial.Rmd`.

* `binomial_joint_alpha`: Logical. If TRUE, use joint alpha in the binomial model. Default is FALSE. More explanation in `catalytic_glm_binomial.Rmd`.

* `binomial_tau_lower`: Lower limit for the tau parameter in the binomial model. Default is 0.05. More explanation in `catalytic_glm_binomial.Rmd`.

* `...`(ellipsis): Denotes additional arguments that can be passed directly to the underlying `rstan::sampling` function used within `cat_glm_bayes` to fit the Bayesian GLM model. These arguments allow for customization of the Bayesian GLM fitting process, such as `control`, `refresh`, or other model-specific settings.

For more details, please refer to `?cat_glm_bayes_joint`.

```{r, eval = FALSE}
cat_glm_bayes_joint_model <- cat_glm_bayes_joint(
  formula = Y ~ ., # Same as `~ .`
  cat_init = cat_init,
  chains = 1, # Default: 4
  iter = 100, # Default: 2000
  warmup = 50, # Default: 1000
  algorithm = "NUTS", # Default: NUTS
  tau_alpha = 2, # Default: 2
  tau_gamma = 1 # Default: 1
)

cat_glm_bayes_joint_model
```

Here shows how users can simplify the input for `cat_glm_bayes_joint`. User do not have to specify `tau_alpha` and other attributes, as `tau_alpha` can derived from `cat_init` , while other attributes have default value, which mentioned above. Here we assign lower value to `chains`, `iter` and `warmup` for quicker processing time.

User can also get the traceplot of the `rstan` model by using `traceplot()` directly into the output from `cat_glm_bayes_joint`.

```{r, eval = FALSE}
traceplot(cat_glm_bayes_joint_model)
```

Like the [`traceplot` shown in the `cat_glm_bayes` function](#cat_glm_bayes) , user can add parameters used in `rstan::traceplot`, like `include` and `inc_warmup`.

```{r, eval = FALSE}
traceplot(cat_glm_bayes_joint_model, inc_warmup = TRUE)
```

Let's check the prediction error. 

```{r, eval = FALSE}
cat_glm_bayes_joint_predicted_y <- predict(
  cat_glm_bayes_joint_model,
  newdata = test_data
)

cat(
  "Catalytic `cat_glm_bayes_joint` - Mean Square Error (Data):",
  mean((cat_glm_bayes_joint_predicted_y - test_data$Y)^2)
)

cat(
  "\nCatalytic `cat_glm_bayes_joint` - Sum Square Error (Coefficients):",
  sum((coef(cat_glm_bayes_joint_model) - true_coefs)^2)
)
```

Let us check the scatter plot of the `predicted_y` from `glm_model` versus the `test_data$Y`, this can be a great way to visually assess the accuracy and performance of the model.

```{r, eval = FALSE}
plot(test_data$Y,
  cat_glm_bayes_joint_predicted_y,
  main = "Scatter Plot of true Y vs Predicted Y (cat_glm_bayes_joint)",
  xlab = "true Y",
  ylab = "Predicted Y (cat_glm_bayes_joint)",
  pch = 19,
  col = "blue"
)

# Add a 45-degree line for reference
abline(a = 0, b = 1, col = "red", lwd = 2)
```

Both `cat_glm_bayes_joint_model` and `cat_glm_bayes_joint_model` objects are outputs from the `cat_glm_bayes_joint` function, providing a list of attributes for further analysis or user inspection.

Here's a breakdown of all attributes except the input parameters:

* `function_name`: The name of the function (`cat_glm_bayes_joint`).

* `tau`: The estimated tau parameter from the MCMC sampling `rstan::sampling`, depending on the model configuration.

* `stan_data`, `stan_model`, `stan_sample_model` and `coefficients` are same in section [Bayesian Posterior Sampling with Fixed Tau](#cat_glm_bayes).

For more details, please refer to `?cat_glm_bayes_joint`.

```{r, eval = FALSE}
names(cat_glm_bayes_joint_model)
```

User can extract items mentioned above from `cat_glm_bayes_joint` object.

```{r, eval = FALSE}
# The estimated tau parameter from the MCMC sampling `rstan::sampling`,
cat_glm_bayes_joint_model$tau

# The mean estimated coefficients from the Bayesian GLM model
cat_glm_bayes_joint_model$coefficients
```

# References