• R语言

  • simulation-interval censoring

一下程序是remon oller给出的一段simulation -interval censoring s-plus的程序。

我实在是不理解right和left是什么意思。我怀疑可能有打印错误？或者s-plus的方程和R不同？

请不吝指教。

function(n,local,scale,rho,p){

#simulate true time t

if(rho==0) t <- rweibull(n,1/scale,exp(local))

else t <- exp(local)*((((1-runif(n,0,1))^-rho)-1)/rho)^scale

t[t>=10] <- 10

t[t<10] <- ceiling(t[t<10])

#simulate 10 bernoulli visit: 0 missing visit 1 not missing

delta <- matrix(c(rbinom(9*n,1,p),rep(1,n)),n,10)

insp <- col(delta)>=t

#calculate lower(left) and upper(right) bound for interval censoring

right<- row(t(delta))[t(pmin(delta,insp)*t(cumsum(t(pmin(delta,insp))))

==c(1,1+cumsum(rowSums(pmin(delta,insp)))[1:(n-1)]))]

delta <- cbind(rep(1,n),delta[,1:9])

insp <- col(delta)<=t

left <- (row(t(delta))[t(pmin(delta,insp)*t(cumsum(t(pmin(delta,insp))))

==cumsum(rowSums(pmin(delta,insp))))])-1

cbind(left,right)

}

步骤大概如此：

In the generation of the simulated data for the lifetime variable T we have considered

a discretization of a variable T under an accelerated lifetime model ...

The censoring mechanism of T mimics a longitudinal study where there is a periodical follow-up with scheduled visits but patients might miss some of the appointments.

Specifically, there are assumed potential monitoring times tj = j for j = 1... 10. The

patients would assist to each of these scheduled visits with probability p. Then, for an

individual i, the observed censoring interval [Li;Ri] is constructed by designing Ri as

the first visit where the event of interest is observed and Li as the previous visit of the patient. ...

The S-plus function gendata() .. implements the generation data process. For instance,

gendata(100,0.9,0.6,1,0.5)

generates a random sample of n = 100 censoring intervals ...

4 年 后
在R package intcox 里面，有类似的生成区间删失的程序，可以参照来看.
The example dataset consists of 200 Weibull distributed random variables with shape=0.75, while scale is derived from (1/λ)^(1/shape) with λ = \exp{β_0+β'X} where β = {0.5,-0.5,0.5,0.5}' and design matrix X which is formed by the four covariates.

Furthermore, an interval was generated where the left interval end is zero and the right interval end is defined by maximum value (T.max) of the which is got by the 0.9-quantile of a Weibull r.v. with shape=0.75 and scale=median(scale).

If the value of event time >=T.max then the event time is right censored (left=T.max and right=NA). Otherwise the interval [0,T.max] was randomly divided into subintervals (grid=10) in order to determine the corresponding interval ends for each event time which is not right censored.

The generating code is given below.

Examples

## Not run:
sim.weibull.intcox.rfc <-function (N=200,beta.0=0.1,beta.cov=c(0.5,-0.5,0.5,0.5),
alpha=0.75,p.cov=c(0.5,0.75),grid=10)
{
x.design<-cbind(rbinom(N,1,p.cov[1]),rbinom(N,1,p.cov[2]),runif(N,-1,1),rnorm(N,0,1))
colnames(x.design)<-paste("x.",1:4,sep="")
lambda<-exp(beta.0+x.design%*%matrix(beta.cov,ncol=1))
scale<-(1/lambda)^(1/alpha)
t.true<-rweibull(N,alpha,scale)
T.max<-max(qweibull(0.9,alpha,median(scale)))
t.left<-NULL
t.right<-NULL
for (i in 1:N) {
tt<-unique(c(0,sort(runif(grid,0,T.max)),T.max))
if (t.true>=T.max) {
x.left<-T.max
x.right<-NA
} else {
x.left<-max(tt[t.true>tt])
x.right<-min(tt[t.true<tt])
}
t.left<-c(t.left,x.left)
t.right<-c(t.right,x.right)
}
return(data.frame(ID=1:N,left=t.left,right=t.right,x.design))
}

## End(Not run)