稳健回归中的 MM 估计

Question

我在 R 中使用不同的线性回归模型。我使用了 DATASET，它有 21263 行和 82 列。

除了使用 R 函数 lmrob 的 MM 估计回归之外，所有回归模型的时间消耗都可以接受。

我等待了 10 多个小时来运行第一个 for 循环 (#Block A)，但它不起作用。 “不起作用”是指两天后它可能会给我一个输出。我用较小的DATASET 尝试了这段代码，它有 9568 行，5 列，它在一分钟内运行。

我使用的是标准笔记本电脑。

我的分析步骤如下

上传和缩放数据集，然后使用 k=30 的 k-folds 拆分，因为我想计算 k 拆分中每个变量的系数方差。

您能给我提供任何指南吗？

wdbc = read.csv("train.csv") #critical_temp is the dependent varaible. 
wdbcc=as.data.frame(scale(wdbc)) # scaling the variables
### k-folds split ###
set.seed(12345)
k = 30
folds <- createFolds(wdbcc$critical_temp, k = k, list = TRUE, returnTrain = TRUE)

############ Start of MM Regression Model #################
#Block A
lmrob = list()
for (i in 1:k) {
    lmrob[[i]] = lmrob(critical_temp~ ., 
                       data = wdbcc[folds[[i]],],setting="KS2014")
}

#Block B
lmrob_coef = list()
lmrob_coef_var = list()

for(j in 1:(lmrob[[1]]$coefficients %>% length())){

    for(i in 1:k){

        lmrob_coef[[i]] = lmrob[[i]]$coefficients[j] 
        lmrob_coef_var[[j]] = lmrob_coef %>% unlist() %>% var()
    }

}

#Block C
lmrob_var = unlist(lmrob_coef_var)
lmrob_df = cbind(coefficients = lmrob[[1]]$coefficients %>% names() %>% as.data.frame()
                 , variance = lmrob_var %>% as.data.frame()) 
colnames(lmrob_df) = c("coefficients", "variance_lmrob")
#Block D
lmrob_var_sum = sum(lmrob_var)

Answer 1

不是答案，而是一些帮助您自己测试的代码。我没有在完整数据集上运行lmrob()，但我在下面显示的所有内容都表明模型的一个完整实现（所有观察结果、所有预测变量）应该在大约 10-20 分钟内运行 [在 10 年的 MacOS 桌面上machine]，这将外推到大约 5 小时进行 30 倍交叉验证。 （看起来时间尺度比观察次数的平方根差一点，并且非线性甚至在对数尺度上与预测变量的数量......）您可以尝试下面的代码看看你的机器上的事情是否慢得多，并预测你认为完成整个问题需要多长时间。其他一般性建议：

• 您是否有可能内存不足？内存限制会使事情的运行速度大大变慢
• 如果问题只是速度太慢，如果您可以访问多个内核，则可以轻松地跨折叠进行并行化（可能不要在笔记本电脑上执行此操作，否则会烧毁它）
• AWS 和其他云服务非常有用

我设置了一个测试函数来记录 lmrob() 在您数据集中的预测变量和观察值的随机子集上运行所花费的时间。

提取数据，加载包：

unzip("superconduct.zip")
xx <- read.csv("train.csv")
library(robustbase)
library(ggplot2); theme_set(theme_bw())
library(cowplot)

定义一个测试函数，用于计时lmrob 针对不同数量的观察和预测变量运行：

nc <- ncol(xx)  ## response vble is last column, "critical_temp"
test <- function(nobs=1000,npred=10,seed=NULL, ...) {
    if (!is.null(seed)) set.seed(seed)
    dd <- xx[sample(nrow(xx),size=nobs),
             c(sample(nc-1,size=npred),nc)]
    tt <- system.time(fit <- lmrob(critical_temp ~ ., data=dd, ...))
    tt[c("user.self","sys.self","elapsed")]
}    

t0 <- test()

这里的最小示例（1000 个观察值，10 个预测变量）非常快（0.2 秒）。 这是我运行的基本循环：

res <- expand.grid(nobs=seq(1000,10000,by=1000), npred=seq(10,30,by=2))
res$user.self <- res$sys.self <- res$elapsed <- NA
for (i in seq(nrow(res))) {
    cat(res$nobs[i],res$npred[i],"\n")
    res[i,-(1:2)] <- test(res$nobs[i],res$npred[i],seed=101)
}

（正如您在下图中看到的那样，我再次使用大量观察值和预测变量进行了此操作，并使用rbind() 将结果组合到单个数据框中。）我还尝试拟合线性模型来预测使用所有预测变量完成完整数据集所需的时间。 （绘图[见下文]表明时间在观察数量上是对数对数线性，但在预测变量数量上是非线性......）

m1 <- lm(log10(elapsed)~poly(log10(npred),2)*log10(nobs), data=resc)
pp <- predict(m1, newdata=data.frame(npred=ncol(xx)-1,nobs=nrow(xx)),
              interval="confidence")  
10^pp  ## convert from log10(predicted seconds) to seconds

测试完整的数据集。

t_all <- test(nobs=nrow(xx),npred=ncol(xx)-1)

然后我意识到您使用的是setting = "KS2014"（如文档中所建议的那样）而不是默认值：这至少慢了 5 倍，如以下比较所示：

test(nobs=10000,npred=30)
test(nobs=10000,npred=30,setting = "KS2014")

我用setting="KS2014" 重新运行了上面的一些内容。对完整数据集进行预测表明运行时间约为 700 秒（CI 从 300 到 2000 秒） - 仍然远没有您建议的那么慢。

gg0 <- ggplot(resc2,aes(x=npred,y=elapsed,colour=nobs,linetype=setting))+
    geom_point()+geom_line(aes(group=interaction(nobs,setting)))+
    scale_x_log10()+scale_y_log10()
gg1 <- ggplot(resc2,aes(x=nobs,y=elapsed,colour=npred, linetype=setting))+
    geom_point()+geom_line(aes(group=interaction(npred,setting)))+
    scale_x_log10()+scale_y_log10()
plot_grid(gg0,gg1,nrow=1)
ggsave("lmrob_times.pdf")