# NOT RUN {
# replicate simulation in Liang et al.
# }
# NOT RUN {
library(data.table)
library(survival)
datasimi <- function(id){
X1 <- runif(1,0,1)
X2 <- rbinom(1,1,0.5)
Z <- rgamma(1,1,1)
Z1 <- rnorm(1,Z-1,1)
gamma <- c(0.5,-0.5)
beta <- c(1,-1)
hazard <- Z*exp(X1/2 - X2/2)
C <- runif(1,0,5.8)
t <- 0
tlast <- t
y <- t + X1-X2 + Z1*X2 + rnorm(1,0,1)
wait <- rexp(1,hazard)
while(tlast+wait<C){
tnew <- tlast+wait
y <- c(y,tnew + X1-X2 + Z1*X2 + rnorm(1,0,1))
t <- c(t,tnew)
tlast <- tnew
wait <- rexp(1,hazard)
}
datai <- list(id=rep(id,length(t)),t=t,y=y,
X1=rep(X1,length(t)),X2=rep(X2,length(t)),C=rep(C,length(t)))
return(datai)
}
sim1 <- function(it,nsubj){
data <- lapply(1:nsubj,datasimi)
data <- as.data.frame(rbindlist(data))
data$event <- 1
C <- tapply(data$C,data$id,mean)
tapply(data$C,data$id,sd)
maxfu <- cbind(1:nsubj,C)
maxfu <- as.data.frame(maxfu)
res <- Liang(data=data, id="id",time="t",Yname="y",
Xnames=c("X1","X2"),
Wnames=c("X2"),Znames=c("X1","X2"), formulaobs=Surv(t.lag,t,event)~X1
+ X2+ frailty(id),invariant=c
("id","X1","X2"),lagvars="t",lagfirst=NA,maxfu=maxfu,
baseline=1)
return(res)
}
# change n to 500 to replicate results of Liang et al.
n <- 10
s <- lapply(1:n,sim1,nsubj=200)
smat <- matrix(unlist(s),byrow=TRUE,ncol=4)
apply(smat,2,mean)
# }
Run the code above in your browser using DataLab