Hi, I’m wondering if fractional polynomials could be achieved in stan.
Following the definition of Patrick Royston et al. (1994), for example, a fractional polynomial of degree 2 can be illustrated as

where x>0, \xi=(\xi_0,\xi_1,\xi_2) is coefficient, p=(p_1,p_2) is the vector of power. These p_j come from a set S={-2,-1,-0.5,0,0.5,1,2,3}. If p_j=0, then x^{(p_j)}=lnx; If p_1=p_2, then x^{p_2}=x^{p_1}lnx.
Maybe doing variable selection by backward elimination is one option, but I don’t know how to achieve it. Any advice would be appreciated.

One thing I am not clear about is what quantities are known and what are parameters to be estimated. If everything is a parameter, then it still should be theoretically possible by marginalizing out the discrete parameters in p with something like (just a sketch, code not tested) - also I am ignoring the special cases and just taking x^{(p_k)} = x^{p_k} (but the special cases should be straightforward to add

EDIT: The following is in fact likely not correct code, sorry for any confusion.

data {
int<lower=1> S_size;
vector[S_size] S;
int N;
vector[N] x;
vector[N] y;
}
parameters {
simplex[S_size * S_size] p_probs; //The probability that p = (p_1,p_2) for all pairs of p_1, p_2
vector[3] ksi[S_size * S_size]; //Separate parameters need to be maintained for each possible polynomial
real<lower=0> sigma;
}
model {
vector[S_size * S_size] partial_lp;
for(i1 in 1:S_size) {
for(i2 in 1:S_size) {
int index = (i1 - 1) * S_size + i2;
real p1 = S[i1];
real p2 = S[i2];
vector[N] mu;
//Compute the likelihood assuming specific values for p1, p2
for(n in 1:N) {
mu[n] = ksi[index, 1] + ksi[index, 2] * x[n] ^ p1 + ksi[index, 3] * x[n] ^ p2;
}
partial_lp[index] = normal_lpdf(y | mu, sigma);
}
}
//Weigh all the partial probabilities by their estimated probability
target += log_sum_exp(log(p_probs) + partial_lp);
//TODO add priors
}

I however fear this will be quite costly and that the model might not be well identified without a lot of data.

Hi, martin, thank you so much. (To be honest I did not expect any reply)
Your code provided a nice example for constructing fractional polynomials.
And as you concerned, I did find that the fracpoly is not a good or efficient choice. With my dataset, simple quadratic or cubic formations could do the same job, with less computation time. I think I will go with simple model. But again, thanks for the example code.