simonoff_tsai {skedastic}

Simonoff-Tsai Tests for Heteroskedasticity in a Linear Regression Model

Description

This function implements the modified profile likelihood ratio test and score test of Simonoff and Tsai (1994) for testing for heteroskedasticity in a linear regression model.

Usage

simonoff_tsai(
  mainlm,
  auxdesign = NA,
  method = c("mlr", "score"),
  hetfun = c("mult", "add", "logmult"),
  basetest = c("koenker", "cook_weisberg"),
  bartlett = TRUE,
  optmethod = "Nelder-Mead",
  statonly = FALSE,
  ...
)

Arguments

mainlm

Either an object of class "lm" (e.g., generated by lm), or a list of two objects: a response vector and a design matrix. The objects are assumed to be in that order, unless they are given the names "X" and "y" to distinguish them. The design matrix passed in a list must begin with a column of ones if an intercept is to be included in the linear model. The design matrix passed in a list should not contain factors, as all columns are treated 'as is'. For tests that use ordinary least squares residuals, one can also pass a vector of residuals in the list, which should either be the third object or be named "e".

auxdesign

A data.frame or matrix representing an auxiliary design matrix of containing exogenous variables that (under alternative hypothesis) are related to error variance, or a character "fitted.values" indicating that the fitted \hat{y}_i values from OLS should be used. If set to NA (the default), the design matrix of the original regression model is used. An intercept is included in the auxiliary regression even if the first column of auxdesign is not a vector of ones.

method

A character specifying which of the tests proposed in Simonoff and Tsai (1994) to implement. "mlr" corresponds to the modified profile likelihood ratio test, and "score" corresponds to the score test.

hetfun

A character describing the form of w(\cdot), the error variance function under the heteroskedastic alternative. Possible values are "mult" (the default), corresponding to w(Z_i,\lambda)=\exp\left\{\sum_{j=1}^{q}\lambda_j Z_{ij}\right\}, "add", corresponding to w(Z_i,\lambda)=\left(1+\sum_{j=1}^{q} \lambda_j Z_{ij}\right)^2, and "logmult", corresponding to w(Z_i,\lambda)=\exp\left\{\sum_{j=1}^{q}\lambda_j \log Z_{ij}\right\}. The multiplicative and log-multiplicative cases are considered in Cook and Weisberg (1983); the additive case is discussed, inter alia, by Griffiths and Surekha (1986). Results for the additive and multiplicative models are identical for this test. Partial matching is used.

basetest

A character specifying the base test statistic which is robustified using the added term described in Details. "koenker" corresponds to the test statistic produced by breusch_pagan with argument koenker set to TRUE, while "cook_weisberg" corresponds to the test statistic produced by cook_weisberg. Partial matching is used. This argument is only used if method is "score".

bartlett

A logical specifying whether a Bartlett correction should be made, as per Ferrari et al. (2004), to improve the fit of the test statistic to the asymptotic null distribution. This argument is only applicable where method is "mlr", and is implemented only where hetfun is "mult" or "logmult".

optmethod

A character specifying the optimisation method to use with optim, if method is "mlr". The default, "Nelder-Mead", corresponds to the default method value in optim. Warnings about Nelder-Mead algorithm being unreliable for one-dimensional optimization have been suppressed, since the algorithm does appear to work for the three implemented choices of hetfun.

statonly

A logical. If TRUE, only the test statistic value is returned, instead of an object of class "htest". Defaults to FALSE.

...

Optional arguments to pass to optim, such as par (initial value of \lambda) and maxit (maximum number of iterations to use in optimisation algorithm), and trace (to provide detailed output on optimisation algorithm). Default initial value of \lambda is rep(1e-3, q).

Details

The Simonoff-Tsai Likelihood Ratio Test involves a modification of the profile likelihood function so that the nuisance parameter will be orthogonal to the parameter of interest. The maximum likelihood estimate of \lambda (called \delta in Simonoff and Tsai (1994)) is computed from the modified profile log-likelihood function using the Nelder-Mead algorithm in optim. Under the null hypothesis of homoskedasticity, the distribution of the test statistic is asymptotically chi-squared with q degrees of freedom. The test is right-tailed.

The Simonoff-Tsai Score Test entails adding a term to either the score statistic of Cook and Weisberg (1983) (a test implemented in cook_weisberg) or to that of Koenker (1981) (a test implemented in breusch_pagan with argument koenker set to TRUE), in order to improve the robustness of these respective tests in the presence of non-normality. This test likewise has a test statistic that is asymptotically \chi^2(q)-distributed and the test is likewise right-tailed.

The assumption of underlying both tests is that \mathrm{Cov}(\epsilon)=\omega W, where W is an n\times n diagonal matrix with ith diagonal element w_i=w(Z_i, \lambda). Here, Z_i is the ith row of an n \times q nonstochastic auxiliary design matrix Z. Note: Z as defined here does not have a column of ones, but is concatenated to a column of ones when used in an auxiliary regression. \lambda is a q-vector of unknown parameters, and w(\cdot) is a real-valued, twice-differentiable function having the property that there exists some \lambda_0 for which w(Z_i,\lambda_0)=0 for all i=1,2,\ldots,n. Thus, the null hypothesis of homoskedasticity may be expressed as \lambda=\lambda_0.

In the score test, the added term in the test statistic is of the form

\sum_{j=1}^{q} \left(\sum_{i=1}^{n} h_{ii} t_{ij}\right) \tau_j

, where t_{ij} is the (i,j)th element of the Jacobian matrix J evaluated at \lambda=\lambda_0:

t_{ij}=\left.\frac{\partial w(Z_i, \lambda)}{\partial \lambda_j}\right|_{\lambda=\lambda_0}

, and \tau=(\bar{J}'\bar{J})^{-1}\bar{J}'d, where d is the n-vector whose ith element is e_i^2\bar{\omega}^{-1}, \bar{\omega}=n^{-1}e'e, and \bar{J}=(I_n-1_n 1_n'/n)J.

Value

An object of class "htest". If object is not assigned, its attributes are displayed in the console as a tibble using tidy.

References

Cook RD, Weisberg S (1983). “Diagnostics for Heteroscedasticity in Regression.” Biometrika, 70(1), 1–10.

Ferrari SL, Cysneiros AH, Cribari-Neto F (2004). “An Improved Test for Heteroskedasticity Using Adjusted Modified Profile Likelihood Inference.” Journal of Statistical Planning and Inference, 124, 423–437.

Griffiths WE, Surekha K (1986). “A Monte Carlo Evaluation of the Power of Some Tests for Heteroscedasticity.” Journal of Econometrics, 31(1), 219–231.

Koenker R (1981). “A Note on Studentizing a Test for Heteroscedasticity.” Journal of Econometrics, 17, 107–112.

Simonoff JS, Tsai C (1994). “Use of Modified Profile Likelihood for Improved Tests of Constancy of Variance in Regression.” Journal of the Royal Statistical Society. Series C (Applied Statistics), 43(2), 357–370.

Examples

mtcars_lm <- lm(mpg ~ wt + qsec + am, data = mtcars)
simonoff_tsai(mtcars_lm, method = "score")
simonoff_tsai(mtcars_lm, method = "score", basetest = "cook_weisberg")
simonoff_tsai(mtcars_lm, method = "mlr")
simonoff_tsai(mtcars_lm, method = "mlr", bartlett = FALSE)
## Not run: simonoff_tsai(mtcars_lm, auxdesign = data.frame(mtcars$wt, mtcars$qsec),
 method = "mlr", hetfun = "logmult")
## End(Not run)