Nothing
### computes the F-bar, Wald(2)-bar and score-bar test statistic ####
conTestF.conRLM <- 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 (!("conRLM" %in% class(object))) {
stop("Restriktor ERROR: object must be of class conRLM.")
}
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
# orignal model call
call.org <- as.list(model.org$call)
# model matrix
X <- model.matrix(model.org)[,,drop=FALSE]
# response variable
y <- as.matrix(model.org$model[, attr(model.org$terms, "response")])
# weights
weights <- model.org$weights
# unconstrained df
df.residual <- object$df.residual
# unconstrained scale
scale <- model.org$s
# 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 restriktions only.")
}
call.org$formula <- call.org$data <- call.org$weights <-
call.org$subset <- call.org$na.action <- call.org$model <-
call.org$x.ret <- call.org$y.ret <- call.org$contrasts <-
call.org$X <- call.org$y <- NULL
# 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 <- c(call.org, list(x = X, y = y, weights = weights,
Amat = AmatG, bvec = bvecG,
meq = nrow(AmatG), tol = tol))
rfit <- do.call("conRLM_fit", CALL)
b.eqrestr <- rfit$coefficients
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
Ts <- robustFm(x = X,
y = y,
b.unrestr = b.unrestr,
b.eqrestr = b.eqrestr,
b.restr = b.restr,
scale = scale,
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))$Ts
} else if (type == "A") {
CALL <- c(call.org, list(x = X, y = y, weights = weights,
Amat = Amat, bvec = bvec,
meq = nrow(Amat), tol = tol))
rfit <- do.call("conRLM_fit", CALL)
b.eqrestr <- rfit$coefficients
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
Ts <- robustFm(x = X,
y = y,
b.unrestr = b.unrestr,
b.eqrestr = b.eqrestr,
b.restr = b.restr,
scale = scale,
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))$Ts
} else if (type == "B") {
if (meq.alt == 0L) {
Ts <- robustFm(x = X,
y = y,
b.unrestr = b.unrestr,
b.eqrestr = b.restr,
b.restr = b.unrestr,
scale = scale,
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))$Ts
} else {
# some equality may be preserved in the alternative hypothesis.
if (meq.alt > 0L && meq.alt <= meq) {
CALL <- c(call.org, list(x = X, y = y, weights = weights,
Amat = Amat[1:meq.alt, , drop = FALSE],
meq = meq.alt, bvec = bvec[1:meq.alt],
tol = tol))
rfit <- do.call("conRLM_fit", CALL)
b.restr.alt <- rfit$coefficients
b.restr.alt[abs(b.restr.alt) < tol] <- 0L
names(b.restr.alt) <- vnames
Ts <- robustFm(x = X,
y = y,
b.unrestr = b.unrestr,
b.eqrestr = b.restr,
b.restr = b.restr.alt,
scale = scale,
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))$Ts
} 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 = "F",
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 = "F",
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 = "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, M.J. (1992) Robust Wald-Type Test of One-Sided Hypotheses in the linear model.
#American Statistical Association, 87, 417.
conTestWald.conRLM <- 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 (!("conRLM" %in% class(object))) {
stop("Restriktor ERROR: object must be of class conRLM.")
}
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
# origianl model call
call.org <- as.list(model.org$call)
# model matrix
X <- model.matrix(model.org)[,,drop=FALSE]
# response variable
y <- as.matrix(model.org$model[, attr(model.org$terms, "response")])
# weights
weights <- model.org$weights
# unconstrained df
df.residual <- object$df.residual
# unconstrained scale
scale <- model.org$s
# 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 restriktions only.")
}
call.org$formula <- call.org$data <- call.org$weights <-
call.org$subset <- call.org$na.action <- call.org$model <-
call.org$x.ret <- call.org$y.ret <- call.org$contrasts <-
call.org$X <- call.org$y <- NULL
# 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 <- c(call.org, list(x = X, y = y, weights = weights,
Amat = AmatG, bvec = bvecG,
meq = nrow(AmatG), tol = tol))
rfit <- do.call("conRLM_fit", CALL)
b.eqrestr <- rfit$coefficients
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
Ts <- robustWald(x = X,
y = y,
b.eqrestr = b.eqrestr,
b.restr = b.restr,
b.unrestr = b.unrestr,
Amat = AmatG,
scale = scale,
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))$Ts
} else if (type == "A") {
CALL <- c(call.org, list(x = X, y = y, weights = weights,
Amat = Amat, bvec = bvec,
meq = nrow(Amat), tol = tol))
rfit <- do.call("conRLM_fit", CALL)
b.eqrestr <- rfit$coefficients
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
Ts <- robustWald(x = X,
y = y,
b.eqrestr = b.eqrestr,
b.restr = b.restr,
b.unrestr = b.unrestr,
Amat = Amat,
scale = scale,
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))$Ts
} else if (type == "B") {
if (meq.alt == 0L) {
Ts <- robustWald(x = X,
y = y,
b.eqrestr = b.restr,
b.restr = b.unrestr,
b.unrestr = b.unrestr,
Amat = Amat,
scale = scale,
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))$Ts
} else {
# some equality may be preserved in the alternative hypothesis.
if (meq.alt > 0L && meq.alt <= meq) {
CALL <- c(call.org, list(x = X, y = y, weights = weights,
Amat = Amat[1:meq.alt, , drop = FALSE],
meq = meq.alt, bvec = bvec[1:meq.alt],
tol = tol))
rfit <- do.call("conRLM_fit", CALL)
b.restr.alt <- rfit$coefficients
b.restr.alt[abs(b.restr.alt) < tol] <- 0L
names(b.restr.alt) <- vnames
Ts <- robustWald(x = X,
y = y,
b.eqrestr = b.restr,
b.restr = b.restr.alt,
b.unrestr = b.restr.alt,
Amat = Amat,
scale = scale,
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))$Ts
} 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 = "Wald",
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 = "Wald",
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 = "Wald",
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
}
conTestScore.conRLM <- 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 (!("conRLM" %in% class(object))) {
stop("Restriktor ERROR: object must be of class conRLM.")
}
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
# original model call
call.org <- as.list(model.org$call)
# model matrix
X <- model.matrix(model.org)[,,drop=FALSE]
# response variable
y <- as.matrix(model.org$model[, attr(model.org$terms, "response")])
# weights
weights <- model.org$weights
# unconstrained df
df.residual <- object$df.residual
# unconstrained scale
scale <- model.org$s
# parameter estimates
b.unrestr <- object$b.unrestr
b.restr <- object$b.restr
b.eqrestr <- NULL
b.restr.alt <- NULL
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 restriktions only.")
}
call.org$formula <- call.org$data <- call.org$weights <-
call.org$subset <- call.org$na.action <- call.org$model <-
call.org$x.ret <- call.org$y.ret <- call.org$contrasts <-
call.org$X <- call.org$y <- NULL
# 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 <- c(call.org, list(x = X, y = y, weights = weights,
Amat = AmatG, bvec = bvecG,
meq = nrow(AmatG), tol = tol))
rfit <- do.call("conRLM_fit", CALL)
b.eqrestr <- rfit$coefficients
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
out0 <- robustScores(x = X,
y = y,
b.eqrestr = b.eqrestr,
b.restr = b.restr,
b.unrestr = b.unrestr,
Amat = AmatG,
scale = scale,
test = "score",
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))
Ts <- out0$Ts
V <- out0$V
} else if (type == "A") {
CALL <- c(call.org, list(x = X, y = y, weights = weights,
Amat = Amat, bvec = bvec,
meq = nrow(Amat), tol = tol))
rfit <- do.call("conRLM_fit", CALL)
b.eqrestr <- rfit$coefficients
b.eqrestr[abs(b.eqrestr) < tol] <- 0L
names(b.eqrestr) <- vnames
out1 <- robustScores(x = X,
y = y,
b.eqrestr = b.eqrestr,
b.restr = b.restr,
b.unrestr = b.unrestr,
Amat = Amat,
scale = scale,
test = "score",
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))
Ts <- out1$Ts
V <- out1$V
} else if (type == "B") {
if (meq.alt == 0L) {
out2 <- robustScores(x = X,
y = y,
b.eqrestr = b.restr,
b.restr = b.unrestr,
b.unrestr = b.unrestr,
scale = scale,
Amat = Amat,
test = "score",
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))
Ts <- out2$Ts
V <- out2$V
} else {
# some equality may be preserved in the alternative hypothesis.
if (meq.alt > 0L && meq.alt <= meq) {
CALL <- c(call.org, list(x = X, y = y, weights = weights,
Amat = Amat[1:meq.alt,,drop=FALSE],
meq = meq.alt,
bvec = bvec[1:meq.alt], tol = tol))
rfit <- do.call("conRLM_fit", CALL)
b.restr.alt <- rfit$coefficients
b.restr.alt[abs(b.restr.alt) < tol] <- 0L
names(b.restr.alt) <- vnames
out3 <- robustScores(x = X,
y = y,
b.eqrestr = b.restr,
b.restr = b.restr.alt,
b.unrestr = b.restr.alt,
scale = scale,
Amat = Amat,
test = "score",
cc = ifelse(is.null(call.org$c), 4.685061, call.org$c))
Ts <- out3$Ts
V <- out3$V
} else {
stop("Restriktor ERROR: neq.alt must not be larger than neq.")
}
}
}
# We need to recalculate the weights based on V.hat = Sigma instead on solve(t(X)%*%X)
# Is this really needed? The differences look very small.
## compute mixing weights
if (!(attr(object$wt.bar, "method") == "none")) {
if (nrow(Amat) == meq) {
# equality constraints only
wt.bar <- rep(0L, ncol(V) + 1)
wt.bar.idx <- ncol(V) - meq + 1
wt.bar[wt.bar.idx] <- 1
} else if (attr(object$wt.bar, "method") == "boot") {
# compute chi-square-bar weights based on Monte Carlo simulation
wt.bar <- con_weights_boot(VCOV = V,
Amat = Amat,
meq = meq,
R = attr(object$wt.bar, "mix_weights_bootstrap_limit"),
seed = seed,
convergence_crit = ifelse(is.null(control$convergence_crit),
1e-03, control$convergence_crit),
chunk_size = ifelse(is.null(control$chunk_size),
5000L, control$chunk_size),
verbose = verbose)
attr(wt.bar, "mix_weights_bootstrap_limit") <- attr(object$wt.bar, "mix_weights_bootstrap_limit")
} else if (attr(object$wt.bar, "method") == "pmvnorm" && (meq < nrow(Amat))) {
# compute chi-square-bar weights based on pmvnorm
wt.bar <- rev(con_weights(Amat %*% V %*% t(Amat), meq = meq,
tolerance = ifelse(is.null(control$tolerance), 1e-15,
control$tolerance),
ridge_constant = ifelse(is.null(control$ridge_constant), 1e-05,
control$ridge_constant),
...))
}
} else {
wt.bar <- as.numeric(NA)
}
attr(wt.bar, "method") <- attr(object$wt.bar, "method")
if (!attr(object$wt.bar, "method") == "none" && boot == "no") {
# compute pvalue based on F-distribution
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,
V = V,
R2.org = object$R2.org,
R2.reduced = object$R2.reduced,
boot = boot,
model.org = 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.