r - 稳健回归中的 MM 估计

Question

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

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

我等了 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)

标签： rmachine-learningregression

Answer 1

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

• 你有没有可能内存不足？内存限制会使事情运行得更慢
• 如果问题只是事情太慢，如果您可以访问多个内核，您可以轻松地跨折叠并行化（可能不要在笔记本电脑上这样做，您会烧毁它）
• 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")