Using splines in RStan

[helen....@gmail.com]

Hi there,

I am trying to create STAN code for a model using spline terms which I have created outside of the stan code:

week.index <- dat$week.index

smooth.week.index <- s(week.index, bs="cr", k=30)

dat$spline_WeekIndex <- smooth.construct(smooth.week.index, dat, knots=NULL)$X

I then create my matrix of predictors:

X <- cbind(Intercept=1, dat$spline_WeekIndex, other predictors..)

and this is my model:

N <- dim(X)[1]

k <- dim(X)[2]

Outcome <- dat$Outcome

data <- list(N=N, X=X, k=k, Outcome=Outcome)

stan_code<-"

data {

    

    int<lower=0> N; 

    int<lower=0> k; //no. of predictors

    real Outcome[N]; //outcome

    matrix[N, k] X; //predictor matrix

     }

parameters {

vector[k] alpha;

real <lower=0> taud;

}

model {

    vector[N] mu;

    

    mu <- X*alpha;

    alpha ~ normal(0,1000);

    taud ~ uniform(0,1000);

DeathRate ~ normal(mu, tau);

}"

I am having real trouble getting this to converge and when it finally does (after many iterations and a long wait to get the Rhat values close to 1) the neff values for the spline terms are all the same and very small which I think means it hasn’t mixed correctly?

Any suggestions how I could make this mix better and run faster?

Thanks in advance

[Andrew Gelman]

Hi, I don’t understand your model as it includes a variable, DeathRate, that doesn’t seem to be defined in your Stan code.  But, in any case, I don’t think normal(0,1000) priors are a good idea.  Better to use something informative.  And uniform(0,1000) is an even worse idea.  Better to use soft constraints than hard constraints.

A

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Helen Bond]

Hi,

Sorry the “DeathRate” should actually just be “Outcome”.

Is there any problem using such highly correlated predictors (the spline terms) in Stan?

Thanks

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/aigQtH8twGk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Andrew Gelman]

I think linear regression should be pretty fast as long as the coefs are of order of magnitude 1, i.e. not around 1000 or .0001.

A

[Helen Bond]

Thanks for the reply.

I have tried rescaling some of my data so that the coefficients are all of magnitude 1 and have tried some more informative priors for alpha and tau but I am still finding it struggles to converge and the neff values are very small (<10) even when using 100000 iterations.

Any other clues as to what could be the problem? This is my code now (I have also rewritten the spline terms inside the stan code):

X <- cbind(Intercept=1, dat$h1n1_rate_ILI, dat$h3n2_rate_ILI, dat$ph1n1_rate_ILI, dat$ph1n1_rate_ILI_Pandemic, dat$B_rate_ILI, dat$RSV_rate_ILI, dat$TEMP, dat$TEMP2, dat$TEMP3, dat$AH, dat$AH2, dat$AH3, dat$NameChange)

N <- dim(X)[1]

k <- dim(X)[2]

DeathRate <- dat$DeathRate

WeekIndex <- dat$week.index

H <- 30 #no. of knots

xi <- seq(1,761, 761/30)

data <- list(N=N, X=X, k=k, DeathRate=DeathRate, WeekIndex=WeekIndex, H=H, xi=xi)

stan_code<-"

data {

    

    int<lower=0> N; //no. of weeks

    int<lower=0> k; //no. of predictors

    real DeathRate[N]; //outcome

    matrix[N, k] X; //predictor matrix

    vector[N] WeekIndex; // vector for spline function

    

    int<lower=0> H;

    vector[H] xi;

    }

transformed data {

    matrix[N, H+4] week_s;

    matrix[N, k+H+4] Z;

    for(n in 1:N) {

        week_s[n,1] <- WeekIndex[n]^0;

        week_s[n,2] <- WeekIndex[n]^1;

        week_s[n,3] <- WeekIndex[n]^2;

        week_s[n,4] <- WeekIndex[n]^3;

        for(h in 1:H)

        week_s[n,h+4] <- fmax(0,(WeekIndex[n] - xi[h])^3);

    }

    Z <- append_col(X, week_s);

  }

parameters {

vector[k+H+4] alpha;

real <lower=0> taud;

}

model {

    vector[N] mud;

    

    mud <- Z*alpha;

    alpha ~ normal(0,10);

  taud ~ cauchy(0,5);

DeathRate ~ normal(mud, taud);

}"

Thanks in advance!

[Jan]

Your original model was simpler and thus easier to debug. Is it possible to post some data,  or simulated data to run it on?

J

[Andrew Gelman]

Hi, I’d suggest centering your predictors before taking them to the first, second, and third power.  Maybe this will make things more stable?  There’s clearly something wrong if you’re running for 100,000 iterations.  100 or 100 should be enough.

A