Mamba413
commented
1 year ago @tomwenseleers , thanks for your questions!
I have reproduced the results based on your pasted code, and I notice that the sample size may be insufficient for recovering such a little bit dense coefficient matrix that has 50 non-zero rows. According to A Polynomial Algorithm for Best-Subset Selection Problem, the recommended maximum support size is $n / (\log(p) * \log(\log(n)))$, i.e., around 44 in your setting. Therefore, recovering 50 non-zero rows is too difficult and exceeds the informatic theoretical limit.
A direct solution is collecting more samples. In the following example, 3 times sample size is used and abess returns an exactly correct result.
library(abess)
# simulate blurred multichannel spike train
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=3*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,
seed = 123
)
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)
# best subset multitask learning using family="mgaussian" using abess
abess_fit <- abess(x = X,
y = Y,
family = "mgaussian",
tune.path = "sequence",
lambda=0,
warm.start = TRUE,
tune.type = "gic") # or cv or ebic or bic or aic
plot(abess_fit, type="tune")
beta_abess = as.matrix(extract(abess_fit)$beta) # coefficient matrix for best subset
sum(rowMax(sim$beta_true)!=0) # 50 true peaks
sum(rowMin(beta_abess)!=0) # 50 peaks detected
real <- as.numeric(rowMax(sim$beta_true)!=0)
est <- as.numeric(rowMin(beta_abess)!=0)
mean(real == est)So, abess can handle nonnegativity regression coefficients. I think your question is still insightful, and I believe the inclusion of nonnegativity (or box) constraints constraint can reduce the sample size for correctly solving this problem. We are now developing a related library to support a nonnegativity constraint. Once this project is mature, we will be glad to inform you if you are interested in this.
As for the question about multitask identity link Poisson models, I conduct a simple experiment to check your idea:
library(abess)
# simulate blurred multichannel spike train
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=3*p,
p=p,
k=round(s*p), # true support size = 0.1*500 = 50
mean_beta = 10000,
sd_logbeta = 1,
family="poisson",
sd_noise = sd_noise,
multichannel=TRUE, sparse=TRUE,
seed = 123
)
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))
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
# best subset multitask learning using family="mgaussian" using abess
abess_fit <- abess(x = X,
y = Y,
weights = W,
family = "mgaussian",
tune.path = "sequence",
# support.size = c(1:50),
lambda=0,
warm.start = TRUE,
tune.type = "gic") # or cv or ebic or bic or aic
plot(abess_fit, type="tune")
beta_abess = as.matrix(extract(abess_fit)$beta) # coefficient matrix for best subset
sum(rowMax(sim$beta_true)!=0) # 50 true peaks
sum(rowMin(beta_abess)!=0) # 51 peaks detected
real <- as.numeric(rowMax(sim$beta_true)!=0)
est <- as.numeric(rowMin(beta_abess)!=0)
mean(real == est) # 0.998The estimated result is very close to the true one.
Feel free to contact us if you still have any questions.
tomwenseleers
bbayukari
Describe the bug I was testing abess on 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 (support size of best subset = 157, while true support size of simulated dataset is 50). Is there a way to resolve this? Does this have to do with the fact that abess currently does not support nonnegativity constraints (all my simulated model coefficients are positive)? Is there any way to allow for nonnegativity or box constraints? Even if not taken into account during fitting it could be handy if the coefficients could be clipped within the allowed range, so that the information criteria would at least indicate the correct support size of the best allowed model. Taking into account the constraints during fitting would of course be even better. I was also interested in fitting a multitask identity link Poisson model - as currently only Poisson log link is incorporated I could potentially still get something close to that by using an observation weights matrix 1/(Y+0.1), which would be approx. 1/Poisson variance weights and then using family="mgaussian". Is that allowed by any chance? Can observation weights be a matrix with family="mgaussian", to allow for different observation weights per outcome channel/task? If not, could this be allowed for?
Code for Reproduction
Paste your code for reproducing the bug:
Expected behavior
Couple of questions here:
Desktop (please complete the following information):
Screenshots
