# Nonlinear Gamma regression Stan error message

I’m trying to fit a simple Gamma regression in a nonlinear framework. The data I am modeling are generated from a function that looks like this:

`Y[i] = X1[i]^a * X2[i]^b * X3[i]^c`

which I want to embed in a Gamma likelihood. I could take logarithms of both sides and fit a normal model instead, but I would like to fit the nonlinear version because it is easier to interpret parameters on the original scale. Therefore, I have started by coding the simple model below using just one predictor in rstan, which produces an error that I do not understand how to resolve. Any input is much appreciated.

Jim

The error concerns the parameter that is an exponent in the model statement and says:

``````SYNTAX ERROR, MESSAGE(S) FROM PARSER:

arguments to ^ must be primitive (real or int); cannot exponentiate vector by real in block=model name
error in 'model240c7f1c5fe_6c57a88a869c8bf7ba65d94d3245c8be' at line 22, column 35
-------------------------------------------------
20:    b ~ normal(0,10);
21:    for(i in 1:N)
22:    y[i] ~ gamma(shape, (shape / x^b));
^
23:
-------------------------------------------------

PARSER EXPECTED: <expression>
Error in stanc(file = file, model_code = model_code, model_name = model_name,  :
failed to parse Stan model '6c57a88a869c8bf7ba65d94d3245c8be' due to the above error.

# R and stan code below here
# First simulate some data assuming a Gamma error

set.seed(999)
N <- 100
x <- runif(N, 0, 1)
b <- 1.7
y_true <- x^b
shape <- 10
y <- rgamma(N, rate = shape / y_true, shape = shape)`

# plot the simulated data

plot(x, y)
lines(sort(x), y_true[order(x)])`

# Now fit the model using Stan

library(rstan)
stan_code <- "

data{

int<lower=0> N;
vector[N] x;
vector[N] y;

}

parameters {

real b;
real<lower=0.001,upper=100> shape;

}

model {

b ~ normal(0,10);
for(i in 1:N)
y[i] ~ gamma(shape, (shape / x^b));

}

generated quantities {

vector[N] y_rep;
for (i in 1:N)

y_rep[i] = gamma_rng(shape, (shape / x^b));

}
"
m_stan <- stan(model_code = stan_code,
data = list(x = x, y = y, N = N))``````

I think you meant `x[i]^b`.

I knew it was something simple. Thank you.