我使用两个分类变量(x1和x2)及其相互作用来进行Cox回归。我需要知道x1,x2以及交互作用的总体意义。
互动的整体效果:
我知道如何使用anova()
找出交互的整体效果:
library(survival)
fit_x1_x2 <- coxph(Surv(time, death) ~ x1 + x2 , data= df)
fit_full <- coxph(Surv(time, death) ~ x1 + x2 + x1:x2, data= df)
anova(fit_x1_x2, fit_full)
但是我们应该如何使用anova()
来找出x1或x2的整体效果?我试过的是:
x1的整体效果
fit_x2_ia <- coxph(Surv(time, death) ~ x2 + x1:x2, data= df)
fit_full <- coxph(Surv(time, death) ~ x1 + x2 + x1:x2, data= df)
anova(fit_x2_ia, fit_full)
x2的整体效果
fit_x1_ia <- coxph(Surv(time, death) ~ x1 + x1:x2, data= df)
fit_full <- coxph(Surv(time, death) ~ x1 + x2 + x1:x2, data= df)
anova(fit_x1_ia, fit_full)
我不确定这是否应该使用anova()
。输出显示自由度为零的事实使我感到怀疑。对于x1和x2的整体效果,这两次都使我感到很困惑,尽管模型的对数似然值相同且Chi值为零,但检验是显着的。
这是我使用的数据
set.seed(1) # make it reproducible
df <- data.frame(x1= rnorm(1000), x2= rnorm(1000)) # generate data
df$death <- rbinom(1000,1, 1/(1+exp(-(1 + 2 * df$x1 + 3 * df$x2 + df$x1 * df$x2)))) # dead or not
library(tidyverse) # for cut_number() function
df$x1 <- cut_number(df$x1, 4); df$x2 <- cut_number(df$x2, 4) # make predictors to groups
df$time <- rnorm(1000); df$time[df$time<0] <- -df$time[df$time<0] # add survival times
您为“总体效果”构建的两个模型实际上似乎并不满足分层(即正确嵌套)的统计属性。具体来说,如果您查看使用该代码构建的实际模型,则应该看到它们实际上是相同的模型,但双向交叉效果的标签不同。在这两种情况下,您都有15个估计系数(因此自由度差为零),并且您不会认为完整模型中的x1
参数与寻找x1的“精简”模型中的x2[-3.2532,-0.6843):x1[-0.6973,-0.0347)
参数具有相同的系数-效果,即0.19729。交叉算子基本上是在所有缺少的像元中填充具有交互结果的主要效果。