[R] simulating 'non-standard' survival data

[vito muggeo]

Dear all,
I'm looking for someone that help me to write an R function to simulate
survival data under complex situations, namely time-varying hazard ratio,
marginal distribution of survival times and covariates. The algorithm is
described in the reference below and it should be not very difficult to
implement it. However I tried but without success....;-(
Below there the code that I used; it works but the relevant results (i.e.
the estimates) are rather far from the true value; I cannot understand where
the error is.

People interested in, could send me an e-mail privately.

best,
vito

@article{mackenzie02,
    author={T. Mackenzie and M. Abrahamowicz},
    title={Marginal and hazard ratio specific random data: Applications to
semi-parametric bootstrapping},
    journal={Statistics and Computing},
    volume={12},
    year={2002},
    pages={245-252}
}

######################à
Code to simulate survival data under complex situations (MacKenzie and
Abrahamowicz, Statistics and Computing 2002)
#######################

#sample size
n<-100
R<-1:n
#explanatory variable
x<-runif(n,3,9)
#survival time and censoring variable:
obsTime<-rexp(n,1)
status<-ifelse(runif(n)>.2,1,0)
obsTime<-obsTime[order(obsTime)]
status<-status[order(obsTime)]
dati<-matrix(-99,n,(1+2))

#hazard ratio as function of the explanatory variable xx, and survival time.
HR<-function(ttime,xx){exp(.5)} #.5 is the log hazard ratio

P<-outer(obsTime,x,HR)
PP<-t(apply(P,1,function(x){sort(cumsum(sort(x,decreasing=T)),decreasing=T)}
))
Pok<-P/PP
for(i in 1:(n-1)){
        j<- if(status[i]==1) sample(R,size=1,prob=Pok[i,]) else
sample(R,size=1)
        ind<-which(R==j)
        dati[i,]<-c(obsTime[j],status[j],x[j])
        filtro<-R!=j
        R<-R[filtro]
        Pok<-Pok[,filtro]
            }
dati[nrow(dati),]<-c(obsTime[nrow(dati)],status[nrow(dati)],x[nrow(dati)])

dati<-data.frame(dati)
names(dati)<-c("SurvTime","cens","x")
library(survival)
coxph(Surv(SurvTime,cens)~x,dati)