用 R 中的交互拟合 Fine-Gray 模型（包 RiskRegression - 函数 FGR）
Fitting a Fine-Gray model with interaction in R (package RiskRegression - function FGR)
我的数据集与 RiskRegression 包中包含的数据集 'Melanoma' 非常相似：3307 名患者，502 名感兴趣的事件（骨折），264 名死亡（竞争风险）。 time 是骨骼检查 (DXA) 后的年数，status 以这种方式编码 O=censored,1=fracture,2=death) .
我正在尝试使用交互拟合 Fine-Gray 模型，但是当我以 var1 * var2) 的形式引入交互项时，我收到一条错误消息：
« 设计错误[pos, , drop = FALSE] : 下标越界 » .
这是我的代码：
fgr<-FGR(Hist(time,status)~age+htot_bmd+tot_bmd+amof+PR+atcdtfam+AlcFR+PR+BMI3C+malchronFR+malchronFR*BMI3C+atcdtfam*PR,data=df2,cause=1)
我尝试了Zhongheng等人论文中提供的代码。 “竞争风险数据的模型验证”与数据集 'Melanoma' 引入交互但出现相同的错误消息。
是否可以引入与 FGR 的相互作用以及如何进行？
谢谢
您可以按如下方式使用 model.matrix 函数。 crr()函数可以做交互。
> library(riskRegression)
> library(survival)
> library(prodlim)
> library(cmprsk)
> data(Melanoma)
> Melanoma$id<-1:nrow(Melanoma)
> set.seed(123)
> ind.split<-sample(1:nrow(Melanoma),
+ round(nrow(Melanoma)*4/5),
+ replace = F)
> dftrain<-Melanoma[ind.split,]
> dftest<-Melanoma[-ind.split,]
> fgr.full<-FGR(Hist(time,status)~age+thick+ici+
+ epicel+ulcer+sex+invasion,
+ data=dftrain,cause=1)
> modMatrix <- model.matrix(~thick+ici+
+ epicel+ulcer*age+invasion,dftrain)[,-1]
>
> fgrMod <- crr(ftime = dftrain$time,
+ fstatus = dftrain$status,
+ cov1 = modMatrix,failcode=2)
> summary(fgrMod)
Competing Risks Regression
Call:
crr(ftime = dftrain$time, fstatus = dftrain$status, cov1 = modMatrix,
failcode = 2)
coef exp(coef) se(coef) z p-value
thick 0.1194 1.127 0.1292 0.924 0.3600
ici1 -0.7607 0.467 1.0721 -0.710 0.4800
ici2 -0.8531 0.426 0.9379 -0.910 0.3600
ici3 -0.1924 0.825 1.0895 -0.177 0.8600
epicelpresent 0.8973 2.453 0.8434 1.064 0.2900
ulcerpresent -0.7101 0.492 1.9776 -0.359 0.7200
age 0.0627 1.065 0.0227 2.766 0.0057
invasionlevel.1 -1.2031 0.300 0.7068 -1.702 0.0890
invasionlevel.2 -2.0365 0.130 1.4121 -1.442 0.1500
ulcerpresent:age 0.0152 1.015 0.0320 0.473 0.6400
exp(coef) exp(-coef) 2.5% 97.5%
thick 1.127 0.887 0.87475 1.45
ici1 0.467 2.140 0.05716 3.82
ici2 0.426 2.347 0.06780 2.68
ici3 0.825 1.212 0.09752 6.98
epicelpresent 2.453 0.408 0.46968 12.81
ulcerpresent 0.492 2.034 0.01019 23.71
age 1.065 0.939 1.01844 1.11
invasionlevel.1 0.300 3.330 0.07515 1.20
invasionlevel.2 0.130 7.664 0.00819 2.08
ulcerpresent:age 1.015 0.985 0.95348 1.08
Num. cases = 164
Pseudo Log-likelihood = -52.3
Pseudo likelihood ratio test = 21.1 on 10 df,
那你可以试试下面的代码：
library(riskRegression)
library(survival)
library(prodlim)
library(cmprsk)
data(Melanoma)
Melanoma$id<-1:nrow(Melanoma)
set.seed(123)
ind.split<-sample(1:nrow(Melanoma),
round(nrow(Melanoma)*4/5),
replace = F)
dftrain<-Melanoma[ind.split,]
dftest<-Melanoma[-ind.split,]
fgr.NoInteraction<-FGR(Hist(time,status)~age+thick+ici+
epicel+ulcer+sex+invasion,
data=dftrain,cause=1)
modMatrix <- model.matrix(~thick+ici+
epicel+ulcer*age+invasion,dftrain)[,-1]
dtInteraction <- cbind(data.frame( modMatrix),status=dftrain$status,
time=dftrain$time)
fgr.Interaction<- FGR(as.formula(paste("Hist(time,status)~",paste(names(dtInteraction[1:9]),collapse = "+"))),
data = dtInteraction,cause = 1)
score.cv<-riskRegression::Score(list("Fine-Gray"= fgr.Interaction),
formula = Hist(time,status)~1,
data=dtInteraction,times = sort(unique(dtInteraction$time)),
cens.method="jackknife",
se.fit=1L,plots="calibration")
plotCalibration(score.cv,times = 3330,cens.method="local")
我们还需要添加与 model.matrix 的交互项。但是，我们只能使用 FGR 中的对象作为 Score 函数的输入。 the paper中的其他数字也可以用类似的技巧完成。
您可以使用以下代码处理您的数据：
> library(riskRegression)
> library(survival)
> library(prodlim)
> library(cmprsk)
> library(readxl)
> df2 <- read_xlsx("/Users/zhang/Downloads/df2.xlsx")
New names:
* `` -> ...1
> df2
# A tibble: 300 x 14
...1 neck_bmd htot_bmd tot_bmd age AlcFR PR atcdtfam malchronFR amof BMI3C time event
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 1 0.960 0.953 1.04 79.1 0 0 0 0 2 3 9.00 Cen
2 2 0.612 0.620 0.988 79.2 0 0 0 0 0 3 4.76 MOF
3 3 0.880 0.990 0.827 78.6 0 0 0 1 1 2 9.14 Cen
4 4 0.869 0.905 0.866 79.0 0 0 0 0 0 2 9.11 Cen
5 5 0.863 0.991 1.17 79.0 1 0 1 0 0 2 10.2 Cen
6 6 0.722 0.902 0.842 78.8 0 0 0 0 0 2 9.09 Cen
7 7 0.853 0.929 1.33 76.9 0 0 0 0 0 3 10.1 Cen
8 8 0.830 0.912 0.947 77.0 0 0 0 1 0 2 8.13 Cen
9 9 0.872 0.968 1.22 77.2 1 0 0 0 0 2 8.12 Cen
10 10 0.639 0.776 0.822 76.7 0 0 0 1 0 2 8.12 Cen
# … with 290 more rows, and 1 more variable: status <dbl>
> modMatrix <- model.matrix(~age+htot_bmd+tot_bmd+amof+atcdtfam+
+ AlcFR+PR+BMI3C*malchronFR+neck_bmd,df2)[,-1]
> dtInteraction <- cbind(data.frame( modMatrix),
+ status=df2$status, time=df2$time)
> fgr.Interaction<- FGR(as.formula(paste("Hist(time,status)~",
+ paste(names(dtInteraction[1:11]),collapse = "+"))),
+ data = dtInteraction,cause = 1)
> score.cv<-riskRegression::Score(list("Fine-Gray"= fgr.Interaction),
+ formula = Hist(time,status)~1,
+ data=dtInteraction,times = sort(unique(dtInteraction$time))[25:200],
+ cens.method="jackknife",
+ se.fit=1L,plots="calibration")
> plotCalibration(score.cv,times = df2$time[11],
+ cens.method="local")
我的数据集与 RiskRegression 包中包含的数据集 'Melanoma' 非常相似：3307 名患者，502 名感兴趣的事件（骨折），264 名死亡（竞争风险）。 time 是骨骼检查 (DXA) 后的年数，status 以这种方式编码 O=censored,1=fracture,2=death) . 我正在尝试使用交互拟合 Fine-Gray 模型，但是当我以 var1 * var2) 的形式引入交互项时，我收到一条错误消息：
« 设计错误[pos, , drop = FALSE] : 下标越界 » .
这是我的代码：
fgr<-FGR(Hist(time,status)~age+htot_bmd+tot_bmd+amof+PR+atcdtfam+AlcFR+PR+BMI3C+malchronFR+malchronFR*BMI3C+atcdtfam*PR,data=df2,cause=1)
我尝试了Zhongheng等人论文中提供的代码。 “竞争风险数据的模型验证”与数据集 'Melanoma' 引入交互但出现相同的错误消息。
是否可以引入与 FGR 的相互作用以及如何进行？
谢谢
您可以按如下方式使用 model.matrix 函数。 crr()函数可以做交互。
> library(riskRegression)
> library(survival)
> library(prodlim)
> library(cmprsk)
> data(Melanoma)
> Melanoma$id<-1:nrow(Melanoma)
> set.seed(123)
> ind.split<-sample(1:nrow(Melanoma),
+ round(nrow(Melanoma)*4/5),
+ replace = F)
> dftrain<-Melanoma[ind.split,]
> dftest<-Melanoma[-ind.split,]
> fgr.full<-FGR(Hist(time,status)~age+thick+ici+
+ epicel+ulcer+sex+invasion,
+ data=dftrain,cause=1)
> modMatrix <- model.matrix(~thick+ici+
+ epicel+ulcer*age+invasion,dftrain)[,-1]
>
> fgrMod <- crr(ftime = dftrain$time,
+ fstatus = dftrain$status,
+ cov1 = modMatrix,failcode=2)
> summary(fgrMod)
Competing Risks Regression
Call:
crr(ftime = dftrain$time, fstatus = dftrain$status, cov1 = modMatrix,
failcode = 2)
coef exp(coef) se(coef) z p-value
thick 0.1194 1.127 0.1292 0.924 0.3600
ici1 -0.7607 0.467 1.0721 -0.710 0.4800
ici2 -0.8531 0.426 0.9379 -0.910 0.3600
ici3 -0.1924 0.825 1.0895 -0.177 0.8600
epicelpresent 0.8973 2.453 0.8434 1.064 0.2900
ulcerpresent -0.7101 0.492 1.9776 -0.359 0.7200
age 0.0627 1.065 0.0227 2.766 0.0057
invasionlevel.1 -1.2031 0.300 0.7068 -1.702 0.0890
invasionlevel.2 -2.0365 0.130 1.4121 -1.442 0.1500
ulcerpresent:age 0.0152 1.015 0.0320 0.473 0.6400
exp(coef) exp(-coef) 2.5% 97.5%
thick 1.127 0.887 0.87475 1.45
ici1 0.467 2.140 0.05716 3.82
ici2 0.426 2.347 0.06780 2.68
ici3 0.825 1.212 0.09752 6.98
epicelpresent 2.453 0.408 0.46968 12.81
ulcerpresent 0.492 2.034 0.01019 23.71
age 1.065 0.939 1.01844 1.11
invasionlevel.1 0.300 3.330 0.07515 1.20
invasionlevel.2 0.130 7.664 0.00819 2.08
ulcerpresent:age 1.015 0.985 0.95348 1.08
Num. cases = 164
Pseudo Log-likelihood = -52.3
Pseudo likelihood ratio test = 21.1 on 10 df,
那你可以试试下面的代码：
library(riskRegression)
library(survival)
library(prodlim)
library(cmprsk)
data(Melanoma)
Melanoma$id<-1:nrow(Melanoma)
set.seed(123)
ind.split<-sample(1:nrow(Melanoma),
round(nrow(Melanoma)*4/5),
replace = F)
dftrain<-Melanoma[ind.split,]
dftest<-Melanoma[-ind.split,]
fgr.NoInteraction<-FGR(Hist(time,status)~age+thick+ici+
epicel+ulcer+sex+invasion,
data=dftrain,cause=1)
modMatrix <- model.matrix(~thick+ici+
epicel+ulcer*age+invasion,dftrain)[,-1]
dtInteraction <- cbind(data.frame( modMatrix),status=dftrain$status,
time=dftrain$time)
fgr.Interaction<- FGR(as.formula(paste("Hist(time,status)~",paste(names(dtInteraction[1:9]),collapse = "+"))),
data = dtInteraction,cause = 1)
score.cv<-riskRegression::Score(list("Fine-Gray"= fgr.Interaction),
formula = Hist(time,status)~1,
data=dtInteraction,times = sort(unique(dtInteraction$time)),
cens.method="jackknife",
se.fit=1L,plots="calibration")
plotCalibration(score.cv,times = 3330,cens.method="local")
我们还需要添加与 model.matrix 的交互项。但是，我们只能使用 FGR 中的对象作为 Score 函数的输入。 the paper中的其他数字也可以用类似的技巧完成。
您可以使用以下代码处理您的数据：
> library(riskRegression)
> library(survival)
> library(prodlim)
> library(cmprsk)
> library(readxl)
> df2 <- read_xlsx("/Users/zhang/Downloads/df2.xlsx")
New names:
* `` -> ...1
> df2
# A tibble: 300 x 14
...1 neck_bmd htot_bmd tot_bmd age AlcFR PR atcdtfam malchronFR amof BMI3C time event
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 1 0.960 0.953 1.04 79.1 0 0 0 0 2 3 9.00 Cen
2 2 0.612 0.620 0.988 79.2 0 0 0 0 0 3 4.76 MOF
3 3 0.880 0.990 0.827 78.6 0 0 0 1 1 2 9.14 Cen
4 4 0.869 0.905 0.866 79.0 0 0 0 0 0 2 9.11 Cen
5 5 0.863 0.991 1.17 79.0 1 0 1 0 0 2 10.2 Cen
6 6 0.722 0.902 0.842 78.8 0 0 0 0 0 2 9.09 Cen
7 7 0.853 0.929 1.33 76.9 0 0 0 0 0 3 10.1 Cen
8 8 0.830 0.912 0.947 77.0 0 0 0 1 0 2 8.13 Cen
9 9 0.872 0.968 1.22 77.2 1 0 0 0 0 2 8.12 Cen
10 10 0.639 0.776 0.822 76.7 0 0 0 1 0 2 8.12 Cen
# … with 290 more rows, and 1 more variable: status <dbl>
> modMatrix <- model.matrix(~age+htot_bmd+tot_bmd+amof+atcdtfam+
+ AlcFR+PR+BMI3C*malchronFR+neck_bmd,df2)[,-1]
> dtInteraction <- cbind(data.frame( modMatrix),
+ status=df2$status, time=df2$time)
> fgr.Interaction<- FGR(as.formula(paste("Hist(time,status)~",
+ paste(names(dtInteraction[1:11]),collapse = "+"))),
+ data = dtInteraction,cause = 1)
> score.cv<-riskRegression::Score(list("Fine-Gray"= fgr.Interaction),
+ formula = Hist(time,status)~1,
+ data=dtInteraction,times = sort(unique(dtInteraction$time))[25:200],
+ cens.method="jackknife",
+ se.fit=1L,plots="calibration")
> plotCalibration(score.cv,times = df2$time[11],
+ cens.method="local")