Update: The siegel Tukey code is now fixed (both on the github page, and the blog's post): https://github.com/talgalili/R-code-snippets/blob/master/siegel.tukey.r http://www.r-statistics.com/2010/02/siegel-tukey-a-non-parametric-test-for-equality-in-variability-r-code/
Best, Tal On Sat, Mar 10, 2012 at 12:05 AM, Tal Galili <tal.gal...@gmail.com> wrote: > 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.
