CODEWITHSUNDEEP
rrcppglmnetmethod-call
mert
mert

Reputation: 371

Call glmnet from Rcpp (Armadillo)

I would like to extract glmnet's coefficient estimates (after cross-validation) in Rcpp Armadillo to use them in another function of Armadillo. I searched for a similar question but I could not find a solution.

I attach my attempt. (Not working) Even if I would get the list result of cv.glmnet, I could not use the coef function to obtain the coefficients.

R code

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

Arguments of 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++ code

#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(...)) ???
}

Upvotes: 2

Views: 466

Answers (1)

Ralf Stubner
Ralf Stubner

Reputation: 26833

Calling functions like cv.glmnet from C++ is complicated (maybe even impossible) since it uses quite a few possibilities offered by R that make the function signature very flexible. However, one can define a wrapper function in R that uses the actually used signature. Instead of getting this function from the (global) environment, I prefer handing it in as a function parameter:

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

Created on 2019-06-12 by the reprex package (v0.3.0)

Of course, you can do more interesting things to the calculated coefficients than returning them to R. The explicit Rcpp::as is necessary since C++ has no way of knowing what type of argument the R function returns. In this case it is a sparse matrix, which can be converted to a arma::sp_mat. BTW, this looses the Dimnames of the matrix, which is why one cannot use all.equal for comparison.

Upvotes: 3

Related Questions