R/bicomprisksim.R In mets: Analysis of Multivariate Event Times

Defines functions sim.F1clustsim.F1F1F1clustilaplapbicomprisksim

```## n <- 1e4; ACE <- c(1,1,1)/3
## logscale <- -4.5; logshape <- .7
## p2 <- .065
## a2 <- -10; b2 <- 0.15
## a1 <- 85; b1 <- 0.1

## pmvn(upper=c(q2,q2),sigma=diag(2)*(1-2/3)+2/3)
## pmvn(c(q2,q2),sigma=diag(2)*(1-2/3)+2/3)

## pnorm(q2,sd=1) ## Marginal / Perfect dependence
## pmvn(upper=c(q2,q2),sigma=diag(2)*(1-2/3)+2/3) ## Concordance
## pnorm(q2,sd=1)^2 ## Independence
## (lambdaR <- pmvn(upper=c(q2,q2),sigma=diag(2)*(1-2/3)+2/3)/pnorm(q2,sd=1)^2)

bicomprisksim <- function(n=1e4,
ACE=c(1/3,1/3,1/3),
logscale=-4.5,logshape=.7,
a1=85,b1=0.1,
a2=-10,b2=0.15,
p2=.065,
tt,
...) {

ACE <- ACE/sum(ACE)
ialpha <- function(v,a,b) -(log(-v)+a)/b
q2 <- qnorm(p2) ## p2: Prostata cancer prevalence
p1 <- 1-p2; q1 <- qnorm(p1) ## Death without cancer
alpha <- function(t,a,b) -exp(-(b*t+a))
F2s <- function(t) pnorm(alpha(t,b=b2,a=a2)+q2) ## Marginal Cumulative Incidence (cancer)
### Random effects
R <- diag(2)*.5+.5
J <- matrix(1,ncol=2,nrow=2)
I <- diag(2)
zyg <- rep(c(0,1),each=n)
A <- rbind(rmvn(n,sigma=R*ACE[1]),rmvn(n,sigma=J*ACE[1]))
C <- rmvn(2*n,sigma=J*ACE[2])
### Random effects 'death'
eta1 <- C
### Random effects 'cancer'
eta2 <- A+C;
### Subject-specific probability of lifetime cancer
probcanc <- pnorm(q2+eta2,sd=ACE[3]^.5)
### Cancer/Death without cancer realizations
cancertrue <- (runif(length(probcanc))<probcanc)*1+1
### Event times given failure and random effect
## inverse P(T<t| eta,cause=1)
iF1 <- function(u,eta,sd,a,b) (qnorm(u,sd=sd)-eta)/b+a
iF2 <- function(u,eta,pr,sd,q,a,b) ialpha(qnorm(u*pr,sd=sd)-eta-q,a=a,b=b)
u <- runif(length(probcanc))
t2 <- iF2(u,eta2,probcanc,sd=sqrt(1-sum(ACE[1:2])),q=q2,a=a2,b=b2)
t1 <- iF1(u,eta1,sd=sqrt(1-sum(ACE[2])),a=a1,b=b1)
t0 <- t1; t0[which(cancertrue==2)] <- t2[which(cancertrue==2)]
### Censoring distributon
cens <- matrix(rep(rweibull(length(probcanc)/2,shape=exp(logshape),scale=1/exp(logscale)),2),
ncol=2)
##hist(cens,xlab=c(0,400),200)
cause <- cancertrue
cause[t0>cens] <- 0
t <- pmin(t0,cens)
(censtab <- table(cause)/length(cause))
### Data.frame
d <- data.frame(t,cause,zyg,cancertrue-1,cens[,1],t0[,1],t0[,2]);
names(d) <- c("time1","time2","cause1","cause2","zyg","cancertrue1","cancertrue2","cens.time","T01","T02")
### Long format
dd <- fast.reshape(d)
dd\$cancer <- (dd\$cause==2)*1

if (missing(tt)) tt <- seq(0,max(dd\$time))
Smz <- J*ACE[1]+J*ACE[2]+I*ACE[3]
Sdz <- R*ACE[1]+J*ACE[2]+I*ACE[3]
rr <- alpha(tt,b=b2,a=a2)+q2
Cmz <- pmvn(upper=cbind(rr,rr),mu=matrix(0,ncol=2,nrow=length(rr)),sigma=Smz)
Cdz <- pmvn(upper=cbind(rr,rr),mu=matrix(0,ncol=2,nrow=length(rr)),sigma=Sdz)

true <- list(p2=p2,
p22mz=pmvn(lower=c(-Inf,-Inf),upper=c(q2,q2),sigma=Smz),
p12mz=pmvn(lower=c(q2,-Inf),upper=c(Inf,q2),sigma=Smz),
p12mz=pmvn(lower=c(q2,-Inf),upper=c(Inf,q2),sigma=Smz),
p11mz=pmvn(lower=c(q2,q2),upper=c(Inf,Inf),sigma=Smz),
p22dz=pmvn(lower=c(-Inf,-Inf),upper=c(q2,q2),sigma=Sdz),
p12dz=pmvn(lower=c(q2,-Inf),upper=c(Inf,q2),sigma=Sdz),
p12dz=pmvn(lower=c(q2,-Inf),upper=c(Inf,q2),sigma=Sdz),
p11dz=pmvn(lower=c(q2,q2),upper=c(Inf,Inf),sigma=Sdz),
time=tt,
F2=F2s(tt),
Cmz=Cmz,
Cdz=Cdz
)
attributes(dd) <- c(attributes(dd),true)
return(dd)
}

##################################################
##################################################
##################################################
##################################################

##################################################
### simulation for gamma distributed cif model
##################################################

lap<-function(theta,t) {
return( (1+t/theta)^(-theta))
}
ilap<-function(theta,t) {
itheta<-1/theta; return((t^(-itheta)-1)/(itheta))
}

F1clust<-function(t,rtheta=1,theta=1,lam0=0.5,beta=0.3,x=0) {
return(1-exp(-rtheta*ilap(theta,exp(-t*lam0-t*x*beta))))
}

###################################################################
return( 1 - exp(-t*lam0-t*x*beta))
}

###F1<-function(t,lam0=0.5,beta=0.3,x=0) # proportional version
###{ return( 1 - exp(-(t*lam0)*exp(x*beta))) }

sim.F1<-function(n,theta=1,lam0=0.5,beta=0.3,crate=2) {
x<-runif(n); tt<-seq(0,1,length=100)
F11x<-F1(1,x=x,beta=beta,lam0=lam0)
cause1<-rbinom(n,1,F11x)

stime<-rep(100,n);
for (i in 1:n)
{
if (cause1[i]==1) {
myhazx<-F1(tt,x=x[i],beta=beta,lam0=lam0)/F11x[i]
stime[i]<-Cpred(cbind(myhazx,tt),runif(1))[1,2]+runif(1,0,0.001)
}
}
ctime<-runif(n)*crate
time<-apply(cbind(ctime,stime),1,min)
status<-(stime<ctime); cause1[status==0]<-0;
data<-data.frame(time=time,status=status,X=x,cause=cause1)
return(data)
}

sim.F1clust<-function(n,theta=1,lam0=0.5,beta=0.3,alpha=0,crate=2,
same.cens=FALSE,fix.cens=FALSE) {
k<-n/2; tt<-seq(0,1,length=100)
rtheta<-rgamma(k,theta,scale=1/theta)
stime<-c();cause1<-c();id<-c();vtheta<-c(); X<-c()
cause1 <- c()

for (i in 1:k) {
x<-runif(2);
X<-c(X,x);
F11x<-F1clust(1,rtheta=rtheta[i],theta=theta,x=x,beta=beta,lam0=lam0)
cause<-rbinom(2,1,F11x);
##cause1<-c(cause1,cause);
id<-c(id,rep(i,2)); vtheta<-c(vtheta,rep(rtheta[i],2))

for (j in 1:2) {
if (cause[j]==1) {
myhazx<-F1clust(tt,x=x[j],rtheta=rtheta[i],theta=theta, beta=beta,lam0=lam0)/F11x[j]
stime<-c(stime,Cpred(cbind(myhazx,tt),runif(1))[1,2]+ runif(1,0,0.001))
cause1 <- c(cause1,1);
} else { stime<-c(stime,runif(1)); cause1 <- c(cause1,2);}
}
}

if (same.cens) ctime <- rep(ctime<-runif(k)*crate,each=2) else ctime<-runif(n)*crate;
if (fix.cens) ctime <- rep(ctime<-rep(crate,each=2))

time<-apply(cbind(ctime,stime),1,min)
status<-(stime<ctime);
cause1[status==0]<-0;
data<-data.frame(time=time,X=X,cause=cause1,stime=stime,ctime=ctime,
id=id,theta=vtheta)
return(data)
}

######################################################
### test ting  #######################################
######################################################

## beta=0; n=200; theta=0.1;
## simdata<-sim.F1clust(n=200,theta=0.1,beta=0,crate=3.0,lam0=0.5)
## table(simdata\$cause)
## simdata\$tv <- rep(1:2,n/2)

## tt <- seq(0,1,by=0.01)
## F1(tt)

## p11<- 1 - (1-F1(tt)) - (1-F1(tt)) + lap(1/theta, 2*ilap(1/theta,1-F1(tt)) )
## plot(tt,F1(tt),type="l")
## lines(tt,p11)
## lines(tt,(F1(tt)^2),col=2)

## beta=0; n=1000; theta=0.1;
## gem <- c()
## for (i in 1:10) {
## simdata<-sim.F1clust(n,theta=theta,beta=beta,crate=3.0,lam0=0.5,same.cens=TRUE)
## table(simdata\$cause)
## simdata\$tv <- rep(1:2,n/2)
## simd <- reshape(simdata,direction="wide",idvar="id",timevar="tv")
## simd\$cancer <- (simd\$cause.1==1 ) + (simd\$cause.2==1 )
## tabd <- table(simd\$cancer)
## casewise <- tabd["2"]/ ( tabd["2"] + tabd["1"]*0.5)
## gem <- c(gem,casewise)
## }
## summary(gem)

## ### plot

## casewiset <- p11/F1(tt)
## plot(tt,casewiset,type="l")

## F1t <- F1(tt)

###################################################################
###################################################################

### increase censoring by crate lower, dependence theta lower

## theta=1.0;
## p11t<- 1 - (1-F1(tt)) - (1-F1(tt)) + lap(theta, 2*ilap(theta,1-F1(tt)) )
## F1t <- F1(tt)
## Gc <- 1 - 2*tt
## Gc[Gc<0] <- 0
## mn11t <- cumsum(diff(p11t)* Gc[-1])
## mndt <- cumsum(diff(F1t)* Gc[-1])
## casewisem <- mn11t[100]/mndt[100]
## casewiset <- p11t/F1(tt)

## gem <- c()
## for (i in 1:10) {
## print(i)
## beta=0; n=10000;
## theta=theta;
## simdata<-sim.F1clust(n,theta=theta,beta=beta,crate=0.5,lam0=0.5,same.cens=TRUE,fix.cens=FALSE)
## table(simdata\$cause)
## simdata\$tv <- rep(1:2,n/2)
## simd <- reshape(simdata,direction="wide",idvar="id",timevar="tv")
## simd\$cancer <- (simd\$cause.1==1 ) + (simd\$cause.2==1 )
## tabd <- table(simd\$cancer)
## casewise <- tabd["2"]/ ( tabd["2"] + tabd["1"]*0.5)
## ###
## outm <- prodlim(Event(time,cause)~+1,data=simdata)
## par(mfrow=c(1,3))
## plot(outm)
## lines(tt,F1t,col=3,lwd=3)
## legend("topleft",c("cif","sand"),lty=1,col=c(1,3))
## ###
## cc <- bicomprisk(Event(time,cause)~+1+id(id),data=simdata,cause=c(1,1),prodlim=TRUE)
## p11t<- 1 - (1-F1(tt)) - (1-F1(tt)) + lap(theta, 2*ilap(theta,1-F1(tt)) )
## plot(cc)
## lines(tt,p11t,col=3,lwd=3)
## lines(tt,F1t^2,col=4,lwd=1)
## legend("topleft",c("conc","sand","indep"),lty=1,col=c(1,3,4))
## ###
## case <- casewise(cc,outm,cause.marg=1)
## ###
## casewiset <- p11t/F1t
## plot(case\$casewise[,1],case\$casewise[,2],type="l")
## lines(tt,casewiset,col=3,lwd=3)
## abline(h=casewise,col=2)
## legend("topleft",c("est","sand","constant"),lty=1,col=c(1,3,2))
## ###
## Gc <- 1 - 2*tt
## Gc[Gc<0] <- 0
## ###
## abline(h=casewisem,col=4)
## gem <- c(gem,casewise)
## }

## abline(h=casewisem,col=4)
## casewisem
## mean(gem)

## pdf("caswise-vs-surv.pdf")
## par(mfrow=c(1,1))

## dev.off()

############################################################################
############################################################################
############################################################################

## library(doMC)
## registerDoMC()
## library(mets)

## onesim <- function(i,n,theta=2,beta=0) {
## 	theta=2; beta=0;
## 	n=400

##         simdata<-sim.F1clust(n,theta=theta,beta=beta,crate=0.5,lam0=0.5,same.cens=TRUE,fix.cens=FALSE)
##         simdata\$tv <- rep(1:2,n/2)
##         simd <- reshape(simdata,direction="wide",idvar="id",timevar="tv")
##         simd\$cancer <- (simd\$cause.1==1 ) + (simd\$cause.2==1 )
##         tabd <- table(simd\$cancer)
##         casewise <- tabd["2"]/ ( tabd["2"] + tabd["1"]*0.5)
## ###
##         outm <- prodlim(Hist(time,cause)~+1,data=simdata)
##         par(mfrow=c(1,3))
##         plot(outm)
##         lines(tt,F1t,col=3,lwd=3)
##         legend("topleft",c("cif","sand"),lty=1,col=c(1,3))
## ###
##         cc <- bicomprisk(Event(time,cause)~+1+id(id),data=simdata,cause=c(1,1),prodlim=TRUE)
##         p11t<- 1 - (1-F1(tt)) - (1-F1(tt)) + lap(theta, 2*ilap(theta,1-F1(tt)) )
##         case <- casewise(cc,outm,cause.marg=1)
##         case\$concordance
##         case\$casewise
##         simdata\$time

##         times <- seq(0,0.4,by=0.01)
##         outm <-comp.risk(Surv(time,cause==0)~+1+cluster(id),data=simdata,simdata\$cause,
## 		                         causeS=1,times=times,max.clust=200,n.sim=0)
##         cifmz <-predict(outm,X=1,uniform=0,resample.iid=1)
##         cc <-bicomprisk(Event(time,status)~+1+id(id),data=simdata,
##                         cause=c(1,1),se.clusters=outm\$clusters)
##         cdz <- cc\$model\$"DZ"
##         cmz <- cc\$model\$"MZ"

##         cdz <- casewise.test(cdz,cifmz,test="case") ## test based on casewise
##         cmz <- casewise.test(cmz,cifmz,test="conc") ## based on concordance:w

##         out <- comp.risk(Event(time,cause)~+1,data=simdata,causeS=1)
##         cc <- bicomprisk(Event(time,cause)~+1+id(id),data=simdata,cause=c(1,1))

##     }
```

Try the mets package in your browser

Any scripts or data that you put into this service are public.

mets documentation built on May 31, 2017, 1:52 a.m.