R中的回归子集算法
Regression Subset Algorithm in R
我想在 R 中为 'beta regression' 模型构建回归子集算法。
R 中有一个适合 beta 回归的包 betareg，我感兴趣的是最大化 'log likelihood'.
的模型
基本上，这是通过选择最佳的 k 因子回归模型来实现的，对于 k = 1,2,...,p，其中 p 是您拥有的变量数。
例如，如果我将 x_1、x_2、x_3 作为我的变量，并将 y 作为我的响应。我想要一些东西：
第 1 步：找到最佳 1 因子模型
mod1 <- betareg(y~x_1, data = test)
mod1.sum <- summary(mod1)
mod2 <- betareg(y~x_2, data = test)
mod2.sum <- summary(mod2)
mod3 <- betareg(y~x_3, data = test)
mod3.sum <- summary(mod3)
现在我已经拟合了所有模型，我想比较每个模型的对数似然：
likelihoods <- c( mod1.sum$loglik, mod2.sum$loglik, mod3.sum$loglik)
which.max(likelihoods)
第 2 步：找到添加到最佳 1 因子模型的最佳因子，我们假设 x_1 是上一步中的最佳因子。然后在这一步中，我们将 x_1 和 x_2 的模型与 x_1 和 x_3 的模型进行比较，选择具有最大对数似然的模型。
第 3 步：将最好的两个变量作为给定，找出对对数似然增加贡献最大的第三个变量。
第 4 步：Return 最佳 1 因子模型，最佳 2 因子模型，...，最佳 p 因子模型，包含的因子及其对应的对数似然。
当 p 很大时，比如 40 左右，我很难有效地做到这一点
这是不使用 betareg 的替代解决方案。输出结果类似，供您考虑。
这是我使用的数据集：
set.seed(12345)
dat <- data.frame(y=runif(50), x_1=runif(50), x_2=runif(50), x_3=runif(50))
使用 leaps 库创建所有可能组合的列表：
library(leaps)
subs<-regsubsets(y~., data=dat, nbest=10, nvmax=100, really.big=T)
subs<-summary(subs)$which[,-1]
all.mods<-lapply(1:nrow(subs), function(x)paste("y", paste(names(which(subs[x,])), collapse="+"), sep="~"))
all.mods
[[1]]
[1] "y~x_2"
[[2]]
[1] "y~x_1"
[[3]]
[1] "y~x_3"
[[4]]
[1] "y~x_2+x_3"
[[5]]
[1] "y~x_1+x_2"
[[6]]
[1] "y~x_1+x_3"
[[7]]
[1] "y~x_1+x_2+x_3"
运行 所有模型的线性回归：
all.lm<-lapply(all.mods, function(x)lm(as.formula(x), data=dat))
检查每个模型的 logLikihood：
lapply(all.lm, logLik)
[[1]]
'log Lik.' -7.051835 (df=3)
[[2]]
'log Lik.' -9.288776 (df=3)
[[3]]
'log Lik.' -9.334048 (df=3)
[[4]]
'log Lik.' -6.904604 (df=4)
[[5]]
'log Lik.' -7.051584 (df=4)
[[6]]
'log Lik.' -9.215915 (df=4)
[[7]]
'log Lik.' -6.888849 (df=5)
据我所知，没有专门有效的最佳子集 selection 用于 beta 回归（在 R 或其他语言中）。然而，有一些通用实现为此提供了近似解决方案，例如，基于 kofnGA 包（在 CRAN and published in JSS 上）等遗传算法。请参见下面的示例。 （要使用 forward 搜索而不是 best-subset selection，请参阅我的其他答案。）
或者，您可以使用近似于 betareg 所做的（广义）线性模型，并为此使用子集 selection。例如，您可以对响应（即 qlogis(y)）进行对数变换，然后通过 leaps（CRAN) or lmSubsets (R-Forge). Or you could use a GLM with family = quasibinomial and use glmulti (CRAN, JSS).然后您可以使用该近似模型的最佳子集结果并将其用于 beta 回归。当然，这不会为您提供 最佳 beta 回归结果，但它可能是进一步分析的有用起点。
因此，回到β回归的直接遗传算法。为了说明如何使用 kofnGA 完成此操作，我们首先加载包和示例数据：
library("betareg")
library("kofnGA")
data("FoodExpenditure", package = "betareg")
然后，我们构建一个包含响应变量 y 和回归矩阵 x 的列表。请注意，我们在这里省略了截距，以便稍后将其强制进入模型（即截距不应受 selection）。
fe_data <- list(
y = with(FoodExpenditure, food/income),
x = model.matrix(~ income + persons, data = FoodExpenditure)[, -1]
)
除了上面设置的两个回归器之外，我们现在向回归器矩阵添加 40 个随机噪声变量
fe_data$x <- cbind(fe_data$x,
matrix(rnorm(40 * nrow(fe_data$x)), ncol = 40))
colnames(fe_data$x)[3:42] <- paste0("x", 1:40)
现在我们可以使用 kofnGA 到 select 可能的 42 个回归变量中有 2 个回归变量的最佳模型（加上始终包含的截距）。由于 kofnGA 最小化目标，我们使用 betareg 提供的负对数似然。主力函数 betareg.fit 而不是 betareg 用于避免不必要的公式解析等
nll <- function(v, data) -betareg.fit(x = cbind(1, data$x[, v]),
y = data$y)$loglik
最后，我们 运行 遗传算法仅用于 100 代，以在这个简短的例子中节省一些计算时间：
set.seed(1)
ga <- kofnGA(n = 42, k = 2, OF = nll, data = fe_data, ngen = 100)
结果输出是
summary(ga)
## Genetic algorithm search, 100 generations
## Number of unique solutions in the final population: 1
##
## Objective function values:
## average minimum
## Initial population -36.56597 -41.74583
## Final population -45.33351 -45.33351
##
## Best solution (found at generation 1):
## 1 2
因此，在这个非常简单的人工设置中，遗传算法确实选择了前 2 个回归变量（来自真实数据），而不是我们添加的任何不相关的随机 40 个回归变量。因此，我们现在可以继续使用回归变量
重新拟合适当的 beta 回归模型
colnames(fe_data$x)[ga$bestsol]
## [1] "income" "persons"
等请注意，此处使用的 beta 回归仅使用固定精度参数（log link）。如果您想要可变色散，则需要相应地修改 nll。
除了我的其他答案显示如何使用 kofnGA 执行 best-subset selection 进行 beta 回归外，我还包括一个例如如何手动 forward selection。
我们再次从加载包和数据开始：
library("betareg")
data("FoodExpenditure", package = "betareg")
我还设置了包含响应和所有回归变量的列表（包括 40 个随机回归变量）。（请注意，与另一个不同，我将截距保留在 x 中，这在这里更方便。）
fe_data <- list(
y = with(FoodExpenditure, food/income),
x = model.matrix(~ income + persons, data = FoodExpenditure)
)
set.seed(123)
fe_data$x <- cbind(fe_data$x,
matrix(rnorm(40 * nrow(fe_data$x)), ncol = 40))
colnames(fe_data$x)[4:43] <- paste0("x", 1:40)
然后我们再次设置负对数似然函数（但不需要手动包括截距，因为它仍在x）。
nll <- function(v, data) -betareg.fit(x = data$x[, v, drop = FALSE],
y = data$y)$loglik
然后我们存储所有可能回归变量的索引 vall 并使用截距 (v <- 1) 和相应的负对数似然 (n) 初始化我们的搜索。
vall <- 1:ncol(fe_data$x)
v <- 1
n <- nll(v, data = fe_data)
然后我们将前向搜索迭代 15 次（以避免在这个小数据集上出现数值不稳定性以获得更多变量）。我们总是 select 最能减少负对数似然的附加变量：
for(i in 1:15) {
vi <- vall[-v]
ni <- sapply(vi, function(vii) nll(v = c(v, vii), data = fe_data))
v <- c(v, vi[which.min(ni)])
n <- c(n, ni[which.min(ni)])
}
选择变量的顺序如下。请注意，首先选择真实回归器，然后选择随机噪声回归器。 （但尝试 set.seed(1) 以上将在真实回归之前包括随机回归...）
colnames(fe_data$x)[v]
## [1] "(Intercept)" "income" "persons" "x28" "x18"
## [6] "x29" "x22" "x11" "x5" "x8"
## [11] "x38" "x24" "x13" "x23" "x36"
## [16] "x16"
负对数似然和相关 BIC 的相应减少可以可视化为：
m <- seq_along(v)
plot(m, n, type = "b",
xlab = "Number of regressors", ylab = "Log-likelihood")
plot(m, n + log(nrow(fe_data$x)) * (m + 1), type = "b",
xlab = "Number of regressors", ylab = "BIC")
所以这确实会选择具有三个真实回归变量的模型作为最佳 BIC 模型。
我想在 R 中为 'beta regression' 模型构建回归子集算法。
R 中有一个适合 beta 回归的包 betareg，我感兴趣的是最大化 'log likelihood'.
基本上，这是通过选择最佳的 k 因子回归模型来实现的，对于 k = 1,2,...,p，其中 p 是您拥有的变量数。
例如，如果我将 x_1、x_2、x_3 作为我的变量，并将 y 作为我的响应。我想要一些东西：
第 1 步：找到最佳 1 因子模型
mod1 <- betareg(y~x_1, data = test)
mod1.sum <- summary(mod1)
mod2 <- betareg(y~x_2, data = test)
mod2.sum <- summary(mod2)
mod3 <- betareg(y~x_3, data = test)
mod3.sum <- summary(mod3)
现在我已经拟合了所有模型，我想比较每个模型的对数似然：
likelihoods <- c( mod1.sum$loglik, mod2.sum$loglik, mod3.sum$loglik)
which.max(likelihoods)
第 2 步：找到添加到最佳 1 因子模型的最佳因子，我们假设 x_1 是上一步中的最佳因子。然后在这一步中，我们将 x_1 和 x_2 的模型与 x_1 和 x_3 的模型进行比较，选择具有最大对数似然的模型。
第 3 步：将最好的两个变量作为给定，找出对对数似然增加贡献最大的第三个变量。
第 4 步：Return 最佳 1 因子模型，最佳 2 因子模型，...，最佳 p 因子模型，包含的因子及其对应的对数似然。
当 p 很大时，比如 40 左右，我很难有效地做到这一点
这是不使用 betareg 的替代解决方案。输出结果类似，供您考虑。
这是我使用的数据集：
set.seed(12345)
dat <- data.frame(y=runif(50), x_1=runif(50), x_2=runif(50), x_3=runif(50))
使用 leaps 库创建所有可能组合的列表：
library(leaps)
subs<-regsubsets(y~., data=dat, nbest=10, nvmax=100, really.big=T)
subs<-summary(subs)$which[,-1]
all.mods<-lapply(1:nrow(subs), function(x)paste("y", paste(names(which(subs[x,])), collapse="+"), sep="~"))
all.mods
[[1]]
[1] "y~x_2"
[[2]]
[1] "y~x_1"
[[3]]
[1] "y~x_3"
[[4]]
[1] "y~x_2+x_3"
[[5]]
[1] "y~x_1+x_2"
[[6]]
[1] "y~x_1+x_3"
[[7]]
[1] "y~x_1+x_2+x_3"
运行 所有模型的线性回归：
all.lm<-lapply(all.mods, function(x)lm(as.formula(x), data=dat))
检查每个模型的 logLikihood：
lapply(all.lm, logLik)
[[1]]
'log Lik.' -7.051835 (df=3)
[[2]]
'log Lik.' -9.288776 (df=3)
[[3]]
'log Lik.' -9.334048 (df=3)
[[4]]
'log Lik.' -6.904604 (df=4)
[[5]]
'log Lik.' -7.051584 (df=4)
[[6]]
'log Lik.' -9.215915 (df=4)
[[7]]
'log Lik.' -6.888849 (df=5)
据我所知，没有专门有效的最佳子集 selection 用于 beta 回归（在 R 或其他语言中）。然而，有一些通用实现为此提供了近似解决方案，例如，基于 kofnGA 包（在 CRAN and published in JSS 上）等遗传算法。请参见下面的示例。 （要使用 forward 搜索而不是 best-subset selection，请参阅我的其他答案。）
或者，您可以使用近似于 betareg 所做的（广义）线性模型，并为此使用子集 selection。例如，您可以对响应（即 qlogis(y)）进行对数变换，然后通过 leaps（CRAN) or lmSubsets (R-Forge). Or you could use a GLM with family = quasibinomial and use glmulti (CRAN, JSS).然后您可以使用该近似模型的最佳子集结果并将其用于 beta 回归。当然，这不会为您提供 最佳 beta 回归结果，但它可能是进一步分析的有用起点。
因此，回到β回归的直接遗传算法。为了说明如何使用 kofnGA 完成此操作，我们首先加载包和示例数据：
library("betareg")
library("kofnGA")
data("FoodExpenditure", package = "betareg")
然后，我们构建一个包含响应变量 y 和回归矩阵 x 的列表。请注意，我们在这里省略了截距，以便稍后将其强制进入模型（即截距不应受 selection）。
fe_data <- list(
y = with(FoodExpenditure, food/income),
x = model.matrix(~ income + persons, data = FoodExpenditure)[, -1]
)
除了上面设置的两个回归器之外，我们现在向回归器矩阵添加 40 个随机噪声变量
fe_data$x <- cbind(fe_data$x,
matrix(rnorm(40 * nrow(fe_data$x)), ncol = 40))
colnames(fe_data$x)[3:42] <- paste0("x", 1:40)
现在我们可以使用 kofnGA 到 select 可能的 42 个回归变量中有 2 个回归变量的最佳模型（加上始终包含的截距）。由于 kofnGA 最小化目标，我们使用 betareg 提供的负对数似然。主力函数 betareg.fit 而不是 betareg 用于避免不必要的公式解析等
nll <- function(v, data) -betareg.fit(x = cbind(1, data$x[, v]),
y = data$y)$loglik
最后，我们 运行 遗传算法仅用于 100 代，以在这个简短的例子中节省一些计算时间：
set.seed(1)
ga <- kofnGA(n = 42, k = 2, OF = nll, data = fe_data, ngen = 100)
结果输出是
summary(ga)
## Genetic algorithm search, 100 generations
## Number of unique solutions in the final population: 1
##
## Objective function values:
## average minimum
## Initial population -36.56597 -41.74583
## Final population -45.33351 -45.33351
##
## Best solution (found at generation 1):
## 1 2
因此，在这个非常简单的人工设置中，遗传算法确实选择了前 2 个回归变量（来自真实数据），而不是我们添加的任何不相关的随机 40 个回归变量。因此，我们现在可以继续使用回归变量
重新拟合适当的 beta 回归模型colnames(fe_data$x)[ga$bestsol]
## [1] "income" "persons"
等请注意，此处使用的 beta 回归仅使用固定精度参数（log link）。如果您想要可变色散，则需要相应地修改 nll。
除了我的其他答案显示如何使用 kofnGA 执行 best-subset selection 进行 beta 回归外，我还包括一个例如如何手动 forward selection。
我们再次从加载包和数据开始：
library("betareg")
data("FoodExpenditure", package = "betareg")
我还设置了包含响应和所有回归变量的列表（包括 40 个随机回归变量）。（请注意，与另一个不同，我将截距保留在 x 中，这在这里更方便。）
fe_data <- list(
y = with(FoodExpenditure, food/income),
x = model.matrix(~ income + persons, data = FoodExpenditure)
)
set.seed(123)
fe_data$x <- cbind(fe_data$x,
matrix(rnorm(40 * nrow(fe_data$x)), ncol = 40))
colnames(fe_data$x)[4:43] <- paste0("x", 1:40)
然后我们再次设置负对数似然函数（但不需要手动包括截距，因为它仍在x）。
nll <- function(v, data) -betareg.fit(x = data$x[, v, drop = FALSE],
y = data$y)$loglik
然后我们存储所有可能回归变量的索引 vall 并使用截距 (v <- 1) 和相应的负对数似然 (n) 初始化我们的搜索。
vall <- 1:ncol(fe_data$x)
v <- 1
n <- nll(v, data = fe_data)
然后我们将前向搜索迭代 15 次（以避免在这个小数据集上出现数值不稳定性以获得更多变量）。我们总是 select 最能减少负对数似然的附加变量：
for(i in 1:15) {
vi <- vall[-v]
ni <- sapply(vi, function(vii) nll(v = c(v, vii), data = fe_data))
v <- c(v, vi[which.min(ni)])
n <- c(n, ni[which.min(ni)])
}
选择变量的顺序如下。请注意，首先选择真实回归器，然后选择随机噪声回归器。 （但尝试 set.seed(1) 以上将在真实回归之前包括随机回归...）
colnames(fe_data$x)[v]
## [1] "(Intercept)" "income" "persons" "x28" "x18"
## [6] "x29" "x22" "x11" "x5" "x8"
## [11] "x38" "x24" "x13" "x23" "x36"
## [16] "x16"
负对数似然和相关 BIC 的相应减少可以可视化为：
m <- seq_along(v)
plot(m, n, type = "b",
xlab = "Number of regressors", ylab = "Log-likelihood")
plot(m, n + log(nrow(fe_data$x)) * (m + 1), type = "b",
xlab = "Number of regressors", ylab = "BIC")
所以这确实会选择具有三个真实回归变量的模型作为最佳 BIC 模型。