SteelConfInt {kSamples}R Documentation

Simultaneous Confidence Bounds Based on Steel's Multiple Comparison Wilcoxon Tests

Description

This function inverts pairwise Wilcoxon tests, comparing a common control sample with each of several treatment samples to provide simultaneous confidence bounds for the respective shift parameters by which the sampled treatment populations may differ from the control population. It is assumed that all samples are independent and that the sampled distributions are continuous to avoid ties. The joint coverage probability for all bounds/intervals is calculated, estimated, or approximated, see Details. For treatment of ties also see Details.

Usage

SteelConfInt(..., data = NULL, conf.level = 0.95, 
	alternative = c("less", "greater", "two.sided"), 
     	method = c("asymptotic", "exact", "simulated"), Nsim = 10000)

Arguments

...

Either several sample vectors, say x1,,xkx_1, \ldots, x_k, with xix_i containing nin_i sample values. ni>4n_i > 4 is recommended for reasonable asymptotic PP-value calculation. The pooled sample size is denoted by N=n1++nkN=n_1+\ldots+n_k. The first vector serves as control sample, the others as treatment samples.

or a list of such sample vectors.

or a formula y ~ g, where y contains the pooled sample values and g (same length as y) is a factor with levels identifying the samples to which the elements of y belong. The lowest factor level corresponds to the control sample, the other levels to treatment samples.

data

= an optional data frame providing the variables in formula y ~ g.

conf.level

= 0.95 (default) the target joint confidence level for all bounds/intervals.

0 < conf.level < 1.

alternative

= c("less", "greater", "two.sided"), where "less" results in simultaneous upper confidence bounds for all shift parameters and any negative upper bound should lead to the rejection of the null hypothesis of all shift parameters being zero or positive in favor of at least one being less than zero.

"greater" results in simultaneous lower confidence bounds for all shift parameters and any positive lower bound should lead to the rejection of the null hypothesis of all shift parameters being zero or negative in favor of at least one being greater than zero.

"two.sided" results in simultaneous confidence intervals for all shift parameters and any interval not straddling zero should lead to the rejection of the null hypothesis of all shift parameters being zero in favor of at least one being different from zero.

method

= c("asymptotic", "exact", "simulated"), where

"asymptotic" uses only an asymptotic normal approximation to approximate the achieved joint coverage probability. This calculation is always done.

"exact" uses full enumeration of all sample splits to obtain the exact achieved joint coverage probability (see Details). It is used only when Nsim is at least as large as the number of full enumerations. Otherwise, method reverts to "simulated" using the given Nsim.

"simulated" uses Nsim simulated random splits of the pooled samples into samples of sizes n1,,nkn_1, \ldots, n_k, to estimate the achieved joint coverage probability.

Nsim

= 10000 (default), number of simulated sample splits to use. It is only used when method = "simulated", or when method = "exact" reverts to method = "simulated", as previously explained.

Details

The first sample is treated as control sample with sample size n1n_1. The remaining s=k1s=k-1 samples are treatment samples. Let W1i,i=2,,kW_{1i}, i=2,\ldots,k denote the respective Wilcoxon statistics comparing the common control sample (index 1) with each of the ss treatment samples (indexed by ii). For each comparison of control and treatment ii sample only the observations of the two samples involved are ranked. By Wi=W1ini(ni+1)/2W_i=W_{1i}-n_i(n_i+1)/2 we denote the corresponding Mann-Whitney test statistic. Furthermore, let Di(j)D_{i(j)} denote the jj-th ordered value (ascending order) of the n1nin_1n_i paired differences between the observations in treatment sample ii and those of the control sample. By simple extension of results in Lehmann (2006), pages 87 and 92, the following equations hold, relating the null distribution of the Mann-Whitney statistics and the joint coverage probabilities of the Di(ji)D_{i(j_i)} for any set of j1,,jsj_1,\ldots,j_s with 1jin1ni1\le j_i \le n_1 n_i.

PΔ(ΔiDi(ji),i=2,,k)=P0(Wiji1,i=2,,k)P_\Delta(\Delta_i \le D_{i(j_i)}, i=2,\ldots,k)=P_0(W_i\le j_i -1, i=2,\ldots,k)

and

PΔ(ΔiDi(ji),i=2,,s)=P0(Win1niji,i=2,,k)P_\Delta(\Delta_i \ge D_{i(j_i)}, i=2,\ldots,s)=P_0(W_{i}\le n_1 n_i -j_i, i=2,\ldots,k)

where PΔP_\Delta refers to the distribution under Δ=(Δ2,,Δk)\Delta=(\Delta_2,\ldots,\Delta_k) and P0P_0 refers to the joint null distribution of the WiW_i when all sampled distributions are the same and continuous. There are k1k-1 indices jij_i that can be manipulated to affect the achieved confidence level. To limit the computational complexity standardized versions of the WiW_i, i.e., (Wiμi)/τi(W_i-\mu_i)/\tau_i with μi\mu_i and τi\tau_i representing mean and standard deviation of WiW_i, are used to choose a common value for (ji1μi)/τi(j_i -1-\mu_i)/\tau_i (satisfying the γ\gamma level) from the multivariate normal approximation for the WiW_i (see Miller (1981) and Scholz (2016)), and reduce that to integer values for jij_i, rounding up, rounding down, and rounding to the nearest integer. These integers jij_i are then used in approximating the actual joint probabilities P0(Wiji1,i=2,,k)P_0(W_i\le j_i -1, i=2,\ldots,k), and from these three coverage probabilities the one that is closest to the nominal confidence level γ\gamma and γ\ge \gamma and also also the one that is closest without the restriction γ\ge \gamma are chosen.

When method = "exact" or = "simulated" is specified, the same process is used, using either the fully enumerated exact distribution of Wi,i=2,,kW_i, i=2,\ldots,k (based on a recursive version of Chase's sequence as presented in Knuth (2011)) for all sample splits, or the simulated distribution of Wi,i=2,,kW_i, i=2,\ldots,k. However, since these distributions are discrete the starting point before rounding up is the smallest quantile such that the proportion of distribution values less or equal to it is at least γ\gamma. The starting point before rounding down is the highest quantile such that the proportion of distribution values less or equal to it is at most γ\gamma. The third option of rounding to the closest integer is performed using the average of the first two.

Confidence intervals are constructed by using upper and lower confidence bounds, each with same confidence level of (1+γ)/2(1+\gamma)/2.

When the original sample data appear to be rounded, and especially when there are ties, one should widen the computed intervals or bounds by the rounding ϵ\epsilon, as illustrated in Lehmann (2006), pages 85 and 94. For example, when all sample values appear to end in one of .0,.2,.4,.6,.8.0, .2, .4, .6, .8, the rounding ϵ\epsilon would be .2.2. Ultimately, this is a judgment call for the user. Such widening of intervals will make the actually achieved confidence level \ge the stated achieved level.

Value

A list of class kSamples with components

test.name

"Steel.bounds"

n1

the control sample size =n1= n_1

ns

vector (n2,,nk)(n_2,\ldots,n_k) of the s=k1s=k-1 treatment sample sizes

N

size of the pooled sample =n1++nk= n_1+\ldots+n_k

n.ties

number of ties in the pooled sample

bounds

a list of data frames. When method = "asymptotic" is specified, the list consists of two data frames named conservative.bounds.asymptotic and closest.bounds.asymptotic. Each data frame consists of ss rows corresponding to the ss shift parameters Δi\Delta_i and three columns, the first column giving the lower bound, the second column the upper bound, while the first row of the third column states the computed confidence level by asymptotic approximation, applying jointly to all ss sets of bounds. For one-sided bounds the corresponding other bound is set to Inf or -Inf, respectively.

In case of conservative.bounds.asymptotic the achieved asymptotic confidence level is targeted to be \ge conf.level, but closest to it among three possible choices (see Details).

In the case of closest.bounds.asymptotic the achieved level is targeted to be closest to conf.level.

When method = "exact" or method = "simulated" is specified the list consists in addition of two further data frames named either

conservative.bounds.exact and closest.bounds.exact or

conservative.bounds.simulated and closest.bounds.simulated.

In either case the structure and meaning of these data frames parallels that of the "asymptotic" case.

method

the method used.

Nsim

the number of simulations used.

j.LU

an ss by 4 matrix giving the indices jij_i used for computing the bounds Di(ji)D_{i(j_i)} for Δi,i=1,,s\Delta_i, i=1,\ldots, s.

warning

method = "exact" should only be used with caution. Computation time is proportional to the number of enumerations. Experiment with system.time and trial values for Nsim to get a sense of the required computing time.

References

Knuth, D.E. (2011), The Art of Computer Programming, Volume 4A Combinatorial Algorithms Part 1, Addison-Wesley

Lehmann, E.L. (2006), Nonparametrics, Statistical Methods Based on Ranks, Revised First Edition, Springer Verlag.

Miller, Rupert G., Jr. (1981), Simultaneous Statistical Inference, Second Edition, Springer Verlag, New York.

Scholz, F.W. (2023), "On Steel's Test with Ties", https://arxiv.org/abs/2308.05873

Examples

z1 <- c(103, 111, 136, 106, 122, 114)
z2 <- c(119, 100,  97,  89, 112,  86)
z3 <- c( 89, 132,  86, 114, 114, 125)
z4 <- c( 92, 114,  86, 119, 131,  94)
set.seed(2627)
SteelConfInt(list(z1,z2,z3,z4),conf.level=0.95,alternative="two.sided", 
   method="simulated",Nsim=10000)
# or with same seed
# SteelConfInt(z1,z2,z3,z4,conf.level=0.95,alternative="two.sided", 
#   method="simulated",Nsim=10000)

[Package kSamples version 1.2-10 Index]