Test statistics of normal distribution-based order-related likelihood ratio tests are often distributed as mixtures of chi-square or beta-distributions with different parameters. These functions determine the mixing weights and the cumulative distribution functions based on these. They can be directly used and are called by function ic.test.
ic.weights(corr, ...)
pchibar(x, df, wt)
pbetabar(x, df1, df2, wt)ic.weights returns the vector of weights,
pchibar and pchibar return the cumulative probability of the
respective distribution.
Function ic.weights relies on package mvtnorm for determining
multivariate normal rectangle probabilities. Note that these calculations
involve Monte Carlo steps so that these weights are not completely repeatable.
corr is the correlation or covariance matrix
(or any multiple thereof) of the data or coefficients
for which weights are to be calculated
... contains further arguments to be given to function
pmvnorm of package mvtnorm for calculating
multivariate normal rectangle probabilities; it is possible
to select an algorithm (default in current version of mvtnorm:
algorithm = GenzBretz()) and
to tune weight accuracy by modifying including additional
parameters into the algorithm specification,
cf. help for GenzBretz
x is the quantile for which the distribution function
is to be calculated
is the vector of the degrees of freedom for the chi-square
distributions that are mixed into the chibar-square-distribution
with the proportions given in wt
each element of wt is the mixing weight of the chi-square
distribution with df as in the corresponding element of df;
such weights can be calculated with function ic.weights
vector of first parameters of the beta-distributions to be mixed into the betabar-distribution
second parameter of the beta-distributions to be
mixed into the betabar-distribution; error degrees of freedom
in the tests implemented for linear models in summary.orlm;
NOTE: see details for the (perhaps unexpected) constancy of df2
Ulrike Groemping, BHT Berlin
Function ic.weights uses results by Kudo (1963)
regarding the calculation of the weights. The weights are the probabilities that
the projection along its covariance onto the non-negative orthant
of a multivariate normal random vector with expectation 0 and
correlation corr lies in faces of dimensions nrow(corr):1
(in this order). It is known that these probabilities coincide with
various other useful probabilities related to order-related hypothesis testing,
cf. e.g. Shapiro (1988). Calculation of the weights involves various calls
to function pmvnorm from package mvtnorm.
Functions pchibar (taken from package ibdreg) and pbetabar
calculate cumulative probabilities from mixtures of chi-square and
beta-distributions, respectively.
IMPORTANT: Contrary to likelihood ratio theory in linear models, the beta
distributions mixed always use the error sum of squares from the unrestricted model,
i.e. the smallest possible error sum of squares with a fixed no. of df. Therefore,
the second df entry is not increased when decreasing the first!
This is appropriate for the test statistics calculated by functions ic.test
or summary.orlm, but not necessarily for test statistics obtained elsewhere.
Kudo, A. (1963) A multivariate analogue of the one-sided test. Biometrika 50, 403--418
Shapiro, A. (1988) Towards a unified theory of inequality-constrained testing in multivariate analysis. International Statistical Review 56, 49--62
Silvapulle, M.J. and Sen, P.K. (2004) Constrained Statistical Inference. Wiley, New York
ic.test, orlm, pmvnorm,
GenzBretz
z <- 0.5
corr <- matrix(c(1,0.9,0.9,1),2,2)
print(wt.plus <- ic.weights(corr))
T <- c(z,z)%*%solve(corr,c(z,z))
1-pchibar(T,2:0,wt.plus)
1-pbetabar(T/(T+10),2:0,10,wt.plus)
corr <- matrix(c(1,0,0,1),2,2)
print(wt.0 <- ic.weights(corr))
T <- c(z,z)%*%solve(corr,c(z,z))
1-pchibar(T,2:0,wt.0)
1-pbetabar(T/(T+10),2:0,10,wt.0)
corr <- matrix(c(1,-0.9,-0.9,1),2,2)
print(wt.minus <- ic.weights(corr))
T <- c(z,z)%*%solve(corr,c(z,z))
1-pchibar(T,2:0,wt.minus)
1-pbetabar(T/(T+10),2:0,10,wt.minus)
Run the code above in your browser using DataLab