r - 从 Rcpp ( Armadillo ) 调用 glmnet

Question

我想提取glmnet Rcpp Armadillo 中的系数估计(交叉验证后)在 Armadillo 的另一个函数中使用它们。 我搜索了类似的问题，但找不到解决方案。

我附上我的尝试。 (不工作) 即使我会得到 cv.glmnet 的列表结果，我无法使用coef函数来获取系数。

R代码

library(glmnet)

set.seed(1)
X = matrix(rnorm(1e3 * 201), 1e3, 201)
beta = -100:100
y = X%*%beta + rnorm(1e3)
cvfit = cv.glmnet(X, y, alpha = 1)
coefs = coef(cvfit, s = "lambda.min")
coefs                                   # get these coefficients from Rcpp

cv.glmnet 的参数

args(cv.glmnet)
> function (x, y, weights, offset = NULL, lambda = NULL, type.measure = c("mse", "deviance", "class", "auc", "mae"), nfolds = 10, foldid, 
    alignment = c("lambda", "fraction"), grouped = TRUE, keep = FALSE, 
    parallel = FALSE, ...) 
NULL

C++代码

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
Rcpp::List f_cpp(const arma::mat &x, const arma::vec &y, 
                 const arma::vec &weights, 
                 const arma::vec &lambda, double alpha, 
                 int nfolds = 10){

  Rcpp::Environment pkg = Rcpp::Environment::namespace_env("glmnet");

  Rcpp::Function f_R = pkg["cv.glmnet"];

  Rcpp::Nullable<arma::vec> offset = pkg["offset"];
  Rcpp::CharacterVector type_measure = pkg["type.measure"];
  arma::vec foldid = pkg["foldid"];
  Rcpp::CharacterVector alignment = pkg["alignment"];
  bool grouped = pkg["grouped"];
  bool keep = pkg["keep"];
  bool parallel = pkg["parallel"];

  return f_R(x, y, weights, offset, lambda, 
             type_measure, nfolds, foldid, 
             alignment, grouped, keep, parallel, alpha = alpha);
  // coef(f_R(...)) ???
}

Answer 1

从 C++ 调用像 cv.glmnet 这样的函数很复杂(甚至可能是不可能的)，因为它使用了 R 提供的相当多的可能性，使函数签名非常灵活。然而，我们可以在 R 中定义一个使用实际使用的签名的包装函数。我更喜欢将其作为函数参数传递，而不是从(全局)环境获取此函数:

library(glmnet)
#> Loading required package: Matrix
#> Loading required package: foreach
#> Loaded glmnet 2.0-16

set.seed(1)
X = matrix(rnorm(1e3 * 201), 1e3, 201)
beta = -100:100
y = X%*%beta + rnorm(1e3)


# set seed since cv.glmnet uses random numbers
set.seed(1)
cvfit = cv.glmnet(X, y, alpha = 1)
coefs = coef(cvfit, s = "lambda.min")

# set seed since cv.glmnet uses random numbers
set.seed(1)
my.glmnet <- function(x, y, alpha) {
    cvfit <- cv.glmnet(x, y, alpha = alpha)
    coef(cvfit, s = "lambda.min")
}
Rcpp::cppFunction(depends = "RcppArmadillo", "
arma::sp_mat f_cpp(const arma::mat &x, const arma::vec &y, double alpha, Rcpp::Function f_R) {
    arma::sp_mat coef = Rcpp::as<arma::sp_mat>(f_R(x, y, alpha));
    return coef;
}")
coefs2 <- f_cpp(X, y, alpha = 1, my.glmnet)

all(coefs - coefs2 == 0)
#> [1] TRUE

^{由reprex package于2019年6月12日创建(v0.3.0)}

当然，你可以对计算出的系数做更多有趣的事情，而不是将它们返回给 R。显式的 Rcpp::as 是必要的，因为 C++ 无法知道 R 的参数是什么类型函数返回。在本例中，它是一个稀疏矩阵，可以转换为 arma::sp_mat。顺便说一句，这会丢失矩阵的 Dimnames，这就是为什么不能使用 all.equal 进行比较。