Nothing
##' Nonparametric multiple contrast tests and simultaneous confidence intervals
##' (repeated measures)
##'
##' The function mctp.rm computes the estimator of nonparametric relative
##' effects based on global rankings, simultaneous confidence intervals for the
##' effects, and adjusted p-values based on contrasts in the setting of a
##' repeated measures design with n independent individuals and d repeated
##' measures. Contrasts include "Tukey", "Dunnett", "Sequen", "Williams",
##' "Changepoint", "AVE", "McDermott", "Marcus", "UmbrellaWilliams",
##' "GrandMean", and "UserDefined". The statistics are computed using
##' multivariate normal distribution, multivariate Satterthwaite
##' t-Approximation, and multivariate transformations (adjusted log odds or
##' Fisher function). The function 'mctp.rm' computes both the one-sided and
##' two-sided simultaneous confidence intervals and adjusted p-values. The
##' confidence intervals can be plotted.
##'
##'
##' @param formula A two-sided 'formula' specifying a numeric response variable
##' and a repeated measures factor with more than two levels. If the factor
##' contains less than 3 levels, an error message will be returned.
##' @param data A dataframe containing the variables specified in formula.
##' @param type Character string defining the type of contrast. It should be
##' one of "Tukey", "Dunnett", "Sequen", "Williams", "Changepoint", "AVE",
##' "McDermott", "Marcus", "UmbrellaWilliams", "GrandMean", "UserDefined".
##' @param conf.level The confidence level for \code{conf.level}-confidence
##' intervals (default is 0.95).
##' @param alternative Character string defining the alternative hypothesis,
##' one of "two.sided", "less", or "greater".
##' @param asy.method Character string defining the asymptotic approximation
##' method, one of "log.odds" (for using the adjusted log odds effect sizes),
##' "mult.t" (for using a multivariate t-distribution with a Satterthwaite
##' Approximation), "fisher" (for using the Fisher transformation function), or
##' "normal" (for using the multivariate normal distribution).
##' @param plot.simci A logical indicating whether you want a plot of the
##' confidence intervals.
##' @param control Character string defining the control group in Dunnett
##' comparisons. By default, it is the first group by definition of the factor
##' variable.
##' @param info A logical whether you want a brief overview with informations
##' about the output.
##' @param rounds Number of rounds for the numeric values of the output
##' (default is 3).
##' @param contrast.matrix User-defined contrast matrix.
##' @param correlation A logical whether the estimated correlation matrix and
##' covariance matrix should be printed.
##' @param const Number used for the adjustment of log odds when the "log.odds"
##' option is chosen.
##' @return
##'
##' \item{Data.Info}{List of samples and sample sizes and estimated effect per
##' repeated measures level.}
##' \item{Contrast}{Contrast matrix.}
##' \item{Analysis}{ Estimator: Estimated relative effect,
##' Lower: Lower limit of the simultaneous confidence intervals,
##' Upper: Upper limit of the simultaneous confidence intervals,
##' Statistic: Test statistic
##' p.Value: Adjusted p-values for the hypothesis by the choosen approximation method.}
##' \item{Analysis.Inf}{The same as \code{Analysis} except that it assumes
##' \code{rounds = Inf}.}
##' \item{Overall}{The critical value and adjusted p-value for the overall hypothesis.}
##' \item{input}{List of input arguments by user.}
##' \item{text.Output}{Character string specifying the alternative hypotheses.}
##' \item{connames}{Character string specifying the contrast names.}
##' \item{AsyMethod}{Character string specifying the approximation method.}
##' @note Estimated relative effects with 0 or 1 are replaced with 0.001 and
##' 0.999.
##'
##' A summary and a graph can be created separately by using the functions
##' \code{\link{summary.mctp.rm}} and \code{\link{plot.mctp.rm}}.
##'
##' For the analysis, the R packages 'multcomp' and 'mvtnorm' are required.
##' @author Marius Placzek, Kimihiro Noguchi
##' @seealso To analyse simple one-way layouts with independent samples use
##' \code{\link{mctp}}.
##' @references F. Konietschke, A.C. Bathke, L.A. Hothorn, E. Brunner: Testing
##' and estimation of purely nonparametric effects in repeated measures
##' designs. Computational Statistics and Data Analysis 54 (2010) 1895-1905.
##'
##' Konietschke, F., Placzek, M., Schaarschmidt, S., Hothorn, L.A. (2015).
##' nparcomp: An R Software Package for Nonparametric Multiple Comparisons and
##' Simultaneous Confidence Intervals. Journal of Statistical Software, 61(10),
##' 1-17.
##'
##' Noguchi, K., Abel, R.S., Marmolejo-Ramos, F., Konietschke, F. (2019).
##' Nonparametric multiple comparisons. Behavior Research Methods,
##' 1-14. DOI:10.3758/s13428-019-01247-9
##' @keywords htest
##' @examples
##'
##' data(panic)
##' a<-mctp.rm(CGI~week, data=panic, type = "Dunnett",
##' alternative = "two.sided",
##' asy.method = "log.odds", plot.simci = FALSE,
##' info = FALSE, contrast.matrix = NULL)
##' summary(a)
##' plot(a)
##'
##' b<-mctp.rm(CGI~week, data=panic, type = "Dunnett",
##' alternative = "two.sided",
##' asy.method = "mult.t", plot.simci = FALSE,
##' info = FALSE, contrast.matrix = NULL)
##' summary(b)
##' plot(b)
##'
##' c<-mctp.rm(CGI~week, data=panic, type = "Dunnett",
##' alternative = "two.sided",
##' asy.method = "fisher", plot.simci = FALSE,
##' info = FALSE, contrast.matrix = NULL)
##' summary(c)
##' plot(c)
##'
##' d<-mctp.rm(CGI~week, data=panic, type = "Tukey",
##' alternative = "two.sided",
##' asy.method = "mult.t", plot.simci = TRUE)
##' summary(d)
##'
##' @export mctp.rm
mctp.rm <- function (formula, data, type = c("Tukey", "Dunnett", "Sequen",
"Williams", "Changepoint", "AVE", "McDermott", "Marcus",
"UmbrellaWilliams", "GrandMean", "UserDefined"), conf.level = 0.95,
alternative = c("two.sided", "less", "greater"), asy.method = c("log.odds",
"fisher", "mult.t", "normal"), plot.simci = FALSE, control = NULL,
info = TRUE, rounds = 3, contrast.matrix = NULL, correlation = FALSE,
const=1/1.702)
{
type <- match.arg(type)
alternative <- match.arg(alternative)
asy.method <- match.arg(asy.method)
input.list <- list(formula = formula, data = data, type = type[1],
conf.level = conf.level, alternative = alternative, asy.method = asy.method,
plot.simci = plot.simci, control = control, info = info,
rounds = rounds, contrast.matrix = contrast.matrix, correlation = correlation, const=const)
conflevel <- conf.level
if (conflevel >= 1 || conflevel <= 0) {
stop("The confidence level must be between 0 and 1!")
if (is.null(alternative)) {
stop("Please declare the alternative! (two.sided, less, greater)")
}
}
if (length(formula) != 3) {
stop("You can only analyse one-way layouts!")
}
dat <- model.frame(formula, data)
if (ncol(dat) != 2) {
stop("Specify one response and only one class variable in the formula")
}
if (is.numeric(dat[, 1]) == FALSE) {
stop("Response variable must be numeric")
}
response <- dat[, 1]
factorx <- as.factor(dat[, 2])
fl <- levels(factorx)
a <- nlevels(factorx)
if (a <= 2) {
stop("You want to perform a two-sample test. Please use the function npar.t.test.paired")
}
samples <- split(response, factorx)
n <- as.numeric(sapply(samples, length))
if (any(n == 1)) {
warn <- paste("The factor level", fl[n == 1], "has got only one observation!")
stop(warn)
}
a <- length(n)
if (type == "UserDefined") {
if (is.null(contrast.matrix)) {
stop("Please enter a contrast matrix!")
}
ch <- contrast.matrix
rownames(ch) <- paste("C", 1:nrow(ch))
colnames(ch) <- fl
}
if (type != "UserDefined") {
if (is.null(control)) {
icon <- 1
}
if (!is.null(control)) {
icon <- which(fl == control)
}
ch <- contrMat(n = n, type, base = icon)
}
nc <- nrow(ch)
connames <- rownames(ch)
Con <- matrix(ch, ncol = a)
rownames(Con) <- connames
colnames(Con) <- colnames(ch)
#Standardizing contrasts
if(asy.method == "log.odds") {
Con <- 2*Con/rowSums(abs(Con))}
tmp1 <- sort(rep(1:a, a))
tmp2 <- rep(1:a, a)
pairRanks <- lapply(1:(a^2), function(arg) rank(c(samples[[tmp1[arg]]],
samples[[tmp2[arg]]])))
intRanks <- lapply(1:a, function(arg) rank(samples[[arg]]))
n <- n[[1]]
plis <- lapply(1:(a^2), function(arg) 1/n * (mean(pairRanks[[arg]][(n +
1):(2 * n)]) - (n + 1)/2))
vec.plis <- as.numeric(as.character(plis))
pd <- c()
for (i in 1:a) {
pd[i] <- 1/a * sum((tmp2 == i) * vec.plis)
}
pd1 <- (pd == 1)
pd0 <- (pd == 0)
pd[pd1] <- 0.999
pd[pd0] <- 0.001
Zlong <- c()
for (i in 1:(a^2)) {
Zlong <- cbind(Zlong, 1/n * (pairRanks[[i]][1:n] - intRanks[[tmp1[i]]] -
pairRanks[[i]][(n + 1):(2 * n)] + intRanks[[tmp2[i]]]))
}
Zquer <- 1/n * colSums(Zlong)
Yd <- matrix(rep(0, a^4), nrow = a^2, ncol = a^2)
for (k in 1:n) {
Yd <- Yd + 1/(n - 1) * (t(Zlong[k, ] - t(Zquer)) %*%
(Zlong[k, ] - t(Zquer)))
}
W <- matrix(rep(1/a * diag(a), a), nrow = a)
sYd <- svd(Yd)
dYd <- diag(sqrt(sYd$d))
hYd <- sYd$u %*% dYd %*% t(sYd$v)
hVd <- W %*% hYd
Vd <- hVd %*% t(hVd)
logit.pd <- log(c(pd/(1 - pd)))
logit.pd.dev <- diag(1/c((pd * (1 - pd))))
hlVd <- logit.pd.dev %*% hVd
lVd <- hlVd %*% t(hlVd)
Lower.logit1 <- logit.pd - qnorm(conf.level)/sqrt(n) * sqrt(c(diag(lVd)))
Upper.logit1 <- logit.pd + qnorm(conf.level)/sqrt(n) * sqrt(c(diag(lVd)))
Lower.logit <- exp(Lower.logit1)/(1 + exp(Lower.logit1))
Upper.logit <- exp(Upper.logit1)/(1 + exp(Upper.logit1))
corr.mat <- function(m, nc) {
rho <- matrix(c(0), ncol = nc, nrow = nc)
for (i in 1:nc) {
for (j in 1:nc) {
rho[i, j] <- m[i, j]/sqrt(m[i, i] * m[j, j])
}
}
return(rho)
}
dfT <- n - 1
ind.plus <- matrix(as.integer(Con > 0), ncol=ncol(Con),nrow=nrow(Con))
ind.minus <- matrix(as.integer(Con < 0), ncol=ncol(Con),nrow=nrow(Con))
Con.plus <- Con*ind.plus
Con.minus<- -Con*ind.minus
Cpd.plus <- Con.plus %*% pd
Cpd.minus <- Con.minus %*% pd
Cpd <- Con %*% pd
hCV <- Con %*% hVd
CV <- hCV %*% t(hCV)
rhobf <- corr.mat(CV, nc)
p.adj <- c()
switch(asy.method, mult.t = {
Tstat <- sqrt(n) * (Cpd)/sqrt(c(diag(CV)))
AsyMethod <- paste("Multi - T with", round(dfT, rounds),
"DF")
switch(alternative, two.sided = {
text.Output <- paste("True differences of relative effects are not equal to 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvt(lower = -abs(Tstat[pp]), upper = abs(Tstat[pp]),
delta = rep(0, nc), df = dfT, corr = rhobf)[1]
}
crit <- qmvt(conflevel, corr = rhobf, tail = "both",
df = dfT)$quantile
Lower <- Cpd - crit/sqrt(n) * sqrt(c(diag(CV)))
Upper <- Cpd + crit/sqrt(n) * sqrt(c(diag(CV)))
}, less = {
text.Output <- paste("True differences of relative effects are less than 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvt(lower = Tstat[pp], upper = Inf,
df = dfT, delta = rep(0, nc), corr = rhobf)
}
crit <- qmvt(conflevel, df = dfT, corr = rhobf, tail = "lower")$quantile
Lower <- rep(-1, nc)
Upper <- Cpd + crit/sqrt(n) * sqrt(c(diag(CV)))
}, greater = {
text.Output <- paste("True differences of relative effects are greater than 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvt(lower = -Inf, upper = Tstat[pp],
df = dfT, delta = rep(0, nc), corr = rhobf)[1]
}
crit <- qmvt(conflevel, corr = rhobf, df = dfT, tail = "lower")$quantile
Lower <- Cpd - crit/sqrt(n) * sqrt(c(diag(CV)))
Upper <- rep(1, nc)
})
}, normal = {
AsyMethod <- "Normal - Approximation"
Tstat <- sqrt(n) * (Cpd)/sqrt(c(diag(CV)))
switch(alternative, two.sided = {
text.Output <- paste("True differences of relative effects are not equal to 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvnorm(lower = -abs(Tstat[pp]),
upper = abs(Tstat[pp]), mean = rep(0, nc), corr = rhobf)[1]
}
crit <- qmvnorm(conflevel, corr = rhobf, tail = "both")$quantile
Lower <- Cpd - crit/sqrt(n) * sqrt(c(diag(CV)))
Upper <- Cpd + crit/sqrt(n) * sqrt(c(diag(CV)))
}, less = {
text.Output <- paste("True differences of relative effects are less than 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvnorm(lower = Tstat[pp], upper = Inf,
mean = rep(0, nc), corr = rhobf)
}
crit <- qmvnorm(conflevel, corr = rhobf, tail = "lower")$quantile
Lower <- rep(-1, nc)
Upper <- Cpd + crit/sqrt(n) * sqrt(c(diag(CV)))
}, greater = {
text.Output <- paste("True differences of relative effects are greater than 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvnorm(lower = -Inf, upper = Tstat[pp],
mean = rep(0, nc), corr = rhobf)
}
crit <- qmvnorm(conflevel, corr = rhobf, tail = "lower")$quantile
Lower <- Cpd - crit/sqrt(n) * sqrt(c(diag(CV)))
Upper <- rep(1, nc)
})
}, fisher = {
AsyMethod <- paste("Fisher with", round(dfT, rounds),
"DF")
Cfisher <- 1/2 * log((1 + Cpd)/(1 - Cpd))
Vfisherdev <- diag(c(1/(1 - Cpd^2)))
hVfisher <- Vfisherdev %*% hCV
Vfisher <- hVfisher %*% t(hVfisher)
Tstat <- sqrt(n) * Cfisher/sqrt(c(diag(Vfisher)))
switch(alternative, two.sided = {
text.Output <- paste("True differences of relative effects are not equal to 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvt(lower = -abs(Tstat[pp]), upper = abs(Tstat[pp]),
delta = rep(0, nc), corr = rhobf, df = dfT)[1]
}
crit <- qmvt(conflevel, corr = rhobf, tail = "both",
df = dfT)$quantile
Lower1 <- Cfisher - crit/sqrt(n) * sqrt(c(diag(Vfisher)))
Upper1 <- Cfisher + crit/sqrt(n) * sqrt(c(diag(Vfisher)))
Lower <- (exp(2 * Lower1) - 1)/(exp(2 * Lower1) +
1)
Upper <- (exp(2 * Upper1) - 1)/(exp(2 * Upper1) +
1)
}, less = {
text.Output <- paste("True differences of relative effects are less than 0")
for (pp in 1:nc) {
p.adj[pp] <- pmvt(lower = Tstat[pp], upper = Inf,
delta = rep(0, nc), df = dfT, corr = rhobf)
}
crit <- qmvt(conflevel, corr = rhobf, tail = "lower",
df = dfT)$quantile
Lower <- rep(-1, nc)
Upper1 <- Cfisher + crit/sqrt(n) * sqrt(c(diag(Vfisher)))
Upper <- (exp(2 * Upper1) - 1)/(exp(2 * Upper1) +
1)
}, greater = {
text.Output <- paste("True differences of relative effects are greater than 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvt(lower = -Inf, upper = Tstat[pp],
delta = rep(0, nc), corr = rhobf, df = dfT)
}
crit <- qmvnorm(conflevel, corr = rhobf, tail = "lower")$quantile
Lower1 <- Cfisher - crit/sqrt(n) * sqrt(c(diag(Vfisher)))
Lower <- (exp(2 * Lower1) - 1)/(exp(2 * Lower1) +
1)
Upper <- rep(1, nc)
})
}, log.odds = {
AsyMethod <- paste("Adjusted log odds with", round(dfT, rounds), "DF")
Clogodds <- const * (log(Cpd.plus/(1-Cpd.plus)) - log(Cpd.minus/(1-Cpd.minus)))
glogodds.plus <- Cpd.plus - Cpd.plus^2
glogodds.minus <- Cpd.minus - Cpd.minus^2
glogodds.plus.mat <- matrix(1,nrow=1,ncol=a) %x% matrix(glogodds.plus, ncol=1)
glogodds.minus.mat <- matrix(1,nrow=1,ncol=a) %x% matrix(glogodds.minus, ncol=1)
Glogodds <- const * (Con.plus/glogodds.plus.mat - Con.minus/glogodds.minus.mat)
hVlogodds <- Glogodds %*% hVd
Vlogodds <- hVlogodds %*% t(hVlogodds)
Tstat <- sqrt(n) * Clogodds/sqrt(c(diag(Vlogodds)))
switch(alternative, two.sided = {
text.Output <- paste("True differences of relative effects are not equal to 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvt(lower = -abs(Tstat[pp]), upper = abs(Tstat[pp]),
delta = rep(0, nc), corr = rhobf, df = dfT)[1]
}
crit <- qmvt(conflevel, corr = rhobf, tail = "both",
df = dfT)$quantile
Lower <- Clogodds - crit/sqrt(n) * sqrt(c(diag(Vlogodds)))
Upper <- Clogodds + crit/sqrt(n) * sqrt(c(diag(Vlogodds)))
LowerPlot <- Lower
UpperPlot <- Upper
}, less = {
text.Output <- paste("True differences of relative effects are less than 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvt(lower = Tstat[pp], upper = Inf,
delta = rep(0, nc), df = dfT, corr = rhobf)
}
crit <- qmvt(conflevel, corr = rhobf, tail = "lower",
df = dfT)$quantile
Lower <- rep(-Inf, nc)
LowerPlot <- Clogodds - crit/sqrt(n) * sqrt(c(diag(Vlogodds)))
Upper <- Clogodds + crit/sqrt(n) * sqrt(c(diag(Vlogodds)))
UpperPlot <- Upper
}, greater = {
text.Output <- paste("True differences of relative effects are greater than 0")
for (pp in 1:nc) {
p.adj[pp] <- 1 - pmvt(lower = -Inf, upper = Tstat[pp],
delta = rep(0, nc), corr = rhobf, df = dfT)
}
crit <- qmvt(conflevel, corr = rhobf, tail = "lower")$quantile
Lower <- Clogodds - crit/sqrt(n) * sqrt(c(diag(Vlogodds)))
LowerPlot <- Lower
Upper <- rep(Inf, nc)
UpperPlot <- Clogodds + crit/sqrt(n) * sqrt(c(diag(Vlogodds)))
})
})
if (plot.simci == TRUE) {
text.Ci <- paste(conflevel * 100, "%", "Simultaneous Confidence Intervals")
Lowerp <- "|"
if(asy.method!="log.odds"){
plot(Cpd, 1:nc, xlim = c(-1, 1), pch = 15, axes = FALSE,
xlab = "", ylab = "")
axis(1, at = seq(-1, 1, 0.1))
points(Lower, 1:nc, pch = Lowerp, font = 2, cex = 2)
points(Upper, 1:nc, pch = Lowerp, font = 2, cex = 2)
abline(v = 0, lty = 3, lwd = 2)
for (ss in 1:nc) {
polygon(x = c(Lower[ss], Upper[ss]), y = c(ss, ss), lwd = 2)
}
}
else{ #log.odds
plot(Clogodds, 1:nc, xlim = c(floor(min(LowerPlot)), ceiling(max(UpperPlot))), pch = 15, axes = FALSE,
xlab = "", ylab = "")
axis(1, at = seq(floor(min(LowerPlot)), ceiling(max(UpperPlot)), 0.1*(ceiling(max(UpperPlot))-floor(min(LowerPlot)))))
hugenumber<-10000000
if(alternative=="two.sided"){
points(LowerPlot, 1:nc, pch = Lowerp, font = 2, cex = 2)
points(UpperPlot, 1:nc, pch = Lowerp, font = 2, cex = 2)
for (ss in 1:nc) {
polygon(x = c(LowerPlot[ss], UpperPlot[ss]), y = c(ss, ss), lwd = 2)
}
}
else if(alternative=="less"){
points(UpperPlot, 1:nc, pch = Lowerp, font = 2, cex = 2)
for (ss in 1:nc) {
polygon(x = c(-hugenumber, UpperPlot[ss]), y = c(ss, ss), lwd = 2)
}
}
else{ #greater
points(LowerPlot, 1:nc, pch = Lowerp, font = 2, cex = 2)
for (ss in 1:nc) {
polygon(x = c(LowerPlot[ss], hugenumber), y = c(ss, ss), lwd = 2)
}
}
abline(v = 0, lty = 3, lwd = 2)
}
axis(2, at = 1:nc, labels = connames)
box()
title(main = c(text.Ci, paste("Type of Contrast:", type),
paste("Method:", AsyMethod)))
}
data.info <- data.frame(Sample = fl, Size = n, Effect = pd,
Lower = Lower.logit, Upper = Upper.logit)
if(asy.method!="log.odds"){
est<-Cpd
}
else{
est<-Clogodds
}
Analysis.of.Relative.Effects <- data.frame(Estimator = round(est,
rounds), Lower = round(Lower, rounds), Upper = round(Upper,
rounds), Statistic = round(Tstat, rounds), p.Value = p.adj)
Analysis.Inf <- data.frame(Estimator = est, Lower = Lower, Upper = Upper,
Statistic = Tstat, p.Value = p.adj)
Overall <- data.frame(Quantile = crit, p.Value = min(p.adj))
result <- list(Data.Info = data.info, Contrast = Con, Analysis = Analysis.of.Relative.Effects,
Analysis.Inf = Analysis.Inf, Overall = Overall)
if (info == TRUE) {
cat("\n", "#----------------Nonparametric Multiple Comparisons for relative effects---------------#",
"\n", "\n", "-", "Alternative Hypothesis: ", text.Output,
"\n", "-", "Estimation Method: Global Pseudo ranks",
"\n", "-", "Type of Contrast", ":", type, "\n", "-",
"Confidence Level:", conflevel * 100, "%", "\n",
"-", "Method", "=", AsyMethod, "\n", "\n", "#--------------------------------------------------------------------------------------#",
"\n", "\n")
}
if (correlation == TRUE) {
result$Covariance <- CV
result$Correlation <- rhobf
}
result$input <- input.list
result$text.Output <- text.Output
result$connames <- connames
result$AsyMethod <- AsyMethod
class(result) <- "mctp.rm"
return(result)
}
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.