R/ad.pval.R

In kSamples: K-Sample Rank Tests and their Combinations

ad.pval<-
function (tx,m,version=1) 
{
# This function "adk.pval" evaluates the p-value of the observed value 
# tx of the T_m statistic in "K-Sample Anderson-Darling Tests" by F.W. Scholz 
# and M.A. Stephens (1987), Journal of the American Statistical Association, 
# Vol 82, No. 399, pp 918-924. Thus this p-value is P(T_m >= tx).
#
# Input: tx = observed value of T_m, tx > 0.
#         m = the index of T_m, m >= 1.
#        version = 1 (default) uses the first version of the AD statistic,
#                    otherwise the second version is used.
#
# Output: a list with components
#         p0 = p-value of tx, i.e., p0 = P(T_m >= tx)
#         extrap = a logical indicator
#                    extrap = TRUE means that linear extrapolation took place
#                    extrap = FALSE means that quadratic interpolation was used.
#
# Computational Details:
#
# This function uses the upper T_m quantiles as obtained via simulation of
# the Anderson-Darling test statistics (Nsim = 2*10^6) with sample sizes n=500
# for each sample, and after standardization, in order to emulate the Table 1 
# values given in the above reference. However, here we estimate p-quantiles
# for p = .00001,.00005,.0001,.0005,.001,.005,.01,.025,.05,.075,
#	.1,.2,.3,.4,.5,.6,.7,.8,.9,.925,.95,.975,.99,.9925,.995,.9975,.999,
# .99925,.9995,.99975,.9999,.999925,.99995,.999975,.99999
# First the appropriate p-quantiles are determined from those simulated
# for ms = 1,2,3,4,6,8,10, Inf, interpolating to the given value of m
# by fitting a smoothing spline in 1/sqrt(ms) to the simulated quantiles.
# Visual inspection of plots of fitted splines versus 1/sqrt(ms) shows
# good agreement for the used smoothing parameter spar = .4.
#
# Next a smoothing spline is used to fit the log((1-p)/p) to
# these interpolated quantiles and the value fitted to tx is
# obtained (extrapolating linearly) beyond p = .00001 and .99999.
#
# The p-values from Table 1 were reproduced mainly with relative error
# bounded by 1%, in 6 cases with respective relative errors of 1.25%, 1.2%,
# 1.3%, 1.5%, 2.7% and 2.8% as can be seen from the relative error table below.
# The columns correspond to Table 1 p-values of .25, .10, .05, .025, .01 and 
# the rows correspond to m = 1, 2 3, 4, 6, 8, 10, Inf.
# 
#         [,1]    [,2]    [,3]    [,4]    [,5]
# [1,] -0.0012  0.0011  0.0038  0.0123  0.0277
# [2,] -0.0105 -0.0109 -0.0022  0.0129  0.0271
# [3,] -0.0077 -0.0052  0.0042  0.0106  0.0152
# [4,] -0.0050 -0.0031  0.0035  0.0067  0.0057
# [5,] -0.0041 -0.0033  0.0055  0.0014 -0.0050
# [6,] -0.0048 -0.0026  0.0002 -0.0021 -0.0082
# [7,] -0.0070 -0.0042 -0.0023 -0.0036 -0.0119
# [8,]  0.0010  0.0005  0.0013  0.0017  0.0023
#
#
#
#  
# Fritz Scholz, March 2012
#=========================================================================================
if(version==1){
table1.adk <- structure(c(1, 2, 3, 4, 6, 8, 10, Inf, -1.1954, -1.5806, -1.8172, 
-2.0032, -2.2526, -2.4204, -2.5283, -4.2649, -1.1786, -1.5394, 
-1.7728, -1.9426, -2.1685, -2.3288, -2.4374, -3.8906, -1.166, 
-1.5193, -1.7462, -1.9067, -2.126, -2.2818, -2.3926, -3.719, 
-1.1407, -1.4659, -1.671, -1.8105, -2.0048, -2.1356, -2.2348, 
-3.2905, -1.1253, -1.4371, -1.6314, -1.7619, -1.9396, -2.0637, 
-2.1521, -3.0902, -1.0777, -1.3503, -1.5102, -1.6177, -1.761, 
-1.8537, -1.9178, -2.5758, -1.0489, -1.2984, -1.4415, -1.5355, 
-1.6625, -1.738, -1.7936, -2.3263, -0.9978, -1.2098, -1.3251, 
-1.4007, -1.4977, -1.5555, -1.5941, -1.96, -0.9417, -1.1187, 
-1.209, -1.2671, -1.3382, -1.379, -1.405, -1.6449, -0.8981, -1.0491, 
-1.1235, -1.1692, -1.2249, -1.2552, -1.2755, -1.4395, -0.8598, 
-0.9904, -1.0513, -1.0879, -1.1317, -1.155, -1.1694, -1.2816, 
-0.7258, -0.7938, -0.8188, -0.8312, -0.8435, -0.8471, -0.8496, 
-0.8416, -0.5966, -0.617, -0.6177, -0.6139, -0.6073, -0.5987, 
-0.5941, -0.5244, -0.4572, -0.4383, -0.419, -0.4033, -0.3834, 
-0.3676, -0.3587, -0.2533, -0.2966, -0.2428, -0.2078, -0.1844, 
-0.1548, -0.1346, -0.1224, 0, -0.1009, -0.0169, 0.0304, 0.0596, 
0.0933, 0.1156, 0.1294, 0.2533, 0.1571, 0.2635, 0.3169, 0.348, 
0.3823, 0.4038, 0.4166, 0.5244, 0.5357, 0.6496, 0.6992, 0.7246, 
0.7528, 0.7683, 0.7771, 0.8416, 1.2255, 1.2989, 1.3202, 1.3254, 
1.3305, 1.3286, 1.3257, 1.2816, 1.5262, 1.5677, 1.5709, 1.5663, 
1.5561, 1.5449, 1.5356, 1.4395, 1.9633, 1.943, 1.919, 1.8975, 
1.8641, 1.8389, 1.8212, 1.6449, 2.7314, 2.5899, 2.5, 2.4451, 
2.3664, 2.3155, 2.2823, 1.96, 3.7825, 3.4425, 3.2582, 3.1423, 
3.0036, 2.9101, 2.8579, 2.3263, 4.1241, 3.716, 3.4984, 3.3651, 
3.2003, 3.0928, 3.0311, 2.4324, 4.6044, 4.0847, 3.8348, 3.6714, 
3.4721, 3.3453, 3.2777, 2.5758, 5.409, 4.7223, 4.4022, 4.1791, 
3.9357, 3.7809, 3.6963, 2.807, 6.4954, 5.5823, 5.1456, 4.8657, 
4.5506, 4.3275, 4.2228, 3.0902, 6.8279, 5.8282, 5.3658, 5.0749, 
4.7318, 4.4923, 4.3642, 3.1747, 7.2755, 6.197, 5.6715, 5.3642, 
4.9991, 4.7135, 4.5945, 3.2905, 8.1885, 6.8537, 6.2077, 5.8499, 
5.4246, 5.1137, 4.9555, 3.4808, 9.3061, 7.6592, 6.85, 6.4806, 
5.9919, 5.6122, 5.5136, 3.719, 9.6132, 7.9234, 7.1025, 6.6731, 
6.1549, 5.8217, 5.7345, 3.7911, 10.0989, 8.2395, 7.4326, 6.9567, 
6.3908, 6.011, 5.9566, 3.8906, 10.8825, 8.8994, 7.8934, 7.4501, 
6.9009, 6.4538, 6.2705, 4.0556, 11.8537, 9.5482, 8.5568, 8.0283, 
7.4418, 6.9524, 6.6195, 4.2649), .Dim = c(8L, 36L), .Dimnames = list(
    NULL, NULL))}else{
table1.adk <- structure(c(1, 2, 3, 4, 6, 8, 10, Inf, -1.1976, -1.5824, -1.8195, 
-2.005, -2.2546, -2.422, -2.5292, -4.2649, -1.1806, -1.5416, 
-1.7747, -1.9434, -2.1687, -2.3301, -2.438, -3.8906, -1.1681, 
-1.5212, -1.7479, -1.9078, -2.1268, -2.2827, -2.3937, -3.719, 
-1.1427, -1.4677, -1.6724, -1.8115, -2.0059, -2.1363, -2.2359, 
-3.2905, -1.1272, -1.4387, -1.6325, -1.7629, -1.9405, -2.0649, 
-2.1527, -3.0902, -1.0794, -1.3518, -1.5112, -1.6187, -1.7617, 
-1.8545, -1.9182, -2.5758, -1.0504, -1.2997, -1.4425, -1.5362, 
-1.6632, -1.7387, -1.7943, -2.3263, -0.999, -1.2109, -1.3259, 
-1.4014, -1.4981, -1.5561, -1.5945, -1.96, -0.9428, -1.1196, 
-1.2098, -1.2677, -1.3386, -1.3795, -1.4054, -1.6449, -0.8991, 
-1.05, -1.1241, -1.1697, -1.2253, -1.2557, -1.2758, -1.4395, 
-0.8607, -0.9911, -1.0518, -1.0883, -1.1321, -1.1555, -1.1698, 
-1.2816, -0.7264, -0.7944, -0.8192, -0.8315, -0.8437, -0.8473, 
-0.8498, -0.8416, -0.597, -0.6173, -0.6179, -0.6141, -0.6074, 
-0.5989, -0.5942, -0.5244, -0.4574, -0.4385, -0.4191, -0.4034, 
-0.3835, -0.3677, -0.3588, -0.2533, -0.2966, -0.2427, -0.2078, 
-0.1844, -0.1548, -0.1347, -0.1224, 0, -0.1007, -0.0168, 0.0305, 
0.0596, 0.0934, 0.1157, 0.1295, 0.2533, 0.1573, 0.2638, 0.3171, 
0.3482, 0.3825, 0.404, 0.4168, 0.5244, 0.5363, 0.6501, 0.6996, 
0.7249, 0.753, 0.7685, 0.7773, 0.8416, 1.2263, 1.2997, 1.3209, 
1.3258, 1.3309, 1.329, 1.326, 1.2816, 1.5274, 1.5686, 1.5716, 
1.5667, 1.5565, 1.5453, 1.536, 1.4395, 1.9644, 1.944, 1.92, 1.8983, 
1.8647, 1.8396, 1.8216, 1.6449, 2.7334, 2.5915, 2.5013, 2.4457, 
2.3671, 2.3162, 2.2827, 1.96, 3.7851, 3.4443, 3.2595, 3.1436, 
3.0046, 2.9111, 2.8585, 2.3263, 4.1255, 3.7175, 3.4997, 3.3661, 
3.2011, 3.0939, 3.0318, 2.4324, 4.6067, 4.0869, 3.8363, 3.6724, 
3.4729, 3.3463, 3.278, 2.5758, 5.4121, 4.7248, 4.4032, 4.1812, 
3.9369, 3.7819, 3.6977, 2.807, 6.5, 5.5856, 5.1469, 4.8683, 4.552, 
4.3284, 4.2229, 3.0902, 6.8324, 5.8302, 5.3678, 5.0769, 4.7332, 
4.4933, 4.3654, 3.1747, 7.278, 6.1999, 5.674, 5.3661, 5.0001, 
4.7147, 4.5956, 3.2905, 8.1926, 6.8586, 6.2082, 5.8524, 5.4265, 
5.115, 4.9571, 3.4808, 9.3096, 7.6673, 6.8522, 6.4825, 5.9934, 
5.6135, 5.5147, 3.719, 9.6207, 7.929, 7.1051, 6.6763, 6.1548, 
5.8229, 5.7358, 3.7911, 10.1076, 8.2437, 7.4349, 6.9593, 6.3923, 
6.0136, 5.9573, 3.8906, 10.8874, 8.9034, 7.8991, 7.4543, 6.9017, 
6.4568, 6.2723, 4.0556, 11.8602, 9.5499, 8.5596, 8.0315, 7.4425, 
6.9537, 6.6213, 4.2649), .Dim = c(8L, 36L), .Dimnames = list(
    NULL, NULL))
}
    extrap <- FALSE
    mt <- table1.adk[, 1]
    sqm1 <- 1/sqrt(mt)
    sqm2 <- sqm1^2
    tm <- NULL
  	p <- 1 - c(.00001,.00005,.0001,.0005,.001,.005,.01,.025,.05,.075,
	.1,.2,.3,.4,.5,.6,.7,.8,.9,.925,.95,.975,.99,.9925,.995,.9975,.999,
	.99925,.9995,.99975,.9999,.999925,.99995,.999975,.99999)
	np <- length(p)
    lp <- log(p/(1 - p))
	sqm0 <- 1/sqrt(m)
    for (i in 1:np) {
    	out <- smooth.spline(sqm1,table1.adk[,i+1],spar=.4)
    	y <- predict(out,sqm0)$y
        tm <- c(tm, y)
    }
	out <- smooth.spline(tm,lp,spar=.25)
    lp0 <- predict(out,tx)$y


    p0 <- exp(lp0)/(1 + exp(lp0))
    names(p0) <- NULL
	p0
}