Couldn't find it anywhere, so for future users who stumble on this thread,
here is some code.
Note: If your data has missing values, delete those observations from the
data first before running this code.
Place these functions in the global environment (ie run the code below).
Then, to obtain the test statistic: H_bisym_stat(x,y) (x and y are the
two vectors you are comparing).
To obtain a p-value: H_bisym_stat(x,y). The p-value is calculated via a
Monte Carlo algorithm.
You can use the permutations argument to H_bisym_stat to increase the
number of Monte Carlo samples.
#Hollander's test of bivariate symmetry
#Calculating test statistic for Hollander's test of bivariate
symmetry
H_bisym_stat <- function(x,y){
minXY <- pmin(x,y);
sortX <- x[order(minXY)];
sortY <- y[order(minXY)];
minXY <- pmin(sortX, sortY);
R <- as.numeric(sortX <= sortY);
maxXY <- pmax(sortX,sortY);
D1 <- matrix(0, length(x), length(x));
D2 <- matrix(0, length(x), length(x));
for (i in 1:length(x)){
for (j in 1:length(x)){
D1[i,j] = as.numeric(minXY[j] <
maxXY[i] & maxXY[i] <= maxXY[j]);
D2[i,j] = as.numeric(minXY[i] <=
minXY[j]);
}
}
D = D1*D2;
S = 2*R - 1;
T = t( t(S) %*% D);
H_squared = (1/length(x)^2)*( t(T) %*% T );
H_squared};
#Create a permutation sample for paired data;
permsamp_paired <- function(x,y){
new <- rbinom(length(x),1,0.5);
x_new <- x*new + y*(1 - new);
y_new <- y*new + x*(1 - new);
list(x_new = x_new, y_new = y_new)};
#Calculate pvalue for Hollander statistic;
H_bisym_pv <- function(x, y, permutations = 1000){
H_squared = H_bisym_stat(x=x, y=y);
numb_greater = 0;
for (i in 1:permutations) {
newvars <- permsamp_paired(x=x, y=y);
x_new <- newvars$x_new;
y_new <- newvars$y_new;
newstat <- H_bisym_stat(x=x_new, y=y_new);
numb_greater = numb_greater +
as.numeric(newstat >= H_squared);
}
pv <- numb_greater/permutations;
pv };
Joe Boyer
Statistical Sciences
GlaxoSmithKline
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
