Hi, I recently ran into the problem that I needed a Siegel-Tukey test for
equal variability based on ranks. Maybe there is a package that has it
implemented, but I could not find it. So I programmed an R function to do
it. The Siegel-Tukey test requires to recode the ranks so that they express
variability rather than ascending order. This is essentially what the code
further below does. After the rank transformation, a regular Mann-Whitney U
test is applied. The "manual" and code are pasted below.
Best,
Daniel
Siegel-Tukey test
Description: Non-parametric Siegel-Tukey test for equality in variability.
The null hypothesis is that the variability of x is equal between two
groups. A rejection of the null indicates that variability differs between
the two groups.
Usage: siegel.tukey(x,y,id.col=FALSE,adjust.median=FALSE,rnd=8, ...)
Arguments:
x: a vector of data
y: Data of the second group (if id.col=FALSE) or group indicator (if
id.col=TRUE). In the latter case, y MUST take 1 or 2 to indicate
observations of group 1 and 2, respectively, and x must contain the data for
both groups.
id.col: If FALSE (default), then x and y are the data columns for group 1
and 2, respectively. If TRUE, the y is the group indicator.
adjust.median: Should between-group differences in medians be leveled
before performing the test? In certain cases, the Siegel-Tukey test is
susceptible to median differences and may indicate significant differences
in variability that, in reality, stem from differences in medians.
rnd: Should the data be rounded and, if so, to which decimal? The default
(-1) uses the data as is. Otherwise, rnd must be a non-negative integer.
Typically, this option is not needed. However, occasionally, differences in
the precision with which certain functions return values cause the merging
of two data frames to fail within the siegel.tukey function. Only then
rounding is necessary. This operation should not be performed if it affects
the ranks of observations.
... arguments passed on to the Wilcoxon test. See ?wilcox.test
Value: Among other output, the function returns rank sums for the two
groups, the associated Wilcoxon's W, and the p-value for a Wilcoxon test on
tie-adjusted Siegel-Tukey ranks (i.e., it performs and returns a
Siegel-Tukey test). If significant, the group with the smaller rank sum has
greater variability.
References: Sidney Siegel and John Wilder Tukey (1960) "A nonparametric sum
of ranks procedure for relative spread in unpaired samples." Journal of the
American Statistical Association. See also, David J. Sheskin (2004)
"Handbook of parametric and nonparametric statistical procedures." 3rd
edition. Chapman and Hall/CRC. Boca Raton, FL.
Notes: The Siegel-Tukey test has relatively low power and may, under certain
conditions, indicate significance due to differences in medians rather than
differences in variabilities (consider using the argument adjust.median).
Output (in this order)
1. Group medians
2. Wilcoxon-test for between-group differences in median (after the median
adjustment if specified)
3. Unique values of x and their tie-adjusted Siegel-Tukey ranks
4. Xs of group 1 and their tie-adjusted Siegel-Tukey ranks
5. Xs of group 2 and their tie-adjusted Siegel-Tukey ranks
6. Siegel-Tukey test (Wilcoxon test on tie-adjusted Siegel-Tukey ranks)
siegel.tukey=function(x,y,id.col=FALSE,adjust.median=F,rnd=-1,alternative="t
wo.sided",mu=0,paired=FALSE,exact=FALSE,correct=TRUE,conf.int=FALSE,conf.lev
el=0.95){
if(id.col==FALSE){
data=data.frame(c(x,y),rep(c(1,2),c(length(x),length(y))))
} else {
data=data.frame(x,y)
}
names(data)=c("x","y")
data=data[order(data$x),]
if(rnd>-1){data$x=round(data$x,rnd)}
if(adjust.median==T){
data$x[data$y==1]=data$x[data$y==1]-(median(data$x[data$y==1])-median(data$x
[data$y==2]))/2
data$x[data$y==2]=data$x[data$y==2]-(median(data$x[data$y==2])-median(data$x
[data$y==1]))/2
}
cat("Median of group 1 = ",median(data$x[data$y==1]),"\n")
cat("Median of group 2 =
",median(data$x[data$y==2]),"\n","\n")
cat("Test of median differences","\n")
print(wilcox.test(data$x[data$y==1],data$x[data$y==y]))
a=rep(seq(ceiling(length(data$x)/4)),each=2)
b=rep(c(0,1),ceiling(length(data$x)/4))
rk.up=c(1,(a*4+b))[1:ceiling(length(data$x)/2)]
rk.down=rev(c(a*4+b-2)[1:floor(length(data$x)/2)])
cat("Performing Siegel-Tukey rank
transformation...","\n","\n")
rks=c(rk.up,rk.down)
unqs=unique(sort(data$x))
corr.rks=tapply(rks,data$x,mean)
cbind(unqs,corr.rks)
rks.data=data.frame(unqs,corr.rks)
names(rks.data)=c("unique values of x","tie-adjusted
Siegel-Tukey rank")
print(rks.data,row.names=F)
names(rks.data)=c("unqs","corr.rks")
data=merge(data,rks.data,by.x="x",by.y="unqs")
rk1=data$corr.rks[data$y==1]
rk2=data$corr.rks[data$y==2]
cat("\n","Tie-adjusted Siegel-Tukey ranks of group
1","\n")
group1=data.frame(data$x[data$y==1],rk1)
names(group1)=c("x","rank")
print(group1,row.names=F)
cat("\n","Tie-adjusted Siegel-Tukey ranks of group
2","\n")
group2=data.frame(data$x[data$y==2],rk2)
names(group2)=c("x","rank")
print(group2,row.names=F)
cat("\n")
cat("Siegel-Tukey test","\n")
cat("Siegel-Tukey rank transformation performed.","Tie adjusted
ranks
computed.","\n")
if(adjust.median==T)cat("Medians adjusted to
equality.","\n") else
cat("Medians not adjusted.","\n")
cat("Rank sum of group 1 =", sum(rk1)," Rank sum of group 2
=",sum(rk2),"\n")
print(wilcox.test(rk1,rk2,alternative=alternative,mu=mu,paired=paired,exactexact,correct=correct,conf.int=conf.int,conf.level=conf.level))
}
#Example:
x=c(4,4,5,5,6,6)
y=c(0,0,1,9,10,10)
siegel.tukey(x,y)
The previously posted code contains bugs. The code below should work:
x: a vector of data
y: Group indicator (if id.col=TRUE); data of the second group (if
id.col=FALSE). If y is the group indicator it MUST take 0 or 1 to indicate
the groups, and x must contain the data for both groups.
id.col: If TRUE (default), then x is the data column and y is the ID column,
indicating the groups. If FALSE, x and y are both data columns. id.col must
be FALSE only if both data columns are of the same length.
adjust.median: Should between-group differences in medians be leveled before
performing the test? In certain cases, the Siegel-Tukey test is susceptible
to median differences and may indicate significant differences in
variability that, in reality, stem from differences in medians.
rnd: Should the data be rounded and, if so, to which decimal? The default
(-1) uses the data as is. Otherwise, rnd must be a non-negative integer.
Typically, this option is not needed. However, occasionally, differences in
the precision with which certain functions return values cause the merging
of two data frames to fail within the siegel.tukey function. Only then
rounding is necessary. This operation should not be performed if it affects
the ranks of observations.
? arguments passed on to the Wilcoxon test. See ?wilcox.test
Value: Among other output, the function returns the data, the Siegel-Tukey
ranks, the associated Wilcoxon?s W and the p-value for a Wilcoxon test on
tie-adjusted Siegel-Tukey ranks (i.e., it performs and returns a
Siegel-Tukey test). If significant, the group with the smaller rank sum has
greater variability.
References: Sidney Siegel and John Wilder Tukey (1960) ?A nonparametric sum
of ranks procedure for relative spread in unpaired samples.? Journal of the
American Statistical Association. See also, David J. Sheskin (2004)
?Handbook of parametric and nonparametric statistical procedures.? 3rd
edition. Chapman and Hall/CRC. Boca Raton, FL.
Notes: The Siegel-Tukey test has relatively low power and may, under certain
conditions, indicate significance due to differences in medians rather than
differences in variabilities (consider using the argument adjust.median).
Output (in this order)
1. Group medians (after median adjustment if specified)
2. Wilcoxon-test for between-group differences in medians (after the median
adjustment if specified)
3. Data, group membership, and the Siegel-Tukey ranks
4. Mean Siegel-Tukey rank by group (smaller values indicate greater
variability)
5. Siegel-Tukey test (Wilcoxon test on tie-adjusted Siegel-Tukey ranks)
siegel.tukey=function(x,y,id.col=TRUE,adjust.median=F,rnd=-1,alternative="two.sided",mu=0,paired=FALSE,exact=FALSE,correct=TRUE,conf.int=FALSE,conf.level=0.95){
if(id.col==FALSE){
data=data.frame(c(x,y),rep(c(0,1),c(length(x),length(y))))
} else {
data=data.frame(x,y)
}
names(data)=c("x","y")
data=data[order(data$x),]
if(rnd>-1){data$x=round(data$x,rnd)}
if(adjust.median==T){
cat("\n","Adjusting
medians...","\n",sep="")
data$x[data$y==0]=data$x[data$y==0]-(median(data$x[data$y==0]))
data$x[data$y==1]=data$x[data$y==1]-(median(data$x[data$y==1]))
}
cat("\n","Median of group 1 =
",median(data$x[data$y==0]),"\n",sep="")
cat("Median of group 2 =
",median(data$x[data$y==1]),"\n","\n",sep="")
cat("Testing median differences...","\n")
print(wilcox.test(data$x[data$y==0],data$x[data$y==1]))
cat("Performing Siegel-Tukey rank
transformation...","\n","\n")
sort.x<-sort(data$x)
sort.id<-data$y[order(data$x)]
data.matrix<-data.frame(sort.x,sort.id)
base1<-c(1,4)
iterator1<-matrix(seq(from=1,to=length(x),by=4))-1
rank1<-apply(iterator1,1,function(x) x+base)
iterator2<-matrix(seq(from=2,to=length(x),by=4))
base2<-c(0,1)
rank2<-apply(iterator2,1,function(x) x+base2)
#print(rank1)
#print(rank2)
if(length(rank1)==length(rank2)){
rank<-c(rank1[1:floor(length(x)/2)],rev(rank2[1:ceiling(length(x)/2)]))
} else{
rank<-c(rank1[1:ceiling(length(x)/2)],rev(rank2[1:floor(length(x)/2)]))
}
unique.ranks<-tapply(rank,sort.x,mean)
unique.x<-as.numeric(as.character(names(unique.ranks)))
rank.matrix<-data.frame(unique.x,unique.ranks)
ST.matrix<-merge(data.matrix,rank.matrix,by.x="sort.x",by.y="unique.x")
print(ST.matrix)
cat("\n","Performing Siegel-Tukey
test...","\n",sep="")
ranks0<-ST.matrix$unique.ranks[ST.matrix$sort.id==0]
ranks1<-ST.matrix$unique.ranks[ST.matrix$sort.id==1]
cat("\n","Mean rank of group 0:
",mean(ranks0),"\n",sep="")
cat("Mean rank of group 1:
",mean(ranks1),"\n",sep="")
print(wilcox.test(ranks0,ranks1,alternative=alternative,mu=mu,paired=paired,exact=exact,correct=correct,conf.int=conf.int,conf.level=conf.level))
}
Examples:
x <- c(33, 62, 84, 85, 88, 93, 97, 4, 16, 48, 51, 66, 98)
id <- c(0,0,0,0,0,0,0,1,1,1,1,1,1)
siegel.tukey(x,id,adjust.median=F,exact=T)
x<-c(0,0,1,4,4,5,5,6,6,9,10,10)
id<-c(0,0,0,1,1,1,1,1,1,0,0,0)
siegel.tukey(x,id)
x <- c(85,106,96, 105, 104, 108, 86)
id<-c(0,0,1,1,1,1,1)
siegel.tukey(x,id)
x<-c(177,200,227,230,232,268,272,297,47,105,126,142,158,172,197,220,225,230,262,270)
id<-c(rep(0,8),rep(1,12))
siegel.tukey(x,id,adjust.median=T)
--
View this message in context:
http://r.789695.n4.nabble.com/Siegel-Tukey-test-for-equal-variability-code-tp1565053p4430608.html
Sent from the R help mailing list archive at Nabble.com.
#The code of rank 1 in the previous post should have read
#rank1<-apply(iterator1,1,function(x) x+base1)
#corrected code below
siegel.tukey=function(x,y,id.col=TRUE,adjust.median=F,rnd=-1,alternative="two.sided",mu=0,paired=FALSE,exact=FALSE,correct=TRUE,conf.int=FALSE,conf.level=0.95){
if(id.col==FALSE){
data=data.frame(c(x,y),rep(c(0,1),c(length(x),length(y))))
} else {
data=data.frame(x,y)
}
names(data)=c("x","y")
data=data[order(data$x),]
if(rnd>-1){data$x=round(data$x,rnd)}
if(adjust.median==T){
cat("\n","Adjusting
medians...","\n",sep="")
data$x[data$y==0]=data$x[data$y==0]-(median(data$x[data$y==0]))
data$x[data$y==1]=data$x[data$y==1]-(median(data$x[data$y==1]))
}
cat("\n","Median of group 1 =
",median(data$x[data$y==0]),"\n",sep="")
cat("Median of group 2 =
",median(data$x[data$y==1]),"\n","\n",sep="")
cat("Testing median differences...","\n")
print(wilcox.test(data$x[data$y==0],data$x[data$y==1]))
cat("Performing Siegel-Tukey rank
transformation...","\n","\n")
sort.x<-sort(data$x)
sort.id<-data$y[order(data$x)]
data.matrix<-data.frame(sort.x,sort.id)
base1<-c(1,4)
iterator1<-matrix(seq(from=1,to=length(x),by=4))-1
rank1<-apply(iterator1,1,function(x) x+base1)
iterator2<-matrix(seq(from=2,to=length(x),by=4))
base2<-c(0,1)
rank2<-apply(iterator2,1,function(x) x+base2)
#print(rank1)
#print(rank2)
if(length(rank1)==length(rank2)){
rank<-c(rank1[1:floor(length(x)/2)],rev(rank2[1:ceiling(length(x)/2)]))
} else{
rank<-c(rank1[1:ceiling(length(x)/2)],rev(rank2[1:floor(length(x)/2)]))
}
unique.ranks<-tapply(rank,sort.x,mean)
unique.x<-as.numeric(as.character(names(unique.ranks)))
rank.matrix<-data.frame(unique.x,unique.ranks)
ST.matrix<-merge(data.matrix,rank.matrix,by.x="sort.x",by.y="unique.x")
print(ST.matrix)
cat("\n","Performing Siegel-Tukey
test...","\n",sep="")
ranks0<-ST.matrix$unique.ranks[ST.matrix$sort.id==0]
ranks1<-ST.matrix$unique.ranks[ST.matrix$sort.id==1]
cat("\n","Mean rank of group 0:
",mean(ranks0),"\n",sep="")
cat("Mean rank of group 1:
",mean(ranks1),"\n",sep="")
print(wilcox.test(ranks0,ranks1,alternative=alternative,mu=mu,paired=paired,exact=exact,correct=correct,conf.int=conf.int,conf.level=conf.level))
}
Examples:
x <- c(33, 62, 84, 85, 88, 93, 97, 4, 16, 48, 51, 66, 98)
id <- c(0,0,0,0,0,0,0,1,1,1,1,1,1)
siegel.tukey(x,id,adjust.median=F,exact=T)
x<-c(0,0,1,4,4,5,5,6,6,9,10,10)
id<-c(0,0,0,1,1,1,1,1,1,0,0,0)
siegel.tukey(x,id)
x <- c(85,106,96, 105, 104, 108, 86)
id<-c(0,0,1,1,1,1,1)
siegel.tukey(x,id)
x<-c(177,200,227,230,232,268,272,297,47,105,126,142,158,172,197,220,225,230,262,270)
id<-c(rep(0,8),rep(1,12))
siegel.tukey(x,id,adjust.median=T)
--
View this message in context:
http://r.789695.n4.nabble.com/Siegel-Tukey-test-for-equal-variability-code-tp1565053p4460705.html
Sent from the R help mailing list archive at Nabble.com.