1
$\begingroup$

I was using a multivariate gaussian (mgaussian) glmnet model to solve the multitask learning problem below (deconvolution of a multichannel signal using a known point spread function / blur kernel by regressing the multichannel signal on shifted copies of the point spread function) but I am experiencing poor performance, where the best subset is inferred to be much larger than it actually is (moderately larger if I set standardize.response to FALSE, much larger if I set standardize.response to TRUE). The solution is also really sensitive to the value I use for thresh (e.g. with 1E-10 I get a completely different, much less sparse & much worse solution). I am not sure in my case if I should set standardize.response to TRUE or FALSE - right now my data of the different channels / tasks can be somewhat different in scale, but despite these differences in scale, the data is simulated to all have the same Gaussian noise SD for the different channels / tasks. Standardizing it would make the scale of the different channels more comparable but would also lead to the noise having different variances for different channels, which in the fitted model would not be taken into account. So I gather I should set standardize.response to FALSE. Doing so, however, leads to model that is not sparse enough & has too large a support size compared to the true support size. Using lambda.1se as lambda value improves things a bit, but the result is still not too great. Is this expected with an L2/L1 block norm penalty (group LASSO penalty) as used here (the real sparser best subset solution I suppose one would be interested in here would use an L0/Linfinity block norm penalty)? The question also is whether I could use some information theoretical criterion perhaps to better identify the lambda that corresponds most closely to the sparsest best subset solution & would recover the true support best (e.g. BIC, EBIC or GIC). Finally, am I also correct that glmnet right now does not allow observation weights to be a matrix to allow different observation weights for different channels / tasks (I would be interested to use 1/(Y+0.1) as approximate 1/Poisson variance weights to do multitask learning for identity link Poisson with nonnegativity constraints on the coefficients)?

Example code:

library(remotes)
remotes::install_github("tomwenseleers/L0glm/L0glm")
library(L0glm)
# simulate blurred multichannel spike train
set.seed(1)
s <- 0.1 # sparsity (% of timepoints where there is a peak)
p <- 500 # 500 variables
# simulate multichannel blurred spike train with Gaussian noise
sd_noise <- 1
sim <- simulate_spike_train(n=p, 
                            p=p, 
                            k=round(s*p), # true support size = 0.1*500 = 50
                            mean_beta = 10000,
                            sd_logbeta = 1,
                            family="gaussian",
                            sd_noise = sd_noise,
                            multichannel=TRUE, sparse=TRUE)
X <- sim$X # covariate matrix with shifted copies of point spread function, n x p matrix
Y <- sim$y # multichannel signal (blurred spike train), n x m matrix  
colnames(X) = paste0("x", 1:ncol(X)) # NOTE: if colnames of X and Y are not set abess gives an error message, maybe fix this?
colnames(Y) = paste0("y", 1:ncol(Y))
true_coefs <- sim$beta_true # true coefficients
m <- ncol(Y) # nr of tasks
n <- nrow(X) # nr of observations
p <- ncol(X) # nr of independent variables (shifted copies of point spread functions)
W <- 1/(Y+0.1) # approx 1/variance Poisson observation weights with family="poisson", n x m matrix

cvfit <- cv.glmnet(X, Y, 
                   alpha = 1, # LASSO 
                   family = "mgaussian", # group LASSO = model with L1/L2 block norm penalty
                   nlambda = 100,
                   nfolds = 5, 
                   standardize = FALSE,
                   standardize.response = FALSE,
                   intercept = FALSE, 
                   relax = FALSE,
                   lower.limits = rep(0, ncol(X)+1)) # impose nonnegativity constraints on coefficients

plot(cvfit)

fit_glmnet <- glmnet(X, 
                     Y,
                     alpha = 1,  # LASSO
                     family = "mgaussian", # group LASSO = model with L1/L2 block norm penalty
                     nlambda = 100,
                     standardize = FALSE,
                     standardize.response = FALSE, 
                     intercept = FALSE, 
                     relax = FALSE, 
                     lower.limits = rep(0, ncol(X)+1)) # for nonnegativity constraints

# best lambda - I am using cvfit$lambda.1se instead of cvfit$lambda.min to get a slightly sparser model, a bit closer to the ground truth in terms of support
best_lambda <- cvfit$lambda

# get the coefficients for each task for the best lambda value
coefs <- coef(fit_glmnet, s = best_lambda)
beta_mgaussian_glmnet <- do.call(cbind, lapply(seq_len(m), function (channel) as.matrix(coefs[[channel]][-1,,drop=F])))
beta_mgaussian_glmnet[abs(beta_mgaussian_glmnet)<0.01] = 0 # slight amount of thresholding of small coefficients

image(x=1:nrow(sim$y), y=1:ncol(sim$y), z=beta_mgaussian_glmnet^0.01, col = topo.colors(255),
      useRaster=TRUE,
      xlab="Time", ylab="Channel", main="nonnegative group Lasso glmnet (red=true support)")
abline(v=(1:nrow(sim$X))[as.vector(rowMax(sim$beta_true)!=0)], col="red")
# abline(v=(1:nrow(sim$X))[as.vector(rowMin(beta_mgaussian_glmnet)!=0)], col="cyan")
sum(rowMax(sim$beta_true)!=0) # 50 true peaks
sum(rowMin(beta_mgaussian_glmnet)!=0) # 72 peaks detected - too large...

enter image description here

enter image description here

$\endgroup$

1 Answer 1

1
$\begingroup$

Seems that doing best subset selection on the variables selected by the nonnegative group LASSO glmnet model using the abess package gives me a very good solution, with near optimal support:

library(abess)
preselvars = which(rowMin(beta_mgaussian_glmnet)!=0) # pre-select variables based on nonnegative mgaussian group LASSO
system.time(abess_fit <- abess(x = X[,preselvars,drop=F], 
                               y = Y, 
                               family = "mgaussian",
                               tune.path = "sequence",
                               support.size = c(1:length(preselvars)),
                               lambda=0, 
                               max.splicing.iter = 20,
                               warm.start = TRUE,
                               tune.type = "ebic") # or cv or bic or aic or gic
            ) # 0.29s
extract(abess_fit)$support.size # 55
plot(abess_fit, type="tune")
            
beta_abess = beta_mgaussian_glmnet
beta_abess[preselvars,] = as.matrix(extract(abess_fit)$beta) # coefficient matrix for best subset
beta_abess[beta_abess < 0] = 0
image(x=1:nrow(sim$y), y=1:ncol(sim$y), z=beta_abess^0.1, col = topo.colors(255),
      useRaster=TRUE,
      xlab="Time", ylab="Channel", main="abess multitask mgaussian (red=true peaks)")
abline(v=(1:nrow(sim$X))[as.vector(rowMax(sim$beta_true)!=0)], col="red")
abline(v=(1:nrow(sim$X))[as.vector(rowMin(beta_abess)!=0)], col="cyan")
sum(rowMax(sim$beta_true)!=0) # 50 true peaks
sum(rowMin(beta_abess)!=0) # 55 peaks detected
length(intersect(which(rowMax(sim$beta_true)!=0), 
                 which(rowMin(beta_abess)!=0))) # 45 out of 50 peaks detected

Using abess mgaussian multitask best subset learning directly on the original data does not work well, however, probably because of the fact that abess does not allow for nonnegativity constraints on the coefficients.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.