How to achieve fractional polynomials in stan

Modeling
1

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

\phi_2(x;\xi,p)=\xi_0+\xi_1 x^{(p_1)}+\xi_2 x^{(p_2)},

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.

2

Sorry for not answering earlier.

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.

Best of luck with your model!

3

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.

1 Like