14.2 异质截距项的线性模型
该节是对A Concave Pairwise Fusion Approach to Subgroup Analysis[32]论文的复现。
\[ y_i = \mu_i+\mathbf{x}_i^T\beta+\varepsilon_i \tag{14.1} \]
该论文在普通线性模型的基础上，考虑了截距项的异质性，即样本内部可能存在亚组，同一亚组共享相同的截距项，不同的亚组之间的截距项不同。该论文利用增广拉格朗日法构造目标函数（引入融合惩罚项\(\mu_i-\mu_j=\eta_{ij}\)），并通过ADMM算法进行求解。
14.2.1 自定义算法
算法逻辑：
传入参数
- y：向量，响应变量
- x：矩阵，预测变量
- scale：是否对x进行标准化，默认为
TRUE - penalty：设置惩罚函数类型(L1、MCP、SCAD)
- \(\theta\)：ADMM算法的惩罚系数
- \(\lambda\)：惩罚函数中的惩罚系数
- \(\gamma\)：惩罚函数中的正则化因子
- tol：收敛精度，默认为0.00001
- max_iter：最大迭代次数，默认为1000
其余符号说明
除了传入参数外，在运算过程中还有其它符号，其含义如下所示。
- \(\upsilon\)：拉格朗日乘子
- \(Q\)：\(X(X^TX)^{-1}X^T\)
- \(Delta\)：\(\{(e_i-e_j), i<j\}^T\)，其中\(e_i\)表示第\(i\)个分量为1，其余分量为0的n维向量
- \(\delta\_\mu\)：\(\mu_i-\mu_j\)
- \(\delta\)：\(\mu_i-\mu_j+\theta^{-1}\upsilon_{ij}\)
- \(\alpha\)：亚组的截距项
初始值
根据ols估计获取\(\beta\)、\(\mu\)、\(\eta\)、\(\upsilon\)的初始值。
迭代
每次循环按照\(\mu, \beta, \delta, \eta, \upsilon\)的顺序进行迭代。
\[ \begin{aligned} \mu^{(m+1)}&=(\theta \Delta^T\Delta+I_n-Q)^{-1}((I_n-Q)y+\theta \Delta^T(\eta^{(m)}-\theta^{-1}\upsilon^{(m)})) \\ \beta^{(m+1)}&=(X^TX)^{-1}X^T(y-\mu^{(m+1)}) \\ \delta_{ij}^{(m+1)}&=\mu_i^{(m+1)}-\mu_j^{(m+1)}+\theta^{-1}\upsilon_{ij}^{(m)} \\ \eta_{ij}^{(m+1)}&=\textrm{Penalty} \\ \upsilon_{ij}^{(m+1)} &= \upsilon_{ij}^{(m)}+\theta (\mu_i^{(m+1)}-\mu_j^{(m+1)}-\eta_{ij}^{(m+1)}) \end{aligned} \]
注意\(\delta_{ij}\)与\(\eta_{ij}\)是针对每一个分量而言的，得将每一个分量代入惩罚函数才能得到对应的结果。
停止
当迭代次数达到最大迭代次数或残差\(r^{(m+1)}=\Delta\mu^{(m+1)}-\eta^{(m+1)}\)的模长足够小时，则停止迭代。
输出
输出每个观测对象的截距项\(\hat \mu_i\)，回归系数\(\hat \beta\)，截距项对应的亚组标签，每个亚组的截距项均值\(\hat \alpha_k\)。
注意，论文涉及到\(i<j\)的排序，自定义算法采取先控制i再迭代j的步骤进行排序，如(1,2)、(1,3)、(1,4)…该顺序影响了后面矩阵的排列方式与lower.tri
由于设备限制，为节省运行时间，并未实现原文中的“\(\lambda\)路径图”及“对\(\lambda\)调参”的功能
自定义算法如下所示
library(tidyverse)
library(igraph)
library(R6)
SubgroupIntercept <- R6Class(
classname <- 'SubgroupIntercept',
public <- list(
# 传入参数
y = NULL, # 响应变量
x = NULL, # 预测变量
n = NULL, # 样本容量
scale = NULL, # 是否标准化
penalty = NULL, # 惩罚函数类型：L1、MCP、SCAD
theta = NULL, # 式(4)二次项的惩罚系数
lambda = NULL, # 惩罚函数中的惩罚系数
gamma = NULL, # 惩罚函数中的正则化因子
tol = NULL, # 收敛精度
max_iter = NULL, # 最大迭代次数
# 初始化
initialize = function(y, x, scale = T, penalty, theta, lambda, gamma, tol = 0.00001, max_iter = 1000){
self$y <- y
self$x <- x
self$n <- length(y)
self$scale <- scale
self$penalty <- penalty
self$theta <- theta
self$lambda <- lambda
self$gamma <- gamma
self$tol <- tol
self$max_iter <- max_iter
},
# 主函数——运行
run = function(){
start_time <- Sys.time()
# 检验输入是否合理
private$validate()
# 标准化协变量
if(self$scale == T) self$x <- apply(self$x, 2, scale)
# 获取初始值
initial_value <- private$initial_value()
private$beta <- initial_value[[1]]
private$mu <- initial_value[[2]]
private$eta <- initial_value[[3]]
private$upsilon <- initial_value[[4]]
# 计算Q矩阵
Q <- private$gen_Q()
# 计算Delta
Delta <- private$gen_Delta()
# 开始迭代
for (i in 1:self$max_iter) {
print(paste0('正在进行第[', i, ']次迭代'))
private$mu <- private$iter_mu(Q, Delta, private$eta, private$upsilon)
private$beta <- private$iter_beta(private$mu)
private$delta_mu <- private$iter_delta_mu(private$mu)
private$delta <- private$iter_delta(private$delta_mu, private$upsilon)
private$eta <- private$iter_eta(private$delta, self$lambda)
private$upsilon <- private$iter_upsilon(private$upsilon, private$delta_mu, private$eta)
# 终止条件
r <- Delta[[1]] %*% private$mu - private$eta
if(sqrt(sum(r^2)) <= self$tol){
print('达到精度要求')
break
}
}
# 获取亚组分类结果
label <- private$cluster_eta(private$eta)
# 根据亚组计算对应的alpha
subgroup_alpha <- tibble(mu = private$mu, label = label)
subgroup_alpha <- subgroup_alpha %>%
group_by(label) %>%
summarise(alpha = mean(mu), size = n())
result <- list(mu = private$mu, beta = private$beta, label = label, alpha = subgroup_alpha)
end_time <- Sys.time()
cost_time <- end_time - start_time
print(paste('花费时间：', round(cost_time, 3), "秒"))
return(result)
}
),
private <- list(
# 迭代参数
mu = NULL, # 截距项
beta = NULL, # 回归系数
delta_mu = NULL, # mu_i-mu_j
delta = NULL, # delta_ij=mu_i-mu_j+1/theta*upsilon
eta = NULL, # 由惩罚函数计算
upsilon =NULL, # 拉格朗日乘子
# 验证输入是否正确
validate = function(){
if(!is.vector(self$y)) stop('Y要求为向量')
if(!is.matrix(self$x)) stop('X要求为矩阵')
if(!is.logical(self$scale)) stop('scale要求为布尔值')
if(!self$penalty %in% c('L1', 'MCP','SCAD')) stop('请选择合适的惩罚函数')
if(self$theta <= 0) stop('请选择合适的theta值')
if(self$lambda <= 0) stop('请选择合适的lambda值')
if(self$gamma <= 0) stop('请选择合适的gamma值')
if(self$tol <= 0) stop('请选择合适的精度要求')
if(self$max_iter <= 0) stop('请选择合适的最大迭代次数')
},
# 获取初始值
initial_value = function(){
df <- cbind(self$y, self$x) %>% as.data.frame()
colnames(df) <- c('y', paste0('x_', 1:dim(x)[2]))
ols <- lm(y~., data = df)
beta_0 <- coef(ols)[-1]
mu_0 <- (self$y - self$x %*% beta_0) %>% as.vector()
eta_0 <- -outer(mu_0, mu_0, '-') # lower.tri按列提取下三角，故添加负号
eta_0 <- eta_0[lower.tri(eta_0)]
upsilon_0 <- rep(0, times = self$n * (self$n - 1)/2)
result <- list(beta_0 = beta_0, mu_0 = mu_0, eta_0 = eta_0, upsilon_0 = upsilon_0)
return(result)
},
# 计算Q矩阵
gen_Q = function(){
Q <- self$x %*% solve(t(self$x) %*% self$x) %*% t(self$x)
return(Q)
},
# 计算Delta矩阵
gen_Delta = function(){
gen_mat <- function(n){
mat_1 <- if(n == self$n-1){
NULL
}else{
matrix(0, nrow = self$n-1-n, ncol = n)
}
mat_2 <- matrix(1, nrow = 1, ncol = n, byrow = T)
mat_3 <- diag(-1, nrow=n)
mat <- rbind(mat_1, mat_2, mat_3)
return(mat)
}
Delta <- as.list(c((self$n-1):1)) %>% map(~gen_mat(.))
Delta <- do.call(cbind, Delta) %>% t()
Delta_Delta <- diag(self$n, nrow = self$n) - matrix(1, nrow = self$n, ncol = self$n)
result <- list(Delta = Delta, Delta_Delta = Delta_Delta)
return(result)
},
# mu迭代式
iter_mu = function(Q, Delta, eta_current, upsilon_current){
term_1 <- solve(self$theta * Delta[[2]] + diag(1, nrow = self$n) - Q)
term_2 <- (diag(1, nrow = self$n) - Q) %*% self$y + self$theta * t(Delta[[1]]) %*% (eta_current - 1/self$theta * upsilon_current)
mu_next <- (term_1 %*% term_2) %>% as.vector()
return(mu_next)
},
# beta迭代式
iter_beta = function(mu_next){
beta_next <- solve(t(self$x) %*% self$x) %*% t(self$x) %*% (self$y - mu_next)
return(beta_next)
},
# eta迭代式
iter_eta = function(delta_current, lambda){
ST <- function(t, lambda) sign(t) * max(abs(t)-lambda,0)
delta_list <- as.list(delta_current)
delta_to_eta <- function(delta_ij){
switch(self$penalty,
'L1' = {
eta_ij <- ST(delta_ij, self$lambda/self$theta)
},
'MCP' = {
if(abs(delta_ij) <= (self$gamma*self$lambda)){
eta_ij <- ST(delta_ij, self$lambda/self$theta)/(1-1/(self$gamma*self$theta))
}else{
eta_ij <- delta_ij
}
},
'SCAD' = {
if(abs(delta_ij) <= (self$lambda+self$lambda/self$theta)){
eta_ij <- ST(delta_ij, self$lambda/self$theta)
}else if(abs(delta_ij) > (self$lambda+self$lambda/self$theta) & abs(delta_ij) <= (self$gamma*self$lambda)){
eta_ij <- ST(delta_ij, self$gamma*self$lambda/((self$gamma-1)*self$theta))/(1-1/((self$gamma-1)*self$theta))
}else{
eta_ij <- delta_ij
}
}
)
return(eta_ij)
}
eta_next <- delta_list %>% map_vec(~delta_to_eta(.))
return(eta_next)
},
# mu_i-mu_j迭代式
iter_delta_mu = function(mu_next){
delta_mu <- -outer(mu_next, mu_next, '-')
delta_mu <- delta_mu[lower.tri(delta_mu)]
return(delta_mu)
},
# delta迭代式
iter_delta = function(delta_mu, upsilon_current){
delta_next <- delta_mu + 1/self$theta * upsilon_current
return(delta_next)
},
# upsilon迭代式
iter_upsilon = function(upsilon_current, delta_mu, eta_next){
upsilon_next <- upsilon_current + self$theta * (delta_mu - eta_next)
return(upsilon_next)
},
# eta归类
# 注意到根据eta为0来分组刚好符合图论中“连通分量”的概念
cluster_eta = function(eta){
mat <- matrix(NA, nrow = self$n, ncol = self$n)
mat[lower.tri(mat)] <- eta # 转化为矩阵方便提取索引
link <- which(mat == 0, arr.ind = T)
node <- data.frame(name = 1:self$n)
graph <- graph_from_data_frame(link, directed = F, vertices = node)
label <- components(graph)$membership
return(label)
}
)
)14.2.2 数据模拟
下面随机生成100条观测，设置3个预测变量，均服从标准正态分布，且两两之间的相关系数为0.3，然后随机生成\(\beta\)值与\(\varepsilon\)，设置截距项以相等的概率随机分为1或-1，根据上述变量生成响应变量。
library(MASS)
set.seed(123)
sigma <- matrix(0.3, ncol = 3, nrow = 3) + diag(0.7, nrow = 3)
x <- mvrnorm(n=100, mu=rep(0,3), Sigma = sigma)
epsilon <- rnorm(100)
mu <- sample(c(1,-1), 100, replace = T, prob = rep(0.5, 2))
beta <- runif(3)
y <- (mu + x %*% beta + epsilon) %>% as.vector()构建模型，并运行。
model <- SubgroupIntercept$new(y, x, penalty = 'MCP', theta = 1, lambda = 0.5, gamma = 3, max_iter = 10000)
result <- model$run()## [1] "正在进行第[1]次迭代"
## [1] "正在进行第[2]次迭代"
## [1] "正在进行第[3]次迭代"
## [1] "正在进行第[4]次迭代"
## [1] "正在进行第[5]次迭代"
## [1] "正在进行第[6]次迭代"
## [1] "正在进行第[7]次迭代"
## [1] "正在进行第[8]次迭代"
## [1] "正在进行第[9]次迭代"
## [1] "正在进行第[10]次迭代"
## [1] "正在进行第[11]次迭代"
## [1] "正在进行第[12]次迭代"
## [1] "正在进行第[13]次迭代"
## [1] "正在进行第[14]次迭代"
## [1] "正在进行第[15]次迭代"
## [1] "正在进行第[16]次迭代"
## [1] "正在进行第[17]次迭代"
## [1] "正在进行第[18]次迭代"
## [1] "正在进行第[19]次迭代"
## [1] "正在进行第[20]次迭代"
## [1] "正在进行第[21]次迭代"
## [1] "正在进行第[22]次迭代"
## [1] "正在进行第[23]次迭代"
## [1] "正在进行第[24]次迭代"
## [1] "正在进行第[25]次迭代"
## [1] "正在进行第[26]次迭代"
## [1] "正在进行第[27]次迭代"
## [1] "正在进行第[28]次迭代"
## [1] "正在进行第[29]次迭代"
## [1] "正在进行第[30]次迭代"
## [1] "正在进行第[31]次迭代"
## [1] "正在进行第[32]次迭代"
## [1] "正在进行第[33]次迭代"
## [1] "正在进行第[34]次迭代"
## [1] "正在进行第[35]次迭代"
## [1] "正在进行第[36]次迭代"
## [1] "正在进行第[37]次迭代"
## [1] "正在进行第[38]次迭代"
## [1] "正在进行第[39]次迭代"
## [1] "正在进行第[40]次迭代"
## [1] "正在进行第[41]次迭代"
## [1] "正在进行第[42]次迭代"
## [1] "正在进行第[43]次迭代"
## [1] "正在进行第[44]次迭代"
## [1] "正在进行第[45]次迭代"
## [1] "正在进行第[46]次迭代"
## [1] "正在进行第[47]次迭代"
## [1] "正在进行第[48]次迭代"
## [1] "正在进行第[49]次迭代"
## [1] "正在进行第[50]次迭代"
## [1] "正在进行第[51]次迭代"
## [1] "正在进行第[52]次迭代"
## [1] "正在进行第[53]次迭代"
## [1] "正在进行第[54]次迭代"
## [1] "正在进行第[55]次迭代"
## [1] "正在进行第[56]次迭代"
## [1] "正在进行第[57]次迭代"
## [1] "正在进行第[58]次迭代"
## [1] "正在进行第[59]次迭代"
## [1] "正在进行第[60]次迭代"
## [1] "正在进行第[61]次迭代"
## [1] "正在进行第[62]次迭代"
## [1] "正在进行第[63]次迭代"
## [1] "正在进行第[64]次迭代"
## [1] "正在进行第[65]次迭代"
## [1] "正在进行第[66]次迭代"
## [1] "正在进行第[67]次迭代"
## [1] "正在进行第[68]次迭代"
## [1] "正在进行第[69]次迭代"
## [1] "正在进行第[70]次迭代"
## [1] "正在进行第[71]次迭代"
## [1] "正在进行第[72]次迭代"
## [1] "正在进行第[73]次迭代"
## [1] "正在进行第[74]次迭代"
## [1] "正在进行第[75]次迭代"
## [1] "正在进行第[76]次迭代"
## [1] "正在进行第[77]次迭代"
## [1] "正在进行第[78]次迭代"
## [1] "正在进行第[79]次迭代"
## [1] "正在进行第[80]次迭代"
## [1] "正在进行第[81]次迭代"
## [1] "正在进行第[82]次迭代"
## [1] "正在进行第[83]次迭代"
## [1] "正在进行第[84]次迭代"
## [1] "正在进行第[85]次迭代"
## [1] "正在进行第[86]次迭代"
## [1] "正在进行第[87]次迭代"
## [1] "正在进行第[88]次迭代"
## [1] "正在进行第[89]次迭代"
## [1] "正在进行第[90]次迭代"
## [1] "正在进行第[91]次迭代"
## [1] "正在进行第[92]次迭代"
## [1] "正在进行第[93]次迭代"
## [1] "正在进行第[94]次迭代"
## [1] "正在进行第[95]次迭代"
## [1] "正在进行第[96]次迭代"
## [1] "正在进行第[97]次迭代"
## [1] "正在进行第[98]次迭代"
## [1] "正在进行第[99]次迭代"
## [1] "正在进行第[100]次迭代"
## [1] "正在进行第[101]次迭代"
## [1] "正在进行第[102]次迭代"
## [1] "正在进行第[103]次迭代"
## [1] "正在进行第[104]次迭代"
## [1] "正在进行第[105]次迭代"
## [1] "正在进行第[106]次迭代"
## [1] "正在进行第[107]次迭代"
## [1] "正在进行第[108]次迭代"
## [1] "正在进行第[109]次迭代"
## [1] "正在进行第[110]次迭代"
## [1] "正在进行第[111]次迭代"
## [1] "正在进行第[112]次迭代"
## [1] "正在进行第[113]次迭代"
## [1] "正在进行第[114]次迭代"
## [1] "正在进行第[115]次迭代"
## [1] "正在进行第[116]次迭代"
## [1] "正在进行第[117]次迭代"
## [1] "正在进行第[118]次迭代"
## [1] "正在进行第[119]次迭代"
## [1] "正在进行第[120]次迭代"
## [1] "正在进行第[121]次迭代"
## [1] "正在进行第[122]次迭代"
## [1] "正在进行第[123]次迭代"
## [1] "正在进行第[124]次迭代"
## [1] "正在进行第[125]次迭代"
## [1] "正在进行第[126]次迭代"
## [1] "正在进行第[127]次迭代"
## [1] "正在进行第[128]次迭代"
## [1] "正在进行第[129]次迭代"
## [1] "正在进行第[130]次迭代"
## [1] "正在进行第[131]次迭代"
## [1] "正在进行第[132]次迭代"
## [1] "正在进行第[133]次迭代"
## [1] "正在进行第[134]次迭代"
## [1] "正在进行第[135]次迭代"
## [1] "正在进行第[136]次迭代"
## [1] "正在进行第[137]次迭代"
## [1] "正在进行第[138]次迭代"
## [1] "正在进行第[139]次迭代"
## [1] "正在进行第[140]次迭代"
## [1] "正在进行第[141]次迭代"
## [1] "正在进行第[142]次迭代"
## [1] "正在进行第[143]次迭代"
## [1] "正在进行第[144]次迭代"
## [1] "正在进行第[145]次迭代"
## [1] "正在进行第[146]次迭代"
## [1] "正在进行第[147]次迭代"
## [1] "正在进行第[148]次迭代"
## [1] "正在进行第[149]次迭代"
## [1] "正在进行第[150]次迭代"
## [1] "正在进行第[151]次迭代"
## [1] "正在进行第[152]次迭代"
## [1] "正在进行第[153]次迭代"
## [1] "正在进行第[154]次迭代"
## [1] "正在进行第[155]次迭代"
## [1] "正在进行第[156]次迭代"
## [1] "正在进行第[157]次迭代"
## [1] "正在进行第[158]次迭代"
## [1] "正在进行第[159]次迭代"
## [1] "正在进行第[160]次迭代"
## [1] "正在进行第[161]次迭代"
## [1] "正在进行第[162]次迭代"
## [1] "正在进行第[163]次迭代"
## [1] "正在进行第[164]次迭代"
## [1] "正在进行第[165]次迭代"
## [1] "正在进行第[166]次迭代"
## [1] "正在进行第[167]次迭代"
## [1] "正在进行第[168]次迭代"
## [1] "正在进行第[169]次迭代"
## [1] "正在进行第[170]次迭代"
## [1] "正在进行第[171]次迭代"
## [1] "正在进行第[172]次迭代"
## [1] "正在进行第[173]次迭代"
## [1] "正在进行第[174]次迭代"
## [1] "正在进行第[175]次迭代"
## [1] "正在进行第[176]次迭代"
## [1] "正在进行第[177]次迭代"
## [1] "正在进行第[178]次迭代"
## [1] "正在进行第[179]次迭代"
## [1] "正在进行第[180]次迭代"
## [1] "正在进行第[181]次迭代"
## [1] "正在进行第[182]次迭代"
## [1] "正在进行第[183]次迭代"
## [1] "正在进行第[184]次迭代"
## [1] "正在进行第[185]次迭代"
## [1] "正在进行第[186]次迭代"
## [1] "正在进行第[187]次迭代"
## [1] "正在进行第[188]次迭代"
## [1] "正在进行第[189]次迭代"
## [1] "正在进行第[190]次迭代"
## [1] "正在进行第[191]次迭代"
## [1] "正在进行第[192]次迭代"
## [1] "正在进行第[193]次迭代"
## [1] "正在进行第[194]次迭代"
## [1] "正在进行第[195]次迭代"
## [1] "正在进行第[196]次迭代"
## [1] "正在进行第[197]次迭代"
## [1] "正在进行第[198]次迭代"
## [1] "正在进行第[199]次迭代"
## [1] "正在进行第[200]次迭代"
## [1] "正在进行第[201]次迭代"
## [1] "正在进行第[202]次迭代"
## [1] "正在进行第[203]次迭代"
## [1] "正在进行第[204]次迭代"
## [1] "正在进行第[205]次迭代"
## [1] "正在进行第[206]次迭代"
## [1] "正在进行第[207]次迭代"
## [1] "正在进行第[208]次迭代"
## [1] "正在进行第[209]次迭代"
## [1] "正在进行第[210]次迭代"
## [1] "正在进行第[211]次迭代"
## [1] "正在进行第[212]次迭代"
## [1] "正在进行第[213]次迭代"
## [1] "正在进行第[214]次迭代"
## [1] "正在进行第[215]次迭代"
## [1] "正在进行第[216]次迭代"
## [1] "正在进行第[217]次迭代"
## [1] "正在进行第[218]次迭代"
## [1] "正在进行第[219]次迭代"
## [1] "正在进行第[220]次迭代"
## [1] "正在进行第[221]次迭代"
## [1] "正在进行第[222]次迭代"
## [1] "正在进行第[223]次迭代"
## [1] "正在进行第[224]次迭代"
## [1] "正在进行第[225]次迭代"
## [1] "正在进行第[226]次迭代"
## [1] "正在进行第[227]次迭代"
## [1] "正在进行第[228]次迭代"
## [1] "正在进行第[229]次迭代"
## [1] "正在进行第[230]次迭代"
## [1] "正在进行第[231]次迭代"
## [1] "正在进行第[232]次迭代"
## [1] "正在进行第[233]次迭代"
## [1] "正在进行第[234]次迭代"
## [1] "正在进行第[235]次迭代"
## [1] "正在进行第[236]次迭代"
## [1] "正在进行第[237]次迭代"
## [1] "正在进行第[238]次迭代"
## [1] "正在进行第[239]次迭代"
## [1] "正在进行第[240]次迭代"
## [1] "正在进行第[241]次迭代"
## [1] "正在进行第[242]次迭代"
## [1] "正在进行第[243]次迭代"
## [1] "正在进行第[244]次迭代"
## [1] "正在进行第[245]次迭代"
## [1] "正在进行第[246]次迭代"
## [1] "正在进行第[247]次迭代"
## [1] "正在进行第[248]次迭代"
## [1] "正在进行第[249]次迭代"
## [1] "正在进行第[250]次迭代"
## [1] "正在进行第[251]次迭代"
## [1] "正在进行第[252]次迭代"
## [1] "正在进行第[253]次迭代"
## [1] "正在进行第[254]次迭代"
## [1] "正在进行第[255]次迭代"
## [1] "正在进行第[256]次迭代"
## [1] "正在进行第[257]次迭代"
## [1] "正在进行第[258]次迭代"
## [1] "正在进行第[259]次迭代"
## [1] "正在进行第[260]次迭代"
## [1] "正在进行第[261]次迭代"
## [1] "正在进行第[262]次迭代"
## [1] "正在进行第[263]次迭代"
## [1] "正在进行第[264]次迭代"
## [1] "正在进行第[265]次迭代"
## [1] "正在进行第[266]次迭代"
## [1] "正在进行第[267]次迭代"
## [1] "正在进行第[268]次迭代"
## [1] "正在进行第[269]次迭代"
## [1] "正在进行第[270]次迭代"
## [1] "正在进行第[271]次迭代"
## [1] "正在进行第[272]次迭代"
## [1] "正在进行第[273]次迭代"
## [1] "正在进行第[274]次迭代"
## [1] "正在进行第[275]次迭代"
## [1] "正在进行第[276]次迭代"
## [1] "正在进行第[277]次迭代"
## [1] "正在进行第[278]次迭代"
## [1] "正在进行第[279]次迭代"
## [1] "正在进行第[280]次迭代"
## [1] "正在进行第[281]次迭代"
## [1] "正在进行第[282]次迭代"
## [1] "正在进行第[283]次迭代"
## [1] "正在进行第[284]次迭代"
## [1] "正在进行第[285]次迭代"
## [1] "正在进行第[286]次迭代"
## [1] "正在进行第[287]次迭代"
## [1] "正在进行第[288]次迭代"
## [1] "正在进行第[289]次迭代"
## [1] "正在进行第[290]次迭代"
## [1] "正在进行第[291]次迭代"
## [1] "正在进行第[292]次迭代"
## [1] "正在进行第[293]次迭代"
## [1] "正在进行第[294]次迭代"
## [1] "正在进行第[295]次迭代"
## [1] "正在进行第[296]次迭代"
## [1] "正在进行第[297]次迭代"
## [1] "正在进行第[298]次迭代"
## [1] "正在进行第[299]次迭代"
## [1] "正在进行第[300]次迭代"
## [1] "正在进行第[301]次迭代"
## [1] "正在进行第[302]次迭代"
## [1] "正在进行第[303]次迭代"
## [1] "正在进行第[304]次迭代"
## [1] "正在进行第[305]次迭代"
## [1] "正在进行第[306]次迭代"
## [1] "正在进行第[307]次迭代"
## [1] "正在进行第[308]次迭代"
## [1] "正在进行第[309]次迭代"
## [1] "正在进行第[310]次迭代"
## [1] "正在进行第[311]次迭代"
## [1] "正在进行第[312]次迭代"
## [1] "正在进行第[313]次迭代"
## [1] "正在进行第[314]次迭代"
## [1] "正在进行第[315]次迭代"
## [1] "正在进行第[316]次迭代"
## [1] "正在进行第[317]次迭代"
## [1] "正在进行第[318]次迭代"
## [1] "正在进行第[319]次迭代"
## [1] "正在进行第[320]次迭代"
## [1] "正在进行第[321]次迭代"
## [1] "正在进行第[322]次迭代"
## [1] "正在进行第[323]次迭代"
## [1] "正在进行第[324]次迭代"
## [1] "正在进行第[325]次迭代"
## [1] "正在进行第[326]次迭代"
## [1] "正在进行第[327]次迭代"
## [1] "正在进行第[328]次迭代"
## [1] "正在进行第[329]次迭代"
## [1] "正在进行第[330]次迭代"
## [1] "正在进行第[331]次迭代"
## [1] "正在进行第[332]次迭代"
## [1] "正在进行第[333]次迭代"
## [1] "正在进行第[334]次迭代"
## [1] "正在进行第[335]次迭代"
## [1] "正在进行第[336]次迭代"
## [1] "正在进行第[337]次迭代"
## [1] "正在进行第[338]次迭代"
## [1] "正在进行第[339]次迭代"
## [1] "正在进行第[340]次迭代"
## [1] "正在进行第[341]次迭代"
## [1] "正在进行第[342]次迭代"
## [1] "正在进行第[343]次迭代"
## [1] "正在进行第[344]次迭代"
## [1] "正在进行第[345]次迭代"
## [1] "正在进行第[346]次迭代"
## [1] "正在进行第[347]次迭代"
## [1] "正在进行第[348]次迭代"
## [1] "正在进行第[349]次迭代"
## [1] "正在进行第[350]次迭代"
## [1] "正在进行第[351]次迭代"
## [1] "正在进行第[352]次迭代"
## [1] "正在进行第[353]次迭代"
## [1] "正在进行第[354]次迭代"
## [1] "正在进行第[355]次迭代"
## [1] "正在进行第[356]次迭代"
## [1] "正在进行第[357]次迭代"
## [1] "正在进行第[358]次迭代"
## [1] "正在进行第[359]次迭代"
## [1] "正在进行第[360]次迭代"
## [1] "正在进行第[361]次迭代"
## [1] "正在进行第[362]次迭代"
## [1] "正在进行第[363]次迭代"
## [1] "正在进行第[364]次迭代"
## [1] "正在进行第[365]次迭代"
## [1] "正在进行第[366]次迭代"
## [1] "正在进行第[367]次迭代"
## [1] "正在进行第[368]次迭代"
## [1] "正在进行第[369]次迭代"
## [1] "正在进行第[370]次迭代"
## [1] "正在进行第[371]次迭代"
## [1] "正在进行第[372]次迭代"
## [1] "正在进行第[373]次迭代"
## [1] "正在进行第[374]次迭代"
## [1] "正在进行第[375]次迭代"
## [1] "正在进行第[376]次迭代"
## [1] "正在进行第[377]次迭代"
## [1] "正在进行第[378]次迭代"
## [1] "正在进行第[379]次迭代"
## [1] "正在进行第[380]次迭代"
## [1] "正在进行第[381]次迭代"
## [1] "正在进行第[382]次迭代"
## [1] "正在进行第[383]次迭代"
## [1] "正在进行第[384]次迭代"
## [1] "正在进行第[385]次迭代"
## [1] "正在进行第[386]次迭代"
## [1] "正在进行第[387]次迭代"
## [1] "正在进行第[388]次迭代"
## [1] "正在进行第[389]次迭代"
## [1] "正在进行第[390]次迭代"
## [1] "正在进行第[391]次迭代"
## [1] "正在进行第[392]次迭代"
## [1] "正在进行第[393]次迭代"
## [1] "正在进行第[394]次迭代"
## [1] "正在进行第[395]次迭代"
## [1] "正在进行第[396]次迭代"
## [1] "正在进行第[397]次迭代"
## [1] "正在进行第[398]次迭代"
## [1] "正在进行第[399]次迭代"
## [1] "正在进行第[400]次迭代"
## [1] "正在进行第[401]次迭代"
## [1] "正在进行第[402]次迭代"
## [1] "正在进行第[403]次迭代"
## [1] "正在进行第[404]次迭代"
## [1] "正在进行第[405]次迭代"
## [1] "正在进行第[406]次迭代"
## [1] "正在进行第[407]次迭代"
## [1] "正在进行第[408]次迭代"
## [1] "正在进行第[409]次迭代"
## [1] "正在进行第[410]次迭代"
## [1] "正在进行第[411]次迭代"
## [1] "正在进行第[412]次迭代"
## [1] "正在进行第[413]次迭代"
## [1] "正在进行第[414]次迭代"
## [1] "正在进行第[415]次迭代"
## [1] "正在进行第[416]次迭代"
## [1] "正在进行第[417]次迭代"
## [1] "正在进行第[418]次迭代"
## [1] "正在进行第[419]次迭代"
## [1] "正在进行第[420]次迭代"
## [1] "正在进行第[421]次迭代"
## [1] "正在进行第[422]次迭代"
## [1] "正在进行第[423]次迭代"
## [1] "正在进行第[424]次迭代"
## [1] "正在进行第[425]次迭代"
## [1] "正在进行第[426]次迭代"
## [1] "正在进行第[427]次迭代"
## [1] "正在进行第[428]次迭代"
## [1] "正在进行第[429]次迭代"
## [1] "正在进行第[430]次迭代"
## [1] "正在进行第[431]次迭代"
## [1] "正在进行第[432]次迭代"
## [1] "正在进行第[433]次迭代"
## [1] "正在进行第[434]次迭代"
## [1] "正在进行第[435]次迭代"
## [1] "正在进行第[436]次迭代"
## [1] "正在进行第[437]次迭代"
## [1] "正在进行第[438]次迭代"
## [1] "正在进行第[439]次迭代"
## [1] "正在进行第[440]次迭代"
## [1] "正在进行第[441]次迭代"
## [1] "正在进行第[442]次迭代"
## [1] "正在进行第[443]次迭代"
## [1] "正在进行第[444]次迭代"
## [1] "正在进行第[445]次迭代"
## [1] "正在进行第[446]次迭代"
## [1] "正在进行第[447]次迭代"
## [1] "正在进行第[448]次迭代"
## [1] "正在进行第[449]次迭代"
## [1] "正在进行第[450]次迭代"
## [1] "正在进行第[451]次迭代"
## [1] "正在进行第[452]次迭代"
## [1] "正在进行第[453]次迭代"
## [1] "正在进行第[454]次迭代"
## [1] "正在进行第[455]次迭代"
## [1] "正在进行第[456]次迭代"
## [1] "正在进行第[457]次迭代"
## [1] "正在进行第[458]次迭代"
## [1] "正在进行第[459]次迭代"
## [1] "正在进行第[460]次迭代"
## [1] "正在进行第[461]次迭代"
## [1] "正在进行第[462]次迭代"
## [1] "正在进行第[463]次迭代"
## [1] "正在进行第[464]次迭代"
## [1] "正在进行第[465]次迭代"
## [1] "正在进行第[466]次迭代"
## [1] "正在进行第[467]次迭代"
## [1] "正在进行第[468]次迭代"
## [1] "正在进行第[469]次迭代"
## [1] "达到精度要求"
## [1] "花费时间： 20.644 秒"截距项的分组结果如下所示
## # A tibble: 4 × 3
## label alpha size
## <dbl> <dbl> <int>
## 1 1 -1.39 44
## 2 2 0.995 51
## 3 3 -3.33 3
## 4 4 2.96 2References
[32]
MA S. HUANG J. A Concave Pairwise Fusion Approach to Subgroup Analysis[J/OL]. Journal of the American Statistical Association, 2017, 112(517): 410-423. DOI:10.1080/01621459.2016.1148039.