Bayesian Adaptive Regression Splines

[Kavir Narsai]

I note nimble is capable of implementing rjmcmc for the purpose of variable selection however, is it possible to use nimble to implement bayesian adaptive regression splines whereby one could randomly vary the number and placement of knots in a regression spline fit?

[Perry de Valpine]

That's an interesting question and idea. The RJMCMC system provided with nimble is pretty basic, set up for one-to-one pairs of indicators (0 or 1) and coefficients. I am not sure but off the top of my head I suspect what you would need is more involved. In that case, you could consider writing your own RJ samplers to work with this specific kind of situation.

On Mon, Jul 22, 2024 at 8:08 AM Kavir Narsai <kavir...@gmail.com> wrote:

I note nimble is capable of implementing rjmcmc for the purpose of variable selection however, is it possible to use nimble to implement bayesian adaptive regression splines whereby one could randomly vary the number and placement of knots in a regression spline fit?

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/d9457b09-d95c-470e-8c2a-ab8349b56b06n%40googlegroups.com.

[Kavir Narsai]

Hi Perry

Noted, thank you for the feedback. Yes, I suppose one of the main elements to cater for is the possibility of the size of the design matrix and associated number of coefficients potentially changing at each iteration. I see that dynamic indexing is not allowed directly in the nimbleCode but is there a way to use a nimbleFunction to enable the indexing of the coefficients and design matrix to change?

Thanks 

[Perry de Valpine]

This is somewhat open-ended and it might help if you can sketch out more about what you would like to do. Dynamic indexing is supported in nimble models, so if you want to say more of what you are seeing about that, please do. For a model that conceptually changes dimension, you will need to implement the maximum sized model and then use model code and/or nimbleFunctions to make use of subsets of it as needed. There is a new feature (as of nimble 1.2.0) that allows nimbleFunctions with setup code to be used in model code. This opens the possibility of user-defined functions and distributions that retain their own internal variables across calls (i.e. they are class objects). That might be useful for example if you need to update internal variables about a design matrix or basis functions as you (have a sampler) change the number or placement of knots. A good way to get support would be if you can abstract your problem down to a toy example that involves the main implementation difficulties you want to solve.

HTH

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/181b3e83-c207-4a3a-9174-c9e24ff0ffe4n%40googlegroups.com.

[Kavir Narsai]

Hi Perry

Below is a simple example in which I would like to fit a bayesian adaptive regression spline using nimble: (Note the nimble code below is meant to illustrate ideally what I would like to achieve but I'm not sure how to specify it correctly using nimbleFunctions:

<R code

## LOAD PACKAGES ---------
library(data.table)

library(splines2)
library(nimble)

## SIMULATE DATA ---------
set.seed(123)
dt <- data.table(x = seq(0.1, 10, 0.1))
betas <- c(20, 30, 50)

dt[, y := nsp(x, knots = 5, intercept = T)%*%betas + rnorm(100, 0, 0.5)
]

## ATTEMPT TO PREPARE BARS IN NIMBLE

## Specify nsp code in R
Build_nsp <- function(x, K) {
  nsp(x, knots = K, intercept = T)
}

## Convert R function to a function that can be called in nimble
Rnsp <- nimbleRcall(function(x = double(1), K = double(1)) {},
                    Rfun = 'Build_nsp',
                    returnType = double(2))

nim_spCode <- nimbleCode({
 
  ## Determine number of coefficients based
  ## on number of knots
  p <- KNum + 2
 
  ## Prep spline design matrix
  B[1:N,1:p] <- Rnsp(x[1:N], K = K)
 
  for (i in 1:N) {
    fit[i] <- inprod(B[i, 1:p], betas[1:p])
  }
 
  ## Simulate the number of knots
  ## such that it is very dependent on
  ## the data
  Alpha[1:Kmax] <- nimRep(0.01, Kmax)
  Theta[1:Kmax] ~ ddirch(Alpha[1:Kmax])
  KNum ~ dcat(Theta[1:Kmax])
 
  ## Attempt to make the knot position random
  for (i in 1:KNum){
    K[i] ~ dunif(0.1, 9.9)
  }
 
  for (i in 1:p) {
    betas[i] ~ dnorm(0, sd = 1)
  }
 
  for (i in 1:N) {
    y[i] ~ dnorm(fit[i], sd = 1)
  }
 
}
)


y_m  <- mean(dt$y)
y_sd <- sd(dt$y)

sp_Const <- list(N = dim(dt)[1], Kmax = 3)
sp_Data <- list(y = (dt$y - y_m) / y_sd,
                x = dt$x)
sp_Inits <- list(betas = rep(0, 3),  KNum = 2 )


sp_mod <- nimbleModel(code      = nim_spCode  ,
                      name      = "sp_mod"    ,
                      constants = sp_Const    ,
                      data      = sp_Data     ,
                      inits     = sp_Inits   )

> 
The error received when running the nimbleModel is as follows: "Defining model

Error in getSymbolicParentNodesRecurse(x, constNames, indexNames, nimbleFunctionNames,  :

  Dynamic indexing found in a vector of indices, 1:p. Only scalar indices, such as 'idx' in 'x[idx]', can be dynamic. One can instead use dynamic indexing in a vector of indices inside a nimbleFunction."

I have been able to bring in a natural cubic spline (provided by the splines2 package) into the nimble code using the nimbleRcall. So it is possible to make a given number of knots random. However, I'm not sure how to enable the number of knots (KNum in the R code above) to vary as well.

Thank you for your assistance Perry.

Best regards

Kavir 

[Daniel Turek]

Kavir, this is great work.  I think you only need to push a few things into nimbleFunctions, and the problem will resolve.

Please take a look at the code below, I'll highlight the changes so you can see them.  I think I got the indexing right, but please double-check that the maths look correct.

Cheers,

Daniel

## Specify nsp code in R
Build_nsp <- function(x, K) {
  nsp(x, knots = K, intercept = T)
}

## Convert R function to a function that can be called in nimble
Rnsp <- nimbleRcall(function(x = double(1), K = double(1)) {},
                    Rfun = 'Build_nsp',
                    returnType = double(2))

calc_fit <- nimbleFunction(
    run = function(x = double(1), K = double(1), betas = double(1), KNum = double(), N = double()) {

        p <- KNum + 2

        B <- Rnsp(x[1:N], K = K[1:KNum])
        fitArray <- B %*% betas[1:p]
        returnType(double(1))
        return(fitArray[1:N,1])
    }
)

nim_spCode <- nimbleCode({
   
    ## Determine number of coefficients based
    ## on number of knots

    ###p <- KNum + 2   ## remove p

   
    ## Prep spline design matrix

    ##B[1:N,1:p] <- Rnsp(x[1:N], K = K)     ## remove B
   
    ###for (i in 1:N) {
    ###    fit[i] <- inprod(B[i, 1:p], betas[1:p])    ## remove
    ###}
    f[1:N] <- calc_fit(x[1:N], K[1:Kmax], betas[1:(Kmax+2)], KNum, N)    ## new

   
    ## Simulate the number of knots
    ## such that it is very dependent on
    ## the data
    Alpha[1:Kmax] <- nimRep(0.01, Kmax)
    Theta[1:Kmax] ~ ddirch(Alpha[1:Kmax])
    KNum ~ dcat(Theta[1:Kmax])
   
    ## Attempt to make the knot position random

    ##for (i in 1:KNum){    ## remove
    for (i in 1:Kmax){      ## instead, define K[i] up to max size of Kmax

        K[i] ~ dunif(0.1, 9.9)
    }
   

    ##for (i in 1:p) {            ## remove
    for (i in 1:(Kmax+2)) {       ## instead, define betas[i] up to max size of Kmax+2

        betas[i] ~ dnorm(0, sd = 1)
    }
   
    for (i in 1:N) {
        y[i] ~ dnorm(fit[i], sd = 1)
    }
   
}
)

y_m  <- mean(dt$y)
y_sd <- sd(dt$y)

sp_Const <- list(N = dim(dt)[1], Kmax = 3)
sp_Data <- list(y = (dt$y - y_m) / y_sd, x = dt$x)

sp_Inits <- list(betas = rep(0,5),  KNum = 2,  ## change initial values for betas to be length = 5
                 K = runif(3, 0.1, 9.9))       ## also add initial values for K[1:Kmax]

sp_mod <- nimbleModel(code      = nim_spCode  ,
                      name      = "sp_mod"    ,
                      constants = sp_Const    ,
                      data      = sp_Data     ,
                      inits     = sp_Inits   )

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/5855562c-618f-443b-83a8-97b706662cadn%40googlegroups.com.

[Perry de Valpine]

And to follow up further:

If you want to be sampling knots by position and also by RJ in and out of the model, you'll have to consider whether the existing RJ is satisfactory or not. You could use it for the beta[i] elements, but then the K[i] would still be mixing and following their prior if not being used in the model. Alternatively you could look at the RJ and make one that jointly samples both beta[i] and K[i] in and out of the model.

And if you want to separate computation of the nsp from computation of the fitArray, you could split those steps once you verify if Daniel's changes do what you want. I am guessing the nsp is a bit costly, and as written that will be done even if one of the beta[i] but not the K[i] is changed by a sampler.

Perry

[Kavir Narsai]

Hi Daniel and Perry

Thank you so much for your assistance with this.

@Daniel, your updates to the code has enabled the model to work. Out of interest I've tested it against running each of the 3 possible models separately (i.e. explicitly specifying 3 different spline fits with each having 1 to 3 knots respectively) and then weighting each fit using the Dirichlet distribution (note that it is not an efficient implementation but simply used for testing). The code used is provided below. The fit of the two implementations are fairly similar. I'd assume the reason for the differences could be due to the Dirichlet implementation (as apposed to the BARS) specifying a separate set of coefficients for each model based on the number of knots.

@Perry your point is well noted and I'll try to consider a RJ implementation of the model. To confirm, my understanding is correct, the main difference between using the 'standard MCMC algorithm' to run the model Vs. a RJ sampler is that the RJ sampler can be more efficient by only sampling the required knot locations (K) and betas depending on the number of knots (KNum) in each iteration?

Thank you.

Best regards

Kavir

### R Code ###

## LOAD PACKAGES ---------
library(data.table)

library(plotly)

library(nimble)


## SIMULATE DATA ---------

set.seed(123)
dt <- data.table(x = seq(0.1, 10, 0.1))
betas <- c(20, 30, 50)

dt[, y := nsp(x, knots = 5, intercept = T)%*%betas + rnorm(100, 0, 0.5)
]

## ATTEMPT TO PREPARE BARS IN NIMBLE -----------------

## Convert R function to a function that can be called in nimble
Rnsp <- nimbleRcall(function(x = double(1), K = double(1)) {},
                    Rfun = 'Build_nsp',
                    returnType = double(2))

calc_fit <- nimbleFunction(

  run = function(x = double(1), K = double(1), betas = double(1), KNum = integer(0)) {
    p <- KNum +2
    B <- Rnsp(x, K[1:KNum])

    fitArray <- B %*% betas[1:p]
    returnType(double(1))

    return(fitArray[, 1])
  }
)


nim_BarsCode <- nimbleCode({
 
  ## Using the nimble function we are able to
  ## randomly subset the coefficients and design matrix
  ## as required to generate the desired fit based on the
  ## randomly selected number of knots
  fit[1:N] <- calc_fit(x[1:N], K[1:Kmax], betas[1:(Kmax+2)], KNum)

 
  ## Simulate the number of knots
  ## such that it is very dependent on
  ## the data
  Alpha[1:Kmax] <- nimRep(0.01, Kmax)
  Theta[1:Kmax] ~ ddirch(Alpha[1:Kmax])
  KNum ~ dcat(Theta[1:Kmax])
 

  ## Note here for the knot locations and
  ## coefficients the model always simulates the maximum possible
  ## number to ensure that the dimension of these simulations do not have
  ## to change. The nimble function then caters for the the number of
  ## knots to be used in the actual fit

  for (i in 1:Kmax){

    K[i] ~ dunif(0.1, 9.9)
  }
 

  for (i in 1:(Kmax+2)) {

    betas[i] ~ dnorm(0, sd = 1)
  }
 
  for (i in 1:N) {
    y[i] ~ dnorm(fit[i], sd = 1)
  }
 
}
)


y_m  <- mean(dt$y)
y_sd <- sd(dt$y)

## Prep constants, data and initial values
Bars_Const <- list(N = dim(dt)[1], Kmax = 3)
Bars_Data <- list(y = (dt$y - y_m) / y_sd,
                x = dt$x)
Bars_Inits <- list(betas = rep(0, 5),  KNum = 2 ,
                 Theta = rep(1/3, 3),

                 K = runif(3, 0.1, 9.9))

## Run MCMC posterior samples
Bars_mcmc.out <- nimbleMCMC(nim_BarsCode       ,
                            Bars_Const         ,
                            Bars_Data          ,
                            inits = Bars_Inits ,
                            nchains = 4        ,
                            niter = 1e4/4      ,
                            summary = T        ,
                            monitors = c("KNum", "Theta", "fit"))


## Summary stats of number of knots simulated
Bars_mcmc.out$summary$all.chains[names(Bars_mcmc.out$summary$all.chains[,1]) %like% 'KNum', ]


## Plot attempt at Bars
plot_ly(dt) %>%
  add_markers(x = ~x, y = ~y, name = "Simulated data") %>%
  add_lines(x = ~x, y = ~Bars_mcmc.out$summary$all.chains[,1][names(Bars_mcmc.out$summary$all.chains[,1]) %like% 'fit']*y_sd+y_m, name = "Nimble Bars fit")



## DIRICHLET IMPLEMENTATION OF BARS ---------

## This model has to separately cater for cases where
## there is a single knot provided as a scalar otherwise
## model does not compile
Rnsp_sc <- nimbleRcall(function(x = double(1), K = double(0)) {},

                       Rfun = 'Build_nsp',
                       returnType = double(2))



nim_spCode <- nimbleCode({
 

  ## Prep all possible spline
  ## design matricies
  B1[1:N,1:p1] <- Rnsp_sc(x[1:N], K = K1)
  B2[1:N,1:p2] <- Rnsp(x[1:N], K = K2[1:KNum2])
  B3[1:N,1:p3] <- Rnsp(x[1:N], K = K3[1:KNum3])
 
  ## Produce each fit separately and then
  ## combine as a weighted average of all
  ## using the dirichlet to weight between
  ## respective number of knot selections

  for (i in 1:N) {

    fit1[i] <- inprod(B1[i, 1:p1], betas1[1:p1])
    fit2[i] <- inprod(B2[i, 1:p2], betas2[1:p2])
    fit3[i] <- inprod(B3[i, 1:p3], betas3[1:p3])
    fit[i] <-  p_f[1]*fit1[i] + p_f[2]*fit2[i] + p_f[3]*fit3[i]
  }
 
  ## Sample weightings of each model
  Alpha[1:3] <- nimRep(1, 3)
  p_f[1:3] ~ ddirch(Alpha[1:3])
 
  ## Sample for each respective number of knots
  K1 ~ dunif(0.1, 9.9)
 
  for (i in 1:KNum2){
    K2[i] ~ dunif(0.1, 9.9)
  }
  for (i in 1:KNum3){
    K3[i] ~ dunif(0.1, 9.9)
  }
 
  ## Sample for each possible beta value
  for (i in 1:p1) {
    betas1[i] ~ dnorm(0, sd = 1)
  }
 
  for (i in 1:p2) {
    betas2[i] ~ dnorm(0, sd = 1)
  }
  for (i in 1:p3) {
    betas3[i] ~ dnorm(0, sd = 1)

  }
 
  for (i in 1:N) {
    y[i] ~ dnorm(fit[i], sd = 1)
  }
 
}
)


y_m  <- mean(dt$y)
y_sd <- sd(dt$y)

sp_Const <- list(N     = dim(dt)[1] ,

                 KNum2 = 2          ,
                 KNum3 = 3          ,
                 p1    = 1 + 2      ,
                 p2    = 2 + 2      ,
                 p3    = 3 + 2

)

sp_Data <- list(y = (dt$y - y_m) / y_sd,
                x = dt$x)

sp_Inits <- list(

  betas1 = rep(0, 3)   ,
  betas2 = rep(0, 4)   ,
  betas3 = rep(0, 5)   ,
  K1     = 1           ,
  K2     = c(1, 2)     ,
  K3     = c(1, 2, 3)  ,
  p_f    = rep(1/3, 3)
)


sp_mcmc.out <- nimbleMCMC(nim_spCode     ,
                          sp_Const    ,
                          sp_Data     ,
                          inits = sp_Inits    ,
                          nchains = 4  ,
                          niter = 1e4/4  ,
                          summary = T  ,
                          monitors = c("p_f", "fit", "K1", "K2", "K3"))


## Obtain samples in a data.table
nim_mcmc_samples <- data.table::rbindlist(lapply(1:length(sp_mcmc.out$samples),
                                                 function(i) {as.data.table(sp_mcmc.out$samples[[i]])[, Chain := i] }))


## Observe the density associated with a single knot
plot(density(nim_mcmc_samples$K1))

## Plots the nimble Bars attempt Vs. the dirichlet weightings of each number of possible knots
plot_ly(dt) %>%
  add_markers(x = ~x, y = ~y, name = "Simulated data") %>%
  add_lines(x = ~x, y = ~Bars_mcmc.out$summary$all.chains[,1][names(Bars_mcmc.out$summary$all.chains[,1]) %like% 'fit']*y_sd+y_m, name = "Nimble Bars fit") %>%
  add_lines(x = ~x, y = ~sp_mcmc.out$summary$all.chains[,1][names(sp_mcmc.out$summary$all.chains[,1]) %like% 'fit']*y_sd+y_m, name = "Nimble Dirichlet fit")

[Daniel Turek]

Kavir, that's great to hear the model is working.  Yes, as best I'm following, it makes sense that the minor differences might result from the BARS implementation, versus the Dirichlet-weighted fit of separate models.

I'll let Perry say more, but there might be more value to the RJMCMC approach than mere efficiency (of only sampling the required knot locations and betas).  It's true, that with a proper RJMCMC implementation, then extraneous (currently not in the model) dimensions of betas and the knots will not be sampled.  But beyond that, an RJMCMC sampler will also inform the "best" number of knots for your model, the number of knots most parsimonious with the model fit.  That is, indicating how many (or how few) knots are most appropriate for your modelling problem, or more accurately, the approximate range of (the number of) knots which is reasonable.

Cheers,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/15bcc9d6-1894-40b7-b995-3ea1511ea0b9n%40googlegroups.com.

[Perry de Valpine]

Yes, that all makes sense. Nice job Kavir getting things working. Regarding RJMCMC: often an alternative implementation is to use indicator variables that take values of 1 (if a term is in the model) and 0 (if it is out of the model), which multiply the term to have their intended effect. In that case, one can sometimes see poor mixing because a parameter being multiplied by a 0 will follow its prior, which may take it on long sojourns away from values that could be reasonable for going into the model, making it hard for the indicator to be sampled back to 1. With RJMCMC, there is a proposal distribution for the value of a parameter to be turned back on in the model (jointly with setting the indicator 1), making for potentially better mixing over which parameters are included in the model. I hope that makes some sense.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAKbe0hpF1iNMFyFPStoqKz0Do9Ho7hTn2-09s5tRkDz6nxWK-A%40mail.gmail.com.

[Kavir Narsai]

Hi Daniel and Perry, thank you both for the feedback.

OK that makes a lot of sense. Will definitely try to get an RJMCMC implementation of the model working to ensure it is more robust.

[Kavir Narsai]

Hi Perry and Daniel

Thank you again for your assistance. To continue where we left off, please see below my attempt at implementing a custom RJMCMC algorithm for running Bayesian Adaptive Regression Splines. (I attempted to adjust the RJ MCMC code provided in the following script https://github.com/nimble-dev/nimble/blob/devel/packages/nimble/R/MCMC_RJ.R to work for this dummy example).

Unfortunately, the code seems to get 'stuck' at the initial number of knots simulated and does not explore different knot values? I also note warnings of LogProb of -Inf for certain knot values which have been set to zero and so should not be in the model in the first place so I'm not sure why the MCMC algorithm is determining its log probability? 

I'm assuming somewhere in the algorithm this is leading to the probability of any proposed changes to the number of knots always being zero but I'm not sure what exactly could be causing this? Please can you assist me in determining how the code below needs to be adjusted to work for this example.

Thank you

Best regards

Kavir

## R CODE:

## ----------------------------------------------------------------------------
## ADAPTIVE BAYESIAN REGRESSION SPLINES CUSTOM SAMPLER
## Attempt to implement BARS using nimble
## ----------------------------------------------------------------------------


## LOAD PACKAGES ----------------


suppressPackageStartupMessages({
  library(data.table)
  library(plotly)
  library(splines2)
  library(ModelMatrixModel)
  library(nimble)
})



## SIM DATA ----------------

set.seed(123)
dt <- data.table(x = seq(0.1, 10, 0.1))
betas <- c(20, 30, 50)

dt[, y := nsp(x, knots = 5, intercept = T)%*%betas + rnorm(100, 0, 0.5)
]

plot_ly(dt) %>%
  add_markers(x = ~x, y = ~y)


## P-BARS MODEL ----------------

## Specify nsp code in R
Build_nsp <- function(x, K) {
  nsp(x, knots = K, intercept = T)
}

Rnsp <- nimbleRcall(function(x = double(1), K = double(1)) {},
                    Rfun = 'Build_nsp',
                    returnType = double(2))

## Produce the fit from the selected spline
calc_fit <- nimbleFunction(
  run = function(x = double(1), K = double(1), betas = double(2), KNum = integer(0)) {

   
    p <- KNum + 2
    B <- Rnsp(x, K[1:KNum])
   

    fitArray <- B %*% betas[1:p, KNum]

    returnType(double(1))
    return(fitArray[, 1])
  }
)

## Specify model code

BARS_code <- nimbleCode({

  Alpha[1:Kmax] <- nimRep(0.1, Kmax)

  Theta[1:Kmax] ~ ddirch(Alpha[1:Kmax])
  KNum ~ dcat(Theta[1:Kmax])
 

  K[1] ~ dunif(x_min, x_max)
  for (i in 2:Kmax){
    K[i] ~ dunif(K[i-1], x_max)
  }
 
  for (j in 1:Kmax) {
    betas[1, j] ~ dnorm(0, sd = 1)
    for (i in 2:(Kmax+2)) {
      betas[i, j] ~ dnorm(betas[i-1, j], sd = 1)
    }
  }
 
  fit[1:N] <- calc_fit(x[1:N], K[1:Kmax], betas[1:(Kmax+2), 1:Kmax], KNum)

 
  for (i in 1:N) {
    y[i] ~ dnorm(fit[i], sd = 1)
  }
 
})

## Prepare required params for the model
BARS_const <- list(N     = dim(dt)[1] ,
                   Kmax  = 5          ,
                   x_min = min(dt$x)  ,
                   x_max = max(dt$x)  
)

BARS_data <- list(x = dt$x,
                  y = dt$y)

BARS_inits <- function(BARS_const) {
 
  KNum  = sample(1:BARS_const$Kmax, 1)
  K     = sort(runif(BARS_const$Kmax, BARS_const$x_min, BARS_const$x_max))
  K[(KNum+1):BARS_const$Kmax] <- 0
  betas = matrix(rnorm((BARS_const$Kmax+2)*BARS_const$Kmax), ncol = BARS_const$Kmax)
  z = matrix(0, ncol = BARS_const$Kmax, nrow = BARS_const$Kmax+2)
  z[, KNum] <- 1
  betas <- betas*z
 
  return(list(Theta = rep(1 / BARS_const$Kmax, BARS_const$Kmax) ,
              KNum  = KNum ,
              K     = K ,
              betas = betas
  ) )
 
}



## Specify the Bars Model
BARS_model <- nimbleModel(code       = BARS_code  ,
                          constants  = BARS_const ,
                          data       = BARS_data  ,
                          inits      = BARS_inits(BARS_const))


## Check all variables have been initialised
BARS_model$initializeInfo()

## Compile model in Cpp
CBARS_model <- compileNimble(BARS_model)

## Configure the base MCMC process
BARS_Conf <- configureMCMC(BARS_model,
                           monitors = c("KNum", "K", "betas", "fit"))


BARS_Conf$getSamplers()


## CUSTOM KNUM SAMPLER ------------------
## Specify a custom sampler of the number of knots which will then require an
## associated proposal distribution for knot locations and coefficients

Sampler_BARSRJ_Knots <- nimbleFunction(
  name = 'Sampler_BARSRJ_Knots',
  contains = sampler_BASE,
  setup = function(model, mvSaved, target, control) {
    x_min <- control$x_min
    x_max <- control$x_max
    Kmax  <- control$Kmax
   
    ## Select all possible parameters that can be effected
    ## for specifying the nodes to be saved in the last step
    calcNodes <- model$getDependencies(c("KNum", "K", "betas"))
   
    ## Due to dependency on the previous knot we need to remove the knot
    ## not being used in the model from the dependency structure
    Current_KNum <- model[[target]]
    KnotNodes <- paste0("K[1:", Current_KNum, "]")
    BetaNodes <- paste0("betas[1:", Kmax+2, ", ",  Current_KNum,"]")
   
    CalcCurrentNodes <- model$getDependencies(c("KNum", KnotNodes, BetaNodes))
    CalcCurrentNodes <- CalcCurrentNodes[CalcCurrentNodes != paste0("K[", Current_KNum+1, "]")]
   
    ## Need to determine required knot selection at this step
    ## cannot make use of r functions in run code like paste0
    if (Current_KNum == 1) {
      ## If we at the minimum number of knots we cannot reduce
      CalcReducedNodes <- CalcCurrentNodes
    } else {
      ## Otherwise we can remove a knot
      ReducedKnotNodes <- paste0("K[1:", Current_KNum-1, "]")
      ReducedBetaNodes <- paste0("betas[1:", Kmax+2, ", ",  Current_KNum-1,"]")
     
      CalcReducedNodes <- model$getDependencies(c("KNum", ReducedKnotNodes, ReducedBetaNodes))
      CalcReducedNodes <- CalcReducedNodes[CalcReducedNodes != paste0("K[", Current_KNum, "]")]
     
    }
   
    if (Current_KNum == Kmax) {
      ## If we at the max number of knots we cannot increase
      CalcAddNodes <- CalcCurrentNodes
    } else {
      ## Otherwise we can remove a knot
      AddKnotNodes <- paste0("K[1:", Current_KNum+1, "]")
      AddBetaNodes <- paste0("betas[1:", Kmax+2, ", ",  Current_KNum+1,"]")
     
      CalcAddNodes <- model$getDependencies(c("KNum", AddKnotNodes, AddBetaNodes))
      CalcAddNodes <- CalcAddNodes[CalcAddNodes != paste0("K[", Current_KNum+2, "]")]
     
    }
   
  },
  run = function() {
    # Current_KNum <- model$KNum
   
    ## Determine proposed change to current num knots
    if (Current_KNum == Kmax) {
      ## If the current number of knots is the maximum allowed
      ## number then propose removing a knot
      K_delta <- -1
    } else {
      if (Current_KNum == 1) {
        ## If the current number of knots is 1
        ##  then propose adding a knot
        K_delta <- 1
      } else {
        ## Randomly add or remove a knot
        if (runif(1)>0.5) {
          K_delta <- 1
        } else {
          K_delta <- -1
        }
      }
    }
   
    ## calculated proposed number of knots
    Proposed_KNum <- Current_KNum+K_delta
   
    if (K_delta == 1) {
      ## Logic for when we adding a new knot
      currentLogProb <- model$getLogProb(CalcCurrentNodes)
     
      ## Update number of knots
      model[[target]] <<- Proposed_KNum
     
      ## In all instances this will be the 2nd or
      ## beyond knot as there is a min of 1 knot
      ## as such we can always use the previous knot
      ## value as a lower bound
      Proposed_knot <- runif(1, model$K[Proposed_KNum-1], x_max)
      model$K[Proposed_KNum] <<- Proposed_knot
     
      ## Propose values for the associated set of
      ## coefficients to be activated.
      Current_Coef <- model$betas[, Current_KNum]
      Proposed_Coef <- rnorm(Kmax+2, 0, 1)
      model$betas[, Proposed_KNum] <<- Proposed_Coef
      model$betas[, Current_KNum] <<- 0
     
      ## Determine associated probability of proposed inclusions
      ## noting that the existing parameter set is being removed
      ## Note at each step we are running the MH algorithm for a specific variable
      ## while keeping all others constant!!
     
      logProbForwardProposal <- dunif(Proposed_knot, model$K[Proposed_KNum-1], x_max, log = TRUE) +
        sum(dnorm(Proposed_Coef, 0, 1, log = TRUE)) - sum(dnorm(Current_Coef, 0, 1, log = TRUE))
      proposalLogProb <- model$calculate(CalcAddNodes)
     
      logAcceptanceProb <- proposalLogProb - currentLogProb - logProbForwardProposal
     
    } else {
      ## Logic for when we remove a knot
      currentLogProb <- model$getLogProb(CalcCurrentNodes)
     
      ## Update number of knots
      model[[target]] <<- Proposed_KNum
     
      ## Define logic to remove the last knot
      Current_knot <- model$K[Current_KNum]
      model$K[Current_KNum] <<- 0
     
      ## Amend coefficient values
      Current_Coef <- model$betas[, Current_KNum]
      Proposed_Coef <- rnorm(Kmax+2, 0, 1)
      model$betas[, Proposed_KNum] <<- Proposed_Coef
      model$betas[, Current_KNum] <<- 0
          
      logProbReverseProposal <- dunif(Current_knot, model$K[Current_KNum-1], x_max, log = TRUE) -
        sum(dnorm(Proposed_Coef, 0, 1, log = TRUE)) + sum(dnorm(Current_Coef, 0, 1, log = TRUE))
     
      proposalLogProb <- model$calculate(CalcReducedNodes)
     
      logAcceptanceProb <- proposalLogProb - currentLogProb + logProbReverseProposal
    }
   
   
    accept <- decide(logAcceptanceProb)
    if(accept) { copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
    } else     { copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE) }
  },
  methods = list(
    reset = function() { }
  )
)

## Next we need to define the toggle logic
## Here we can try to leverage the existing rj method of only
## sampling if current value != fixed value (0)
## Note that each target node of the knots and betas will need to be looped through this logic
BARS_RJ_toggled <- nimbleFunction(
  name = 'BARS_RJ_toggled',
  contains = sampler_BASE,
  setup = function(model, mvSaved, target, control) {
    ## control list extraction
    fixedValue  <- 0
    samplerType <- if(!is.null(control$samplerType)) control$samplerType else stop('must provide \'samplerType\' control list element to RJ_toggled sampler')
    ## nested function and function list definitions
    samplerFL <- nimbleFunctionList(sampler_BASE)
    samplerFL[[1]] <- samplerType$buildSampler(model = model, mvSaved = mvSaved)
  },
  run = function() {
    if(model[[target]] != fixedValue)
      samplerFL[[1]]$run()
  },
  methods = list(
    reset = function() {
      samplerFL[[1]]$reset()
    }
  )
)

## Configure number of knot sampler

BARS_Conf$removeSamplers("KNum")

BARS_Conf$addSampler(type = "Sampler_BARSRJ_Knots",
                     target = "KNum",
                     control = list(x_min = BARS_const$x_min,
                                    x_max = BARS_const$x_max,
                                    Kmax  = BARS_const$Kmax
                     ))

## Check samplers
BARS_Conf$getSamplers()

## Configuration function for knots and coefficients
## to be toggled in and out of the model
configureBARSRJ_toggled <- function(conf, targetNodes) {
  model <- conf$model
  nNodes <- length(targetNodes)
 
  for(i in 1:nNodes) {
    nodeAsScalar <- model$expandNodeNames(targetNodes[i], returnScalarComponents = TRUE)
   
    for(j in 1:length(nodeAsScalar)) {
      currentConf <- conf$getSamplers(nodeAsScalar[j])
     
      ## Add sampler for the coefficient variable (when is in the model)
      conf$removeSamplers(nodeAsScalar[j])
      conf$addSampler(type = BARS_RJ_toggled,
                      target = nodeAsScalar[j],
                      control = list(samplerType = currentConf[[1]]))
    }
  }
 
  return(invisible(NULL))
}


## Configure the parameters to be toggled in and out
## of the sampling process
configureBARSRJ_toggled(BARS_Conf, c("K", "betas"))

## Check samplers
BARS_Conf$getSamplers()

## Build the MCMC object
nimBarsRJMCMC <- buildMCMC(BARS_Conf)


## Compile the Sampler
CnimBarsRJMCMC <- compileNimble(nimBarsRJMCMC)

BARS_model$KNum
BARS_model$K

RJBarssamples <- runMCMC(CnimBarsRJMCMC,
                         niter = 1e4,
                         nchains = 4
)


## NB! PROBLEM: WARNINGS PROVIDED OF -INF FOR KNOT LOCATIONS > CURRENT NUMBER OF KNOTS
## THE SAMPLING IS NOT MOVING BETWEEN DIFFERENT NUMBER OF KNOTS??

nim_rjmcmc_samples <- data.table::rbindlist(lapply(1:length(RJBarssamples),
                                                   function(i) {as.data.table(RJBarssamples[[i]])[, Chain := i] }))


summary(nim_rjmcmc_samples$KNum)

Bars_fit <- colMeans(nim_rjmcmc_samples[, .SD, .SDcols = patterns('fit')])

plot_ly(dt) %>%
  add_markers(x = ~x, y = ~y, name = "Simulated data") %>%

  add_lines(x = ~x, y = ~Bars_fit, name = "Nimble Bars fit")

[Perry de Valpine]

Hi Kavir,

Thanks for sending the code. I've tried to follow, and it's not entirely clear to me but here are a few thoughts. The setup code of Sampler_BARSRJ_Knots creates sets of nodes for calculations based on the value of KNum in the model at the time the sampler is built. This is run once and only once in R. The important thing to realize is that sets of nodes like CalcCurrentNodes and others will then never change, so I think the calculations in the run function won't work as you want. I am also not clear why the betas is set up as a matrix, with column indexed by KNum. That will mean there is a separate set (column) of betas for each possible value of KNum. But when you are trying to increase or decrease KNum (if you are able to fix this), I think you could just use the same vector of betas. The "toggled" sampler as adapted here is not making sense to me. The target node the way this was written in nimble should be an indicator variable (taking values of 0 or 1) in the model, but I think you are using it for each knot (K) and each beta, which doesn't make sense to me. I'm not sure you even really need a "toggle" scheme because the role of that in the variable selection RJMCMC system is to toggle indicators that represent taking a coefficient in and/or out of a model, separately from each other. With the knots situation, it looks like you want to sample KNum higher or lower, but you don't have (or need) it set up with a set of independent indicator variables. (A more minor point is that from 1 or Kmax, the proposal for KNum would not be symmetric, so that may need accounting in the acceptance probability.)

It may be that the pieces of RJMCMC and what you are setting up are not entirely clear, making it hard to get things right, so one suggestion is that you can run the MCMC uncompiled. Then you can either put debuggers on functions or put browser() in functions to step through and follow more closely what is happening (e.g. put browser() at the start of the sampler's run method). You can do this for nimble's samplers as well. If you want to watch the entire MCMC stepped through, you can set debug(CnimBarsRJMCMC$run). 

I hope that helps.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/7f92385d-95be-4eef-aa72-4be69a390dc2n%40googlegroups.com.

[Perry de Valpine]

P.S. I meant "debug(nimBarsRJMCMC$run)" (using the uncompiled version, not the compiled one).

[Kavir Narsai]

Hi Perry

Noted thank you so much for the feedback. 

I will attempt to review and debug the code and I'll let you know the outcome. 

With regard to the matrix of coefficients, my understanding is that for a given spline fit, a change in the number of knots considered completely changes the curved representation of the covariate and as such would require a different set of coefficients. 

In terms of the toggled sampler, my intention is to add/remove knot locations and coefficients as required depending on the number of knots(KNum) selected in a specific sampling step. I will also review this step to see how I can include an indicator variable to make this step work as required. 

Thank you Perry.

Best regards

Kavir 

[Kavir Narsai]

Hi Perry

As part of my attempt to dynamically change the nodes referenced when calculating the log probability of the proposals in the Sampler_BARSRJ_Knots sampling process, how do I change the nodes being calculated in each step based on the number of knots (KNum) proposed in each iteration. I have tried the following and these have not worked:

1. It appears that the call model$getDependencies does not work in the run step of the sampling procedure and as such, I kept getting the following error when using it to determine the required nodes: "Cannot handle this expression (looks like a model with an invalid member function call?): model$getDependencies(cnim(..."
2. Attempted to build a nested list in the setup step to store all possible node values that can be indexed by KNum however, this produced an error when compiling the model regarding the type of list being used
3. Attempted to make use of a nimbleList to store all the required node values based on the number of knots but this produced a very unusual error when compiling which is as follows when running printErrors():

"C:\Users\XKV\AppData\Local\Programs\R\R-4.2.0patched\library\nimble\include/nimble/RcppNimbleUtils.h:229:6: note: candidate: 'SEXPREC* NimArr_2_SEXP(const NimArr<ndim, double>&) [with int ndim = 1; SEXP = SEXPREC*]'
  229 | SEXP NimArr_2_SEXP(const NimArr<ndim, double> &val);
      |      ^~~~~~~~~~~~~
C:\Users\XKV\AppData\Local\Programs\R\R-4.2.0patched\library\nimble\include/nimble/RcppNimbleUtils.h:231:6: note: candidate: 'SEXPREC* NimArr_2_SEXP(const NimArr<ndim, int>&) [with int ndim = 1; SEXP = SEXPREC*]'
  231 | SEXP NimArr_2_SEXP(const NimArr<ndim, int> &val);
      |      ^~~~~~~~~~~~~
make: *** [C:/Users/XKV/AppData/Local/Programs/R/R-42~1.0PA/etc/x64/Makeconf:260: P_23_nfRefClass_R_GlobalEnv98.o] Error 1

"

Ideally the nodes included in the calculation of the log probability of the proposal should include/exclude knots and associated coefficients as required, however, I'm not sure how else to try to adjust the nodes to be calculated in the run portion of the sampler?

Thank you Perry.

Best regards

Kavir

[Perry de Valpine]

Hi Kavir,

For this step of what you are doing, I think the nimbleFunctionList will be a good approach. It is correct that sets of nodes, such as determined by model$getDependencies, can only be established in setup code. In run code they cannot be modified. A solution is a list of nimbleFunction objects, one for each set of nodes your sampler might need to handle. (And of course, the nimble compiler system is limited in handling some of the other things you tried because it is just not that general for many reasons.)

Here is a toy example that I hope illustrates the key ideas for you.

library(nimble)

set.seed(1)
data <- list(y = rnorm(10))
constants <- list(p = rep(.1, 10))
inits <- list(x = rnorm(10), k = 4)

# Here is a model where k will be the largest index that we want to include in calculations.

# e.g. is k==4, we will want to use x[1:4] and y[1:4]
m <- nimbleModel(
  nimbleCode({
    for(i in 1:10) x[i] ~ dnorm(0, 1)
    for(i in 1:10) y[i] ~ dnorm(x[i], 1)
    k ~ dcat(p[1:10])
  }),
  data = data, constants = constants, inits = inits
)

# we need a base class
calc_fixed_k_base <- nimbleFunctionVirtual(
  methods = list(
    calc = function() {
      returnType(double())
    }
  )
)

# Here is a nimbleFunction for one set of nodes, based on one value of k.
calc_fixed_k <- nimbleFunction(
  contains = calc_fixed_k_base,
  setup = function(model, nodes) {
    calcNodes <- model$getDependencies(nodes)
  },
  methods = list(
    calc = function() {
      ans <- model$calculate(calcNodes)
      returnType(double())
      return(ans)
    }
  )
)

# Here is a nimbleFunction that makes a nimbleFunctionList of the above, with one for each possible value of k.

# Some things are simplified, like hard-coding the names "x" and "k", and could be made more general.

# You could use this kind of approach in a sampler, but here I am just showing a simpler nimbleFunction.
calc_any_k <- nimbleFunction(
  setup = function(model) {
    mycalc_list <- nimbleFunctionList(calc_fixed_k_base)
    xnodes <- model$expandNodeNames("x")
    for(k in 1:10) {
      mycalc_list[[k]] <- calc_fixed_k(model, xnodes[1:k])
    }
  },
  run = function() {
    k <- model[["k"]]
    ans <- mycalc_list[[k]]$calc()
    return(ans)
    returnType(double())
  }
)

# Make one instance:
my_calc_any_k <- calc_any_k(m)

# Verify that it works when running uncompiled
my_calc_any_k$run()
m$calculate(c("x[1:4]", "y[1:4]")) # This shows the previous answer is correct (with k == 4)

# Compile
cm <- compileNimble(m)
c_my_calc_any_k <- compileNimble(my_calc_any_k, project = m)

# verify that it works when running compiled, with a new value of k for good measure.
cm$k <- 6
c_my_calc_any_k$run()
cm$calculate(c("x[1:6]", "y[1:6]"))# This shows the previous answer is correct (with k == 6)

Finally, the issues of whether you want to re-use a vector with different extents for different numbers of knots or to have a matrix with columns corresponding to different numbers of knots is something you will have to work out based on the statistical ideas you want to implement. I was just trying to point out that either one is an option in terms of software mechanics.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/03dd15ab-13f0-4e42-a93e-7c569b151e24n%40googlegroups.com.

[Kavir Narsai]

Hi Perry

Thank you so much for your feedback. The nimbleFunctionList is working as required and is dynamically adjusting the nodes as it samples.

I really appreciate your assistance with this.

Best regards

Kavir

[Colton Padilla]

Hey Perry,

I wanted to bring this back to the top as I am currently working on implementing a similar sampler within nimble but I am having trouble getting the model to mix correctly. This is my first run at attempting to write my own sampler in nimble, so I'm certain there are things I don't quite understand about this. I have added in a reproducible example here that unfortunately does not have data along with it as the data is sensitive, so if I would need to find a way to simulate this data, let me know. The trouble I am having is the number of knots is not being sampled very well and I am assuming it has something to do with how I am sampling the locations and the coefficients in the model. My direct question as to not go too far in multiple directions at this point is: how should I be sampling the locs and b_lag variables in this context given that sometimes coefficients are being added, removed, or minorly perturbed within the sampler? Additionally, I wrote this sampler following two papers (Dennison et al. 1998 and Dimatteo et al. 2001) that are quite frankly a small amount over my head when it comes to mathematics so it is possible that I have just written the logic incorrectly. Thank you for any help you can provide here.

[Perry de Valpine]

Hi Colton,

This looks interesting. I wasn't able to look up the papers -- not sure if you meant to provide links for them, but in any case I was primarily looking at the code and use of nimble and not at the correctness of the math.

One minor detail that is not related to your mixing outcome is that you should be able to do "currentLogProb <- model$getLogProb(calcNodes)". That will use the cached log probability values rather than re-calculating them. It's part of the MCMC protocol that upon entry you can assume the model is in a fully calculated state (and upon exit you must ensure that is the case for the next sampler, which you do with the nimCopy calls). 

As far as I can see the code looks mechanically reasonable. I can't check the math or say anything about what to expect for mixing. I wondered if it would help to draw the new beta value (when adding a knot) from a distribution centered around the neighboring beta values to make it more likely to be accepted. But I don't know whether or how that would affect the acceptance probability for the reversible jump situation.

I could see that `lag_effect` might potentially be made more efficient if the new basis function calculations could be done only for a new knot location, but that is an efficiency consideration and again wouldn't be related to mixing.

You probably know that many people now prefer to establish a fixed number of knots, with plenty of them, and then manage smoothness vs wiggliness with a variance parameter. But I'm taking at face value that you want to have the knots and their placement be estimated.

Sorry I can't be more helpful. If you want to see what is going on internally, since you already have a nimbleRcall, you could insert a browser() into `lag_effect` and then from that browser you could inspect the state of the compiled model to see what proposal is being considered. Or you could make a new nimbleRcall solely for debugging purposes and call it from the run code of the new sampler (passing in any values you want to inspect).

HTH

Perry

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/f1d7d4ae-8324-4224-a893-88baba8824f4n%40googlegroups.com.

[Colton Padilla]

Hey Perry,

Thanks for the insights. I would under normal circumstances use a fixed basis matrix for these splines, but right now we are trying to establish some type of statistical metric to indicate when the trend in the data is stabilizing/changing. We can get the model to fit the data pretty well with the fixed basis matrix, but it lacks some type of statistical metric to indicate which knots are showing where the distribution is stabilizing. Part of this difficulty is due to autocorrelation between individuals as these animals are highly social and it results in dense knot placement for the true trend. Basically, it seems as though the animals are following the exact same pattern where the mean value of the metric is declining and then stabilizing overtime, but the autocorrelation between individuals ends up causing the spline to wiggle way more than we would want it to. I have attached an example where I fit a spline using a dense knot placement and ddexp shrinkage priors where it still allows a lot of these wiggles to occur purely because the individuals are extremely autocorrelated. This may be more of a data processing thing than it is a spline knot selection thing, but we are still trying to find some type of way to show statistical significance or something relevant on the knot time periods themselves.

The two papers I am referencing are here (https://academic.oup.com/biomet/article/88/4/1055/225982?guestAccessKey=) and here (https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00128). Now that I look at them more, I'm not sure that I do have this parameterized correctly. My first issue that I am dealing with is that when there are only a certain amount of knots proposed, the way I have this coded, it truncates where the knots can be placed based on the where the nth knot is placed. So if I place 67 locations along the line but the initial value is 32, only about 1/2 of the entire distribution of TST values is actually sampled. I am going to rework my sampler and see if I can get around that behavior here in a minute. If I am able to walk my way through it and find a way to implement this sampler, I will share it with the group.

[Perry de Valpine]

Thanks for the explanations and good luck.

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/212096db-896f-4e37-8772-5990751cf310n%40googlegroups.com.