Nothing
R/ltrend.test.R
In lawstat: Tools for Biostatistics, Public Policy, and Law
#' Test for a Linear Trend in Variances
#'
#' Test for a linear trend in variances.
#'
#' @inherit lnested.test details
#'
#'
#' @inheritParams lnested.test
#' @param score weights to be used in testing an increasing/decreasing trend in
#' group variances, \code{score} coincides by default with \code{group};
#' it can be chosen as a linear, quadratic or any other monotone function.
#'
#'
#' @return A list of class \code{"htest"} containing the following components:
#' \item{statistic}{the value of the test statistic expressed in terms of correlation
#' (Pearson, Kendall, or Spearman).}
#' \item{p.value}{the \eqn{p}-value of the test.}
#' \item{method}{type of test performed.}
#' \item{data.name}{a character string giving the name of the data.}
#' \item{t.statistic}{the value of the test statistic from Student's t-test.}
#' \item{non.bootstrap.p.value}{the \eqn{p}-value of the test without bootstrap method.}
#' \item{log.p.value}{the log of the \eqn{p}-value}
#' \item{log.q.value}{the log of the (one minus the \eqn{p}-value).}
#'
#' @references
#' \insertAllCited{}
#'
#' @seealso \code{\link{neuhauser.hothorn.test}}, \code{\link{levene.test}},
#' \code{\link{lnested.test}}, \code{\link{mma.test}}, \code{\link{robust.mmm.test}}
#'
#' @keywords htest robust variability
#'
#' @author Kimihiro Noguchi, W. Wallace Hui, Yulia R. Gel, Joseph L. Gastwirth, Weiwen Miao
#'
#' @export
#' @examples
#' data(pot)
#' ltrend.test(pot[, "obs"], pot[, "type"], location = "median", tail = "left",
#' correction.method = "zero.correction")
#'
#' ## Bootstrap version of the test. The calculation may take up a few minutes
#' ## depending on the number of bootstrap samples.
#' ltrend.test(pot[, "obs"], pot[, "type"], location = "median", tail = "left",
#' correction.method = "zero.correction",
#' bootstrap = TRUE, num.bootstrap = 500)
#'
ltrend.test <-
function (y,
group,
score = NULL,
location = c("median", "mean", "trim.mean"),
tail = c("right", "left", "both"),
trim.alpha = 0.25,
bootstrap = FALSE,
num.bootstrap = 1000,
correction.method = c("none",
"correction.factor",
"zero.removal",
"zero.correction"),
correlation.method = c("pearson", "kendall", "spearman"))
{
### assign score to each group ###
if (is.null(score))
{
score <- group
}
### stop the code if the length of y does not match the length of group ###
if (length(y) != length(group))
{
stop("the length of the data (y) does not match the length of the group")
}
### assign location, tail, and a correction method ###
location <- match.arg(location)
tail <- match.arg(tail)
correlation.method <- match.arg(correlation.method)
correction.method <- match.arg(correction.method)
DNAME <- deparse(substitute(y))
y <- y[!is.na(y)]
score <- score[!is.na(y)]
group <- group[!is.na(y)]
### stop the code if the location "trim.mean" is selected and trim.alpha is too large ###
if ((location == "trim.mean") & (trim.alpha > 0.5))
{
stop("trim.alpha value of 0 to 0.5 should be provided for the trim.mean location")
}
### sort the order just in case the input is not sorted by group ###
reorder <- order(group)
group <- group[reorder]
y <- y[reorder]
score <- score[reorder]
gr <- score
group <- as.factor(group)
### define the measure of central tendency (mean, median, trimmed mean) ###
if (location == "mean")
{
means <- tapply(y, group, mean)
METHOD <-
"ltrend test based on classical Levene's procedure using the group means"
}
else if (location == "median")
{
means <- tapply(y, group, median)
METHOD <-
"ltrend test based on the modified Brown-Forsythe Levene-type procedure using the group medians"
} else {
location <- "trim.mean"
means <- tapply(y, group, mean, trim = trim.alpha)
METHOD <-
"ltrend test based on the modified Brown-Forsythe Levene-type procedure using the group trimmed means"
}
### calculate the sample size of each group and absolute deviation from center ###
n <- tapply(y, group, length)
ngroup <- n[group]
resp.mean <- abs(y - means[group])
### assign no correction technique if the central tendency is median, and ###
### any technique other than "correction.factor" is chosen ###
if (location != "median" &&
correction.method != "correction.factor")
{
METHOD <-
paste(
METHOD,
"(",
correction.method,
"not applied because the location is not set to median",
")"
)
correction.method <- "none"
}
### multiply the correction factor to each observation if "correction.factor" is chosen ###
if (correction.method == "correction.factor")
{
METHOD <- paste(METHOD, "with correction factor applied")
correction <- 1 / sqrt(1 - 1 / ngroup)
resp.mean <- resp.mean * correction
}
### perform correction techniques for "zero.removal" (Hines and Hines, 2000) ###
### or "zero.correction" (Noguchi and Gel, 2009). ###
if (correction.method == "zero.removal" ||
correction.method == "zero.correction")
{
if (correction.method == "zero.removal")
{
METHOD <-
paste(METHOD,
"with Hines-Hines structural zero removal method")
}
if (correction.method == "zero.correction")
{
METHOD <-
paste(
METHOD,
"with modified structural zero removal method and correction factor"
)
}
### set up variables for calculating the deviation from center ###
resp.mean <- y - means[group]
k <- length(n)
temp <- double()
endpos <- double()
startpos <- double()
### calculate the absolute deviation from mean and remove zeros ###
for (i in 1:k)
{
group.size <- n[i]
j <- i - 1
### calculate the starting and ending index of each group ###
if (i == 1)
start <- 1
else
start <- sum(n[1:j]) + 1
startpos <- c(startpos, start)
end <- sum(n[1:i])
endpos <- c(endpos, end)
### extract the deviation from center for the ith group ###
sub.resp.mean <- resp.mean[start:end]
sub.resp.mean <- sub.resp.mean[order(sub.resp.mean)]
### remove structural zero for the odd-sized group ###
if (group.size %% 2 == 1)
{
mid <- (group.size + 1) / 2
temp2 <- sub.resp.mean[-mid]
### multiply the correction factor for the "zero.correction" option ###
if (correction.method == "zero.correction")
{
ntemp <- length(temp2) + 1
correction <- sqrt((ntemp - 1) / ntemp)
temp2 <- correction * temp2
}
}
### remove structural zero for the even-sized group ###
if (group.size %% 2 == 0)
{
mid <- group.size / 2
### set up the denominator value for the transformation ###
### set sqrt(2) for the "zero.removal" option ###
if (correction.method == "zero.removal")
{
denom <- sqrt(2)
}
### set 1 for the "zero.correction" option ###
else
{
denom <- 1
}
### perform the orthogonal transformation ###
replace1 <- (sub.resp.mean[mid + 1] - sub.resp.mean[mid]) / denom
temp2 <- sub.resp.mean[c(-mid, -mid - 1)]
temp2 <- c(temp2, replace1)
### multiply the correction factor for the "zero.correction" option ###
if (correction.method == "zero.correction")
{
ntemp <- length(temp2) + 1
correction <- sqrt((ntemp - 1) / ntemp)
temp2 <- correction * temp2
}
}
### collect the transformed variables into the vector ###
temp <- c(temp, temp2)
}
### calculate the absolute deviation from center with modifications ###
ngroup2 <- ngroup[-endpos] - 1
resp.mean <- abs(temp)
zero.removal.gr <- gr[-endpos]
}
### set correction.method to be "none" if specified other than those in the option ###
else
{
correction.method = "none"
}
### calculate z ###
mu <- mean(resp.mean)
z <- as.vector(resp.mean - mu)
### set zero.removal.gr as d for methods with structural zero removal ###
if (correction.method == "zero.removal" ||
correction.method == "zero.correction")
{
d <- as.numeric(zero.removal.gr)
}
### set the original gr as d otherwise ###
else
{
d <- as.numeric(gr)
}
### calculate the t-statistic using a simple linear regression ###
t.statistic <- summary(lm(z ~ d))$coefficients[2, 3]
df <- summary(lm(z ~ d))$df[2]
### calculate if the tail is set to "left" ###
if (correlation.method == "pearson")
{
### calculate the correlation between d and z ###
correlation <- cor(d, z, method = "pearson")
if (tail == "left")
{
METHOD <-
paste(METHOD,
"(left-tailed with Pearson correlation coefficient)")
p.value <- pt(t.statistic, df, lower.tail = TRUE)
log.p.value <- pt(t.statistic,
df,
lower.tail = TRUE,
log.p = TRUE)
log.q.value <- pt(t.statistic,
df,
lower.tail = FALSE,
log.p = TRUE)
}
### calculate if the tail is set to "right" ###
else if (tail == "right")
{
METHOD <-
paste(METHOD,
"(right-tailed with Pearson correlation coefficient)")
p.value <- pt(t.statistic, df, lower.tail = FALSE)
log.p.value <- pt(t.statistic,
df,
lower.tail = FALSE,
log.p = TRUE)
log.q.value <- pt(t.statistic,
df,
lower.tail = TRUE,
log.p = TRUE)
}
### calculate if the tail is set to "both" ###
else
{
tail <- "both"
METHOD <-
paste(METHOD,
"(two-tailed with Pearson correlation coefficient)")
p.value <- pt(t.statistic, df, lower.tail = TRUE)
log.p.value <- pt(t.statistic,
df,
lower.tail = TRUE,
log.p = TRUE)
log.q.value <- pt(t.statistic,
df,
lower.tail = FALSE,
log.p = TRUE)
if (p.value >= 0.5)
{
p.value <- pt(t.statistic, df, lower.tail = FALSE)
log.p.value <- pt(t.statistic,
df,
lower.tail = FALSE,
log.p = TRUE)
log.q.value <- pt(t.statistic,
df,
lower.tail = TRUE,
log.p = TRUE)
}
p.value <- p.value * 2
log.p.value <- log.p.value + log(2)
log.q.value <- log(1 - p.value)
}
}
else if (correlation.method == "kendall")
{
### calculate the correlation between d and z ###
correlation <- cor(d, z, method = "kendall")
if (tail == "left")
{
METHOD <-
paste(METHOD,
"(left-tailed with Kendall correlation coefficient)")
p.value.temp <- Kendall(d, z)$sl
if (correlation < 0)
{
p.value <- p.value.temp / 2
}
else
{
p.value <- 1 - p.value.temp / 2
}
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
if (tail == "right")
{
METHOD <-
paste(METHOD,
"(right-tailed with Kendall correlation coefficient)")
p.value.temp <- Kendall(d, z)$sl
if (correlation > 0)
{
p.value <- p.value.temp / 2
}
else
{
p.value <- 1 - p.value.temp / 2
}
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
if (tail == "both")
{
METHOD <-
paste(METHOD,
"(two-tailed with Kendall correlation coefficient)")
p.value <- Kendall(d, z)$sl
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
}
else
{
### calculate the correlation between d and z ###
correlation <- cor(d, z, method = "spearman")
if (tail == "left")
{
METHOD <-
paste(METHOD,
"(left-tailed with Spearman correlation coefficient)")
p.value.temp <- cor.test(d, z, method = "spearman")$p.value
if (correlation < 0)
{
p.value <- p.value.temp / 2
}
else
{
p.value <- 1 - p.value.temp / 2
}
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
if (tail == "right")
{
METHOD <-
paste(METHOD,
"(right-tailed with Spearman correlation coefficient)")
p.value.temp <- cor.test(d, z, method = "spearman")$p.value
if (correlation > 0)
{
p.value <- p.value.temp / 2
}
else
{
p.value <- 1 - p.value.temp / 2
}
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
if (tail == "both")
{
METHOD <-
paste(METHOD,
"(two-tailed with Spearman correlation coefficient)")
p.value <- cor.test(d, z, method = "spearman")$p.value
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
}
### store the non-boostrap p-value ###
non.bootstrap.p.value <- p.value
### perform bootstrapping (followed Lim and Loh(1996)) ###
if (bootstrap == TRUE)
{
METHOD = paste("bootstrap", METHOD)
### step 2 of Lim and Loh (1996): initialize variables ###
R <- 0
N <- length(y)
### step 3 of Lim and Loh (1996): calculate the fractional trimmed mean ###
frac.trim.alpha = 0.2
b.trimmed.mean <- function(y)
{
nn <- length(y)
wt <- rep(0, nn)
y2 <- y[order(y)]
lower <- ceiling(nn * frac.trim.alpha) + 1
upper <- floor(nn * (1 - frac.trim.alpha))
if (lower > upper)
stop("frac.trim.alpha value is too large")
m <- upper - lower + 1
frac <- (nn * (1 - 2 * frac.trim.alpha) - m) / 2
wt[lower - 1] <- frac
wt[upper + 1] <- frac
wt[lower:upper] <- 1
return(weighted.mean(y2, wt))
}
b.trim.means <- tapply(y, group, b.trimmed.mean)
rm <- y - b.trim.means[group]
### step 7 of Lim and Loh (1996): enter a loop ###
for (j in 1:num.bootstrap)
{
### step 4 of Lim and Loh (1996): obtain a bootstrap sample ###
sam <- sample(rm, replace = TRUE)
boot.sample <- sam
### step 5 of Lim and Loh (1996): smooth the variables if n_i < 10 for at least one sample size ###
if (min(n) < 10)
{
U <- runif(1) - 0.5
means <- tapply(y, group, mean)
v <- sqrt(sum((y - means[group]) ^ 2) / N)
boot.sample <- ((12 / 13) ^ (0.5)) * (sam + v * U)
}
### step 6 of Lim and Loh (1996): compute the bootstrap statistic, and increment R to R + 1 if necessary ###
if (location == "mean") {
boot.means <- tapply(boot.sample, group, mean)
} else if (location == "median") {
boot.means <- tapply(boot.sample, group, median)
} else {
location <- "trim.mean"
boot.means <- tapply(boot.sample, group, mean, trim = trim.alpha)
}
### calculate bootstrap statistic ###
resp.boot.mean <- abs(boot.sample - boot.means[group])
if (correction.method == "correction.factor")
{
correction <- 1 / sqrt(1 - 1 / ngroup)
resp.boot.mean <- resp.boot.mean * correction
}
### perform correction techniques for "zero.removal" (Hines and Hines, 2000) ###
### or "zero.correction" (Noguchi and Gel, 2009). ###
if (correction.method == "zero.removal" ||
correction.method == "zero.correction")
{
### set up variables for calculating the deviation from center ###
resp.mean <- boot.sample - boot.means[group]
k <- length(n)
temp <- double()
endpos <- double()
startpos <- double()
### calculate the absolute deviation from mean and remove zeros ###
for (i in 1:k)
{
group.size <- n[i]
j <- i - 1
### calculate the starting and ending index of each group ###
if (i == 1)
start <- 1
else
start <- sum(n[1:j]) + 1
startpos <- c(startpos, start)
end <- sum(n[1:i])
endpos <- c(endpos, end)
### extract the deviation from center for the ith group ###
sub.resp.mean <- resp.mean[start:end]
sub.resp.mean <- sub.resp.mean[order(sub.resp.mean)]
### remove structural zero for the odd-sized group ###
if (group.size %% 2 == 1)
{
mid <- (group.size + 1) / 2
temp2 <- sub.resp.mean[-mid]
### multiply the correction factor for the "zero.correction" option ###
if (correction.method == "zero.correction")
{
ntemp <- length(temp2) + 1
correction <- sqrt((ntemp - 1) / ntemp)
temp2 <- correction * temp2
}
}
### remove structural zero for the even-sized group ###
if (group.size %% 2 == 0)
{
mid <- group.size / 2
### set up the denominator value for the transformation ###
### set sqrt(2) for the "zero.removal" option ###
if (correction.method == "zero.removal")
{
denom <- sqrt(2)
}
### set 1 for the "zero.correction" option ###
else
{
denom <- 1
}
### perform the orthogonal transformation ###
replace1 <- (sub.resp.mean[mid + 1] - sub.resp.mean[mid]) / denom
temp2 <- sub.resp.mean[c(-mid, -mid - 1)]
temp2 <- c(temp2, replace1)
### multiply the correction factor for the "zero.correction" option ###
if (correction.method == "zero.correction")
{
ntemp <- length(temp2) + 1
correction <- sqrt((ntemp - 1) / ntemp)
temp2 <- correction * temp2
}
}
### collect the transformed variables into the vector ###
temp <- c(temp, temp2)
}
### calculate the absolute deviation from center with modifications ###
ngroup2 <- ngroup[-endpos] - 1
resp.boot.mean <- abs(temp)
zero.removal.gr <- gr[-endpos]
}
### set zero.removal.gr as d for methods with structural zero removal ###
if (correction.method == "zero.removal" ||
correction.method == "zero.correction")
{
d <- as.numeric(zero.removal.gr)
}
### set the original gr as d otherwise ###
else
{
d <- as.numeric(gr)
}
boot.mu <- mean(resp.boot.mean)
boot.z <- as.vector(resp.boot.mean - boot.mu)
correlation2 <- cor(boot.z, d, method = correlation.method)
if (tail == "right")
{
if (correlation2 > correlation)
R <- R + 1
}
else if (tail == "left")
{
if (correlation2 < correlation)
R <- R + 1
}
else
{
tail = "both"
if (abs(correlation2) > abs(correlation))
R <- R + 1
}
}
### step 8 of Lim and Loh (1996): calculate the bootstrap p-value ###
p.value <- R / num.bootstrap
}
### display output ###
STATISTIC = correlation
names(STATISTIC) = "Test Statistic (Correlation)"
structure(
list(
statistic = STATISTIC,
p.value = p.value,
method = METHOD,
data.name = DNAME,
t.statistic = t.statistic,
non.bootstrap.p.value = non.bootstrap.p.value,
log.p.value = log.p.value,
log.q.value = log.q.value
),
class = "htest"
)
}
Try the lawstat package in your browser
Any scripts or data that you put into this service are public.
lawstat documentation built on April 6, 2023, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' Test for a Linear Trend in Variances
#'
#' Test for a linear trend in variances.
#'
#' @inherit lnested.test details
#'
#'
#' @inheritParams lnested.test
#' @param score weights to be used in testing an increasing/decreasing trend in
#' group variances, \code{score} coincides by default with \code{group};
#' it can be chosen as a linear, quadratic or any other monotone function.
#'
#'
#' @return A list of class \code{"htest"} containing the following components:
#' \item{statistic}{the value of the test statistic expressed in terms of correlation
#' (Pearson, Kendall, or Spearman).}
#' \item{p.value}{the \eqn{p}-value of the test.}
#' \item{method}{type of test performed.}
#' \item{data.name}{a character string giving the name of the data.}
#' \item{t.statistic}{the value of the test statistic from Student's t-test.}
#' \item{non.bootstrap.p.value}{the \eqn{p}-value of the test without bootstrap method.}
#' \item{log.p.value}{the log of the \eqn{p}-value}
#' \item{log.q.value}{the log of the (one minus the \eqn{p}-value).}
#'
#' @references
#' \insertAllCited{}
#'
#' @seealso \code{\link{neuhauser.hothorn.test}}, \code{\link{levene.test}},
#' \code{\link{lnested.test}}, \code{\link{mma.test}}, \code{\link{robust.mmm.test}}
#'
#' @keywords htest robust variability
#'
#' @author Kimihiro Noguchi, W. Wallace Hui, Yulia R. Gel, Joseph L. Gastwirth, Weiwen Miao
#'
#' @export
#' @examples
#' data(pot)
#' ltrend.test(pot[, "obs"], pot[, "type"], location = "median", tail = "left",
#' correction.method = "zero.correction")
#'
#' ## Bootstrap version of the test. The calculation may take up a few minutes
#' ## depending on the number of bootstrap samples.
#' ltrend.test(pot[, "obs"], pot[, "type"], location = "median", tail = "left",
#' correction.method = "zero.correction",
#' bootstrap = TRUE, num.bootstrap = 500)
#'
ltrend.test <-
function (y,
group,
score = NULL,
location = c("median", "mean", "trim.mean"),
tail = c("right", "left", "both"),
trim.alpha = 0.25,
bootstrap = FALSE,
num.bootstrap = 1000,
correction.method = c("none",
"correction.factor",
"zero.removal",
"zero.correction"),
correlation.method = c("pearson", "kendall", "spearman"))
{
### assign score to each group ###
if (is.null(score))
{
score <- group
}
### stop the code if the length of y does not match the length of group ###
if (length(y) != length(group))
{
stop("the length of the data (y) does not match the length of the group")
}
### assign location, tail, and a correction method ###
location <- match.arg(location)
tail <- match.arg(tail)
correlation.method <- match.arg(correlation.method)
correction.method <- match.arg(correction.method)
DNAME <- deparse(substitute(y))
y <- y[!is.na(y)]
score <- score[!is.na(y)]
group <- group[!is.na(y)]
### stop the code if the location "trim.mean" is selected and trim.alpha is too large ###
if ((location == "trim.mean") & (trim.alpha > 0.5))
{
stop("trim.alpha value of 0 to 0.5 should be provided for the trim.mean location")
}
### sort the order just in case the input is not sorted by group ###
reorder <- order(group)
group <- group[reorder]
y <- y[reorder]
score <- score[reorder]
gr <- score
group <- as.factor(group)
### define the measure of central tendency (mean, median, trimmed mean) ###
if (location == "mean")
{
means <- tapply(y, group, mean)
METHOD <-
"ltrend test based on classical Levene's procedure using the group means"
}
else if (location == "median")
{
means <- tapply(y, group, median)
METHOD <-
"ltrend test based on the modified Brown-Forsythe Levene-type procedure using the group medians"
} else {
location <- "trim.mean"
means <- tapply(y, group, mean, trim = trim.alpha)
METHOD <-
"ltrend test based on the modified Brown-Forsythe Levene-type procedure using the group trimmed means"
}
### calculate the sample size of each group and absolute deviation from center ###
n <- tapply(y, group, length)
ngroup <- n[group]
resp.mean <- abs(y - means[group])
### assign no correction technique if the central tendency is median, and ###
### any technique other than "correction.factor" is chosen ###
if (location != "median" &&
correction.method != "correction.factor")
{
METHOD <-
paste(
METHOD,
"(",
correction.method,
"not applied because the location is not set to median",
")"
)
correction.method <- "none"
}
### multiply the correction factor to each observation if "correction.factor" is chosen ###
if (correction.method == "correction.factor")
{
METHOD <- paste(METHOD, "with correction factor applied")
correction <- 1 / sqrt(1 - 1 / ngroup)
resp.mean <- resp.mean * correction
}
### perform correction techniques for "zero.removal" (Hines and Hines, 2000) ###
### or "zero.correction" (Noguchi and Gel, 2009). ###
if (correction.method == "zero.removal" ||
correction.method == "zero.correction")
{
if (correction.method == "zero.removal")
{
METHOD <-
paste(METHOD,
"with Hines-Hines structural zero removal method")
}
if (correction.method == "zero.correction")
{
METHOD <-
paste(
METHOD,
"with modified structural zero removal method and correction factor"
)
}
### set up variables for calculating the deviation from center ###
resp.mean <- y - means[group]
k <- length(n)
temp <- double()
endpos <- double()
startpos <- double()
### calculate the absolute deviation from mean and remove zeros ###
for (i in 1:k)
{
group.size <- n[i]
j <- i - 1
### calculate the starting and ending index of each group ###
if (i == 1)
start <- 1
else
start <- sum(n[1:j]) + 1
startpos <- c(startpos, start)
end <- sum(n[1:i])
endpos <- c(endpos, end)
### extract the deviation from center for the ith group ###
sub.resp.mean <- resp.mean[start:end]
sub.resp.mean <- sub.resp.mean[order(sub.resp.mean)]
### remove structural zero for the odd-sized group ###
if (group.size %% 2 == 1)
{
mid <- (group.size + 1) / 2
temp2 <- sub.resp.mean[-mid]
### multiply the correction factor for the "zero.correction" option ###
if (correction.method == "zero.correction")
{
ntemp <- length(temp2) + 1
correction <- sqrt((ntemp - 1) / ntemp)
temp2 <- correction * temp2
}
}
### remove structural zero for the even-sized group ###
if (group.size %% 2 == 0)
{
mid <- group.size / 2
### set up the denominator value for the transformation ###
### set sqrt(2) for the "zero.removal" option ###
if (correction.method == "zero.removal")
{
denom <- sqrt(2)
}
### set 1 for the "zero.correction" option ###
else
{
denom <- 1
}
### perform the orthogonal transformation ###
replace1 <- (sub.resp.mean[mid + 1] - sub.resp.mean[mid]) / denom
temp2 <- sub.resp.mean[c(-mid, -mid - 1)]
temp2 <- c(temp2, replace1)
### multiply the correction factor for the "zero.correction" option ###
if (correction.method == "zero.correction")
{
ntemp <- length(temp2) + 1
correction <- sqrt((ntemp - 1) / ntemp)
temp2 <- correction * temp2
}
}
### collect the transformed variables into the vector ###
temp <- c(temp, temp2)
}
### calculate the absolute deviation from center with modifications ###
ngroup2 <- ngroup[-endpos] - 1
resp.mean <- abs(temp)
zero.removal.gr <- gr[-endpos]
}
### set correction.method to be "none" if specified other than those in the option ###
else
{
correction.method = "none"
}
### calculate z ###
mu <- mean(resp.mean)
z <- as.vector(resp.mean - mu)
### set zero.removal.gr as d for methods with structural zero removal ###
if (correction.method == "zero.removal" ||
correction.method == "zero.correction")
{
d <- as.numeric(zero.removal.gr)
}
### set the original gr as d otherwise ###
else
{
d <- as.numeric(gr)
}
### calculate the t-statistic using a simple linear regression ###
t.statistic <- summary(lm(z ~ d))$coefficients[2, 3]
df <- summary(lm(z ~ d))$df[2]
### calculate if the tail is set to "left" ###
if (correlation.method == "pearson")
{
### calculate the correlation between d and z ###
correlation <- cor(d, z, method = "pearson")
if (tail == "left")
{
METHOD <-
paste(METHOD,
"(left-tailed with Pearson correlation coefficient)")
p.value <- pt(t.statistic, df, lower.tail = TRUE)
log.p.value <- pt(t.statistic,
df,
lower.tail = TRUE,
log.p = TRUE)
log.q.value <- pt(t.statistic,
df,
lower.tail = FALSE,
log.p = TRUE)
}
### calculate if the tail is set to "right" ###
else if (tail == "right")
{
METHOD <-
paste(METHOD,
"(right-tailed with Pearson correlation coefficient)")
p.value <- pt(t.statistic, df, lower.tail = FALSE)
log.p.value <- pt(t.statistic,
df,
lower.tail = FALSE,
log.p = TRUE)
log.q.value <- pt(t.statistic,
df,
lower.tail = TRUE,
log.p = TRUE)
}
### calculate if the tail is set to "both" ###
else
{
tail <- "both"
METHOD <-
paste(METHOD,
"(two-tailed with Pearson correlation coefficient)")
p.value <- pt(t.statistic, df, lower.tail = TRUE)
log.p.value <- pt(t.statistic,
df,
lower.tail = TRUE,
log.p = TRUE)
log.q.value <- pt(t.statistic,
df,
lower.tail = FALSE,
log.p = TRUE)
if (p.value >= 0.5)
{
p.value <- pt(t.statistic, df, lower.tail = FALSE)
log.p.value <- pt(t.statistic,
df,
lower.tail = FALSE,
log.p = TRUE)
log.q.value <- pt(t.statistic,
df,
lower.tail = TRUE,
log.p = TRUE)
}
p.value <- p.value * 2
log.p.value <- log.p.value + log(2)
log.q.value <- log(1 - p.value)
}
}
else if (correlation.method == "kendall")
{
### calculate the correlation between d and z ###
correlation <- cor(d, z, method = "kendall")
if (tail == "left")
{
METHOD <-
paste(METHOD,
"(left-tailed with Kendall correlation coefficient)")
p.value.temp <- Kendall(d, z)$sl
if (correlation < 0)
{
p.value <- p.value.temp / 2
}
else
{
p.value <- 1 - p.value.temp / 2
}
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
if (tail == "right")
{
METHOD <-
paste(METHOD,
"(right-tailed with Kendall correlation coefficient)")
p.value.temp <- Kendall(d, z)$sl
if (correlation > 0)
{
p.value <- p.value.temp / 2
}
else
{
p.value <- 1 - p.value.temp / 2
}
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
if (tail == "both")
{
METHOD <-
paste(METHOD,
"(two-tailed with Kendall correlation coefficient)")
p.value <- Kendall(d, z)$sl
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
}
else
{
### calculate the correlation between d and z ###
correlation <- cor(d, z, method = "spearman")
if (tail == "left")
{
METHOD <-
paste(METHOD,
"(left-tailed with Spearman correlation coefficient)")
p.value.temp <- cor.test(d, z, method = "spearman")$p.value
if (correlation < 0)
{
p.value <- p.value.temp / 2
}
else
{
p.value <- 1 - p.value.temp / 2
}
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
if (tail == "right")
{
METHOD <-
paste(METHOD,
"(right-tailed with Spearman correlation coefficient)")
p.value.temp <- cor.test(d, z, method = "spearman")$p.value
if (correlation > 0)
{
p.value <- p.value.temp / 2
}
else
{
p.value <- 1 - p.value.temp / 2
}
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
if (tail == "both")
{
METHOD <-
paste(METHOD,
"(two-tailed with Spearman correlation coefficient)")
p.value <- cor.test(d, z, method = "spearman")$p.value
q.value <- 1 - p.value
log.p.value <- log(p.value)
log.q.value <- log(q.value)
}
}
### store the non-boostrap p-value ###
non.bootstrap.p.value <- p.value
### perform bootstrapping (followed Lim and Loh(1996)) ###
if (bootstrap == TRUE)
{
METHOD = paste("bootstrap", METHOD)
### step 2 of Lim and Loh (1996): initialize variables ###
R <- 0
N <- length(y)
### step 3 of Lim and Loh (1996): calculate the fractional trimmed mean ###
frac.trim.alpha = 0.2
b.trimmed.mean <- function(y)
{
nn <- length(y)
wt <- rep(0, nn)
y2 <- y[order(y)]
lower <- ceiling(nn * frac.trim.alpha) + 1
upper <- floor(nn * (1 - frac.trim.alpha))
if (lower > upper)
stop("frac.trim.alpha value is too large")
m <- upper - lower + 1
frac <- (nn * (1 - 2 * frac.trim.alpha) - m) / 2
wt[lower - 1] <- frac
wt[upper + 1] <- frac
wt[lower:upper] <- 1
return(weighted.mean(y2, wt))
}
b.trim.means <- tapply(y, group, b.trimmed.mean)
rm <- y - b.trim.means[group]
### step 7 of Lim and Loh (1996): enter a loop ###
for (j in 1:num.bootstrap)
{
### step 4 of Lim and Loh (1996): obtain a bootstrap sample ###
sam <- sample(rm, replace = TRUE)
boot.sample <- sam
### step 5 of Lim and Loh (1996): smooth the variables if n_i < 10 for at least one sample size ###
if (min(n) < 10)
{
U <- runif(1) - 0.5
means <- tapply(y, group, mean)
v <- sqrt(sum((y - means[group]) ^ 2) / N)
boot.sample <- ((12 / 13) ^ (0.5)) * (sam + v * U)
}
### step 6 of Lim and Loh (1996): compute the bootstrap statistic, and increment R to R + 1 if necessary ###
if (location == "mean") {
boot.means <- tapply(boot.sample, group, mean)
} else if (location == "median") {
boot.means <- tapply(boot.sample, group, median)
} else {
location <- "trim.mean"
boot.means <- tapply(boot.sample, group, mean, trim = trim.alpha)
}
### calculate bootstrap statistic ###
resp.boot.mean <- abs(boot.sample - boot.means[group])
if (correction.method == "correction.factor")
{
correction <- 1 / sqrt(1 - 1 / ngroup)
resp.boot.mean <- resp.boot.mean * correction
}
### perform correction techniques for "zero.removal" (Hines and Hines, 2000) ###
### or "zero.correction" (Noguchi and Gel, 2009). ###
if (correction.method == "zero.removal" ||
correction.method == "zero.correction")
{
### set up variables for calculating the deviation from center ###
resp.mean <- boot.sample - boot.means[group]
k <- length(n)
temp <- double()
endpos <- double()
startpos <- double()
### calculate the absolute deviation from mean and remove zeros ###
for (i in 1:k)
{
group.size <- n[i]
j <- i - 1
### calculate the starting and ending index of each group ###
if (i == 1)
start <- 1
else
start <- sum(n[1:j]) + 1
startpos <- c(startpos, start)
end <- sum(n[1:i])
endpos <- c(endpos, end)
### extract the deviation from center for the ith group ###
sub.resp.mean <- resp.mean[start:end]
sub.resp.mean <- sub.resp.mean[order(sub.resp.mean)]
### remove structural zero for the odd-sized group ###
if (group.size %% 2 == 1)
{
mid <- (group.size + 1) / 2
temp2 <- sub.resp.mean[-mid]
### multiply the correction factor for the "zero.correction" option ###
if (correction.method == "zero.correction")
{
ntemp <- length(temp2) + 1
correction <- sqrt((ntemp - 1) / ntemp)
temp2 <- correction * temp2
}
}
### remove structural zero for the even-sized group ###
if (group.size %% 2 == 0)
{
mid <- group.size / 2
### set up the denominator value for the transformation ###
### set sqrt(2) for the "zero.removal" option ###
if (correction.method == "zero.removal")
{
denom <- sqrt(2)
}
### set 1 for the "zero.correction" option ###
else
{
denom <- 1
}
### perform the orthogonal transformation ###
replace1 <- (sub.resp.mean[mid + 1] - sub.resp.mean[mid]) / denom
temp2 <- sub.resp.mean[c(-mid, -mid - 1)]
temp2 <- c(temp2, replace1)
### multiply the correction factor for the "zero.correction" option ###
if (correction.method == "zero.correction")
{
ntemp <- length(temp2) + 1
correction <- sqrt((ntemp - 1) / ntemp)
temp2 <- correction * temp2
}
}
### collect the transformed variables into the vector ###
temp <- c(temp, temp2)
}
### calculate the absolute deviation from center with modifications ###
ngroup2 <- ngroup[-endpos] - 1
resp.boot.mean <- abs(temp)
zero.removal.gr <- gr[-endpos]
}
### set zero.removal.gr as d for methods with structural zero removal ###
if (correction.method == "zero.removal" ||
correction.method == "zero.correction")
{
d <- as.numeric(zero.removal.gr)
}
### set the original gr as d otherwise ###
else
{
d <- as.numeric(gr)
}
boot.mu <- mean(resp.boot.mean)
boot.z <- as.vector(resp.boot.mean - boot.mu)
correlation2 <- cor(boot.z, d, method = correlation.method)
if (tail == "right")
{
if (correlation2 > correlation)
R <- R + 1
}
else if (tail == "left")
{
if (correlation2 < correlation)
R <- R + 1
}
else
{
tail = "both"
if (abs(correlation2) > abs(correlation))
R <- R + 1
}
}
### step 8 of Lim and Loh (1996): calculate the bootstrap p-value ###
p.value <- R / num.bootstrap
}
### display output ###
STATISTIC = correlation
names(STATISTIC) = "Test Statistic (Correlation)"
structure(
list(
statistic = STATISTIC,
p.value = p.value,
method = METHOD,
data.name = DNAME,
t.statistic = t.statistic,
non.bootstrap.p.value = non.bootstrap.p.value,
log.p.value = log.p.value,
log.q.value = log.q.value
),
class = "htest"
)
}
Try the lawstat package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.