如何使用 anova() 查找虚拟交互中主要效果的总体意义
How to find overall significance for main effects in a dummy interaction using anova()
I 运行 具有两个分类变量（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
您为 "overall effect" 构建的两个模型似乎确实不满足分层（即正确嵌套）的统计 属性。具体来说，如果您查看使用该代码构建的实际模型，您应该会发现它们实际上是同一个模型，但双向交叉效应的标签不同。在这两种情况下，您都有 15 个估计系数（因此自由度差异为零），并且您不会认为完整模型中的 x1 参数与 "reduced" 中的 x2[-3.2532,-0.6843):x1[-0.6973,-0.0347) 参数具有相同的系数] 寻找 x1 效应的模型，即 0.19729。交叉运算符基本上是用交互结果填充主效应的所有缺失单元格。
如果您想保持在普遍接受的统计实践范围内，那么在没有所有主要影响的情况下查看交互模型确实没有什么价值。
如果您键入：
fit_full
...您应该得到模型的摘要，该模型具有 x1 水平、x2 水平和交互水平的 p 值。因为您选择按四个任意分割点对它们进行分类，所以您最终将得到总共 15 个参数估计值。相反，如果您不进行任何削减并对线性效应和线性之间的交互作用建模，则可以直接获得三个 p 值。我猜有人怀疑效果不是线性的，如果是这样的话，我认为三次样条模型可能更简约，并且比离散化成 4 个不相交的水平更能扭曲生物现实。如果您认为效果可能是非线性的但是有序的，则可以使用因子分类变量的有序版本，但结果通常会让外行感到困惑。
42- 的回答很有参考意义，但在阅读之后我仍然不知道如何确定三个 p 值，或者这是否可能。因此，我与我大学的生物统计学教授进行了交谈。他的回答很简单，我分享一下，以防其他人有类似的问题。
在此设计中，无法确定 x1、x2 及其相互作用的整体效果的三个 p 值。如果我们想知道三个整体效应的 p 值，我们需要保持连续变量不变。但是将变量分组回答了一个不同的问题，因此无论我们使用哪种统计模型，我们都无法检验整体效果的假设。
I 运行 具有两个分类变量（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
您为 "overall effect" 构建的两个模型似乎确实不满足分层（即正确嵌套）的统计 属性。具体来说，如果您查看使用该代码构建的实际模型，您应该会发现它们实际上是同一个模型，但双向交叉效应的标签不同。在这两种情况下，您都有 15 个估计系数（因此自由度差异为零），并且您不会认为完整模型中的 x1 参数与 "reduced" 中的 x2[-3.2532,-0.6843):x1[-0.6973,-0.0347) 参数具有相同的系数] 寻找 x1 效应的模型，即 0.19729。交叉运算符基本上是用交互结果填充主效应的所有缺失单元格。
如果您想保持在普遍接受的统计实践范围内，那么在没有所有主要影响的情况下查看交互模型确实没有什么价值。
如果您键入：
fit_full
...您应该得到模型的摘要，该模型具有 x1 水平、x2 水平和交互水平的 p 值。因为您选择按四个任意分割点对它们进行分类，所以您最终将得到总共 15 个参数估计值。相反，如果您不进行任何削减并对线性效应和线性之间的交互作用建模，则可以直接获得三个 p 值。我猜有人怀疑效果不是线性的，如果是这样的话，我认为三次样条模型可能更简约，并且比离散化成 4 个不相交的水平更能扭曲生物现实。如果您认为效果可能是非线性的但是有序的，则可以使用因子分类变量的有序版本，但结果通常会让外行感到困惑。
42- 的回答很有参考意义，但在阅读之后我仍然不知道如何确定三个 p 值，或者这是否可能。因此，我与我大学的生物统计学教授进行了交谈。他的回答很简单，我分享一下，以防其他人有类似的问题。
在此设计中，无法确定 x1、x2 及其相互作用的整体效果的三个 p 值。如果我们想知道三个整体效应的 p 值，我们需要保持连续变量不变。但是将变量分组回答了一个不同的问题，因此无论我们使用哪种统计模型，我们都无法检验整体效果的假设。