Setting one element of parameter vector to zero

[Ted]

Hi Stan users,

I'm attempting to convert a dynamic time series JAGS program to Stan, and I have a question.  We are predicting mine safety (measured by injury rate) from lagged mine-quarter level data.  The unit of observation is a mine-quarter.  To allow seasonal and regional effects we include an indicator variable for each district-season (eg NewEnglandQ1, NewEnglandQ2, ... NorthwestQ1, NorthwestQ2, etc.).  I have coded this up using the parameter districtseason_eff, a vector with D elements, where D is the number of district-seasons.  I want to set districseason_eff[1] <- 0 as the base case and and estimate the other elements of the vector, but I'm not exactly sure in which section to make the assignment of the first element.  Currently I have it in transformed parameters, as below:

data {
    int<lower=1> N; // number of observations
    int<lower=1> Q; // number of quarters
    int<lower=1> J; // number of variables
    int<lower=1> D; // number of district-seasons
    ...
}


parameters {
    real beta[J, Q + 1];
    real intercept;
    real districtseason_eff[D];
    real<lower=.0001, upper=10000> sig[J];
}


transformed parameters {
    districtseason_eff[1] <- 0;
}


model {
    ...
    for(k in 2:D)
        districtseason_eff[k] ~ normal(0,10000)T(-10,10);
    ...
}

Is this correct?  Is there any way to vectorize setting the prior of the remaining elements of the vector? Such as:

districtseason_eff[2:D] ~ normal(0,10000)T(-10,10);

The full model file is attached.

Thanks,

Ted

[Marcus Brubaker]

What you want is something like the following:

parameters {
...
real<lower=-10,upper=10> districtseason_eff_params[D-1];
...
}

transformed parameters {
real districtseason_eff[D];
districtseason_eff[0] <- 0;
for (i in 1:(D-1)) districtseason_eff[i+1] <- districtseason_eff_params[i];
}

model {
...
districtseason_eff_params ~ normal(0,10000) T(-10,10);
...
}

Note that you don't set the parameters, you create a new vector with
one extra (non-parameter) field and copy the real parameters in.
Also, note the upper and lower bounds on the parameters. They are
needed because you use a truncated normal distribution.

Cheers,
Marcus

> --
>
>

[Marcus Brubaker]

Whoops,

districtseason_eff[0] <- 0;

should be

districtseason_eff[1] <- 0;

Cheers,
Marcus

[Ben Goodrich]

On Tuesday, October 30, 2012 4:34:23 PM UTC-4, Ted wrote: 

I want to set districseason_eff[1] <- 0 as the base case and and estimate the other elements of the vector, but I'm not exactly sure in which section to make the assignment of the first element.  Currently I have it in transformed parameters

I think the best thing to do in situations like these is to copy free parameters into a larger container in the model{} block and set the additional elements to fixed values. That would then solve your problem with the prior because you can just specify the prior on the free parameters. In other words, in your case,

parameters {

    real beta[J, Q + 1];
    real intercept;


    real<lower=-10,upper=10> districtseason_eff[D-1]; // note added constraints
    real<lower=.0001, upper=10000> sig[J];
}
model {
    real Districtseason_eff[D];
    Districtseason_eff[1] <- 0.0;
    for (d in 2:D)
        Districtseason_eff[d] <- districtseason_eff[d-1];

    for(d in 1:(D-1))
        districtseason_eff[d] ~ normal(0,10000)T(-10,10);
    ...
}

That is legal, but it doesn't make a lot of sense to have a normal with a standard deviation of 10000 truncated to the (-10,10) interval. Why not just do an implicit uniform prior over the (-10,10) interval by omitting the explicit prior on districtseason_eff?

Ben

[Ted]

Thanks Marcus!  That helped.

[Ted]

Thanks Ben.  Regarding the prior, another guy put in those priors before I started on the project.  I first want to check that the results from Stan are roughly the same as the results from JAGS before I start tweaking things.  But that's a good point about the truncated normal prior.

[Ted]

Marcus, with your method, can I use districtseason_eff later in the model?  I also have

    traum_injuries ~ poisson(exp(log_total_hours_worked 

                            + intercept

                            + beta[1,qtr] * any_union

                            +  ... // other covariates

                            + beta[31,qtr] * mean_bed_thickness_A4_dum

                            + districtseason_eff[districtseason]);

districtseason is the data vector containing the district-season for each observation as an integer 1,...,D.  Will this still work with your method?  I'm not familiar with the order that Stan executes assignments and parameter transformations.

On Tuesday, October 30, 2012 1:51:49 PM UTC-7, Marcus wrote:

[Ben Goodrich]

On Tuesday, October 30, 2012 5:50:16 PM UTC-4, Ted wrote:

Marcus, with your method, can I use districtseason_eff later in the model?  I also have

    traum_injuries ~ poisson(exp(log_total_hours_worked 

                            + intercept

                            + beta[1,qtr] * any_union

                            +  ... // other covariates

                            + beta[31,qtr] * mean_bed_thickness_A4_dum

                            + districtseason_eff[districtseason]);

districtseason is the data vector containing the district-season for each observation as an integer 1,...,D.  Will this still work with your method?  I'm not familiar with the order that Stan executes assignments and parameter transformations.

See Appendix C of the Stan manual. This is one of the major differences between Stan and BUGS. In Stan, things happen in the order that you write them. So, yes you can do something like the above.

Ben