This function implements the modified profile likelihood ratio test and score test of Simonoff94;textualskedastic for testing for heteroskedasticity in a linear regression model.
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,
...
)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".
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.
A character specifying which of the tests proposed in
Simonoff94;textualskedastic to implement. "mlr"
corresponds to the modified profile likelihood ratio test, and
"score" corresponds to the score test.
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
Cook83;textualskedastic; the additive case is discussed,
inter alia, by Griffiths86;textualskedastic.
Results for the additive and multiplicative models are identical for this
test. Partial matching is used.
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".
A logical specifying whether a Bartlett correction should be
made, as per Ferrari04;textualskedastic, 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".
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.
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).
An object of class "htest". If object is
not assigned, its attributes are displayed in the console as a
tibble using tidy.
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
Simonoff94;textualskedastic) is computed from the modified
profile log-likelihood function using the Nelder-Mead algorithm in
maxLik. 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 Cook83;textualskedastic (a test implemented in
cook_weisberg) or to that of
Koenker81;textualskedastic (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)=\sigma^2 W\), where \(W\) is an \(n\times n\) diagonal matrix with \(i\)th diagonal element \(w_i=w(Z_i, \lambda)\). Here, \(Z_i\) is the \(i\)th 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 \(i\)th element is \(e_i^2\bar{\sigma}^{-2}\), \(\bar{\sigma}^2=n^{-1}e'e\), and \(\bar{J}=(I_n-1_n 1_n'/n)J\).
# NOT RUN {
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")
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab