为 gbm 交互强度实现零分布

Posted 2023-03-12

【中文标题】为 gbm 交互强度实现零分布

【英文标题】：Implement null distribution for gbm interaction strength

【发布时间】：2019-06-25 03:55:08

【问题描述】：

我正在尝试使用 Friedman 和 Popescu 2008 https://projecteuclid.org/euclid.aoas/1223908046 中描述的方法来确定 gbm 模型中的哪些交互是重要的。我的 gbm 是一个有 9 个不同类别的分类模型。 我正在努力解决如何将第 8.3 节转换为在 R 中运行的代码。

我认为整个过程是：

使用 max.depth = 1 训练模型版本 模拟来自该模型的响应数据 在此数据上训练一个新模型，其 max.depth 与真实模型相同 获取此模型的交互强度 重复步骤 1-4 以创建交互强度的零分布

我发现最令人困惑的部分是实现方程式 48 和 49。（您必须查看链接的文章，因为我无法在此处复制它们）

这是我认为我理解的，但如果我错了，请纠正我：

y_i 是一个新的响应向量，我们将使用它来训练一个新模型，该模型将提供交互统计的零分布。

F_A(x_i) 是使用 max.depth = 1 训练的 gbm 模型版本的预测

b_i 是基于加性模型 F_A(x_i)

预测的介于 0 和 1 之间的概率

问题

什么是下标i？是引导程序中的迭代次数吗？ 每个人工数据集与其他数据集有何不同？ 我们是否将 Pr(b_i = 1) 代入方程 48？ 如何通过多项分类来做到这一点？ 如何在 R 中实现这一点？最好使用 gbm 包。

欢迎任何想法或参考！

【问题讨论】：

您在stats.stackexchange.com 上发布此问题可能会更幸运，至少对于您在整个过程和解释符号方面的问题。 R 代码很有可能会包含在您获得的任何答案中，一旦您有一些代码可以开始，您可以随时在此处再次发布以寻求实施帮助。

【参考方案1】：

总体而言，该过程是一种通过置换/重新分配建模对交互的额外贡献来抵消y 中交互效果的优雅方式。额外的贡献可以通过完整模型和加法模型之间的边际来捕获。

什么是下标i？是引导程序中的迭代次数吗？

它是样本的索引。每次迭代都有N个样本。

每个人工数据集与其他数据集有何不同？

预测变量 X 在数据集中是相同的。由于equation 47 中边距的随机排列和equation 48 中的随机实现（仅用于分类结果），响应值Y~ 不同。

我们是否将 Pr(b_i = 1) 替换为 equation 48？

是的，如果结果 Y 是二进制的。

如何通过多项分类来做到这一点？

一种方法是随机排列每个类别的对数赔率中的边距。然后根据加法模型的概率随机实现。

如何在 R 中实现这一点？最好使用gbm 包。

我尝试按照您的整体流程实施它。

首先，模拟训练数据集X1,X2,Y，大小为N=200，其中Y具有三个类别（Y1，Y2，Y3），由X1确定的概率实现， X2。交互部分X1*X2在Y1，而附加部分在Y2,Y3。

set.seed(1)
N <- 200
X1 <- rnorm(N) # 2 predictors
X2 <- rnorm(N)

#log-odds for 3 categories
Y1 <- 2*X1*X2 + rnorm(N, sd=1/10) # interaction
Y2 <- X1^2 + rnorm(N, sd=1/10) #additive
Y3 <- X2^2 + rnorm(N, sd=1/10) #additive
Y <- rep(NA, N) # Multinomial outcome with 3 categories
for (i in 1:N)
  
    prob <- 1 / (1 + exp(-c(Y1[i],Y2[i],Y3[i]))) #logistic regression
    Y[i] <- which.max(rmultinom(1, 10000, prob=prob)) #realisation from prob
  
Y <- factor(Y)
levels(Y) <- c('Y1','Y2','Y3')
table(Y)
#Y1 Y2 Y3 
#38 75 87
dat = data.frame(Y, X1, X2)
head(dat)
# Y         X1         X2
# 2 -0.6264538  0.4094018
# 3  0.1836433  1.6888733
# 3 -0.8356286  1.5865884
# 2  1.5952808 -0.3309078
# 3  0.3295078 -2.2852355
# 3 -0.8204684  2.4976616

分别使用 max.depth = 2 和 1 训练完整模型和加法模型。

library(gbm)
n.trees <- 100
F_full <- gbm(Y ~ ., data=dat, distribution='multinomial', n.trees=n.trees, cv.folds=3,
            interaction.depth=2) # consider interactions
F_additive <- gbm(Y ~ ., data=dat, distribution='multinomial', n.trees=n.trees, cv.folds=3,
            interaction.depth=1) # ignore interactions

#use improved prediction as interaction strength
interaction_strength_original <- min(F_additive$cv.error) - min(F_full$cv.error)
> 0.1937891

模拟来自该模型的响应数据。

#randomly permute margins (residuals) of log-odds to remove any interaction effects
margin <- predict(F_full, n.trees=gbm.perf(F_full, plot.it=FALSE), type='link')[,,1] -
         predict(F_additive, n.trees=gbm.perf(F_additive, plot.it=FALSE), type='link')[,,1]
margin <- apply(margin, 2, sample) #independent permutation for each category (Y1, Y2, Y3)

Y_art <- rep(NA, N) #response values of an artificial dataset
for (i in 1:N)
  
    prob <- predict(F_additive, n.trees=gbm.perf(F_additive, plot.it=FALSE), type='link',
                  newdata=dat[i,])
    prob <- prob + margin[i,] # equation (47)
    prob <- 1 / (1 + exp(-prob))
    Y_art[i] <- which.max(rmultinom(1, 1000, prob=prob)) #Similar to random realisation in equation (49)
  
Y_art <- factor(Y_art)
levels(Y_art) = c('Y1','Y2','Y3')

table(Y_art)
#Y1 Y2 Y3 
#21 88 91

使用max.depth（2）在此人工数据上训练一个新模型，与真实模型相同

F_full_art = gbm(Y_art ~ ., distribution='multinomial', n.trees=n.trees, cv.folds=3,
            data=data.frame(Y_art, X1, X2),
            interaction.depth=2)
F_additive_art = gbm(Y_art ~ ., distribution='multinomial', n.trees=n.trees, cv.folds=3,
            data=data.frame(Y_art, X1, X2),
            interaction.depth=1)

获取此模型的交互强度

interaction_strength_art = min(F_additive_art$cv.error) - min(F_full_art$cv.error)
> 0.01323959 # much smaller than interaction_strength_original in step 1.

重复步骤 2-4 以创建交互强度的零分布。正如预期的那样，中和数据集中的交互效应（-0.0527 到 0.0421）远低于原始训练数据集（0.1938）。

interaction_strength_art <- NULL
for (j in 1:10)
   
    print(j)
    interaction_strength_art <- c(interaction_strength_art, step_2_to_4())
  

summary(interaction_strength_art)
#     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#-0.052648 -0.019415  0.001124 -0.004310  0.012759  0.042058

interaction_strength_original
> 0.1937891

【讨论】：

感谢您的回复！不幸的是，自从我发布以来已经很长时间了，我不再为此工作了，所以我无法检查答案。不过，我会尽快恢复使用 gbms。