### simPBC.R ---
#----------------------------------------------------------------------
## Author: Thomas Alexander Gerds
## Created: Nov 1 2022 (18:03)
## Version:
## Last-Updated: Nov 21 2022 (18:10)
## By: Thomas Alexander Gerds
## Update #: 8
#----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
#----------------------------------------------------------------------
##
### Code:
##' This function can be used to simulate data alike the pbc data from the survival package.
##'
##' using lava to synthesize data
##' @title simulating data alike the pbc data
##' @param n Sample size
##' @return The simulated data.
##' @export simPBC
##' @examples
##' library(survival)
##' library(lava)
##' # simulate data alike pbc data
##' set.seed(98)
##' d=simPBC(847)
##' d$protimegrp1 <- d$protimegrp=="10-11"
##' d$protimegrp2 <- d$protimegrp==">11"
##' d$sex <- factor(d$sex,levels=0:1,labels=c("m","f"))
##' sF1 <- survreg(Surv(time,status==1)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4,
##' data=d)
##' coxF1 <- coxph(Surv(time,status==1)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4,
##' data=d)
##' # load real pbc data
##' data(pbc,package="survival")
##' pbc <- na.omit(pbc[,c("time","status","age","sex","stage","bili","protime","trt")])
##' pbc$stage <- factor(pbc$stage)
##' levels(pbc$stage) <- list("1/2"=c(1,2),"3"=3,"4"=4)
##' pbc$logbili <- log(pbc$bili)
##' pbc$logprotime <- log(pbc$protime)
##' pbc$stage3 <- 1*(pbc$stage=="3")
##' pbc$stage4 <- 1*(pbc$stage=="4")
##' pbc$protimegrp <- cut(pbc$protime,c(-Inf,10,11,Inf),labels=c("<=10","10-11",">11"))
##' pbc$protimegrp1 <- pbc$protimegrp=="10-11"
##' pbc$protimegrp2 <- pbc$protimegrp==">11"
##' form1=Surv(time,status==1)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4
##' F1 <- survival::survreg(form1,data=pbc)
##' form2=Surv(time,status==2)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4
##' F2 <- survival::survreg(form1,data=pbc)
##' sF2 <- survreg(Surv(time,status==2)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4,
##' data=d)
##' G <- survreg(Surv(time,status==0)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4,
##' data=pbc)
##' sG <- survreg(Surv(time,status==0)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4,
##' data=d)
##' # compare fits in real and simulated pbc data
##' cbind(coef(F1),coef(sF1))
##' cbind(coef(F2),coef(sF2))
##' cbind(coef(G),coef(sG))
##' cbind(coef(glm(protimegrp1~age+sex+logbili,data=pbc,family="binomial")),
##' coef(glm(protimegrp1~age+sex+logbili,data=d,family="binomial")))
##' cbind(coef(lm(logbili~age+sex,data=pbc)),coef(lm(logbili~age+sex,data=d)))
##' @author Thomas A. Gerds <tag@@biostat.ku.dk>
simPBC <- function(n){
requireNamespace("survival")
requireNamespace("lava")
## library(survival)
## library(prodlim)
## library(lava)
utils::data(pbc,package = "survival")
pbc$stage <- factor(pbc$stage)
levels(pbc$stage) <- list("1/2"=c(1,2),"3"=3,"4"=4)
pbc <- na.omit(pbc[,c("time","status","age","sex","stage","bili","protime","trt")])
pbc$logbili <- log(pbc$bili)
pbc$logprotime <- log(pbc$protime)
pbc$stage3 <- 1*(pbc$stage=="3")
pbc$stage4 <- 1*(pbc$stage=="4")
pbc$protimegrp <- cut(pbc$protime,c(-Inf,10,11,Inf),labels=c("<=10","10-11",">11"))
pbc$protimegrp1 <- pbc$protimegrp=="10-11"
pbc$protimegrp2 <- pbc$protimegrp==">11"
## f1 <- coxph(Surv(time,status==1)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4 ,data=pbc)
## f2 <- coxph(Surv(time,status==2)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4 ,data=pbc)
## g <- coxph(Surv(time,status==0)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4 ,data=pbc)
## parametric survival
F1 <- survival::survreg(Surv(time,status==1)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4,data=pbc)
## F1a <- survival::survreg(Surv(time,status==1)~sex+age+logbili+protimegrp+stage,data=pbc)
F2 <- survival::survreg(Surv(time,status==2)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4,data=pbc)
G <- survreg(Surv(time,status==0)~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4,data=pbc)
## G <- survival::survreg(Surv(time,status==0)~1,data=pbc)
## logistic lava::regression
M <- lava::lvm()
lava::distribution(M,~sex) <- lava::binomial.lvm(p=mean(pbc$sex=="f"))
lava::distribution(M,~age) <- lava::normal.lvm(mean=50,sd=10)
lava::distribution(M,~trt) <- lava::binomial.lvm(p=.5)
lava::distribution(M,~logbili) <- lava::normal.lvm(mean=mean(pbc$logbili),sd=sd(pbc$logbili))
lava::distribution(M,~logprotime) <- lava::normal.lvm(mean=mean(pbc$logprotime),sd=sd(pbc$logprotime))
lava::distribution(M,~time.transplant) <- lava::coxWeibull.lvm(scale=exp(-F1$coef["(Intercept)"]/F1$scale),shape=1/F1$scale)
lava::distribution(M,~time.death) <- lava::coxWeibull.lvm(scale=exp(-F2$coef["(Intercept)"]/F2$scale),shape=1/F2$scale)
lava::distribution(M,~time.cens) <- lava::coxWeibull.lvm(scale=exp(-G$coef["(Intercept)"]/G$scale),shape=1/G$scale)
## bili coefficients
lava::regression(M,logbili~age+sex) <- coef(lm(logbili~age+sex,data=pbc))[-1]
## protime coefficients
lava::regression(M,protimegrp1~age+sex+logbili) <- coef(glm(protimegrp1~age+sex+logbili,data=pbc,family="binomial"))[-1]
lava::regression(M,protimegrp2~age+sex+logbili) <- coef(glm(protimegrp2~age+sex+logbili,data=pbc,family="binomial"))[-1]
## stage coefficients
lava::regression(M,stage3~age+sex+protimegrp1+protimegrp2+logbili) <- coef(glm(stage3~age+sex+protimegrp1+protimegrp2+logbili,data=pbc,family="binomial"))[-1]
lava::regression(M,stage4~age+sex+protimegrp1+protimegrp2+logbili) <- coef(glm(stage4~age+sex+protimegrp1+protimegrp2+logbili,data=pbc,family="binomial"))[-1]
## survival coefficients
lava::regression(M,time.transplant~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4) <- -coef(F1)[-1]/F1$scale
lava::regression(M,time.death~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4) <- -coef(F2)[-1]/F2$scale
lava::regression(M,time.cens~sex+age+logbili+protimegrp1+protimegrp2+stage3+stage4) <- -coef(G)[-1]/G$scale
M <- lava::categorical(M,~stage,labels=c("1/2","3","4"),K=3,p=c(mean(pbc$stage=="3"),mean(pbc$stage=="4")))
lava::transform(M,stage3~stage) <- function(x){1*(x==1)}
lava::transform(M,stage4~stage) <- function(x){1*(x==2)}
lava::transform(M,protime~logprotime) <- function(x){exp(x)}
lava::transform(M,protimegrp~protime) <- function(x){cut(x[,"protime"], c(-Inf,10,11,Inf), labels=c("<=10","10-11",">11"))}
lava::transform(M,protimegrp1~protimegrp) <- function(x){1*(x=="10-11")}
lava::transform(M,protimegrp2~protimegrp) <- function(x){1*(x==">11")}
M <- lava::eventTime(M,time~min(time.cens=0,time.transplant=1,time.death=2),"status")
d <- lava::sim(M,n)
d
}
######################################################################
### simPBC.R ends here
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.