I sent two e-mails to maintainer("hal9001"), on Sep 15 and Oct 7, but did not receive a response. Is GitHub preferred? The gist is that hal9001 fails its checks under forthcoming version 1.6-2 of Matrix. Matrix 1.6-2 is due to be released in November, so a patched version of hal9001 should be submitted to CRAN "soon" to avoid unnecessary archival. The relevant part of the first e-mail is pasted below.

Matrix 1.6-2 is more strict about conformation of x and y in crossprod(x=<vector>, y=<Matrix>) and tcrossprod(x=<Matrix>, y=<vector>), throwing an error in exactly those cases where base R would throw an error if the Matrix were replaced with a traditional matrix. Your packages seem to have relied on the old base-incompatible behaviour, as they now fail with errors "non-conformable arguments" tracing back to crossprod or tcrossprod calls.

Please revise your packages and test that they pass their checks under the latest Matrix-devel, which you can install with

install.packages("Matrix", repos = "http://R-Forge.R-project.org")

IIRC only minimal changes are required. Do let me know if you are unsure how to adjust your code. Then I can investigate and try to provide a patch.

Request to add feature to reduce sparse basis, in terms of both the proportion of zeros and the proportion of ones

I am getting weird answers when inputting the family argument

devtools::install_github("jlstiles/Simulations")
library(hal9001)
library(Simulations)

n=200
data =gendata(n, g0=g0_linear,Q0=Q0_trig1)
X=data
X$Y=NULL
Y=data$Y
X0=X1=X
X0$A=0
X1$A=1
newdata = rbind(X,X1,X0)

time = proc.time()
halres9001 <- fit_hal(Y = Y, X = X, family = "binomial")
timehal9001 = proc.time() - time

Qk = predict(halres9001, new_data = newdata)[1:n]
Q1k = predict(halres9001, new_data = newdata)[n+1:n]
Q0k = predict(halres9001, new_data = newdata)[2*n+1:n]
esthal9001 = var(Q1k-Q0k)
esthal9001ATE = mean(Q1k-Q0k)

esthalATE
esthal9001ATE

When I tried to install hal9001 this error came up
clang: error: unsupported option '-fopenmp'
make: *** [assign.o] Error 1
ERROR: compilation failed for package ‘data.table’

• removing ‘/Library/Frameworks/R.framework/Versions/3.5/Resources/library/data.table’
• restoring previous ‘/Library/Frameworks/R.framework/Versions/3.5/Resources/library/data.table’
  Error in i.p(...) :
  (converted from warning) installation of package ‘/var/folders/ph/l_06hx0d75x66h6f7dr3s72w0000gn/T//RtmpOicTJA/file5eb3aa883a7/data.table_1.12.3.tar.gz’ had non-zero exit status

R version 3.5.3 (2019-03-11) -- "Great Truth"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)

Hi,
Am trying to install the package and am having trouble to do it. Any help ?
Asma

We need a custom implementation of iteratively re-weighted least squares (IRLS) in order to properly support the custom ("lassi") implementation of HAL when fitting GLMs other than with the Gaussian error family.

Currently, the custom implementation (fit_type = "origami") supports family = "gaussian" while cv.glmnet must be used when calling family = "binomial".

Getting predictions from a hal fit that uses the "reduce_basis" argument fails. The error message is: "Error in .local(x, y, ...) : Cholmod error 'X and/or Y have wrong dimensions' at file ../MatrixOps/cholmod_sdmult.c, line 90".
Here is a reproducible example:

set.seed(385971)
n = 100
p = 3
x <- xmat <- matrix(rbinom(n * p, 1, .3), n, p)
y <- x[, 1] * sin(x[, 2]) + rnorm(n, mean = 0, sd = 0.2)

\# fit the HAL regression two ways
hal_fit1 <- fit_hal(X = x, Y = y)
pred1 <- predict(hal_fit1, new_data = x)  # this works
pct <- 10^-12
hal_fit2 <- fit_hal(X = x, Y = y, reduce_basis = pct)
pred2 <- predict(hal_fit2, new_data = x)  # this fails for many values of pct, e.g. .1, .2

clever things that might help us be fast:

• lassi_fit_cd should accept an index and perform fitting on that subset
• make things into a core class as per #20
• benchmark further cases of n, p for how lassi is performing vs glmnet
• coordinate descent wrt cv-risk -- only do things that minimize cv-risk
• use of R code for cv_lasso might be (probably is) causing significant delays
• could consider moving CV steps into RcppParallel instead of origami
• Wainwright

clever things that might help us be fast:

• lassi_fit_cd should accept an index and perform fitting on that subset
• make things into a core class as per #20
• benchmark further cases of n, p for how lassi is performing vs glmnet
• coordinate descent wrt cv-risk -- only do things that minimize cv-risk
• use of R code for cv_lasso might be (probably is) causing significant delays
• could consider moving CV steps into RcppParallel instead of origami
• Wainwright http://jmlr.csail.mit.edu/papers/volume17/15-605/15-605.pdf

Probably not a disastrous error but would be nice if this didn't happen:

n_obs <- 1000
w <- cbind(rbinom(n_obs, 1, 0.8), rbinom(n_obs, 1, 0.5),
           rbinom(n_obs, 1, 0.4), rbinom(n_obs, 1, 0.6))
colnames(w) <-  paste("W", seq_len(ncol(w)), sep = "_")
prob_a <- plogis(rowMeans(w[,-1]))
a <- rbinom(n_obs, 1, prob_a)

fit_gn <- fit_hal(X = as.matrix(w), Y = as.numeric(a),
                  fit_type = "glmnet", family = "binomial",
                  standardize = FALSE, cv_select = FALSE,
                  return_x_basis = TRUE, return_lasso = TRUE,
                  lambda = exp(seq(-1, -10, length = 1e4)),
                  yolo = FALSE)
est_gn <- predict(fit_gn, new_data = as.matrix(w))

where the call to predict fails with a glmnet error

Error in .local(x, y, ...) :
  Cholmod error 'X and/or Y have wrong dimensions' at file ../MatrixOps/cholmod_sdmult.c, line 90

Note that this does not happen if cv_select = TRUE in the above call to fit_hal.

The current custom implementation of the Lasso only supports continuous outcomes (that is, the case of the Gaussian error family). Some degree of revision will be necessary to support other types of outcomes, the most common case being that of binary outcomes.

For the time being, those interested in using Lasso with the binary outcomes can simply use cv.glmnet and pass it the option family = "binomial".

A further step in the custom implementation of the lasso would be to more tightly integrate the CV and lasso steps so that time is not wasted in evaluating the full sequence of lambdas. This would be done by performing cross-validation on each lambda in the initial sequence one at a time.

By evaluating the CV-MSE at each lambda, we would be able to observe shifts in the CV-MSE and safely stop at a value of lambda where the CV-risk is minimized but before the full sequence of lambdas is evaluated. This saves time by avoiding fitting the model on very small values of lambda for which a disproportionate amount of time is required (since there are more coefficients to fit).

library(hal9001)
set.seed(10)

# defining a propensity score
g0 = function(W1) plogis(-.1-.5*sin(W1) - .4*(abs(W1)>1)*W1^2)

# generating pscores and outcomes for sample of 1000
n=1000
W = rnorm(n)
ps = g0(W)
A = rbinom(n, 1, ps)

# going to fit over 1/2 the data
df_train = data.frame(W = W[1:(n/2)])

# defining validation over other half, one df including A
df_valA = data.frame(A = A[(n/2+1):n], W = W[(n/2+1):n])
df_val = data.frame(W = W[(n/2+1):n])

fitg_hal = fit_hal(X = df_train, Y = A[1:(n/2)], degrees = NULL, fit_type = "glmnet",
                   n_folds = 10, use_min = TRUE, family = "binomial",
                   return_lasso = FALSE, yolo = TRUE)

# predicting with dataframe including A should give an error or better yet, 
# the right answer, instead of wrong answer
pred_halA = predict(fitg_hal, new_data = df_valA, type = 'response')
# prediction with dataframe not including A, gives right answer
pred_hal = predict(fitg_hal, new_data = df_val, type = 'response')

# big difference
pred_hal - pred_halA

# correct risk
risk_hal = mean(with(df_valA, -A*log(pred_hal)-(1-A)*log(1-pred_hal)))
risk_hal

# ridiculous risk for the wrong preds
risk_hal_wrong = mean(with(df_valA, -A*log(pred_halA)-(1-A)*log(1-pred_halA)))
risk_hal_wrong

# glm fcn gives correct preds either way
fit_glm = glm(A[1:(n/2)] ~ W, data = df_train, family = binomial)
pred_glmA = predict(fit_glm, newdata = df_valA, type = 'response')
pred_glm = predict(fit_glm, newdata = df_val, type = 'response')
pred_glmA - pred_glm

In David's hal, there is an error given when inputting df_valA. I think it should give the right answer like glm does, if the data.frame or matrix contains the right names.

Scaling of the design matrix (of basis functions) provided considerable speedups in our custom implementation of the LASSO -- a ~2x speedup, in fact.

Based on the above observation, it appears that centering of the design matrix might provide further (hopefully considerable) speedups / efficiency gains. That said, it is not immediately clear how to go about efficiently centering sparse matrices. This will require some investigation.

Currently, we're implementing the LASSO by hand, including standard gradient descent. Theory has shown that faster convergence is possible using methods that incorporate momentum. Numerical experiments have shown that Nesterov momentum can be applied to the LASSO fairly easily, and with good results (see slides 27-29 of this deck for details). @jeremyrcoyle Would it be worthwhile for us to add in this method as an option and assess whether this would improve the performance of HAL?

Ideally, example usage demonstrates the utility of the software in solving real-world problems. The examples using simulated data are good and helpful for the user, but it would be good to have examples that use this software with real, large datasets. If these examples do exist in the repo, then they should be more clearly labeled as they are not obvious. Maybe making an /examples directory would help.

Thanks for all the work on this package.

If I am reading the documentation correctly, the new default values of num_knots uses a fixed number of knots for specifying a basis irrespective of sample size. Is that correct? If so, is there a motivation for this?

The paper does not talk about competing work at all. If this paper is the first that has been written on software implementation of HAL, this should be clearly stated and then my concern about this will be ameliorated. Otherwise, the authors should more clearly highlight previous implementation of HAL. I see the van der Laan and Dudoit paper series, but it is unclear if these papers are the only prior work, if they address the same questions as this paper, etc.

Where I met the issue:

I have the following:
X: (W, A), a n*2 numeric matrix
Y: binary outcome, numeric vector with length n.

Then I fit HAL with the codes:
fit_obj <- fit_hal(X = X, Y = Y, family = "binomial",
return_x_basis = TRUE,
num_knots = hal9001:::num_knots_generator(
max_degree = 1,
smoothness_orders = 1
))

I still have interactions. First few rows of the output:

summary(fit_obj)
Summary of non-zero coefficients is based on lambda of 0.002024853

        coef                                                                term

-5.575820e+00 (Intercept)
1.710380e+00 [ I(A >= 0)(A - 0)^1 ]
-3.430264e-01 [ I(W >= -0.804)(W - -0.804)^1 ] * [ I(A >= 0.396)*(A - 0.396)^1 ]

Expand formula interface to allow nonpenalization of some terms/variables (e.g. ~ u(t) + h(.,.))

the package has been pulled from CRAN due to failure of the solaris builds, see https://cran-archive.r-project.org/web/checks/2021/2021-11-03_check_results_hal9001.html for the test error log. To resolve this, we can debug the error after install a Solaris VM, as documented at https://github.com/r-hub/solarischeck/tree/master/packer

Hi - just fitting a grouped binomial which works great by setting Y to be a 2 columned matrix of negatives and positives. To align with other model fitting functions (GLM, GAM etc.) you might want to stick to the convention of positives and negatives, i.e. columns the other way around.

Thanks for the package!

Hi,

I am tryiny to add the "SL.hal9001" function to SuperLearner() function, and there is no problem to fit such model. But when I try to do prediction with new data set, by using predict(mod, newdataset, onlyST = T), it returns me "Error in predict.hal9001(object$object, new_data = newX) : argument "newX" is missing, with no default" if the weight for hal9001 is not 0. I ma wondering how to deal with this case?

Best,

Hao

https://github.com/jeremyrcoyle/hal9001/blob/9ab56bb964dbeda31641bc4b0c5380cad3613434/tests/testthat/test-hal.R#L48-L54

Currently (as of ebf0084), the LASSO fitting step of the HAL algorithm is a call to the cv.glmnet function from the R package glmnet. The results of simple benchmarks indicate that the vast majority of time spent on execution of a call to hal9001::fit_hal is spent in the cv.glmnet step (in fact, up to 95% of the total execution time).

Re-writing the LASSO in C++ with a Rcpp wrapper is a simple solution to this problem.

It should be a priority to implement unit tests for simple/easy statistical problems, examples including

• does hal_fit return the same thing as cv.glmnet for an easy problem?
• ...
• ...

Incredibly grateful for you releasing this.

Is there any plan for a python implementation?

Regards,

We should have one, partial implementation available here.

The Eddelbuettel references should have "R" and "C++" capitalized.

Part of this JOSS review.

Asymptotically, it appears that a call to fit_hal, wrapping cv.glmnet, will lead to severely constrained fits that fail to capture a rich enough set of basis functions asymptotically. In a previous implementation, this was ameliorated by explicitly setting the lambda.min.ratio argument of glmnet to a low value (e.g., 0.01); we should consider something similar for fit_hal.

#79 introduced a regression in which support for providing predictions along a grid of lambda values instead returns a vector of apparently identical values. This is clear from the README example, which used to return:

r$> set.seed(385971)
r$> packageVersion("hal9001")
[1] ‘0.2.7’
r$> n <- 100
r$> p <- 3
r$> x <- matrix(rnorm(n * p), n, p)
r$> y <- x[, 1] * sin(x[, 2]) + rnorm(n, mean = 0, sd = 0.2)
r$> hal_fit <- fit_hal(X = x, Y = y, cv_select = FALSE)
[1] "I'm sorry, Dave. I'm afraid I can't do that."
r$> hal_fit$times
                  user.self sys.self elapsed user.child sys.child
enumerate_basis       0.003        0   0.003          0         0
design_matrix         0.006        0   0.006          0         0
reduce_basis          0.000        0   0.000          0         0
remove_duplicates     0.004        0   0.004          0         0
lasso                 0.028        0   0.028          0         0
total                 0.041        0   0.042          0         0
r$> preds <- predict(hal_fit, new_data = x)

r$> preds
               s0          s1          s2          s3          s4          s5         s6
  [1,] 0.06356236  0.06478106  0.06559040  0.06310132  0.06072537  0.05844328  0.0562756
                s7          s8          s9         s10          s11         s12
  [1,]  0.05463778  0.05627511  0.05339973  0.05226898  0.051860713  0.05198880
                s13          s14         s15         s16         s17         s18
  [1,]  0.052159151  0.052403705  0.05170359  0.05074108  0.04971674  0.04872899
               s19         s20         s21         s22          s23         s24
  [1,]  0.04718589  0.04224513  0.03713774  0.03313715  0.029578804  0.02611531
               s25         s26          s27          s28          s29          s30
  [1,]  0.02220966  0.01861265  0.016745656  0.014962968  0.014016019  0.013150862
                 s31           s32          s33          s34          s35          s36
  [1,]  0.0084326033  0.0016781764 -0.004379224 -0.009756481 -0.018208534 -0.027420644
                s37          s38          s39          s40          s41           s42
  [1,] -0.035692904 -0.042651339 -0.047731797 -0.057171332 -0.068810925 -0.0814902028
                s43           s44          s45           s46          s47          s48
  [1,] -0.094145457 -0.1063458006 -0.112258682 -0.1169236800 -0.122157452 -0.126232513
                s49          s50          s51          s52          s53          s54
  [1,] -0.128004295 -0.129774524 -0.129584798 -0.129109598 -0.128860797 -0.128431573
                s55          s56           s57           s58          s59          s60
  [1,] -0.128002734 -0.127860160 -1.276471e-01 -1.263560e-01 -0.124513722 -0.123259610
                s61           s62          s63          s64          s65          s66
  [1,] -0.122217608 -0.1212981615 -0.119173537 -0.117948634 -0.118493536 -0.118916254
               s67          s68          s69          s70          s71          s72
  [1,] -0.11921098 -0.119095745 -0.118599351 -0.116668209 -0.115855176 -0.115227862
                s73          s74          s75          s76          s77          s78
  [1,] -0.114295608 -0.113795088 -0.113272632 -0.113057443 -0.113127302 -0.112956354
               s79          s80          s81          s82          s83          s84
  [1,] -0.11290637 -0.114645294 -0.116617272 -0.118362758 -0.119785603 -0.121354928
                 s85          s86           s87          s88          s89          s90
  [1,] -0.1219591307 -0.121072044 -0.1201841224 -0.119764167 -0.119247460 -0.118394374
                s91          s92          s93          s94          s95          s96
  [1,] -0.117574887 -0.117310636 -0.116898537 -0.116803168 -0.117258712 -0.117751223
                s97         s98          s99
  [1,] -0.118306951 -0.11876637 -0.119148360
 [ reached getOption("max.print") -- omitted 99 rows ]

...but now returns

r$> set.seed(385971)
r$> packageVersion("hal9001")
[1] ‘0.3.0’
r$> n <- 100
r$> p <- 3
r$> x <- matrix(rnorm(n * p), n, p)
r$> y <- x[, 1] * sin(x[, 2]) + rnorm(n, mean = 0, sd = 0.2)
r$> hal_fit <- fit_hal(X = x, Y = y, yolo = TRUE, cv_select = FALSE)
[1] "I'm sorry, Dave. I'm afraid I can't do that."
r$> hal_fit$times
                  user.self sys.self elapsed user.child sys.child
enumerate_basis       0.188    0.000   0.187          0         0
design_matrix         0.010    0.000   0.010          0         0
reduce_basis          0.000    0.000   0.000          0         0
remove_duplicates     0.000    0.000   0.000          0         0
lasso                 0.054    0.001   0.054          0         0
total                 0.253    0.001   0.253          0         0
r$> preds <- predict(hal_fit, new_data = x)

r$> preds
  [1] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
  [8] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [15] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [22] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [29] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [36] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [43] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [50] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [57] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [64] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [71] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [78] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [85] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [92] 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236 0.06356236
 [99] 0.06356236 0.06356236
 [ reached getOption("max.print") -- omitted 9900 entries ]

The implementation of Lasso available in the lasso2 package (GitHub | CRAN) appears to allow options for fitting the Lasso model while imposing bounds on the L1 norm. This could be a highly desirable feature for HAL regressions. Writing a backend/wrapper for hal9001 to incorporate the lasso2 package --- similar to the wrapper for glmnet --- would likely be worth the effort. It's possible that lasso2 does not directly allow for cross-validation, meaning some of the functionality of CV-Lasso (available in glmnet) will have to be replicated by hand. Idea originally from @benkeser.

It would be very cool to add a formula interface for fitting HAL additive models.

So if X has two columns x1, x2 you could write

fit_hal(..., formula = "x1 + x2")

and it would fit a 1-degree HAL

fit_hal(..., formula = "x1*x2")

would fit the 2-degree HAL, etc...

The number of interaction terms increases dramatically, when the number of variables increases. This makes it almost impossible to produce model results in a reasonable time when ncol(X) is large. Is it possible to change this default value ncol(X) >= 20 to something smaller?

Is it possible to fit poisson models using this package? It is not listed as a supported family in the documentation, and it doesn't work when I try it, but I wonder how hard it would be to add support for it. Is this something that requires just R code (in which case I might be interested in working on a PR), or would it require changes to Rcpp code underneath?

Thanks in advance. This package is awesome.

We're passing the same pile of variables between a pile of functions. This is why classes exist.

In this paper, Hastie et al. report significant gains in efficiency of their algorithmic implementation of the Lasso using covariance updates rather than naive updates of the coefficients. By incorporating this process into our implementation of the Lasso -- at least for the Gaussian error family case -- we ought to be able to obtain further gains in efficiency.

To achieve the quartic convergence rate dictated by theory, it is necessary to have the HAL regression function reach a saturated model asymptotically; however, this is currently not possible in the family = "binomial" implementation of glmnet(), as a call to glmnet() will ignore any specification of lambda.min.ratio in this case. This can be avoided in two ways: (1) setting family = "gaussian" and specifying a very low value of lambda.min.ratio (essentially a linear probability model in terms of HAL basis functions); or (2) directly passing in a sequence of L1 constraint parameters that includes a very low value through the lambda argument of glmnet(). Issue originally reported and confirmed by @idiazst.

Consider filling in lines 634, 635, and 655 (title, author names, year) of the license.

Part of this JOSS review.

Currently, by default, fit_hal fits the lasso on the full set of interactions between all covariates presented to it. This is extremely slow and/or impractical for even modestly large data sets. Perhaps it ought to be changed?

In their Lasso implementation in the glmnet package, Hastie et al. note that

Considerable speedup is obtained by organizing the iterations around the
active set of features—those with nonzero coefficients. After a complete cycle
through all the variables, we iterate on only the active set till convergence. If
another complete cycle does not change the active set, we are done, otherwise
the process is repeated. Active-set convergence is also mentioned in Meier
et al. [2008] and Krishnapuram and Hartemink [2005].

We could very likely speed up our implementation by incorporating this procedure.

Issue
predict.hal9001 ignores the lambda parameter and always shows predictions for all lambdas specified during fitting. Having had a cursory look at the source code, lambda does not seem to be used at all in predict.hal9001 but is used in predict.lassi.

This is not clear from the documentation, and led to confusions in my case since it makes the lambda parameter redundant and diverges from glmnet's behaviour. There might be theoretical reasons behind this but it would good to remove this parameter in this case, or state this behaviour clearly in the manual if the parameter is needed for compatibility with predict.lassi.

Reproducible example

# Simulated data taken from glmnet help page ?glmnet::glmnet
x = matrix(rnorm(100 * 20), 100, 20)
y = sample(1:2, 100, replace = TRUE)

# Fit a elastic net and predict with different lambda 
net_fit = glmnet(x, y, family = "binomial", lambda = c(0.001, 0.01, 0.1))
predict(net_fit, head(x), type = "response")  # Use lambda used during fitting
predict(net_fit, head(x), s = c(0.1, 0.05), type = "response") # Use custom lambdas
# --> glmnet correctly shows the specified lambdas, interpolating or refitting for the previously unseen lambda

# Fit a highly adaptive lasso model and try to predict with different lambda
hal_fit = fit_hal(X = x, Y = y, family = "binomial", lambda = c(0.001, 0.01, 0.1), cv_select = FALSE, fit_type = "glmnet")
predict(hal_fit, new_data = head(x), type = "response") # fitted lambdas 
predict(hal_fit, new_data = head(x), lambda = c(0.1, 0.05), type = "response") # custom lambdas
# --> hal9001 ignores the lambda parameter provided to predict, showing the same result in each case

Hi Nima,

No hurry, but if you could implement family="quasibinomial" as an option in HAL now that cvglmnet allows it, that would be great!

Thanks,
Lauren