r - 为 gbm 交互强度实现零分布
问题描述
我正在尝试使用 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 包。
欢迎任何想法或参考!
解决方案
y
总体而言,该过程是一种通过置换/重新分配建模对交互的额外贡献来中和交互效果的优雅方式。额外的贡献可以通过完整模型和加法模型之间的边际来捕获。
- 什么是下标
i
?它是引导程序中的迭代次数吗?
它是样本的索引。每次迭代都有N
样本。
- 每个人工数据集与其他数据集有何不同?
X
跨数据集的预测变量是相同的。由于边缘的随机排列和随机实现(仅用于分类结果),响应值Y~
不同。equation 47
equation 48
- 我们是否将 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
推荐阅读
- python - 将特定于类别的值分配为新列中的观察值
- firebase - 使用 Spring Boot 将图像上传到 Firebase
- redux - 使用(Material-ui + React Redux + SSR 中的 Next.js)重新渲染未更新传递的属性
- windows - 显示内存内容时,windbg 是否会自行处理内核中的分页代码?
- python - 当特定语句需要更多时间执行时,如何在 Python for 循环中跳过迭代?
- csv - 用于搜索(大)未排序数组的索引(类型)
- python - Scrapy CrawlSpider 无法抓取首页
- azure-devops - 如何创建 ADO 工作项,将迭代路径设置为 REST API 中的当前迭代?
- c++ - 使用递归的c ++中的分段错误
- google-app-engine - 使用服务帐户的 GCP 服务器到服务器身份验证