Adding a new mean for each quarter

Let’s review the Stan code from before:

data {
  int ndates;
  vector[ndates] flights;
}

parameters {
  real mu;
  real<lower=0> sigma;
}

model {
  flights ~ normal(mu,sigma);
  mu ~ normal(32,10);
  sigma ~ uniform(0,20);
}

The model statement, which connects our parameters to our number of flights, will clearly need to be extended. We start with the flights ~ normal(mu,sigma) statement. We are now saying that the number of flights is affected by more things, like day of week and quarter. So let’s keep it simple for now and only add quarter; say that flights ~ normal(qtr_effect,sigma). Note we eliminate mu altogether. I had actually kept it in a previous iteration, but got something that didn’t make a lot of sense. Why do you think that might be?

Now there are 4 quarters, so qtr_effect really needs to be a vector with 4 entries. So let’s try again, on the whole Stan code:

data {
  int ndates;
  vector[ndates] flights;
}

parameters {
  real<lower=0> sigma;
  vector[4] qtr_effect;
}

model {
  flights ~ normal(qtr_effect[qtr],sigma);
  sigma ~ uniform(0,20);
}

We’re almost there. We defined a vector of length 4 with our quarter effect in our parameters block, and we implemented an index for qtr_effect in the model block. We need to put a prior on our quarter effects, which we do with a for statement. For ease, we will assume the same prior on each, a normal with mean 0 and sigma of 10. The last thing we need to add is this qtr bit we added, which now has to map the day corresponding to the number of flights to an integer between 1 and 4 (so it can index qtr_effect). As we pass in number of flights as data, so must we pass in this mapping as data. So our final Stan code for number of flights accounting for quarter looks as follows:

data {
  int ndates;
  vector[ndates] flights;
  int qtr[ndates];
}

parameters {
  real<lower=0> sigma;
  vector[4] qtr_effect;
}

model {
  flights ~ normal(qtr_effect[qtr],sigma);
  sigma ~ uniform(0,20);
  for (i in 1:4)
    qtr_effect[i] ~ normal(32,10);
}

Why we represent continuous vectors as vector[length] varname and integer vectors as integer varname[length] I don’t quite get yet. But we’re done with the Stan code. I opened a text file and saved the above as flights_qtr.stan, to be called below.

Now we can proceed as the first post in this series, except now we need to define this qtr vector as the number of the quarter (from 1 to 4) that corresponds to each date and pass it to Stan. The lubridate package can make quick work of this. Remember to update data_to_pass and the file name of the Stan file.

ndates <- nrow(counts_depart)
flights <- counts_depart$n
qtr <- quarter(counts_depart$date)

data_to_pass <- c("ndates","flights","qtr")
qtr_model_fit <- stan("flights_qtr.stan",data=data_to_pass)
save(qtr_model_fit,file="qtr_model_fit.RData")

Analysis

As before, we examine the printout from qtr_model_fit:

qtr_model_fit

## Inference for Stan model: flights_qtr.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##                  mean se_mean   sd    2.5%     25%     50%     75%   97.5%
## sigma            3.13    0.00 0.12    2.90    3.04    3.12    3.21    3.37
## qtr_effect[1]   13.86    0.00 0.33   13.21   13.64   13.86   14.08   14.48
## qtr_effect[2]   20.20    0.00 0.32   19.59   19.99   20.20   20.42   20.84
## qtr_effect[3]   17.71    0.01 0.40   16.94   17.44   17.71   17.98   18.50
## qtr_effect[4]   17.52    0.00 0.33   16.88   17.30   17.52   17.73   18.16
## lp__          -553.91    0.03 1.55 -557.67 -554.72 -553.60 -552.79 -551.84
##               n_eff Rhat
## sigma          4972    1
## qtr_effect[1]  4671    1
## qtr_effect[2]  4959    1
## qtr_effect[3]  4811    1
## qtr_effect[4]  5140    1
## lp__           2223    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Jan 08 16:28:06 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

As we did last time, we look at traceplots before saying too much more. (Note also that Rhat=1 for each of the parameters above, so that looks pretty good.) Because we have a vector, and it would be a pain in the hindside to type out all of the qtr_effect[1] ... qtr_effect[4] parameters, we use the regex_pars option of mcmc_trace to examine the trace plots.

p <- mcmc_trace(as.array(qtr_model_fit),pars="sigma",regex_pars = "qtr_effect",
                facet_args = list(nrow = 5, labeller = label_parsed))
p

Things look pretty good, so we can proceed. We look at density plots as well, with the same bit of trickery with regex_pars as above:

p <- mcmc_areas(as.matrix(qtr_model_fit),pars="sigma",regex_pars = "qtr_effect",prob = 0.5) +
  ggtitle("Posterior distributions with medians and 50% intervals")
p

Now we can plainly see that in Quarter 2, the number of flights is much higher, averaging around 20 departures per day (perhaps due to all the summer vacations), with a pullback in Quarters 3 and 4, and with a lot less in Quarter 1.

Discussion

After having introduced the Greenville flight data, a simple model for analysis, and how to analyze output, we discussed how to exploit some of the structure in the data, specifically the fact that airlines tend to run flights by seasons and day of week. We saw how this manifested in the Stan code, working from the model block back to the data block. But there are storm clouds on the horizon, as we started getting some curious results and some more difficulty with convergence.