##########################################################################
####This program computes the distribution-free rectangular reference ####
####region by using Tukey's equivalence blocks. The region is a 3-D   ####
####rectangle and the sample size nn must be large enough in order    ####
####for the rectangle to be constructed.                              ####
####nn: the sample size;                                              ####
####rp: the content;                                                  ####
####rgamma: the confidence level;                                     ####
####xx: the sample available.                                         ####
##########################################################################

###############################
### Read in the data first ####
###############################

####Clear the workspace first####
ls()
rm(list=ls())

rp <- 0.95
rgamma <- 0.95

####Read data into R####
#setwd("C:/Users/wl/Documents/??")
getwd()
list.files()
Kidney.data <- read.table("KidneyData.txt", header=T)   #Read in the data
#mean(Kidney.data$UACR)  #Check that the data are read in correctly
#sd(Kidney.data$UACR)
#mean(Kidney.data$gender)
#sd(Kidney.data$gender)

####Check it is a data.frame####
str(Kidney.data)
which(is.na(Kidney.data$gender))    #Any missing values: 0
which(is.na(Kidney.data$UACR))    #Any missing values: 0
which(is.na(Kidney.data$UA))   #Any missing values: 0
which(is.na(Kidney.data$SC))   #Any missing values: 0
summary(Kidney.data)

Kidney.male.data <- subset(Kidney.data,gender==1)  #Data for male
n.male <- length(Kidney.male.data$gender)  #sample size
Kidney.female.data <- subset(Kidney.data,gender==2)  #Data for female
n.female <- length(Kidney.female.data$gender)  #sample size

mean.male.UACR <- mean(Kidney.male.data$UACR)
sd.male.UACR <- sd(Kidney.male.data$UACR)
hist(Kidney.male.data$UACR)
mean.male.UA <- mean(Kidney.male.data$UA)
sd.male.UA <- sd(Kidney.male.data$UA)
hist(Kidney.male.data$UA)
mean.male.SC <- mean(Kidney.male.data$SC)
sd.male.SC <- sd(Kidney.male.data$SC)
hist(Kidney.male.data$SC)

mean.female.UACR <- mean(Kidney.female.data$UACR)
sd.female.UACR <- sd(Kidney.female.data$UACR)
hist(Kidney.female.data$UACR)
mean.female.UA <- mean(Kidney.female.data$UA)
sd.female.UA <- sd(Kidney.female.data$UA)
hist(Kidney.female.data$UA)
mean.female.SC <- mean(Kidney.female.data$SC)
sd.female.SC <- sd(Kidney.female.data$SC)
hist(Kidney.female.data$SC)


##Define the function that finds the smallest sample size nn so that ##
##a (rp, rgamma) 1-dimension 2-sided distribution-free tolerance     ##
##interval can be constructed as (xx_[1], xx_[nn]).                  ##
#######################################################################
minSampSize.1d.2s <- function(rp, rconf){
  n0 = ceiling( log(1-rconf)/log(rp) ) #this is the min sample size for 1-d 1-s
  
  nk=n0
  while ( pbeta(rp, nk-1, 2, lower.tail=FALSE) < rconf ){
    nk <- nk + 1
  }
  
  return(nk)  #the minimum sample size required
}
#nn<-minSampSize.1d.2s(rp, rgamma)
#pbeta(rp, nn-1, 2, lower.tail=FALSE) #this should be >= rgamma
#pbeta(rp, nn-2, 2, lower.tail=FALSE) #this should be < rgamma


##Define the function that finds k so that the union of k##
##equivalence blocks is a (rp, rgamma) tolerance region. ##
##This works irrespective the dimension of an observation##
###########################################################
numEquiBloc<-function(rp,rconf,nn){
  if(nn<minSampSize.1d.2s(rp,rconf)){
    print(c("The sample size is too small"))
    return("The sample size is too small")
  }else{
    nk=nn
    while ( pbeta(rp, nk, nn+1-nk, lower.tail=FALSE) > rconf ){
      nk <- nk - 1
    }
    k = nk + 1   #The  k  required
    return(k)
  } 
}
nn=2529 #nn=2726
nk=numEquiBloc(rp,rgamma,nn) #The number of blocks required to form the tolerance region
pbeta(rp, nk, nn+1-nk, lower.tail=FALSE) #the true confidence level; this sould be > rgamma
pbeta(rp, nk-1, nn+2-nk, lower.tail=FALSE) #this sould be < rgamma


##Define the function that finds the 1-dimensional distribution-free##
##2-sided and upper/lower (rp, rgamma) tolerance intervals with     ##
##given 1-d sample xx                                               ##
######################################################################
DistFreeTolInts.1d<-function(rp,rgamma,xx){
  nn <- length(xx)
  if(nn<minSampSize.1d.2s(rp,rgamma)){
    print(c("The sample size is too small"))
    return("The sample size is too small")
  }else{
    nk=numEquiBloc(rp,rgamma,nn) #The number of equi blocks required
    xx.order <- sort(xx)
    upper.limit <- xx.order[nk]
    lower.limit <- xx.order[nn-nk+1]
    k1 <- floor((nn-nk+1)/2)
    if( floor((nn-nk+1)/2)==floor((nn-nk+1)/2 +0.6) ){ 
      two.s.A.upp.limit <- xx.order[k1+nk]
      two.s.A.low.limit <- xx.order[k1]
      two.s.B.upp.limit <- two.s.A.upp.limit
      two.s.B.low.limit <- two.s.A.low.limit
    }else{
      two.s.A.upp.limit <- xx.order[k1+nk]
      two.s.A.low.limit <- xx.order[k1]
      two.s.B.upp.limit <- xx.order[k1+nk+1]
      two.s.B.low.limit <- xx.order[k1+1]
    }
    return(list("Lower.limit",lower.limit,"Upper.limit",upper.limit,"two.sided",
     c(two.s.A.low.limit,two.s.A.upp.limit,two.s.B.low.limit,two.s.B.upp.limit) ))
  } 
}

##1-d tolerance intervals for the Kidney data##
DistFreeTolInts.1d(rp,rgamma,Kidney.male.data$UACR)
DistFreeTolInts.1d(rp,rgamma,Kidney.male.data$UA)
DistFreeTolInts.1d(rp,rgamma,Kidney.male.data$SC)
DistFreeTolInts.1d(rp,rgamma,Kidney.female.data$UACR)
DistFreeTolInts.1d(rp,rgamma,Kidney.female.data$UA)
DistFreeTolInts.1d(rp,rgamma,Kidney.female.data$SC)


#################################
###Now 3-dimensional situation###
#################################

##Define the function that finds the smallest sample size nn so that ##
##a (rp, rgamma) distribution-free 3-d tolerance region can be       ##
##constructed via the union of several equivalence blocks            ##
#######################################################################
minSampSize.3d.2s <- function(rp, rgamma){
  n0 = ceiling( log(1-rgamma)/log(rp) )
  
  nk=n0
  while ( pbeta(rp, nk-5, 6, lower.tail=FALSE) < rgamma ){
    nk <- nk + 1
  }
  
  return(nk)  #the minimum sample size required
}
nn<-minSampSize.3d.2s(rp, rgamma)
pbeta(rp, nn-5, 6, lower.tail=FALSE) #this should be >= rgamma
pbeta(rp, nn-6, 6, lower.tail=FALSE) #this should be < rgamma




#########################################################################
##Define the main function that finds the (rp,rgamma) distribution-free##
##tolerance box by deleting nn+1-nk equivalence blocks, in the         ## 
## order of R(ight),T(op),F(ront),L(eft),B(ottom),H(ind).              ##
##rp: the content;                                                     ##
##rgamma: the confidence level;                                        ##
##xx: the sample available in the form of a nn*3 matrix.               ##
#########################################################################
DistFreeTolRect.3d <- function(xx,rp,rgamma){
  nn <- dim(xx)[1]
  nk <- numEquiBloc(rp,rgamma,nn)
  nkD <- nn+1-nk #number of equiv blocks to delete
  nkDCycl <- floor(nkD/6) #number of cycles for striping 6 blocks R,T,F,L,B,H
  nkDRem <- nkD - 6*nkDCycl #number of blocks (<6) to strip further
  
  for (i in 1:nkDCycl) { #Strip away nkDCycl cycles of  6 blocks R,T,F,L,B,H
    
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample; [1] is used
                            # in case there may be ties
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
    
    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
    
    minLind <- which(xx[,1]==min(xx[,1])) #Indices of the obs having the min x_1 value
    minLval <- xx[minLind[1],1] #Record the min x_1 value
    xx <- xx[-minLind[1],]  #Delete this observation from the sample
    
    minBind <- which(xx[,2]==min(xx[,2])) #Indices of the obs having the min x_2 value
    minBval <- xx[minBind[1],2] #Record the min x_2 value
    xx <- xx[-minBind[1],]  #Delete this observation from the sample
    
    minHind <- which(xx[,3]==min(xx[,3])) #Indices of the obs having the min x_3 value
    minHval <- xx[minHind[1],3] #Record the min x_3 value
    xx <- xx[-minHind[1],]  #Delete this observation from the sample
  }
  
  if(nkDRem == 0){ #Strp the remaining (<6) block(s)
    "Do nothing"
  }else if(nkDRem == 1){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
  }else if(nkDRem == 2){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
  }else if(nkDRem == 3){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
    
    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
  }else if(nkDRem == 4){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
    
    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
    
    minLind <- which(xx[,1]==min(xx[,1])) #Indices of the obs having the min x_1 value
    minLval <- xx[minLind[1],1] #Record the min x_1 value
    xx <- xx[-minLind[1],]  #Delete this observation from the sample
  }else if(nkDRem == 5){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
    
    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
    
    minLind <- which(xx[,1]==min(xx[,1])) #Indices of the obs having the min x_1 value
    minLval <- xx[minLind[1],1] #Record the min x_1 value
    xx <- xx[-minLind[1],]  #Delete this observation from the sample
    
    minBind <- which(xx[,2]==min(xx[,2])) #Indices of the obs having the min x_2 value
    minBval <- xx[minBind[1],2] #Record the min x_1 value
    xx <- xx[-minBind[1],]  #Delete this observation from the sample
  }    
  return(list( c(maxRval,maxTval,maxFval,minLval,minBval,minHval), c(nn,nk, nkD) ))
}


##Now for the Kidney data##
#is.matrix(cbind(Kidney.male.data$UACR,Kidney.male.data$UA,Kidney.male.data$SC))
bounds.male <- DistFreeTolRect.3d(cbind(Kidney.male.data$UACR,Kidney.male.data$UA,Kidney.male.data$SC),rp,rgamma)[[1]]
bounds.female <- DistFreeTolRect.3d(cbind(Kidney.female.data$UACR,Kidney.female.data$UA,Kidney.female.data$SC),rp,rgamma)[[1]]

#Plot the figure tolerance boxes#
#dev.off()
opar <- par()  # Record the old setting
par()$mfrow
par(mfrow = c(3,2))

plot(Kidney.male.data$UACR,Kidney.male.data$UA, col="red", xlim=c(0.0, 40), ylim=c(0,16), 
     xlab = expression(UACA (mg/g)),ylab = expression(UA (mg/dL)))
xxx <- c(bounds.male[4],bounds.male[1],bounds.male[1],bounds.male[4])
yyy <- c(bounds.male[5],bounds.male[5],bounds.male[2],bounds.male[2])
polygon(xxx, yyy, col = NA, border = "blue")
title("UA vs. UACA reference region (Male)")

plot(Kidney.male.data$UACR,Kidney.male.data$SC, col="red", xlim=c(0.0, 40), ylim=c(0,5), 
     xlab = expression(UACA (mg/g)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.male[4],bounds.male[1],bounds.male[1],bounds.male[4])
yyy <- c(bounds.male[6],bounds.male[6],bounds.male[3],bounds.male[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UACA reference region (Male)")

plot(Kidney.male.data$UA,Kidney.male.data$SC, col="red", xlim=c(0.0, 16), ylim=c(0,5), 
     xlab = expression(UA (mg/dL)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.male[5],bounds.male[2],bounds.male[2],bounds.male[5])
yyy <- c(bounds.male[6],bounds.male[6],bounds.male[3],bounds.male[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UA reference region (Male)")


plot(Kidney.female.data$UACR,Kidney.female.data$UA, col="red", xlim=c(0.0, 40), ylim=c(0,16), 
     xlab = expression(UACA (mg/g)),ylab = expression(UA (mg/dL)))
xxx <- c(bounds.female[4],bounds.female[1],bounds.female[1],bounds.female[4])
yyy <- c(bounds.female[5],bounds.female[5],bounds.female[2],bounds.female[2])
polygon(xxx, yyy, col = NA, border = "blue")
title("UA vs. UACA reference region (Female)")

plot(Kidney.female.data$UACR,Kidney.female.data$SC, col="red", xlim=c(0.0, 40), ylim=c(0,5), 
     xlab = expression(UACA (mg/g)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.female[4],bounds.female[1],bounds.female[1],bounds.female[4])
yyy <- c(bounds.female[6],bounds.female[6],bounds.female[3],bounds.female[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UACA reference region (Female)")

plot(Kidney.female.data$UA,Kidney.female.data$SC, col="red", xlim=c(0.0, 16), ylim=c(0,5), 
     xlab = expression(UA (mg/dL)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.female[5],bounds.female[2],bounds.female[2],bounds.female[5])
yyy <- c(bounds.female[6],bounds.female[6],bounds.female[3],bounds.female[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UA reference region (Female)")

par(opar)  #Return to the old setting





#########################################################################
##Define the main function that finds the (rp,rgamma) distribution-free##
##tolerance mixed-rectangle, upper-limit for x_1=UACA and two-sided    ##
##intervals for x_2=UA and x_3=SC, by deleting nn+1-nk equivalence     ## 
##blocks, in the order of R(ight),T(op),F(ront),B(ottom),H(ind).       ##
##rp: the content;                                                     ##
##rgamma: the confidence level;                                        ##
##xx: the sample available in the form of a nn*3 matrix.               ##
#########################################################################
DistFreeTolMixRect.3d <- function(xx,rp,rgamma){
  nn <- dim(xx)[1]
  nk <- numEquiBloc(rp,rgamma,nn)
  nkD <- nn+1-nk #number of equiv blocks to delete
  nkDCycl <- floor(nkD/5) #number of cycles for striping 5 blocks R,T,F,B,H
  nkDRem <- nkD - 5*nkDCycl #number of blocks (<5) to strip off further
  
  for (i in 1:nkDCycl) { #Strip away nkDCycl cycles of  5 blocks R,T,F,B,H
    
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value; [1] is used in case of ties
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample 
                                                      
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
    
    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
    
#    minLind <- which(xx[,1]==min(xx[,1])) #Indices of the obs having the min x_1 value
#    minLval <- xx[minLind[1],1] #Record the min x_1 value
#    xx <- xx[-minLind[1],]  #Delete this observation from the sample
    
    minBind <- which(xx[,2]==min(xx[,2])) #Indices of the obs having the min x_2 value
    minBval <- xx[minBind[1],2] #Record the min x_2 value
    xx <- xx[-minBind[1],]  #Delete this observation from the sample
    
    minHind <- which(xx[,3]==min(xx[,3])) #Indices of the obs having the min x_3 value
    minHval <- xx[minHind[1],3] #Record the min x_3 value
    xx <- xx[-minHind[1],]  #Delete this observation from the sample
  }
  
  if(nkDRem == 0){ #Strp the remaining (<5) block(s)
    "Do nothing"
  }else if(nkDRem == 1){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
  }else if(nkDRem == 2){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
  }else if(nkDRem == 3){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
    
    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
#  }else if(nkDRem == 4){
#    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
#    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
#    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
#    
#    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
#    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
#    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
#    
#    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
#    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
#    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
#    
#    minLind <- which(xx[,1]==min(xx[,1])) #Indices of the obs having the min x_1 value
#    minLval <- xx[minLind[1],1] #Record the min x_1 value
#    xx <- xx[-minLind[1],]  #Delete this observation from the sample
#  }else if(nkDRem == 5){
#    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
#    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
#    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
#    
#    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
#    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
#    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
#    
#    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
#    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
#    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
#    
#    minLind <- which(xx[,1]==min(xx[,1])) #Indices of the obs having the min x_1 value
#    minLval <- xx[minLind[1],1] #Record the min x_1 value
#    xx <- xx[-minLind[1],]  #Delete this observation from the sample
#    
#    minBind <- which(xx[,2]==min(xx[,2])) #Indices of the obs having the min x_2 value
#    minBval <- xx[minBind[1],2] #Record the min x_1 value
#    xx <- xx[-minBind[1],]  #Delete this observation from the sample
#  }    
  }else if(nkDRem == 4){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
        
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
      
    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
        
    minBind <- which(xx[,2]==min(xx[,2])) #Indices of the obs having the min x_2 value
    minBval <- xx[minBind[1],2] #Record the min x_1 value
    xx <- xx[-minBind[1],]  #Delete this observation from the sample
  }    
  return(list( c(maxRval,maxTval,maxFval,0,minBval,minHval), c(nn,nk, nkD) ))
}


##Now for the Kidney data##
#is.matrix(cbind(Kidney.male.data$UACR,Kidney.male.data$UA,Kidney.male.data$SC))
bounds.male <- DistFreeTolMixRect.3d(cbind(Kidney.male.data$UACR,Kidney.male.data$UA,Kidney.male.data$SC),rp,rgamma)[[1]]
bounds.female <- DistFreeTolMixRect.3d(cbind(Kidney.female.data$UACR,Kidney.female.data$UA,Kidney.female.data$SC),rp,rgamma)[[1]]

#Plot the figure of mix-sided boxes: figure3 in the paper#
#dev.off()
opar <- par()  # Record the old setting
par()$mfrow
par(mfrow = c(3,2))

plot(Kidney.male.data$UACR,Kidney.male.data$UA, col="red", xlim=c(0.0, 40), ylim=c(0,16), 
     xlab = expression(UACA (mg/g)),ylab = expression(UA (mg/dL)))
xxx <- c(bounds.male[4],bounds.male[1],bounds.male[1],bounds.male[4])
yyy <- c(bounds.male[5],bounds.male[5],bounds.male[2],bounds.male[2])
polygon(xxx, yyy, col = NA, border = "blue")
title("UA vs. UACA reference region (Male)")

plot(Kidney.male.data$UACR,Kidney.male.data$SC, col="red", xlim=c(0.0, 40), ylim=c(0,5), 
     xlab = expression(UACA (mg/g)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.male[4],bounds.male[1],bounds.male[1],bounds.male[4])
yyy <- c(bounds.male[6],bounds.male[6],bounds.male[3],bounds.male[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UACA reference region (Male)")

plot(Kidney.male.data$UA,Kidney.male.data$SC, col="red", xlim=c(0.0, 16), ylim=c(0,5), 
     xlab = expression(UA (mg/dL)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.male[5],bounds.male[2],bounds.male[2],bounds.male[5])
yyy <- c(bounds.male[6],bounds.male[6],bounds.male[3],bounds.male[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UA reference region (Male)")


plot(Kidney.female.data$UACR,Kidney.female.data$UA, col="red", xlim=c(0.0, 40), ylim=c(0,16), 
     xlab = expression(UACA (mg/g)),ylab = expression(UA (mg/dL)))
xxx <- c(bounds.female[4],bounds.female[1],bounds.female[1],bounds.female[4])
yyy <- c(bounds.female[5],bounds.female[5],bounds.female[2],bounds.female[2])
polygon(xxx, yyy, col = NA, border = "blue")
title("UA vs. UACA reference region (Female)")

plot(Kidney.female.data$UACR,Kidney.female.data$SC, col="red", xlim=c(0.0, 40), ylim=c(0,5), 
     xlab = expression(UACA (mg/g)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.female[4],bounds.female[1],bounds.female[1],bounds.female[4])
yyy <- c(bounds.female[6],bounds.female[6],bounds.female[3],bounds.female[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UACA reference region (Female)")

plot(Kidney.female.data$UA,Kidney.female.data$SC, col="red", xlim=c(0.0, 16), ylim=c(0,5), 
     xlab = expression(UA (mg/dL)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.female[5],bounds.female[2],bounds.female[2],bounds.female[5])
yyy <- c(bounds.female[6],bounds.female[6],bounds.female[3],bounds.female[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UA reference region (Female)")

par(opar)  #Return to the old setting






#########################################################################
##Define the main function that finds the (rp,rgamma) distribution-free##
##tolerance upper-rectangle, upper-limits for x_1=UACA, x_2=UA and     ##
##x_3=SC, by deleting nn+1-nk equivalence                              ## 
##blocks, in the order of R(ight),T(op),F(ront).                       ##
##For lower-limit for x_i, simply use -x_i as the i-th variable.       ##
##rp: the content;                                                     ##
##rgamma: the confidence level;                                        ##
##xx: the sample available in the form of a nn*3 matrix.               ##
#########################################################################
DistFreeTolUppRect.3d <- function(xx,rp,rgamma){
  nn <- dim(xx)[1]
  nk <- numEquiBloc(rp,rgamma,nn)
  nkD <- nn+1-nk #number of equiv blocks to delete
  nkDCycl <- floor(nkD/3) #number of cycles for striping 3 blocks R,T,F
  nkDRem <- nkD - 3*nkDCycl #number of blocks (<2) to strip off further
  
  for (i in 1:nkDCycl) { #Strip away nkDCycl cycles of  3 blocks R,T,F
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value; [1] is used in case of ties
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample 
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
    
    maxFind <- which(xx[,3]==max(xx[,3])) #Indices of the obs having the max x_3 value
    maxFval <- xx[maxFind[1],3] #Record the max x_3 value
    xx <- xx[-maxFind[1],]  #Delete this observation from the sample
    
    #    minLind <- which(xx[,1]==min(xx[,1])) #Indices of the obs having the min x_1 value
    #    minLval <- xx[minLind[1],1] #Record the min x_1 value
    #    xx <- xx[-minLind[1],]  #Delete this observation from the sample
    
    #minBind <- which(xx[,2]==min(xx[,2])) #Indices of the obs having the min x_2 value
    #minBval <- xx[minBind[1],2] #Record the min x_2 value
    #xx <- xx[-minBind[1],]  #Delete this observation from the sample
    
    #minHind <- which(xx[,3]==min(xx[,3])) #Indices of the obs having the min x_3 value
    #minHval <- xx[minHind[1],3] #Record the min x_3 value
    #xx <- xx[-minHind[1],]  #Delete this observation from the sample
  }
  
  if(nkDRem == 0){ #Strip off the remaining (<3) block(s)
    "Do nothing"
  }else if(nkDRem == 1){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
  }else if(nkDRem == 2){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
  }    
  return(list( c(maxRval,maxTval,maxFval,0,0,0), c(nn,nk, nkD) ))
}


##Now for the Kidney data##
#is.matrix(cbind(Kidney.male.data$UACR,Kidney.male.data$UA,Kidney.male.data$SC))
bounds.male <- DistFreeTolUppRect.3d(cbind(Kidney.male.data$UACR,Kidney.male.data$UA,Kidney.male.data$SC),rp,rgamma)[[1]]
bounds.female <- DistFreeTolUppRect.3d(cbind(Kidney.female.data$UACR,Kidney.female.data$UA,Kidney.female.data$SC),rp,rgamma)[[1]]

#Plot the figure of upper-sided boxes#
#dev.off()
opar <- par()  # Record the old setting
par()$mfrow
par(mfrow = c(3,2))

plot(Kidney.male.data$UACR,Kidney.male.data$UA, col="red", xlim=c(0.0, 40), ylim=c(0,16), 
     xlab = expression(UACA (mg/g)),ylab = expression(UA (mg/dL)))
xxx <- c(bounds.male[4],bounds.male[1],bounds.male[1],bounds.male[4])
yyy <- c(bounds.male[5],bounds.male[5],bounds.male[2],bounds.male[2])
polygon(xxx, yyy, col = NA, border = "blue")
title("UA vs. UACA reference region (Male)")

plot(Kidney.male.data$UACR,Kidney.male.data$SC, col="red", xlim=c(0.0, 40), ylim=c(0,5), 
     xlab = expression(UACA (mg/g)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.male[4],bounds.male[1],bounds.male[1],bounds.male[4])
yyy <- c(bounds.male[6],bounds.male[6],bounds.male[3],bounds.male[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UACA reference region (Male)")

plot(Kidney.male.data$UA,Kidney.male.data$SC, col="red", xlim=c(0.0, 16), ylim=c(0,5), 
     xlab = expression(UA (mg/dL)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.male[5],bounds.male[2],bounds.male[2],bounds.male[5])
yyy <- c(bounds.male[6],bounds.male[6],bounds.male[3],bounds.male[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UA reference region (Male)")


plot(Kidney.female.data$UACR,Kidney.female.data$UA, col="red", xlim=c(0.0, 40), ylim=c(0,16), 
     xlab = expression(UACA (mg/g)),ylab = expression(UA (mg/dL)))
xxx <- c(bounds.female[4],bounds.female[1],bounds.female[1],bounds.female[4])
yyy <- c(bounds.female[5],bounds.female[5],bounds.female[2],bounds.female[2])
polygon(xxx, yyy, col = NA, border = "blue")
title("UA vs. UACA reference region (Female)")

plot(Kidney.female.data$UACR,Kidney.female.data$SC, col="red", xlim=c(0.0, 40), ylim=c(0,5), 
     xlab = expression(UACA (mg/g)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.female[4],bounds.female[1],bounds.female[1],bounds.female[4])
yyy <- c(bounds.female[6],bounds.female[6],bounds.female[3],bounds.female[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UACA reference region (Female)")

plot(Kidney.female.data$UA,Kidney.female.data$SC, col="red", xlim=c(0.0, 16), ylim=c(0,5), 
     xlab = expression(UA (mg/dL)),ylab = expression(SC (mg/dL)))
xxx <- c(bounds.female[5],bounds.female[2],bounds.female[2],bounds.female[5])
yyy <- c(bounds.female[6],bounds.female[6],bounds.female[3],bounds.female[3])
polygon(xxx, yyy, col = NA, border = "blue")
title("SC vs. UA reference region (Female)")


par(opar)  #Return to the old setting