Nothing
### computes the F-bar, LRT-bar and score-bar test statistic ####
# REFs:
# Silvapulle and Sen (2005). Constrained statistical inference. Chapter 2.
# Wolak, F. An exact test for multiple inequality and equality
# constraints in the linear regression model Journal of the American
# statistical association, 1987, 82, 782-793
conTestF.conLM <- function(object, type = "A", neq.alt = 0, boot = "no", R = 9999,
p.distr = rnorm, parallel = "no", ncpus = 1L,
cl = NULL, seed = 1234, verbose = FALSE,
control = NULL, ...) {
# rename for internal use
meq.alt <- neq.alt
# checks
if (!inherits(object, "conLM")) {
stop("Restriktor ERROR: object must be of class conLM.")
}
if (type != "global") {
type <- toupper(type)
}
if(!(type %in% c("A","B","global"))) {
stop("Restriktor ERROR: type must be \"A\", \"B\", or \"global\"")
}
if(!(boot %in% c("no","residual","model.based","parametric"))) {
stop("Restriktor ERROR: boot method unknown.")
}
if (boot == "residual") {
boot <- "model.based"
}
# original model
model.org <- object$model.org
# model matrix
X <- model.matrix(object)[,,drop=FALSE]
# response variable
y <- as.matrix(model.org$model[, attr(model.org$terms, "response")])
# sample size
#n <- dim(X)[1]
# weights
w <- weights(model.org)
# unconstrained df
df.residual <- object$df.residual
# unconstrained covariance matrix
Sigma <- vcov(model.org)
# parameter estimates
b.unrestr <- object$b.unrestr
b.restr <- object$b.restr
b.eqrestr <- NULL
b.restr.alt <- NULL
# length parameter vector
p <- length(b.unrestr)
# variable names
vnames <- names(b.unrestr)
# constraints stuff
Amat <- object$constraints
bvec <- object$rhs
meq <- object$neq
control <- c(object$control, control)
# remove duplicated elements from control list
control <- control[!duplicated(control)]
# get tolerance for control if exists
tol <- ifelse(is.null(control$tol), sqrt(.Machine$double.eps), control$tol)
# check for equalities only
if (meq == nrow(Amat)) {
stop("Restriktor ERROR: test not applicable for object with equality restrictions only.")
}
# check for intercept
intercept <- any(attr(terms(model.org), "intercept"))
if (type == "global") {
if (intercept) {
AmatG <- cbind(rep(0, (p - 1)), diag(rep(1, p - 1)))
} else {
AmatG <- diag(1, p)
for (i in 1:p) {
AmatG[i,i-1] <- -1
}
AmatG <- AmatG[-1,]
}
AmatX <- AmatG %*% (diag(rep(1, p)) - t(Amat) %*%
MASS::ginv(Amat %*% t(Amat)) %*% Amat)
if (all(abs(AmatX) < tol)) {
type <- "A"
attr(type, "type.org") <- "global"
} else {
# remove all rows with only zeros
AmatX <- AmatX[!rowSums(abs(AmatX) < tol) == p,, drop = FALSE]
rAmatX <- GaussianElimination(t(AmatX), tol = tol)
AmatX <- AmatX[rAmatX$pivot,, drop = FALSE]
}
AmatG <- rbind(AmatX, Amat)
bvecG <- c(rep(0, nrow(AmatX)), bvec)
attr(Amat, "Amat.global") <- AmatG
attr(bvec, "bvec.global") <- bvecG
}
if (type == "global") {
# call quadprog
b.eqrestr <- con_solver_lm(X = X,
y = y,
#b.unrestr = b.unrestr,
w = w,
Amat = AmatG,
bvec = bvecG,
meq = nrow(AmatG),
absval = ifelse(is.null(control$absval), 1e-09,
control$absval),
maxit = ifelse(is.null(control$maxit), 1e04,
control$maxit))$qp$solution
# fix estimates < tol to zero
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
# compute global test statistic
Ts <- c( t(b.restr - b.eqrestr) %*% solve(Sigma, b.restr - b.eqrestr) )
} else if (type == "A") {
b.eqrestr <- con_solver_lm(X = X,
y = y,
#b.unrestr = b.unrestr,
w = w,
Amat = Amat,
bvec = bvec,
meq = nrow(Amat),
absval = ifelse(is.null(control$absval), 1e-09,
control$absval),
maxit = ifelse(is.null(control$maxit), 1e04,
control$maxit))$qp$solution
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
# compute test statistic for hypothesis test type A
Ts <- c( t(b.restr - b.eqrestr) %*% solve(Sigma, b.restr - b.eqrestr) )
} else if (type == "B") {
if (meq.alt == 0L) {
# compute test statistic for hypothesis test type B when no equalities are
# preserved in the alternative hypothesis.
Ts <- c( t(b.unrestr - b.restr) %*% solve(Sigma, b.unrestr - b.restr) )
} else {
if (meq.alt > 0L && meq.alt <= meq) {
b.restr.alt <- con_solver_lm(X = X,
y = y,
#b.unrestr = b.unrestr,
w = w,
Amat = Amat[1:meq.alt,,drop = FALSE],
bvec = bvec[1:meq.alt],
meq = meq.alt,
absval = ifelse(is.null(control$absval), 1e-09,
control$absval),
maxit = ifelse(is.null(control$maxit), 1e04,
control$maxit))$qp$solution
b.restr.alt[abs(b.restr.alt) < tol] <- 0L
names(b.restr.alt) <- vnames
# compute test statistic for hypothesis test type B when some equalities may
# be preserved in the alternative hypothesis.
Ts <- c( t(b.restr - b.restr.alt) %*% solve(Sigma, b.restr - b.restr.alt) )
} else {
stop("Restriktor ERROR: neq.alt must not be larger than neq.")
}
}
}
# The test statistics based on inequality constraints are often
# distributed as mixtures of chi-squares. These mixing weights can be computed
# using the multivariate normal distribution with additional Monte Carlo steps
# or via bootstrapping. The pvalue can also be computed directly via
# the parametric bootstrap or model based bootstrap, without fist computing
# the mixing weights.
if (!(attr(object$wt.bar, "method") == "none") && boot == "no") {
wt.bar <- object$wt.bar
pvalue <- con_pvalue_Fbar(wt.bar = wt.bar,
Ts.org = Ts,
df.residual = df.residual,
type = type,
Amat = Amat,
bvec = bvec,
meq = meq,
meq.alt = meq.alt)
attr(pvalue, "wt.bar" ) <- as.numeric(wt.bar)
attr(pvalue, "wt.bar.method" ) <- attr(wt.bar, "method")
attr(pvalue, "converged" ) <- attr(wt.bar, "converged")
attr(pvalue, "convergence_crit" ) <- attr(wt.bar, "convergence_crit")
attr(pvalue, "total_bootstrap_draws" ) <- attr(wt.bar, "total_bootstrap_draws")
attr(pvalue, "error.idx" ) <- attr(wt.bar, "error.idx")
attr(pvalue, "mix_weights_bootstrap_limit") <- attr(wt.bar, "mix_weights_bootstrap_limit")
attr(pvalue, "wt_bar_chunk") <- attr(wt.bar, "wt_bar_chunk")
attr(pvalue, "chunk_size" ) <- attr(wt.bar, "chunk_size_org")
attr(pvalue, "total_chunks") <- attr(wt.bar, "total_chunks")
attr(pvalue, "chunk_iter" ) <- attr(wt.bar, "chunk_iter")
} else if (boot == "parametric") {
if (!is.function(p.distr)) {
p.distr <- get(p.distr, mode = "function")
}
arguments <- list(...)
pnames <- names(formals(p.distr))
pm <- pmatch(names(arguments), pnames, nomatch = 0L)
pm <- names(arguments)[pm > 0L]
formals(p.distr)[pm] <- unlist(arguments[pm])
pvalue <- con_pvalue_boot_parametric(object,
Ts.org = Ts,
type = type,
test = "F",
meq.alt = meq.alt,
R = R,
p.distr = p.distr,
parallel = parallel,
ncpus = ncpus,
cl = cl,
seed = seed,
control = control,
verbose = verbose)
} else if (boot == "model.based") {
pvalue <- con_pvalue_boot_model_based(object,
Ts.org = Ts,
type = type,
test = "F",
meq.alt = meq.alt,
R = R,
parallel = parallel,
ncpus = ncpus,
cl = cl,
seed = seed,
control = control,
verbose = verbose)
} else {
pvalue <- as.numeric(NA)
}
# necessary for the print function
if (!is.null(attr(type, "type.org"))) {
type <- "global"
}
OUT <- list(CON = object$CON,
Amat = Amat,
bvec = bvec,
meq = meq,
meq.alt = meq.alt,
iact = object$iact,
type = type,
test = "F",
Ts = Ts,
df.residual = df.residual,
pvalue = pvalue,
b.eqrestr = b.eqrestr,
b.unrestr = b.unrestr,
b.restr = b.restr,
b.restr.alt = b.restr.alt,
R2.org = object$R2.org,
R2.reduced = object$R2.reduced,
boot = boot,
model.org = model.org)
class(OUT) <- "conTest"
OUT
}
# REF: Silvapulle and Sen (2005). Constrained statistical inference. Chapter 3.
# Wolak, F. An exact test for multiple inequality and equality constraints in the
# linear regression model Journal of the American statistical association, 1987, 82, 782-793
conTestLRT.conLM <- function(object, type = "A", neq.alt = 0, boot = "no", R = 9999,
p.distr = rnorm, parallel = "no", ncpus = 1L,
cl = NULL, seed = 1234, verbose = FALSE,
control = NULL, ...) {
# rename for internal use
meq.alt <- neq.alt
# checks
if (!inherits(object, "conLM")) {
stop("Restriktor ERROR: object must be of class conLM.")
}
if (type != "global") {
type <- toupper(type)
}
if(!(type %in% c("A","B","global"))) {
stop("Restriktor ERROR: type must be \"A\", \"B\", or \"global\"")
}
if(!(boot %in% c("no", "residual", "model.based", "parametric"))) {
stop("Restriktor ERROR: boot method unknown.")
}
if (boot == "residual") {
boot <- "model.based"
}
# original model
model.org <- object$model.org
# model matrix
X <- model.matrix(object)[,,drop=FALSE]
# response variable
y <- as.matrix(model.org$model[, attr(model.org$terms, "response")])
# weights
w <- weights(model.org)
# unconstrained df
df.residual <- object$df.residual
# unconstrained covariance matrix
#Sigma <- vcov(model.org)
# parameter estimates
b.unrestr <- object$b.unrestr
b.restr <- object$b.restr
b.eqrestr <- NULL
b.restr.alt <- NULL
# length parameter vector
p <- length(b.unrestr)
# variable names
vnames <- names(b.unrestr)
# constraints stuff
Amat <- object$constraints
bvec <- object$rhs
meq <- object$neq
control <- c(object$control, control)
# remove duplicated elements from control list
control <- control[!duplicated(control)]
# get tolerance for control if exists
tol <- ifelse(is.null(control$tol), sqrt(.Machine$double.eps), control$tol)
# check for equalities only
if (meq == nrow(Amat)) {
stop("Restriktor ERROR: test not applicable for object with equality restrictions only.")
}
# check for intercept
intercept <- any(attr(terms(model.org), "intercept"))
if (type == "global") {
if (intercept) {
AmatG <- cbind(rep(0, (p - 1)), diag(rep(1, p - 1)))
} else {
AmatG <- diag(1, p)
for (i in 1:p) {
AmatG[i,i-1] <- -1
}
AmatG <- AmatG[-1,]
}
AmatX <- AmatG %*% (diag(rep(1, p)) - t(Amat) %*%
MASS::ginv(Amat %*% t(Amat)) %*% Amat)
if (all(abs(AmatX) < tol)) {
type <- "A"
attr(type, "type.org") <- "global"
} else {
# remove all rows with only zeros
AmatX <- AmatX[!rowSums(abs(AmatX) < tol) == p,, drop = FALSE]
rAmatX <- GaussianElimination(t(AmatX), tol = tol)
AmatX <- AmatX[rAmatX$pivot,, drop = FALSE]
}
AmatG <- rbind(AmatX, Amat)
bvecG <- c(rep(0, nrow(AmatX)), bvec)
attr(Amat, "Amat.global") <- AmatG
attr(bvec, "bvec.global") <- bvecG
}
if (type == "global") {
b.eqrestr <- con_solver_lm(X = X,
y = y,
#b.unrestr = b.unrestr,
w = w,
Amat = AmatG,
bvec = bvecG,
meq = nrow(AmatG),
absval = ifelse(is.null(control$absval), 1e-09,
control$absval),
maxit = ifelse(is.null(control$maxit), 1e04,
control$maxit))$qp$solution
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
#b.eqrestr <- as.vector(b.eqrestr)
fitted0 <- X %*% b.eqrestr
residuals0 <- y - fitted0
object.eqrestr <- list(residuals = residuals0)
object.eqrestr$weights <- object$weights
ll0 <- con_loglik_lm(object.eqrestr)
ll1 <- object$loglik
Ts <- -2*(ll0 - ll1)
} else if (type == "A") {
b.eqrestr <- con_solver_lm(X = X,
y = y,
#b.unrestr = b.unrestr,
w = w,
Amat = Amat,
bvec = bvec,
meq = nrow(Amat),
absval = ifelse(is.null(control$absval), 1e-09,
control$absval),
maxit = ifelse(is.null(control$maxit), 1e04,
control$maxit))$qp$solution
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
#b.eqrestr <- as.vector(b.eqrestr)
fitted0 <- X %*% b.eqrestr
residuals0 <- y - fitted0
object.eqrestr <- list(residuals = residuals0)
object.eqrestr$weights <- object$weights
ll0 <- con_loglik_lm(object.eqrestr)
ll1 <- object$loglik
Ts <- -2*(ll0 - ll1)
} else if (type == "B") {
if (meq.alt == 0L) {
ll0 <- object$loglik
fitted1 <- X %*% object$b.unrestr
residuals1 <- y - fitted1
object.unrestr <- list(residuals = residuals1)
object.unrestr$weights <- object$weights
ll1 <- con_loglik_lm(object.unrestr)
Ts <- -2*(ll0 - ll1)
}
else {
# some equality may be preserved in the alternative hypothesis.
if (meq.alt > 0L && meq.alt <= meq) {
b.restr.alt <- con_solver_lm(X = X,
y = y,
#b.unrestr = b.unrestr,
w = w,
Amat = Amat[1:meq.alt,,drop=FALSE],
bvec = bvec[1:meq.alt],
meq = meq.alt,
absval = ifelse(is.null(control$absval), 1e-09,
control$absval),
maxit = ifelse(is.null(control$maxit), 1e04,
control$maxit))$qp$solution
b.restr.alt[abs(b.restr.alt) < tol] <- 0L
names(b.restr.alt) <- vnames
ll0 <- object$loglik
fitted1 <- X %*% b.restr.alt
residuals1 <- y - fitted1
object.restr.alt <- list(residuals = residuals1)
object.restr.alt$weights <- object$weights
ll1 <- con_loglik_lm(object.restr.alt)
Ts <- -2*(ll0 - ll1)
}
else {
stop("Restriktor ERROR: neq.alt must not be larger than neq.")
}
}
}
if (!(attr(object$wt.bar, "method") == "none") && boot == "no") {
wt.bar <- object$wt.bar
pvalue <- con_pvalue_Fbar(wt.bar = wt.bar,
Ts.org = Ts,
df.residual = df.residual,
type = type,
Amat = Amat,
bvec = bvec,
meq = meq,
meq.alt = meq.alt)
attr(pvalue, "wt.bar" ) <- as.numeric(wt.bar)
attr(pvalue, "wt.bar.method" ) <- attr(wt.bar, "method")
attr(pvalue, "converged" ) <- attr(wt.bar, "converged")
attr(pvalue, "convergence_crit" ) <- attr(wt.bar, "convergence_crit")
attr(pvalue, "total_bootstrap_draws" ) <- attr(wt.bar, "total_bootstrap_draws")
attr(pvalue, "error.idx" ) <- attr(wt.bar, "error.idx")
attr(pvalue, "mix_weights_bootstrap_limit") <- attr(wt.bar, "mix_weights_bootstrap_limit")
attr(pvalue, "wt_bar_chunk") <- attr(wt.bar, "wt_bar_chunk")
attr(pvalue, "chunk_size" ) <- attr(wt.bar, "chunk_size_org")
attr(pvalue, "total_chunks") <- attr(wt.bar, "total_chunks")
attr(pvalue, "chunk_iter" ) <- attr(wt.bar, "chunk_iter")
} else if (boot == "parametric") {
if (!is.function(p.distr)) {
p.distr <- get(p.distr, mode = "function")
}
arguments <- list(...)
pnames <- names(formals(p.distr))
pm <- pmatch(names(arguments), pnames, nomatch = 0L)
pm <- names(arguments)[pm > 0L]
formals(p.distr)[pm] <- unlist(arguments[pm])
pvalue <- con_pvalue_boot_parametric(object,
Ts.org = Ts,
type = type,
test = "LRT",
meq.alt = meq.alt,
R = R,
p.distr = p.distr,
parallel = parallel,
ncpus = ncpus,
cl = cl,
seed = seed,
verbose = verbose)
} else if (boot == "model.based") {
pvalue <- con_pvalue_boot_model_based(object,
Ts.org = Ts,
type = type,
test = "LRT",
meq.alt = meq.alt,
R = R,
parallel = parallel,
ncpus = ncpus,
cl = cl,
seed = seed,
verbose = verbose)
} else {
pvalue <- as.numeric(NA)
}
# necessary for the print function
if (!is.null(attr(type, "type.org"))) {
type <- "global"
}
OUT <- list(CON = object$CON,
Amat = Amat,
bvec = bvec,
meq = meq,
meq.alt = meq.alt,
iact = object$iact,
type = type,
test = "LRT",
Ts = Ts,
df.residual = df.residual,
pvalue = pvalue,
b.eqrestr = b.eqrestr,
b.unrestr = b.unrestr,
b.restr = b.restr,
b.restr.alt = b.restr.alt,
R2.org = object$R2.org,
R2.reduced = object$R2.reduced,
boot = boot,
model.org = model.org)
class(OUT) <- "conTest"
OUT
}
# REF: Robertson, Wright and Dykstra (1988, p. 321). Order constrained statistical inference.
# Global score
conTestScore.conLM <- function(object, type = "A", neq.alt = 0, boot = "no", R = 9999,
p.distr = rnorm, parallel = "no", ncpus = 1L,
cl = NULL, seed = 1234, verbose = FALSE,
control = NULL, ...) {
# rename for internal use
meq.alt <- neq.alt
# checks
if (!inherits(object, "conLM")) {
stop("Restriktor ERROR: object must be of class conLM.")
}
if (type != "global") {
type <- toupper(type)
}
if(!(type %in% c("A","B","global"))) {
stop("Restriktor ERROR: type must be \"A\", \"B\", or \"global\"")
}
if(!(boot %in% c("no", "residual", "model.based", "parametric"))) {
stop("Restriktor ERROR: boot method unknown.")
}
if (boot == "residual") {
boot <- "model.based"
}
# original model
model.org <- object$model.org
# model matrix
X <- model.matrix(object)[,,drop=FALSE]
# response variable
y <- as.matrix(model.org$model[, attr(model.org$terms, "response")])
# unconstrained df
df.residual <- object$df.residual
# unconstrained covariance matrix
#Sigma <- vcov(model.org)
# sample size
n <- dim(X)[1]
# number of parameters
p <- dim(X)[2]
# weights
w <- weights(model.org)
if (is.null(w)) {
w <- rep(1, n)
}
#W <- diag(w)
# parameter estimates
b.unrestr <- object$b.unrestr
b.restr <- object$b.restr
b.eqrestr <- NULL
b.restr.alt <- NULL
# length parameter vector
p <- length(b.unrestr)
# variable names
vnames <- names(b.unrestr)
# restraints stuff
Amat <- object$constraints
bvec <- object$rhs
meq <- object$neq
control <- c(object$control, control)
# remove duplicated elements from control list
control <- control[!duplicated(control)]
# get tolerance for control if exists
tol <- ifelse(is.null(control$tol), sqrt(.Machine$double.eps), control$tol)
# check for equalities only
if (meq == nrow(Amat)) {
stop("Restriktor ERROR: test not applicable for object with equality restrictions only.")
}
# check for intercept
intercept <- any(attr(terms(model.org), "intercept"))
if (type == "global") {
if (intercept) {
AmatG <- cbind(rep(0, (p - 1)), diag(rep(1, p - 1)))
} else {
AmatG <- diag(1, p)
for (i in 1:p) {
AmatG[i,i-1] <- -1
}
AmatG <- AmatG[-1,]
}
AmatX <- AmatG %*% (diag(rep(1, p)) - t(Amat) %*%
MASS::ginv(Amat %*% t(Amat)) %*% Amat)
if (all(abs(AmatX) < tol)) {
type <- "A"
attr(type, "type.org") <- "global"
} else {
# remove all rows with only zeros
AmatX <- AmatX[!rowSums(abs(AmatX) < tol) == p,, drop = FALSE]
rAmatX <- GaussianElimination(t(AmatX), tol = tol)
AmatX <- AmatX[rAmatX$pivot,, drop = FALSE]
}
AmatG <- rbind(AmatX, Amat)
bvecG <- c(rep(0, nrow(AmatX)), bvec)
attr(Amat, "Amat.global") <- AmatG
attr(bvec, "bvec.global") <- bvecG
}
if (type == "global") {
b.eqrestr <- con_solver_lm(X = X,
y = y,
#b.unrestr = b.unrestr,
w = w,
Amat = AmatG,
bvec = bvecG,
meq = nrow(AmatG),
absval = ifelse(is.null(control$absval), 1e-09,
control$absval),
maxit = ifelse(is.null(control$maxit), 1e04,
control$maxit))$qp$solution
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
res0 <- y - X %*% b.eqrestr
res1 <- y - X %*% object$b.restr
# score vector
df0 <- n - (p - nrow(AmatG))
s20 <- sum(res0^2) / df0
df1 <- n - (p - qr(Amat[0:meq,])$rank)
s21 <- sum(res1^2) / df1
G0 <- colSums(as.vector(res0) * w * X) / s20
G1 <- colSums(as.vector(res1) * w * X) / s21
# information matrix under the null-hypothesis
I0 <- 1 / s20 * crossprod(X)#(t(X) %*% X)
# score test-statistic
Ts <- c((G0 - G1) %*% solve(I0, (G0 - G1)))
###############################################
# df0 <- n - (p - nrow(AmatG))
# df1 <- n - (p - qr(Amat[0:meq,])$rank)
# s20 <- sum((y - X %*% b.eqrestr)^2) / df0
# s21 <- sum((y - X %*% b.restr)^2) / df1
# d0 <- 1/s20 * (t(X) %*% (w*(y - X %*% b.eqrestr)))
# d1 <- 1/s21 * (t(X) %*% (w*(y - X %*% b.restr)))
# I0 <- 1/s20 * (t(X) %*% X)
# Ts2 <- t(c(d0 - d1)) %*% solve(I0) %*% c(d0 - d1)
###############################################
} else if (type == "A") {
b.eqrestr <- con_solver_lm(X = X,
y = y,
#b.unrestr = b.unrestr,
w = w,
Amat = Amat,
bvec = bvec,
meq = nrow(Amat),
absval = ifelse(is.null(control$absval), 1e-09,
control$absval),
maxit = ifelse(is.null(control$maxit), 1e04,
control$maxit))$qp$solution
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
res0 <- y - X %*% b.eqrestr
res1 <- y - X %*% object$b.restr
# score vector
df0 <- n - (p - nrow(Amat))
s20 <- sum(res0^2) / df0
df1 <- n - (p - qr(Amat[0:meq,])$rank)
s21 <- sum(res1^2) / df1
G0 <- colSums(as.vector(res0) * w * X) / s20
G1 <- colSums(as.vector(res1) * w * X) / s21
# information matrix under the null-hypothesis
I0 <- 1 / s20 * crossprod(X)#(t(X) %*% X)
# score test-statistic
Ts <- c((G0 - G1) %*% solve(I0, (G0 - G1)))
#############################################
# df0 <- n - (p - nrow(Amat))
# df1 <- n - (p - qr(Amat[0:meq,])$rank)
# s20 <- sum((y - X %*% b.eqrestr)^2) / df0
# s21 <- sum((y - X %*% b.restr)^2) / df1
# d0 <- 1/s20 * (t(X) %*% (w*(y - X %*% b.eqrestr)))
# d1 <- 1/s21 * (t(X) %*% (w*(y - X %*% b.restr)))
# I0 <- 1/s20 * (t(X) %*% X)
# Ts2 <- t(c(d0 - d1)) %*% solve(I0) %*% c(d0 - d1)
###############################################
} else if (type == "B") {
if (meq.alt == 0L) {
res0 <- y - X %*% object$b.restr
res1 <- y - X %*% object$b.unrestr
#res0 <- residuals(object)
#res1 <- residuals(model.org)
# score vector
df0 <- n - (p - qr(Amat[0:meq,])$rank)
s20 <- sum(res0^2) / df0
df1 <- n - p
s21 <- sum(res1^2) / df1
G0 <- colSums(as.vector(res0) * w * X) / s20
G1 <- colSums(as.vector(res1) * w * X) / s21
# information matrix under the null-hypothesis
I0 <- 1 / s20 * crossprod(X)#(t(X) %*% X)
# score test-statistic
Ts <- c((G0 - G1) %*% solve(I0, (G0 - G1)))
#########
# df0 <- n - (p - qr(Amat[0:meq,])$rank)
# df1 <- n - p
# s20 <- sum((y - X %*% b.restr)^2) / df0
# s21 <- sum((y - X %*% b.unrestr)^2) / df1
# d0 <- 1/s20 * (t(X) %*% (w*(y - X %*% b.restr)))
# d1 <- 1/s21 * (t(X) %*% (w*(y - X %*% b.unrestr)))
# I0 <- 1/s20 * (t(X) %*% X)
# Ts2 <- t(c(d0 - d1)) %*% solve(I0) %*% c(d0 - d1)
###############################################
} else {
# some equality may be preserved in the alternative hypothesis.
if (meq.alt > 0L && meq.alt <= meq) {
b.restr.alt <- con_solver_lm(X = X,
y = y,
#b.unrestr = b.unrestr,
w = w,
Amat = Amat[1:meq.alt,,drop=FALSE],
bvec = bvec[1:meq.alt],
meq = meq.alt,
absval = ifelse(is.null(control$absval), 1e-09,
control$absval),
maxit = ifelse(is.null(control$maxit), 1e04,
control$maxit))$qp$solution
b.restr.alt[abs(b.restr.alt) < tol] <- 0L
names(b.restr.alt) <- vnames
res0 <- y - X %*% object$b.restr
res1 <- y - X %*% b.restr.alt
# res0 <- residuals(object)
# res1 <- y - X %*% b.restr.alt
# score vector
df0 <- n - (p - qr(Amat[0:meq,])$rank)
s20 <- sum(res0^2) / df0
df1 <- n - (p - qr(Amat[0:meq,])$rank)
s21 <- sum(res1^2) / df1
G0 <- colSums(as.vector(res0) * w * X) / s20
G1 <- colSums(as.vector(res1) * w * X) / s21
# information matrix
I0 <- 1 / s20 * crossprod(X)#(t(X) %*% X)
# score test-statistic
Ts <- c((G0 - G1) %*% solve(I0, (G0 - G1)))
########################
# df1 <- df0 <- n - (p - qr(Amat[0:meq,])$rank)
# s20 <- sum((y - X %*% b.restr)^2) / df0
# s21 <- sum((y - X %*% b.restr.alt)^2) / df1
# d0 <- 1/s20 * (t(X) %*% (w*(y - X %*% b.restr)))
# d1 <- 1/s21 * (t(X) %*% (w*(y - X %*% b.restr.alt)))
# I0 <- 1/s20 * (t(X) %*% X)
# Ts2 <- t(c(d0 - d1)) %*% solve(I0) %*% c(d0 - d1)
###############################################
}
else {
stop("Restriktor ERROR: neq.alt must not be larger than neq.")
}
}
}
if (!(attr(object$wt.bar, "method") == "none") && boot == "no") {
wt.bar <- object$wt.bar
pvalue <- con_pvalue_Fbar(wt.bar = wt.bar,
Ts.org = Ts,
df.residual = df.residual,
type = type,
Amat = Amat,
bvec = bvec,
meq = meq,
meq.alt = meq.alt)
attr(pvalue, "wt.bar" ) <- as.numeric(wt.bar)
attr(pvalue, "wt.bar.method" ) <- attr(wt.bar, "method")
attr(pvalue, "converged" ) <- attr(wt.bar, "converged")
attr(pvalue, "convergence_crit" ) <- attr(wt.bar, "convergence_crit")
attr(pvalue, "total_bootstrap_draws" ) <- attr(wt.bar, "total_bootstrap_draws")
attr(pvalue, "error.idx" ) <- attr(wt.bar, "error.idx")
attr(pvalue, "mix_weights_bootstrap_limit") <- attr(wt.bar, "mix_weights_bootstrap_limit")
attr(pvalue, "wt_bar_chunk") <- attr(wt.bar, "wt_bar_chunk")
attr(pvalue, "chunk_size" ) <- attr(wt.bar, "chunk_size_org")
attr(pvalue, "total_chunks") <- attr(wt.bar, "total_chunks")
attr(pvalue, "chunk_iter" ) <- attr(wt.bar, "chunk_iter")
} else if (boot == "parametric") {
if (!is.function(p.distr)) {
p.distr <- get(p.distr, mode = "function")
}
arguments <- list(...)
pnames <- names(formals(p.distr))
pm <- pmatch(names(arguments), pnames, nomatch = 0L)
pm <- names(arguments)[pm > 0L]
formals(p.distr)[pm] <- unlist(arguments[pm])
pvalue <- con_pvalue_boot_parametric(object,
Ts.org = Ts,
type = type,
test = "score",
meq.alt = meq.alt,
R = R,
p.distr = p.distr,
parallel = parallel,
ncpus = ncpus,
cl = cl,
seed = seed,
verbose = verbose)
} else if (boot == "model.based") {
pvalue <- con_pvalue_boot_model_based(object,
Ts.org = Ts,
type = type,
test = "score",
meq.alt = meq.alt,
R = R,
parallel = parallel,
ncpus = ncpus,
cl = cl,
seed = seed,
verbose = verbose)
} else {
pvalue <- as.numeric(NA)
}
# necessary for the print function
if (!is.null(attr(type, "type.org"))) {
type <- "global"
}
OUT <- list(CON = object$CON,
Amat = Amat,
bvec = bvec,
meq = meq,
meq.alt = meq.alt,
iact = object$iact,
type = type,
test = "Score",
Ts = Ts,
df.residual = df.residual,
pvalue = pvalue,
b.eqrestr = b.eqrestr,
b.unrestr = b.unrestr,
b.restr = b.restr,
b.restr.alt = b.restr.alt,
R2.org = object$R2.org,
R2.reduced = object$R2.reduced,
boot = boot,
model.org = model.org)
class(OUT) <- "conTest"
OUT
}
## hypothesis test type C is based on a t-distribution.
## intersection-union test
## REF: S. Sasabuchi (1980). A Test of a Multivariate Normal Mean with
## Composite Hypotheses Determined by Linear Inequalities. Biometrika Trust, 67 (2), 429-439.
conTestC.restriktor <- function(object, ...) {
if (!(inherits(object, "restriktor"))) {
stop("Restriktor ERROR: object must be of class restriktor.")
}
Amat <- object$constraints
bvec <- object$rhs
meq <- object$neq
if (meq > 0L) {
stop("Restriktor ERROR: test not applicable with equality restrictions.")
}
b.unrestr <- object$b.unrestr
Sigma <- object$Sigma
df.residual <- object$df.residual
Ts <- as.vector(min((Amat %*% b.unrestr - bvec) /
sqrt(diag(Amat %*% Sigma %*% t(Amat)))))
pvalue <- 1 - pt(Ts, df.residual)
OUT <- list(CON = object$CON,
Amat = Amat,
bvec = bvec,
meq = meq,
meq.alt = 0L,
iact = object$iact,
type = "C",
test = "t",
Ts = Ts,
df.residual = df.residual,
pvalue = pvalue,
b.unrestr = b.unrestr,
b.restr = object$b.restr,
Sigma = Sigma,
R2.org = ifelse(!is.null(object$R2.org), object$R2.org, as.numeric(NA)),
R2.reduced = ifelse(!is.null(object$R2.reduced), object$R2.reduced, as.numeric(NA)),
boot = "no",
model.org = object$model.org)
class(OUT) <- "conTest"
OUT
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.