Unable to convert projected model to matrix in develop about projpred HOT 15 CLOSED

On 4 Jul 2020, 1:39 AM +0300, andrewdigby ***@***.***>, wrote: Hi, I've been using the develop branch of projpred with some success, but have hit a snag trying to plot the projected posteriors for my own logistic regression model, with interactions and random effects. I'm unsure if this is a bug, or something wrong with my model. The command as.matrix(proj) fails when applying to my own model (fit - see below): priorh <- set_prior("horseshoe(scale_global = tau0, scale_slab=1)", class="b") cfull_formula_fixed <- as.formula( y ~ v1_cts * f1 + v2_cts * f2 + f1:f2 + f3 + v3_cts + v4_cts ) fit <- brm(cfull_formula_fixed, data=data, family="bernoulli", iter=5000, prior = priorh, sample_prior = T, control=list(adapt_delta=0.999,max_treedepth = 15), save_all_pars = F, refresh=0) # Data: data (Number of observations: 186) vs <- varsel(fit, method='forward') # fails without method='forward' (using method='L1') proj <- project(vs, nterms = 3, ns = 500) mcmc_areas(as.matrix(proj)) + coord_cartesian(xlim = c(-2, 2)) Error: Error in dimnames(x) <- dn : length of 'dimnames' [2] not equal to array extent I've been unable to diagnose the problem by applying as.matrix.projection directly, since the line: res <- t(do.call(cbind, lapply(proj$sub_fit, as.matrix.lm))) returns Error in t.default(do.call(cbind, lapply(proj$sub_fit, as.matrix.lm))) : argument is not a matrix since proj$sub_fit is a list. Apologies for not providing a reproducible example; I can upload the projected (proj) or fitted model (fit) if necessary. The logistic regression example at the bottom of this post (adapted from the vignettes) works as expected, so this must be an issue with my model. I suspect it may be because my model contains multi-level factors and interactions? Have you seen similar behaviour elsewhere? Any suggestions of how to diagnose this issue would be very much appreciated. Many thanks in advance for your help and the great package. This works: library(projpred) library(brms) library(tidyr) library(dplyr) library(ggplot2) library(bayesplot) # Binomial example data('df_binom', package = 'projpred') split_structure <- break_up_matrix_term(y ~ x, data = df_binom) df_binom <- split_structure$data formula <- split_structure$formula d <- df_binom n <- nrow(df_binom) # 100 D <- ncol(df_binom[, -1]) # 20 p0 <- 5 # prior guess for the number of relevant variables tau0 <- p0/(D-p0) * 1/sqrt(n) # scale for tau priorh <- set_prior("horseshoe(scale_global = tau0, scale_slab=1)", class="b") fit2 <- brm(formula, data=df_binom, family=bernoulli(), iter=500, chains=2, prior = priorh, sample_prior = T, control=list(adapt_delta=0.999,max_treedepth = 15), save_all_pars = TRUE, refresh=0) vs2 <- cv_varsel(fit2) plot(vs2, stats = c('elpd', 'rmse')) proj2 <- project(vs2, nterms = 3, ns = 500) mcmc_areas(as.matrix(proj2)) + coord_cartesian(xlim = c(-2, 2)) — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

On 5. Jul 2020, 3.01 +0300, andrewdigby ***@***.***>, wrote: Here's the original model object (brms_cfull in the below): this is the smallest part I can upload! Hope that's OK. Thank you. projpred_model.Rdata.zip reffull <- get_refmodel(brms_cfull) vs <- varsel(reffull, method='forward') proj <- project(vs, nterms = 3, ns = 500) — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

On 9 Jul 2020, 8:07 AM +0300, andrewdigby ***@***.***>, wrote: Hi @AlejandroCatalina. Thanks very much for the fix, and for the explanation regarding the varsel plots. I really appreciate your help with this. The varsel(..., method='L1' fix works for me, but I still get the same error with as.matrix(proj): packageVersion("projpred") [1] ‘2.0.0’ > mcmc_areas(as.matrix(proj)) + coord_cartesian(xlim = c(-2, 2)) Error in dimnames(x) <- dn : length of 'dimnames' [2] not equal to array extent I'm still trying the varsel(…, ndraws_pred = 400) fix for the model with the random effects - it's taking a very long time to run....! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

On 12 Jul 2020, 6:06 AM +0300, andrewdigby ***@***.***>, wrote: Thanks, Alejandro. I haven't got 400 draws to run yet, but have achieved 50 :-) I'm still getting a flat varsel plot which doesn't get near the reference when I add the random terms. Is that what you'd expect with insufficient draws? The order of solution_terms is also not what I expect from other model comparisons. priorh <- set_prior("horseshoe(scale_global = tau0, scale_slab=1)", class="b") cfull_formula_rand <- as.formula( y ~ v1_cts * f1 + v2_cts * f2 + f1:f2 + f3 + v3_cts + v4_cts + v5_cts + (1|r1) + (1|r2) + (1|r3) ) fit <- brm(cfull_formula_rand, data=data, family="bernoulli", iter=5000, prior = priorh, sample_prior = T, control=list(adapt_delta=0.999,max_treedepth = 15), save_all_pars = F, refresh=0) vs <- varsel(brms_cfull, method='forward', ndraws_pred=50)#ndraws_pred=400) # ndraws_pred=400 crashes/doesn’t finish. solution_terms(vs) "(1 | r1)" "(1 | r2)" "(1 | r3)” "v1_cts" “f3” "v2_cts” "v5_cts” "f2" “v2_cts:f2” "f1” “v1_cts:f1” “f1:f2” “v4_cts” plot(vs, stats=c('elpd', 'rmse'), baseline='ref') When I remove the random terms, the varsel results are as expected. solution_terms(vs) “f2” “f3” “v2_cts” “f1” “f1:f2” “v4_cts” “v1_cts” “f1:vt_cts” “f2:v2_cts” “v5_cts" I'll try to get varsel(..., ndraw_pred=400) to complete, but in the meantime if you think there's another issue here, I'd be grateful for any suggestions! Thanks again for your help. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

On 13 Jul 2020, 8:19 AM +0300, andrewdigby ***@***.***>, wrote: OK, thanks for explaining. I'll work off the fixed effect model for now, and will wait and hope that you can tweak projpred to constrain the projections. I'm guessing that's a longer-term fix? Please let me know if you need more information about this example to help with that. Thanks, Andrew — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

On 13. Jul 2020, 8.19 +0300, andrewdigby ***@***.***>, wrote: OK, thanks for explaining. I'll work off the fixed effect model for now, and will wait and hope that you can tweak projpred to constrain the projections. I'm guessing that's a longer-term fix? Please let me know if you need more information about this example to help with that. Thanks, Andrew — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

On 13 Jul 2020, 12:19 +0300, andrewdigby ***@***.***>, wrote: Thanks for the update about the plan. So I'll anticipate it being a few months off, and will find a workaround for now. Regarding the as.matrix(proj) issue, I get this error for method="L1": vs <- varsel(brms_cfull, method='L1', ndraws_pred=400) Error in stop("l1 search is not supported for multilevel models", ) : argument is missing, with no default Using method="forward" works OK. I'm using the develop version: packageVersion("projpred") [1] ‘2.0.0’ — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.