如何利用anova()找到虚拟交互中主要影响的总体意义

Question

我用两个分类变量(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()。产出显示自由度为零这一事实使我感到怀疑。更令我困惑的是，尽管模型的对数似然值相同，Chi值为零，但对于x1和x2的总体效果而言，两者的检验都是显著的。

这里是我使用的数据

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

Answer 1

答案从42-是信息，但阅读后，我仍然不知道如何确定三个p值，或如果这是可能的。于是我和我们大学的生物统计学教授谈了起来。他的回答很简单，如果其他人有类似的问题，我也同意。

在此设计中，不可能确定x1、x2及其相互作用的总体效果的三个p值。如果我们想知道这三个整体效应的p值，我们需要保持连续变量的原样。但是，把变量分解成组来回答一个不同的问题，因此无论我们使用哪种统计模型，我们都无法检验总体效应的假设。

Answer 2

您为“总体效应”构建的两个模型实际上并不满足分层的统计特性，即正确嵌套。具体来说，如果您查看使用该代码构造的实际模型，您应该会发现它们实际上是同一个模型，具有不同的标签，用于双向交叉效果。在这两种情况下，您都有15个估计系数(因此是零自由度差)，而且您将不会认为完整模型中的x1参数与寻找x1效应的“缩减”模型中的x2[-3.2532,-0.6843):x1[-0.6973,-0.0347)参数具有相同的系数，即0.19729。交叉算子基本上是填补所有缺失细胞的主要影响与相互作用的结果。

如果你想保持在普遍接受的统计实践的范围内，那么如果你想要保持在普遍接受的统计实践的范围内，那么在没有所有主要影响的情况下，研究交互模型真的没有什么价值。

如果您键入：

fit_full

..。您应该得到模型的摘要，其中包含x1级别、x2级别和交互级别的p值。因为您选择按四个任意切点对这些数据进行分类，因此最终将得到总共15个参数估计。如果你不做任何切割，并模拟线性效应和线性的线性交互作用，你可以直接得到三个p值。我猜人们怀疑这种效应不是线性的，如果是这样的话，我认为三次样条模型可能更简练，比离散化更能扭曲生物现实，变成4个不相交的层次。如果你认为影响可能是非线性的，但顺序的，有一个因素分类变量的序号版本，但结果通常是混乱的不熟悉的人。