# An Introduction to the distfreereg Package

# 1 Introduction

The distfreereg package has two main functions: `distfreereg()` and `compare()`. The former is introduced in this vignette; `compare()` is introduced here.

The function `distfreereg()` implements the distribution-free regression testing procedure introduced by Khmaladze (2021). Its purpose is to test whether or not a specified function \(Y=f(X;\theta)\) is the conditional mean, \({\rm E}(Y\mathrel|X)\). Specifically, the test is of the following type: \[\begin{equation} H_0{:}\ \exists\theta\in\Theta\,|\,{\rm E}(Y\mathrel|X)=f(X;\theta) \quad\hbox{against}\quad H_1{:}\ \forall\theta\in\Theta\,|\,{\rm E}(Y\mathrel|X)\neq f(X;\theta). \label{eqn:hypo1} \end{equation}\]

This vignette starts by discussing the testing of a linear model created using `lm()`. Following this is an analogous discussion of testing a non-linear model created using `nls()`. Both scenarios can be replicated using the formulas directly. In the most general implementation, the mean function is specified by an R function. To use `distfreereg()` with a mean function implemented in software other than R, see the section on the default method.

# 2 Testing a Linear Model

Suppose that a linear model object is created using `lm()`, and we want to test whether or not a linear model is appropriate; specifically, whether or not the true mean function follows a specified form that is linear in its parameters.

In some cases, a residual plot suffices to show that the specified form is wrong. Below, the true mean function is \({\rm E}(y\mathrel|x)=x^{2.1}\), but the formula used to build the model is linear in \(x\). The residuals plot clearly indicates model mis-specification.

    set.seed(20240304)
    n <- 3e2
    form_lm <- y ~ x
    data_lm <- data.frame(x = runif(n, min = 0, max = 3))
    data_lm$y <- data_lm$x^2.1 + rnorm(n, sd = 0.3)
    m_1 <- lm(form_lm, data = data_lm)
    plot(m_1, which = 1)

If, however, we use a form that is quadratic in \(x\), the residuals plot is suggestive, but it leaves us without a clear conclusion. How certain can we be that the specified form is wrong?

    form_lm_2 <- y ~ I(x^2)
    m_2 <- lm(form_lm_2, data = data_lm)
    plot(m_2, which = 1)

In this case, there appears to be evidence of a mis-specified mean, but we need a formal test to quantify our (un)certainty. The function `distfreereg()` implements such a test.

    set.seed(20240304)
    (dfr_lm_2 <- distfreereg(test_mean = m_2))

    ## Number of observations: 300
    ## Length of EPSP: 300
    ## Monte Carlo simulations: 10000
    ##
    ## Estimated parameter values:
    ## (Intercept)      I(x^2)
    ##  -7.718e-02   1.101e+00
    ##
    ## Observed statistics:
    ## Stat      Value  Pr(>Value)       MCSE
    ##   KS  1.237e+00   7.270e-02  2.596e-03
    ##  CvM  4.506e-01   5.260e-02  2.232e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

In fact, tests based on two statistics are included by default: a Kolmogorov–Smirnov (KS) statistic and a Cramér–von Mises (CvM) statistic, defined as \(\max_i|w_i|\) and \({1\over n}\sum_iw_i^2\), respectively, where \(w\mathrel{\raise 0.4pt\hbox{$:$}\hskip-2pt=}(w_1,\ldots,w_n)\) is the partial sum process calculated by `distfreereg()`. The preceding output shows a table summarizing the tests’ results, including p-values and Monte Carlo standard errors thereof. Both of the p-values are marginal, supplying a measure of our intuitive conclusion based on the residuals plot.

# 3 Testing a Non-Linear Model

The procedure just illustrated applies to objects created with `nls()`, as well. Here, the data generating function has the form being tested, so we should not expect the null to be rejected. Note that the `weights` argument in `nls()` is set equal to the inverse of the variance vector.

    set.seed(20240304)
    n <- 3e2
    sds <- runif(n, min = 0.5, max = 5)
    data_nls <- data.frame(x = rnorm(n), y = rnorm(n))
    data_nls$z <- exp(3*data_nls$x) - 2*data_nls$y^2 + rnorm(n, sd = sds)
    form_nls <- z ~ exp(a*x) - b*y^2
    m_2 <- nls(form_nls, data = data_nls, weights = 1/sds^2, start = c(a = 1, b = 1))

Testing this form of the mean structure is done just as with an `lm` object:

    set.seed(20240304)
    (dfr_nls <- distfreereg(test_mean = m_2))

    ## Number of observations: 300
    ## Length of EPSP: 300
    ## Monte Carlo simulations: 10000
    ##
    ## Estimated parameter values:
    ##         a         b
    ## 3.000e+00 1.957e+00
    ##
    ## Observed statistics:
    ## Stat      Value  Pr(>Value)       MCSE
    ##   KS  6.821e-01   6.835e-01  4.651e-03
    ##  CvM  9.343e-02   6.155e-01  4.865e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

Both tests correctly fail to reject the null.

# 4 Means Specified by Formulas

If no `lm` or `nls` object has been created yet, the formula and data can be supplied directly to `distfreereg()`. Below, we recreate the examples from above without first creating the `lm` and `nls` objects.

The most notable practical difference in implementing the tests in this context is that three arguments are required rather than just one:

1. `test_mean`: a `formula` object specifying the mean function;

2. `data`: a data frame containing the data;

3. `method`: specifies the function to use to create the model.

## 4.1 Formulas with `lm()`

The following call conducts the same test as in `dfr_lm_2` above.

    set.seed(20240304)
    (dfr_form_lm_2 <- distfreereg(test_mean = form_lm_2, data = data_lm, method = "lm"))

    ## Number of observations: 300
    ## Length of EPSP: 300
    ## Monte Carlo simulations: 10000
    ##
    ## Estimated parameter values:
    ## (Intercept)      I(x^2)
    ##  -7.718e-02   1.101e+00
    ##
    ## Observed statistics:
    ## Stat      Value  Pr(>Value)       MCSE
    ##   KS  1.237e+00   7.270e-02  2.596e-03
    ##  CvM  4.506e-01   5.260e-02  2.232e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

The key pieces of the output are identical to those of `dfr_lm`:

    all(identical(dfr_lm_2$observed_stats, dfr_form_lm_2$observed_stats),
        identical(dfr_lm_2$p, dfr_form_lm_2$p))

    ## [1] TRUE

## 4.2 Formulas with `nls()`

An additional element of `method_args` is used here to specify starting values for parameter estimation. (This is optional, since `nls()` will use default starting values with a warning if not supplied. See `nls()` documentation for further details.) The following call conducts the same test as in `dfr_nls`.

    set.seed(20240304)
    (dfr_form_nls <- distfreereg(test_mean = form_nls, data = data_nls, method = "nls",
                                 method_args = list(weights = 1/sds^2,
                                                    start = c(a = 1, b = 1))))

    ## Number of observations: 300
    ## Length of EPSP: 300
    ## Monte Carlo simulations: 10000
    ##
    ## Estimated parameter values:
    ##         a         b
    ## 3.000e+00 1.957e+00
    ##
    ## Observed statistics:
    ## Stat      Value  Pr(>Value)       MCSE
    ##   KS  6.821e-01   6.835e-01  4.651e-03
    ##  CvM  9.343e-02   6.155e-01  4.865e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

The key pieces of the output are identical to those of `dfr_nls`:

    all(identical(dfr_nls$observed_stats, dfr_form_nls$observed_stats),
        identical(dfr_nls$p, dfr_form_nls$p))

    ## [1] TRUE

# 5 Means Specified by Functions

The most general case that can be handled entirely within R is that in which the mean regression function is specified by an R function.

Suppose that we want to replicate the situation in the `nls` examples above, but we want to account for errors that have a non-diagonal covariance matrix. The example below illustrates how to do this.

The function being tested is specified by `test_mean`, which is a function with two arguments, `X` and `theta`. Here, `X` represents the entire matrix of covariates. Therefore, to replicate the references to `x` and `y` in this example, we use `X[,1]` and `X[,2]`, respectively. The arguments `X` and `Y` supply the covariate matrix and the response vector, respectively, in lieu of the `data` argument.

    set.seed(20240304)
    n <- 3e2
    true_mean <- function(X, theta) exp(theta[1]*X[,1]) - theta[2]*X[,2]^2
    test_mean <- true_mean
    theta <- c(3,-2)
    Sigma <- rWishart(1, df = n, Sigma = diag(n))[,,1]
    X <- matrix(rnorm(2*n), nrow = n)
    Y <- distfreereg:::f2ftheta(true_mean, X)(theta) +
      distfreereg:::rmvnorm(n = n, reps = 1, SqrtSigma = distfreereg:::matsqrt(Sigma))

    (dfr_1 <- distfreereg(test_mean = test_mean, Y = Y, X = X,
                          covariance = list(Sigma = Sigma),
                          theta_init = rep(1, length(theta))))

    ## Number of observations: 300
    ## Length of EPSP: 300
    ## Monte Carlo simulations: 10000
    ##
    ## Estimated parameter values:
    ##     theta1     theta2
    ##  3.000e+00 -2.037e+00
    ##
    ## Observed statistics:
    ## Stat      Value  Pr(>Value)       MCSE
    ##   KS  1.016e+00   2.220e-01  4.156e-03
    ##  CvM  1.332e-01   4.413e-01  4.965e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

As with the `nls` example, testing the mean function used to generate the data produces test results that correctly fail to reject the null hypothesis.

# 6 Using the Algorithm with External Functions

All of the examples above illustrate how `distfreereg()` can be used when the mean function being tested is defined in R. If the mean function is not defined in R, `distfreereg()` can still be used as long as certain objects can be imported into R. Specifically, `distfreereg()` needs the model’s fitted values and the Jacobian of the mean function, evaluated at the estimated parameter vector for each observation’s covariate values.

To illustrate this, we extract certain elements from the example above, and pretend that these were created using external software, and then imported into R.

    J <- dfr_1$J
    fitted_values <- fitted(dfr_1)

These can be supplied to `distfreereg()` in the `override` argument to implement the test. Note that `test_mean` is set to `NULL`, and since no parameter estimation is done here, no value for `theta_init` is specified.

    distfreereg(test_mean = NULL, Y = Y, X = X, covariance = list(Sigma = Sigma),
                override = list(fitted_values = fitted_values, J = J))

    ## Number of observations: 300
    ## Length of EPSP: 300
    ## Monte Carlo simulations: 10000
    ##
    ## Observed statistics:
    ## Stat      Value  Pr(>Value)       MCSE
    ##   KS  1.016e+00   2.238e-01  4.168e-03
    ##  CvM  1.332e-01   4.478e-01  4.973e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

# 7 Summary and Diagnostic Plots

Below is an introduction to the three types of plots available for `distfreereg` objects. See the plotting vignette for more details regarding plot customization.

## 7.1 A Summary Plot

A summary of the results of the test can be plotted easily.

    plot(dfr_1)

The density of the simulated statistic is plotted, as is a vertical line showing the observed statistic. The upper tail is shaded, and the p-value (the area of that shaded region) is shown. Finally, a 95% simultaneous confidence band of the density function is plotted, and the region is shaded in grey.

## 7.2 Diagnostic Plots

Two diagnostic plots are also available. The first produces a time-series-like plot of the residuals in the order specified by `res_order` in the `distfreereg` object.

    plot(dfr_1, which = "residuals")

The second produces a plot of the empirical partial sum process.

    plot(dfr_1, which = "epsp")

# 8 Specifying the Covariance Structure

The covariance structure of the errors, conditional on the covariates, must be specified for the testing algorithm. When `test_mean` is not a `function`, the covariance structure is estimated using the supplied/created model object. When `test_mean` is a `function`, this structure must be specified by the user using the `covariance` argument. The value of `covariance` is a list with at least one named element that specifies one of four matrices:

1. `Sigma`, the covariance matrix;

2. `SqrtSigma`, the square root of `Sigma`;

3. `P`, the precision matrix (that is, the inverse of `Sigma`); and

4. `Q`, the square root of `P`.

Internally, the algorithm only needs `Q`, so some efficiency can be gained by supplying this directly if it is known. Supplying more than one of the four matrices is not forbidden, but is strongly discouraged. In this case, no verification is done that the supplied matrices are consistent with each other, and `Q` will be calculated using the most convenient supplied element.

Each of the four named elements can be one of four types of object:

1. A positive-definite matrix, the most general way to specify a covariance structure;

2. A list of positive-definite matrices;

3. a numeric vector of positive values whose length is the sample size; and

4. a length-1 numeric vector specifying a positive number.

Specifying a list of matrices is theoretically equivalent to specifying a block diagonal matrix with the listed matrices as the blocks. Specifying a numeric vector is theoretically equivalent to specifying a diagonal matrix with the given vector along the diagonal. Specifying a single number is theoretically equivalent to specifying a diagonal matrix with that value along the diagonal. Using the simplest possible expression in each case is recommended for conceptual and computational simplicity.

The following code shows that supplying `Q` produces observed statistics identical to those in the previous example.

    Q <- distfreereg:::matinv(distfreereg:::matsqrt(Sigma), tol = .Machine$double.eps)
    dfr_2 <- distfreereg(Y = Y, X = X, test_mean = test_mean,
                         covariance = list(Q = Q),
                         theta_init = rep(1, length(theta)))
    identical(dfr_1$observed_stats, dfr_2$observed_stats)

    ## [1] TRUE

# 9 Methods

Several common generic functions have `distfreereg` methods, including `coef()`, `predict()`, and `update()`. See the documentation for the complete list.

# 10 Adding New Statistics

The `stat` argument of `distfreereg()` allows the user to specify other statistics by specifying any function whose input is a numeric vector and whose output is a numeric vector of length one.

    new_stat <- function(x) sum(abs(x))
    update(dfr_1, stat = "new_stat")

    ## Number of observations: 300
    ## Length of EPSP: 300
    ## Monte Carlo simulations: 10000
    ##
    ## Estimated parameter values:
    ##     theta1     theta2
    ##  3.000e+00 -2.037e+00
    ##
    ## Observed statistics:
    ##     Stat      Value  Pr(>Value)       MCSE
    ## new_stat  8.082e+01   5.373e-01  4.986e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

If necessary, the default statistics can be selected using “`KS`” and “`CvM`” as explicit values in the character vector supplied to `stat`.

    update(dfr_1, stat = c("new_stat", "KS", "CvM"))

    ## Number of observations: 300
    ## Length of EPSP: 300
    ## Monte Carlo simulations: 10000
    ##
    ## Estimated parameter values:
    ##     theta1     theta2
    ##  3.000e+00 -2.037e+00
    ##
    ## Observed statistics:
    ##     Stat      Value  Pr(>Value)       MCSE
    ## new_stat  8.082e+01   5.358e-01  4.987e-03
    ##       KS  1.016e+00   2.202e-01  4.144e-03
    ##      CvM  1.332e-01   4.437e-01  4.968e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

Checks are done to verify that the functions provided calculate legitimate statistics; that is, that they take a numeric vector as input and output a single number.

    bad_stat <- function(x) x[1,2]
    try(update(dfr_1, stat = c("KS", "CvM", "bad_stat")))

    ## Error in value[[3L]](cond) :
    ##   Unable to evaluate stat[i](rnorm(100)): Error in x[1, 2]: incorrect number of dimensions

    try(update(dfr_1, stat = c("KS", "CvM", "undefined_stat")))

    ## Error in value[[3L]](cond) :
    ##   Unable to evaluate stat[i](rnorm(100)): Error in get(stat[i]): object 'undefined_stat' not found

# 11 Non-Normal Errors

One of the strengths of the test implemented by `distfreereg()` is that it makes only weak assumptions about the distribution of the errors. Specifically, it only requires that the error covariance matrix be finite and that the error distribution be strongly mixing (also known as \(\alpha\)-mixing). Below is an example illustrating that we get the expected results when the error distribution is the \(t\) distribution. Note that the \(t\) distribution with three degrees of freedom has variance 3.

The first call to `distfreereg()` below tests the true mean function, and the test correctly fails to reject the null hypothesis. The second call tests a wrong mean function, and the test correctly rejects the null hypothesis.

    set.seed(20241003)
    n <- 1e2
    func <- function(X, theta) theta[1] + theta[2]*X[,1] + 0.5*X[,1]^2
    theta <- c(2,5)
    X <- matrix(rexp(n))
    Y <- theta[1] + theta[2]*X[,1] + 0.5*X[,1]^2 + rt(n, df = 3)

    dfr <- distfreereg(Y = Y, X = X, test_mean = func, theta_init = c(1,1),
                       covariance = list(Sigma = 3))
    dfr

    ## Number of observations: 100
    ## Length of EPSP: 100
    ## Monte Carlo simulations: 10000
    ##
    ## Estimated parameter values:
    ##    theta1    theta2
    ## 1.805e+00 4.934e+00
    ##
    ## Observed statistics:
    ## Stat      Value  Pr(>Value)       MCSE
    ##   KS  6.079e-01   7.623e-01  4.257e-03
    ##  CvM  7.765e-02   7.062e-01  4.555e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

    alt_func <- function(X, theta) theta[1] + theta[2]*X[,1]
    update(dfr, test_mean = alt_func)

    ## Number of observations: 100
    ## Length of EPSP: 100
    ## Monte Carlo simulations: 10000
    ##
    ## Estimated parameter values:
    ##    theta1    theta2
    ## 1.252e+00 6.400e+00
    ##
    ## Observed statistics:
    ## Stat      Value  Pr(>Value)       MCSE
    ##   KS  1.466e+00   1.990e-02  1.397e-03
    ##  CvM  6.932e-01   1.280e-02  1.124e-03
    ## ---
    ## `MCSE' is the Monte Carlo standard error of the estimated p-value.

# References

Khmaladze, E. V. (2021), “Distribution-free testing in linear and parametric regression,” Annals of the Institute of Statistical Mathematics, Springer Science; Business Media LLC, 73, 1063–1087. https://doi.org/10.1007/s10463-021-00786-3.