0

我正在对生存数据执行审查分位数回归，以解释我的队列中某些百分位数的生存时间差异。我正在使用quantregCRAN 的软件包。

这是一个可作为可重现示例的代码

id<-seq(1,2000, by=1)
set.seed(123)
months <- runif(2000, min=0, max=36) #Create a time variable
set.seed(123)
event <- sample(c(0,1), 2000, replace=TRUE, prob=c(0.9,0.2)) #Create a censor/event variable
set.seed(123)
sex <- sample(c(0,1), 2000, replace=TRUE, prob=c(0.55,0.45)) #Create a sex variable
df <- as.data.frame(cbind(id,months,event,sex))

#Generate some random NA to simulate real-world dataset
set.seed(123)
navector1 <- sample(seq(1,2000, by=1), 120, replace=FALSE) 
set.seed(123)
navector2 <- sample(seq(1,2000, by=1), 150, replace=FALSE)
set.seed(123)
navector3 <- sample(seq(1,2000, by=1), 200, replace=FALSE)
df$months <- replace(df$months, navector1, NA)
df$event <- replace(df$event, navector2, NA)
df$sex <- replace(df$sex, navector3, NA)

require(survival)
require(quantreg)
fit <- crq(Surv(months,event) ~ sex, method="Portnoy", data=df)
fit
summary(fit)

请注意，在我的df数据集中，您有一个时间变量 ( months)、一个事件变量 ( event) 和一个分组变量 ( sex)。我随机添加了一些 NA 值来模拟观测的真实数据。

如您所见，fit存储crq函数的结果。结果如下：

Call:
crq(formula = Surv(months, event) ~ sex, data = df, method = "Portnoy")

Coefficients:
            tau= 0.2 tau= 0.4 tau= 0.6 tau= 0.8
(Intercept)       NA       NA       NA       NA
sex               NA       NA       NA       NA

这些NA可能与事件在我的数据库中非常罕见有关，因此两组的总生存率可能高于 20%（即 0.2 个百分点）。但是，我对较低的百分位数感兴趣（应该报告非 NA 值）。当我运行summary(fit)命令时，结果如下：

Error in B$A[1, , , drop = TRUE] : incorrect number of dimensions

我没有尝试过任何可以帮助我的代码工作的方法。请注意，我尝试在 crq 的调用中指定tau,taus和grid属性，但这些似乎不会影响输出（它一直在发布“NA”表）。

谁能帮我？

4

0 回答 0