With coordination with the code's author (Daniel),
The updated code has been uploaded to github here:
https://github.com/talgalili/R-code-snippets/blob/master/siegel.tukey.r
And also the following post was updated with the code:
http://www.r-statistics.com/2010/02/siegel-tukey-a-non-parametric-test-for-equality-in-variability-r-code/
I suspect that the code still needs some tweaks so it will be able to take
care of two vectors of different lengths.
Tal
----------------Contact
Details:-------------------------------------------------------
Contact me: tal.gal...@gmail.com | 972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
----------------------------------------------------------------------------------------------
On Fri, Mar 9, 2012 at 11:11 PM, Daniel Malter <dan...@umd.edu> wrote:
> #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.
>
> ______________________________________________
> 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.
>
[[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.
