Unable to convert projected model to matrix in develop

[andrewdigby]

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))

[AlejandroCatalina]

Oh! Just to be clear, this fails when fitting the model with l1 search? I believe I know the issue. The model is correct and the search is correct too, the issue simply is transforming the fit object to matrix and that may happen because the underlying class of the su model is not “lm” but other that I didn’t implement a method for. If you upload your projection I can debug it myself as it can be tricky! Thank you for your feedback! Best, Alejandro

…

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.

[andrewdigby]

Thanks, Alejandro.

L1 search in varsel fails with a different error:

> vs <- varsel(reffull, method='L1')
Error in Reduce(function(current, pattern) { : subscript out of bounds

The forward search in varsel works, but then fails at as.matrix(proj):

vs <- varsel(reffull, method='forward')
proj <- project(vs, nterms = 3


, ns = 500)
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 think you're right about the model type, as when I add random terms as.matrix(proj) works (but this model gives poor varsel plots, but that's for another post!).

Thanks very much for your help. I've tried to upload the proj object but it's huge (~120MB). Will find out how to reduce it.

[andrewdigby]

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)

[AlejandroCatalina]

I see! I have now pushed a hotfix for both issues (projecting with `method = “forward”` and calling `varsel(…, method = “l1”)`). We can now talk about those weird varsel plots when including group terms. An important thing to consider when running `projpred` for multilevel models is that posterior identifiability becomes more of an issue that can mess up with the projections a bit. I’ll try to explain myself. When you have a somehow vague posterior in a linear model (GLM), `projpred` will try to deal with it and it might return a suboptimal projection (the ELPD of the projected model gets close but not exactly to the reference performance). No worries, running `varsel(…, ndraws_pred = 400)` will almost surely fix it. The extra argument is indicating `projpred` to basically project more draws from the posterior and account for the extra uncertainty that the default 5 draws can’t capture. While the same principle applies to multilevel models (running more draws will definitely help), now it’s easier to have weird geometrical shapes in the posterior, and some combination of such parameters will more easily end up in a singular fit, which, for our purposes, basically means a possibly suboptimal projection. With this in mind, it’s important to work a bit harder on the reference model so that the posterior is more easily identifiable, i.e. narrower marginal posteriors for the parameters, so that singular fits are less likely. This can happen even if the standard MCMC diagnostics are fine for the reference model, it does not need to be super bad, just vague enough so that the projection is singular. This is, for instance, very common in correlated posteriors, nothing wrong with them, but identifiability is much harder. Of course, the issue in your case may be different, but most of the “issues” we have had with multilevel models somehow fit in this kind of situation, so that might help you already! Best, Alejandro

…

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.

[andrewdigby]

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....!

[AlejandroCatalina]

Yes projecting 400 draws is usually a lot, depends on the specific model, and even more expensive in non gaussian cases. I’ll take a look I must have missed that case! Best, Alejandro

…

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.

[andrewdigby]

Thanks, Alejandro.

Even with 400 draws, I'm still getting a flat varsel plot which doesn't get near the reference when I add the random terms. 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(fit, method='forward', ndraws_pred=400)  

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'd be grateful for any suggestions as to what's going wrong here, or tips for diagnosis.

It'd be great if you can fix the as.matrix(proj) issue, so I can look at the projected posterior.

Thanks again for your help.

[AlejandroCatalina]

Okay that’s interesting. What’s happening is that the random effects are capturing the whole variance of the outcome and projpred gets stuck there. We are aware of this phenomenon and have it on the pipeline next! It’s not that projpred is wrong, it’s that we have to constrain the projections slightly different to enforce them to consider the population structure as well. I hope this clarifies the issue!1 Best, Alejandro

…

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.

[andrewdigby]

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

[AlejandroCatalina]

It’s not inmediate. The short term goal is to merge develop into master when I’ve fixed the l1 projection bug you reported before August. Then, we’ll work on setting the next achievable set of goals for a minor release, as well as next major goals! This issue is among the next highest priority things, but we don’t have a clear priority list for after the merge yet. Best, Alejandro

…

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.

[AlejandroCatalina]

By the way, I just tested that running `as.matrix(proj)`, with `proj` being the projection from a `varsel(…, method=“l1”)` works fine! Did you test the latest develop? Best, Alejandro

…

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.

[andrewdigby]

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’

[AlejandroCatalina]

Yes of course l1 search is not available for multilevel models. I was referring to the previous fixed effects model. Best, Alejandro

…

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.