If-else statement with estimated parameter

[ugolin...@gmail.com]

Dear all,

I am trying to code a model in which there is a if-else statement

However, it implies parameters that are not in the "constant" section.

It is written as follow : 

if (X[obs] < bp) { ...

With X an explanatory variable contained in the ModelConst and bp a parameter to estimate.

(see this page to see the original code with stan: https://janhove.github.io/analysis/2018/07/04/bayesian-breakpoint-model )

I am getting this error : "Condition must be able to be evaluated based on values in 'constants'."

As I understand it, it is because I am not allowed to use bp in the statement because it is not a constant. Why is that ? And is there  is something planned or known to get around this issue ?

Thank you in advance.

Ugoline

[ugolin...@gmail.com]

For information I also tried with this alternative : 

        conditional_mean[obs][which(X[obs] < bp)]<- intercept + slope_before * (X[obs] - bp)

        conditional_mean[obs][which(X[obs] >= bp)]<- intercept + slope_after * (X[obs] - bp)

But apparently "which" function is not available in Nimble either.

Error: 

defining model...

Error in genReplacementsAndCodeRecurse(x, constAndIndexNames, nimbleFunctionNames,  : 

  Error, R function "which" has non-replaceable node values as arguments.  Must be a nimble function.

Thank you

[Perry de Valpine]

Dear Ugoline,

Thanks for the message.

The if-then-else provided for models is evaluated at the time of model definition.  It provides a way to create multiple model variants starting from the same code.  It does not provide if-then-else logic during model calculations.  Once it is built, the model is a graphical model with fixed relationships of what depends on what.  For this reason, if-then-else logic that would change what depends on what during model calculations is not part of nimble.  That is similar to the reason that the second idea you tried would not make sense for nimble.

I think either of the following two ideas can work:

conditional_mean[obs] <- intercept + slope_before_or_after[ (X[obs] < bp) + 1] * (X[obs] - bp)

This would rely on a logical result of (X[obs] < bp)  being 0 for FALSE and 1 for TRUE, so slope_before_or_after would have length of 2.  This kind of model flexibility is supported because nimble will see that the index for slope_before_or_after depends on a model variable, and internally it will handle conditional_mean[obs]  as depending on both elements of slope_before_or_after, even though it only depends on one element at a time.

Another option, which would be more general, such as for more slope categories, would be:

conditional_mean[obs] <- intercept +  choose_slope(slopes[1:max_slopes], X[obs], bp) * (X[obs] - bp)

choose_slope <- nimbleFunction(

 run = function(slopes = double(1), xobs = double(), bp = double()) {

   if(xobs < bp) return(slopes[1])

   return(slopes[2])

   returnType(double())

})

Again nimble will in this case handle conditional_mean[obs] as depending on all of slopes[1:max_slopes].  (The if-then inside of choose_slope is not in the model code itself and so is not used for deciding what depends on what in the model graph.)

I hope it helps.

-Perry

--
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/da61cf01-e0d7-4c63-9af6-0f7bd2d462c9n%40googlegroups.com.

[ugolin...@gmail.com]

Dear Perry,

Indeed the first one would not be good because later I will use more than two slopes.

However I don't know why the second solution is not working (I think I do something wrong but don't know what).

Thank you for your help.

He is my script as it is now (I tested the choose_slope nimbleFunction outside the model and it works)  : 

choose_slope <- nimbleFunction(

  run = function(slopes = double(1), xobs = double(), bp = double()) {

    if(xobs < bp) return(slopes[1])

    return(slopes[2])

    returnType(double())

  })

modelCodeSLR <- nimbleCode(

  {

    intercept ~ dnorm(150, 25)  # Average Y at breakpoint

    slope_before ~ dnorm(0, 5)  # Slope before breakpoint

    slope_after ~ dnorm(0, 5)   # Slope after breakpoint

    bp ~ dnorm(12, 6)           # Breakpoint X

    error ~ dnorm(0, 20)        # Residual error

    both_slopes <-  c(slope_before,slope_after) 

      for (obs in 1:Nobs) {

        conditional_mean[obs] <- intercept +  choose_slope(both_slopes, X[obs], bp) * (X[obs] - bp)

      }

    for (obs in 1:Nobs) {

      Y[obs] ~ dnorm(conditional_mean[obs], error)

    }

  })

Here is the error: 

> Model <- nimbleModel(code = modelCodeSLR, name = 'Nimble', constants = ModelConsts, data = ModelData)

defining model...

building model...

setting data and initial values...

running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... Error in if (xobs < bp) return(slopes[1]) : 

  valeur manquante là où TRUE / FALSE est requis

checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization).

model building finished.

Warning message:

In model$both_slopes[1] <<- nimC(model$slope_before[1], model$slope_after[1]) :

  le nombre d'objets à remplacer n'est pas multiple de la taille du remplacement

> CModel <- compileNimble(Model)

compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.

Warning, in eigenizing model_both_slopes[1] the [ is still there but nDim is not 0 (not a scalar).

Erreur : Failed to create the shared library. Run 'printErrors()' to see the compilation errors.

[ugolin...@gmail.com]

It seems to work just like that :  

        conditional_mean[obs] <- intercept +  choose_slope(slopes=c(slope_before,slope_after), xobs=X[obs], bp=bp) * (X[obs] - bp)

Not using "both_slopes"

Thank you a lot for your help.

[Perry de Valpine]

In the case that didn't work, I think you would need explicit indexing in this declaration:

both_slopes[1:2] <-  c(slope_before,slope_after) 

The version that did work for you should be essentially doing the same thing internally.

Glad you've got something that works.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/19c458b9-9c34-4952-b433-79b21137f5a6n%40googlegroups.com.

[ugolin...@gmail.com]

Indeed, it's been a while since I used Nimble. I forgot that each position needed to be specified. Thank you for the reminder.

Thank you again for your precious help.

Ugoline

[ugolin...@gmail.com]

Dear Perry,

The model is working but I have an error message (that does not prevent to compile, run or even estimate properly the parameters and that seems to be linked to the nimbleFunction), but I wanted to know if this error should be adressed or safely ignored : 

> Model <- nimbleModel(code = modelCodeSLR, name = 'Nimble', constants = ModelConsts, data = ModelData)

defining model...

building model...

setting data and initial values...

running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... Error in if (xobs < bp) { : 

  valeur manquante là où TRUE / FALSE est requis

checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization).

model building finished.

Is it also possible to use the "else if" statement in nimbleFunction ? I could not find an anwer for that. 

Like this one : 

choose_intercept <- nimbleFunction(

  run = function(intercepts = double(1), xobs = double(),Xbp=double(1)) {

    if(xobs < Xbp[1]) {return(intercepts[1])}

    else if (xobs >= Xbp[1] & xobs < Xbp[2]) {return(intercepts[2])}

    else if (xobs >= Xbp[2] & xobs < Xbp[3]) {return(intercepts[3])}

    else {return(intercepts[4])}

    returnType(double())

  })

Thank you for all your help.

[Perry de Valpine]

Hi Ugoline,

If the MCMC seems to be working, it is probably ok.  The message is telling you something is not fully initialized when the model is created, but it may be initialized later (by the MCMC, for example) and that is why the MCMC might run fine (or might not, if its attempt at initialization doesn't go well).  My guess is that the input to xobs or Xbp is not initialized in some call from the model to choose_intercept, and that's why you're getting an error from choose_intercept.  It's like a lot of error messages in R: there might be something hard to make sense of from an inner level of the call stack and some higher-level information from an outer level.  The outer call layer doesn't necessarily stop the inner call layer from emitting its confusing message, in case it is helpful.

Yes, the nested ifs should work. If you test a compiled version outside of choose_intercept (which is a good idea, simply from Cchoose_intercept <- compileNimble(choose_intercept) ), and it doesn't work, please let me know!

Glad this is making progress for you.

-Perry

 

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/1e85983e-70a6-46f3-9d88-9c5ab4c10302n%40googlegroups.com.