library(survival)
## --------------------------------------------------------------
## 1 censored-covariate and 2 non-censored covariates
## no censoring, to compare result with survival::survreg
## modify prop.cens to introduce censoring of covariate
## --------------------------------------------------------------
set.seed(158)
size <- 100
lambda <- exp(-2)
gamma <- 1.5
## vector of regression parameters: the last entry is the one for the censored covariate
ParamBeta <- c(0.3, -0.2, 0.25)
true <- c(lambda, gamma, ParamBeta)
## simulate from a Weibull regression model
Censdata <- rnorm(size, mean = 5, sd = 0.5)
min.value <- min(Censdata)
prop.cens <- 0
LOD <- qnorm(prop.cens, mean = 5, sd = 0.5)
Censdata[Censdata <= LOD] <- LOD
index.noncens <- which(Censdata > LOD)
InfoCens <- matrix(nrow = size, ncol = 1)
InfoCens[index.noncens] <- Censdata[index.noncens]
CovariateCens <- cbind(InfoCens, Censdata)
NonCensdata <- data.frame("var1" = rnorm(size, mean = 4, sd = 0.5),
"var2" = rnorm(size, mean = 4, sd = 0.5))
index.cens <- which(Censdata <= LOD)
index.cens <- index.cens[index.cens >= 250]
InfoCens[index.cens] <- min.value
time <- TimeSampleWeibull(covariate_noncens = NonCensdata, covariate_cens = Censdata,
lambda = lambda, gamma = gamma, beta=ParamBeta)
IndicTime <- matrix(1, nrow = size, ncol = 1)
TimeCensVector <- rweibull(size, 1.5, exp(3))
IndicTime[time >= TimeCensVector] <- 0
index.cens.time <- which(IndicTime >= TimeCensVector)
time[index.cens.time] <- TimeCensVector[index.cens.time]
## specify the true density for the censored covariate:
DensityCens <- function(value){return(dnorm(value, mean = 5 , sd = 0.5))}
## use Weibull regression where each censored covariate value is set
## to LOD ("naive" method)
naive <- survreg(Surv(time, IndicTime) ~ NonCensdata[,1] +
NonCensdata[,2] + Censdata, dist = "weibull")
initial <- as.vector(ConvertWeibull(naive)$vars[, 1])
## use new method that takes into account the left-censoring of one covariate
cens1 <- SurvRegCens(time = time, event = IndicTime, CovariateNonCens = NonCensdata,
CovariateCens = CovariateCens, Density = DensityCens,
initial = initial, namCens = "biomarker")
## compare estimates
tab <- data.frame(cbind(true, initial, cens1$coeff[, 1]))
colnames(tab) <- c("true", "naive", "Weibull MLE")
tab
## compare confidence intervals
ConvertWeibull(naive)$HR[, 2:3]
cens1$coeff[, 7:8]
## --------------------------------------------------------------
## model without the non-censored covariates
## --------------------------------------------------------------
naive2 <- survreg(Surv(time, IndicTime) ~ Censdata, dist = "weibull")
initial2 <- as.vector(ConvertWeibull(naive2)$vars[, 1])
## use new method that takes into account the left-censoring of one covariate
cens2 <- SurvRegCens(time = time, event = IndicTime, CovariateNonCens = NULL,
CovariateCens = CovariateCens, Density = DensityCens,
initial = initial2, namCens = "biomarker")
## compare estimates
tab <- data.frame(cbind(true[c(1, 2, 5)], initial2, cens2$coeff[, 1]))
colnames(tab) <- c("true", "naive", "Weibull MLE")
tab
## compare confidence intervals
ConvertWeibull(naive2)$HR[, 2:3]
cens2$coeff[, 7:8]
Run the code above in your browser using DataLab