R/hochberg.R

hochberg <-
function(x,x2=NA,cil=NA,crit=NA,con=0,tr=.2,alpha=.05,iter=10000,SEED=TRUE){
#
# A generalization of Hochberg's method
# method to trimmed mean.
#
# x contains first stage data
# x2 contains second stage data
#
# cil is the desired length of the confidence intervals.
# That is, cil is the distance between the upper and lower
# ends of the confidence intervals.
#
x3<-x2
if(is.matrix(x))x<-listm(x)
if(!is.list(x))stop("Data must be stored in list mode or in matrix mode.")
J<-length(x)
tempn<-0
svec<-NA
for(j in 1:J){
temp<-x[[j]]
temp<-temp[!is.na(temp)] # Remove missing values.
tempn[j]<-length(temp)
x[[j]]<-temp
svec[j]<-winvar(temp,tr=tr)/(1-2*tr)^2
}
tempt<-floor((1-2*tr)*tempn)
A<-sum(1/(tempt-1))
df<-J/A
print(paste("If using the tables of Studentized range distribution,"))
print(paste("the degrees of freedom are:",df))
if(!is.list(x2) && !is.matrix(x2)){
x2<-list()
for(j in 1:J)x2[[j]]<-NA
}
if(is.na(cil))stop("To proceed, you must specify the maximum length of the confidence intervals.")
if(is.na(crit)){
print("Approximating critical value")
crit<-trange(tempn-1,alpha=alpha,iter=iter,SEED=SEED)
print(paste("The critical value is ",crit))
}
#
if(con[1] == 0){
               Jm<-J-1
               ncon <- (J^2 - J)/2
                con <- matrix(0, J, ncon)
                id <- 0
                for(j in 1:Jm) {
                        jp <- j + 1
                        for(k in jp:J) {
                                id <- id + 1
                                con[j, id] <- 1
                                con[k, id] <- 0 - 1
                        }
                }
        }
        ncon <- ncol(con)
avec<-NA
for(i in 1:ncon){
temp<-con[,i]
avec[i]<-sum(temp[temp>0])
}
dvec<-(cil/(2*crit*avec))^2
d<-max(dvec)
n.vec<-NA
for(j in 1:J){
n.vec[j]<-max(tempn[j],floor(svec[j]/d)+1)
print(paste("Need an additional ", n.vec[j]-tempn[j],
" observations for group", j))
}
#
# Do second stage if data are supplied
#
if(is.matrix(x2))x2<-listm(x2)
temp2<-n.vec-tempn
if(!is.list(x3) && !is.matrix(x3) && sum(temp2)>0)stop("No second stage data supplied; this function is terminating")
if(length(x) != length(x2))warning("Number of groups in first stage data does not match the number in the second stage.")
ci.mat<-NA
if(!is.na(x2[1]) || sum(temp2)==0){
xtil<-NA
nvec2<-NA
for(j in 1:J){
nvec2[j]<-0
temp<-x2[[j]]
if(!is.na(temp[1]))nvec2[j]<-length(x2[[j]])
if(nvec2[j] <n.vec[j]-tempn[j])warning(paste("The required number of observations for group",j," in the second stage is ",n.vec[j]-tempn[j]," but only ",nvec2[j]," are available"))
xtil[j]<-mean(c(x[[j]],x2[[j]]),tr=tr,na.rm=TRUE)
}
ci.mat<-matrix(0,ncol=3,nrow=ncon)
dimnames(ci.mat)<-list(NULL,c("con.num","ci.low","ci.high"))
for(ic in 1:ncon){
ci.mat[ic,1]<-ic
bvec<-con[,ic]*sqrt(svec/(nvec2+tempn))
A<-sum(bvec[bvec>0])
C<-0-sum(bvec[bvec<0])
D<-max(A,C)
ci.mat[ic,2]<-sum(con[,ic]*xtil)-crit*D
ci.mat[ic,3]<-sum(con[,ic]*xtil)+crit*D
}}
list(ci.mat=ci.mat,con=con)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.