---
title: "The gofreg package: Perform goodness-of-fit tests for parametric
regression"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{gofreg}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
comment = "#>",
fig.width = 6
)
```
```{r setup, echo = FALSE, message = FALSE}
set.seed(123)
library(gofreg)
library(ggplot2)
```
## Typical usage
1. Generate or load a dataset
2. Fit a parametric regression model to the given data
3. Compute the p-value of the fitted model using one of the available test
statistics
## Example
In this example, we will fit a generalized linear model (GLM) to an artificially
created dataset. It consists of two covariates, one of them normally and the
other one uniformly distributed, and the response variable following a classical
linear model with normal distribution.
```{r}
set.seed(123)
n <- 100
x <- cbind(rnorm(n, mean = 3), runif(n, min = 1, max = 10))
model_true <- GLM.new(distr = "normal", linkinv = identity)
params_true <- list(beta = c(2, 6), sd = 1)
y <- model_true$sample_yx(x, params_true)
data <- dplyr::tibble(x = x, y = y)
```
First, we fit the correct model to the data.
``` {r}
model_test <- GLM.new(distr = "normal", linkinv = identity)
model_test$fit(data, params_init = list(beta = c(1,1), sd = 5), inplace = TRUE)
print(model_test$get_params())
```
The parameters estimates are very close to the true values. To assess whether
the fitted model fits to the given data, we perform a bootstrap-based
goodness-of-fit test using the conditional Kolmogorov test statistic for the
marginal distribution of Y.
``` {r}
gt <- GOFTest$new(data = data, model_fitted = model_test,
test_stat = CondKolmY$new(), nboot = 100)
print(gt$get_pvalue())
```
As we would expect, the p-value is rather high, so the model is not rejected.
Next, we will fit a wrong model to the data. In particular, we exclude the
second covariate.
``` {r}
model_test <- GLM.new(distr = "normal", linkinv = identity)
data_miss <- dplyr::tibble(x = data$x[,1], y = data$y)
model_test$fit(data_miss, params_init = list(beta = c(2), sd = 2),
inplace = TRUE)
print(model_test$get_params())
```
It can be seen that the variance was estimated to be rather high which is
reasonable as it includes the part of the variance that could be explained when
taking the second covariate into account. The corresponding
p-value is computed in the following code chunk.
``` {r}
gt2 <- GOFTest$new(data = data_miss, model_fitted = model_test,
test_stat = CondKolmY$new(), nboot = 100)
print(gt2$get_pvalue())
```
As the p-value is very low, the model hypothesis should be rejected. So the test
reveals the mistake in the model assumption. To further investigate the
discrepancy, we could look at the plots of the processes underlying the test
statistic.
```{r}
gt2$plot_procs()
```
It can be seen that the original process (red line) behaves very differently
than its bootstrap versions (gray lines). For comparison purposes, we also plot
the processes in case the correct model was fitted.
```{r}
gt$plot_procs()
```
This time, the original process (red line) behaves very similar to its bootstrap
versions (gray lines). It does not show any more extreme behavior.
## Parametric Regression Models
Here is a list of parametric regression models that are available in the
`gofreg` package:
* `NormalGLM`: Generalized linear model with normal distribution
* `GammaGLM`: Generalized linear model with gamma distribution
* `ExpGLM`: Generalized linear model with exponential distribution
* `WeibullGLM`: Generalized linear model with Weibull distribution
* `NegBinomGLM`: Generalized linear model with negative binomial distribution
The package also offers the option to use other user-defined models. For
instructions on how to implement new models see `vignette("New-Models")`.
## Test Statistics
Here is a list of test statistics that are available in the `gofreg` package:
* `CondKolmXY`: Conditional Kolmogorov of the joint distribution of $(X,Y)$
defined in Andrews (1997) [doi:10.2307/2171880](https://doi.org/10.2307/2171880)
* `SICM`: Simulated integrated conditional moment test defined in Bierens &
Wang (2012) [doi:10.1017/S0266466611000168](https://doi.org/10.1017/S0266466611000168)
* `MEP`: Marked Empirical Process defined in Dikta & Scheer (2021)
[doi:10.1007/978-3-030-73480-0](https://doi.org/10.1007/978-3-030-73480-0)
* `CondKolmY`: Conditional Kolmogorov of the marginal distribution of $Y$
defined in Kremling & Dikta (2024) [arXiv:2409.20262](https://arxiv.org/abs/2409.20262)
<!-- * `CondKolmbXY`: Conditional Kolmogorov of the joint distribution of $(\beta^T X, Y)$ -->
The package also offers the option to use other user-defined test statistics.
For instructions on how to implement new test statistics see
`vignette("New-TestStatistics")`.
## Censored data
The package can also be used to fit parametric regression models and perform
goodness-of-fit tests for randomly right-censored survival times $Y$. In this
case, the `loglik` and `resample` arguments in the `ParamRegrModel$fit()` and
`GOFTest$new()` methods have to be specified. Moreover, the `data` object needs
to be a `data.frame()` with tags "x", "z" and "delta" with $X$ representing the
covariates, $Z = \min(Y, C)$ the censored times and $\delta = 1_{\{Y \le C\}}$
the censoring indicators. A test statistic for the censored setting is given by
`CondKolmY_RCM`.
Here is an example with artificial data generated from a normal GLM with
normally distributed censoring times.
``` {r, echo = FALSE, message = FALSE}
set.seed(123)
```
``` {r}
n <- 100
x <- cbind(runif(n), rbinom(n, 1, 0.5))
model <- NormalGLM$new()
y <- model$sample_yx(x, params = list(beta = c(2, 3), sd = 1))
c <- rnorm(n, mean(y) * 1.2, sd(y) * 0.5)
data <- dplyr::tibble(x = x, z = pmin(y, c), delta = as.numeric(y <= c))
model$fit(data, params_init = list(beta = c(1, 1), sd = 3), inplace = TRUE,
loglik = loglik_xzd)
print(model$get_params())
```
It can be seen that the estimated parameters are close to the true parameters
$\beta = (2,3)$ and $\sigma = 1$. Now, we compute the corresponding p-value
using the Conditional Kolmogorov test statistic for the marginal distribution
of $Y$ under random censorship.
``` {r}
gt <- GOFTest$new(
data = data, model_fitted = model, test_stat = CondKolmY_RCM$new(),
nboot = 100, resample = resample_param_cens, loglik = loglik_xzd
)
print(gt$get_pvalue())
```
The p-value is rather high and the model is not rejected which is expected since
we fitted the correct model.