Pooled Nonlinear Regression

Modeling
1

I’m trying to fit a simple three parameter nonlinear logistic function to a pooled dataset consisting of dose/response measurements. I’ve specified the logistic function as follows:
y=β1/((1+x/β2 )^(β3 ))

The data to be modeled is as follows:
Analysis Data.csv (426 Bytes)

The Stan code is as follows:

data {
  int<lower=1> nY;  // Number samples
  vector[nY] x;        // Dose
  vector[nY] y;       // Response
} 

parameters {
  real b1; 
  real b2;
  real b3;
  real<lower=0> tau;
} 

transformed parameters { 
  real<lower=0> sigma; 
  sigma = 1 / sqrt(tau); 

} 
model {

  // priors
  b1 ~ uniform(1,2000000);
  b2 ~ uniform(0.1,100);
  b3 ~ normal(0, 1);
  tau ~ gamma(0.1, 0.1);

  // likelihood
  y ~ normal(b1/(1+(x/b2)**b3), sigma);
}

The stan code generated the following error with b2 underlined in red:
PARSER EXPECTED: expression

As I’m new to stan, I’d greatly appreciate insight regarding my coding mistake. Thanks in advance for your help.

2

The error is because you’re using ** for exponents, the Stan syntax is ^, so you need to change your likelihood to:

y ~ normal(b1/(1+(x/b2)^b3), sigma);
1 Like
3

Modified the code which generated the following error message:

Arguments to ^ must be primitive (real or int); cannot exponentiate vector by real in block=model name.

Thank you very much for your quick response to my question.

4

Ah sorry, missed that it was a vector. There is an element-wise exponentiation operator: .^, but it’s not available in the Stan version used by RStan.

There are two options here. You can use the cmdstanR interface which uses the latest Stan version: https://mc-stan.org/cmdstanr/, and then write your likelihood as:

y ~ normal(b1/(1+(x/b2) .^ b3), sigma);

Alternatively, you can code around this in the current RStan by using:

y ~ normal(b1/(1+exp(b3*log(x/b2))), sigma);
5

I tried option 1 by installing the cmdstanR interface and updating the likelihood. Unfortunately, still throwing an error:

No matches for:
Available argument signatures for operator /:
Expression is ill-formed.

Thank you very much for your time and help. Anxious to get this right.

6

Right, so when you mix arithmetic operations of real and vector/matrix types, you need to specify that they’ll be applied elementwise. So because b1 is a real and 1+(x/b2) .^ b3 is a vector, you need to request elementwise divison (./):

y ~ normal(b1 ./ (1+(x/b2) .^ b3), sigma);
7

Apologies, for being such a pain as I’m still a relatively new R user. Now throwing the following error:

PARSER EXPECTED: “)”

And should I be developing this model using cmdstanR? Do I need to set the path for the cmdstanR package?

Thank you again for your time and expertise.

8

PARSER EXPECTED: “)”

For this I’ll need to see your full model to see where that’s coming from

And should I be developing this model using cmdstanR?

If you use the .^ operator, then you’ll need cmdstanR, otherwise RStan is fine

Do I need to set the path for the cmdstanR package?

The setup is all pretty automated for cmdstanR these days, as long as you follow the install instructions it generally all ‘just works’

9

The full model is a simple 3 parameter logistic function given in the initial post. The analysis is a nonlinear regression of response against dose ignoring day. That is, each day is a replicate for each dose. Please let me know what additional information you need.

10

I’m able to compile your model in cmdstanR using:

stan_code = "data {
  int<lower=1> nY;  // Number samples
  vector[nY] x;        // Dose
  vector[nY] y;       // Response
} 

parameters {
  real b1; 
  real b2;
  real b3;
  real<lower=0> tau;
} 

transformed parameters { 
  real<lower=0> sigma; 
  sigma = 1 / sqrt(tau); 

} 
model {

  // priors
  b1 ~ uniform(1,2000000);
  b2 ~ uniform(0.1,100);
  b3 ~ normal(0, 1);
  tau ~ gamma(0.1, 0.1);

  // likelihood
  y ~ normal(b1 ./ (1 + (x/b2) .^ b3), sigma);
}"

library(cmdstanr)
mod = cmdstan_model(write_stan_file(stan_code))

Note that while I’m just using a string in R for the Stan code, we recommend saving the code to an external file. This lets you use the line numbers from error messages to track down where in the code an issue is.

Let me know if you have any trouble compiling that model

11

The model compiles and executes without throwing errors. The compiler seems to ignore the "PARSER EXPECTED: “)” " error appearing on the stan code file. Is there documentation/tutorial explaining the commands needed to generate typical diagnostic summaries and plots derived from the fitted cmdstan model?

Thank you vey much for your help.

12

Great!

cmdstanr is generally used in combination with the posterior and bayesplot packages for statistical and graphical summaries, there’s a basic overview in this section of the getting started guide: https://mc-stan.org/cmdstanr/articles/cmdstanr.html#posterior-summary-statistics

Alternatively, I’m a *huge* fan of the shinystan package for quickly summarising and exploring results. For that, you’ll need to read the cmdstan results into an rstan object that the package can interface with:

library(cmdstanr)
library(shinystan)

# Sample model in cmdstanR
mod = cmdstan_model(write_stan_file(stan_code))
samp = mod$sample(
  data = your_data,
  parallel_chains = 4,
  # Don't read the results into R yet
  validate_csv = FALSE
)

# Read output into rstan object
rstan_samp = rstan::read_stan_csv(samp$output_files())

# Open ShinyStan interface
launch_shinystan(rstan_samp)
13

This is awesome. Thank you.

1 Like