gehan.obj<-function(beta,y,delta,x, wt)
{
x<-as.matrix(x)
n<-dim(x)[1]
p<-length(beta)

residual=y-as.vector(x%*%beta)

index<-order(residual)
residual.order<-residual[index]
delta.order<-delta[index]*wt[index]
 
s1<-cumsum(delta.order)
s2<-(n:1)*delta.order
return(sum((s1-s2)*residual.order))
}



gehan.fit=function(y, delta, x, wt, beta.ini=NULL)
  {if(length(beta.ini)==0)
     {beta1=lm(y~x)$coef[-1]
      beta2=lm(y~x, subset=(delta==1))$coef[-1]
      beta.ini=(beta1+beta2)/2
      }

   fit=optim(beta.ini, gehan.obj, y=y, delta=delta, x=x, wt=wt)
   beta=fit$par
   return(beta)
   }


regaft=function(y, delta, x, B=400, seed=100)
   {x<-as.matrix(x)
    n<-dim(x)[1]
    p<-dim(x)[2]

    wt0=rep(1, n)
    beta0=gehan.fit(y, delta, x, wt0)

    beta.bt=matrix(0, B, p)
    set.seed(seed)
    for(b in 1:B)
       {wt=rexp(n)
        beta.bt[b,]=gehan.fit(y, delta, x, wt=wt, beta.ini=beta0)
        }
    beta.sd=apply(beta.bt,2,sd)
    z=beta0/beta.sd
    pv=1-pchisq(z^2,1)
    result=matrix(0, p, 4)
    result[,1]=beta0
    result[,2]=beta.sd
    result[,3]=z
    result[,4]=pv
    colnames(result)=c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
    rownames(result)=colnames(x)
    return(result)
    }


#####################################################
#####################################################
####@@@@@@@@@########################################
### y:     log(observed time)   #####################
### delta: censoring indicator  #####################
### x:     covariates matrix    #####################
### B:     number of resampling #####################
### seed:  initial seed         #####################
#####################################################

### Example #########################################

p=2
n=100
x=matrix(rnorm(n*p), n, p)
lt=x[,1]+rnorm(n)
lc=rnorm(n)
y=pmin(lt, lc)
delta=1*(lt<lc)

regaft(y, delta, x)