11 Feb 2021

Overview of today’s material

  • Using STAN for MAP estimation
  • Multivariate time series models
  • DFA models
  • Writing our own Stan code

MAP (maximum a posteriori) estimation

  • Point estimate of unknown parameters
  • Optimization algorithms: BFGS, Newton, etc
  • Estimates similar to maximum likelihood, but also incorporate prior
  • Quick for checking models, etc

MAP estimation in Stan

  • We left off with a Bayesian DLM with random walk in intercept
bayes_fit = fit_stan(y = SalmonSurvCUI$logit.s, 
               model_name="dlm-intercept")

MAP estimation in Stan

  • Switching this to MAP estimation
map_fit = fit_stan(y = SalmonSurvCUI$logit.s, 
               model_name="dlm-intercept", map_estimation=TRUE)

MAP estimation in Stan

print(map_fit)
## $par
##            x0       beta[1]    pro_dev[1]    pro_dev[2]    pro_dev[3] 
## -3.460000e+00 -1.503290e+00  1.046740e-01 -1.943951e-01  4.112205e-01 
##    pro_dev[4]    pro_dev[5]    pro_dev[6]    pro_dev[7]    pro_dev[8] 
## -4.336506e-01  1.943952e-01 -4.336507e-01 -1.943952e-01 -4.710349e-01 
##    pro_dev[9]   pro_dev[10]   pro_dev[11]   pro_dev[12]   pro_dev[13] 
## -6.205692e-01  7.850578e-01  8.074876e-01 -1.659836e+00  4.186974e-01 
##   pro_dev[14]   pro_dev[15]   pro_dev[16]   pro_dev[17]   pro_dev[18] 
##  4.261744e-01 -2.542092e-01 -1.719652e-01  1.345815e-01  5.607551e-01 
##   pro_dev[19]   pro_dev[20]   pro_dev[21]   pro_dev[22]   pro_dev[23] 
## -2.317783e-01 -6.729076e-02  3.514065e-01 -8.074876e-01  8.972066e-02 
##   pro_dev[24]   pro_dev[25]   pro_dev[26]   pro_dev[27]   pro_dev[28] 
##  1.271049e-01 -1.869187e-01  6.729066e-01 -7.476740e-01 -6.729063e-01 
##   pro_dev[29]   pro_dev[30]   pro_dev[31]   pro_dev[32]   pro_dev[33] 
##  6.729059e-02  8.972085e-01 -5.233691e-02 -5.308487e-01  1.375720e+00 
##   pro_dev[34]   pro_dev[35]   pro_dev[36]   pro_dev[37]   pro_dev[38] 
## -8.859995e-08  4.261745e-01 -1.046745e-01 -2.841160e-01 -4.486054e-02 
##   pro_dev[39]   pro_dev[40]   pro_dev[41] sigma_process     sigma_obs 
## -8.897319e-01  2.542093e-01 -5.607554e-01  1.337481e+00  5.537342e-07 
##       pred[1]       pred[2]       pred[3]       pred[4]       pred[5] 
## -3.460000e+00 -3.320000e+00 -3.580000e+00 -3.030000e+00 -3.610000e+00 
##       pred[6]       pred[7]       pred[8]       pred[9]      pred[10] 
## -3.350000e+00 -3.930000e+00 -4.190000e+00 -4.820000e+00 -5.650000e+00 
##      pred[11]      pred[12]      pred[13]      pred[14]      pred[15] 
## -4.600000e+00 -3.520000e+00 -5.740000e+00 -5.180000e+00 -4.610000e+00 
##      pred[16]      pred[17]      pred[18]      pred[19]      pred[20] 
## -4.950000e+00 -5.180000e+00 -5.000000e+00 -4.250000e+00 -4.560000e+00 
##      pred[21]      pred[22]      pred[23]      pred[24]      pred[25] 
## -4.650000e+00 -4.180000e+00 -5.260000e+00 -5.140000e+00 -4.970000e+00 
##      pred[26]      pred[27]      pred[28]      pred[29]      pred[30] 
## -5.220000e+00 -4.320000e+00 -5.320000e+00 -6.220000e+00 -6.130000e+00 
##      pred[31]      pred[32]      pred[33]      pred[34]      pred[35] 
## -4.930000e+00 -5.000000e+00 -5.710000e+00 -3.870000e+00 -3.870000e+00 
##      pred[36]      pred[37]      pred[38]      pred[39]      pred[40] 
## -3.300000e+00 -3.440000e+00 -3.820000e+00 -3.880000e+00 -5.070000e+00 
##      pred[41]      pred[42]  intercept[1]  intercept[2]  intercept[3] 
## -4.730000e+00 -5.480000e+00 -3.460000e+00 -3.320000e+00 -3.580000e+00 
##  intercept[4]  intercept[5]  intercept[6]  intercept[7]  intercept[8] 
## -3.030000e+00 -3.610000e+00 -3.350000e+00 -3.930000e+00 -4.190000e+00 
##  intercept[9] intercept[10] intercept[11] intercept[12] intercept[13] 
## -4.820000e+00 -5.650000e+00 -4.600000e+00 -3.520000e+00 -5.740000e+00 
## intercept[14] intercept[15] intercept[16] intercept[17] intercept[18] 
## -5.180000e+00 -4.610000e+00 -4.950000e+00 -5.180000e+00 -5.000000e+00 
## intercept[19] intercept[20] intercept[21] intercept[22] intercept[23] 
## -4.250000e+00 -4.560000e+00 -4.650000e+00 -4.180000e+00 -5.260000e+00 
## intercept[24] intercept[25] intercept[26] intercept[27] intercept[28] 
## -5.140000e+00 -4.970000e+00 -5.220000e+00 -4.320000e+00 -5.320000e+00 
## intercept[29] intercept[30] intercept[31] intercept[32] intercept[33] 
## -6.220000e+00 -6.130000e+00 -4.930000e+00 -5.000000e+00 -5.710000e+00 
## intercept[34] intercept[35] intercept[36] intercept[37] intercept[38] 
## -3.870000e+00 -3.870000e+00 -3.300000e+00 -3.440000e+00 -3.820000e+00 
## intercept[39] intercept[40] intercept[41] intercept[42]    log_lik[1] 
## -3.880000e+00 -5.070000e+00 -4.730000e+00 -5.480000e+00  1.326919e+01 
##    log_lik[2]    log_lik[3]    log_lik[4]    log_lik[5]    log_lik[6] 
##  1.348558e+01  1.348143e+01  1.343850e+01  1.346689e+01  1.348198e+01 
##    log_lik[7]    log_lik[8]    log_lik[9]   log_lik[10]   log_lik[11] 
##  1.340442e+01  1.338355e+01  1.341164e+01  1.348663e+01  1.337676e+01 
##   log_lik[12]   log_lik[13]   log_lik[14]   log_lik[15]   log_lik[16] 
##  1.348383e+01  1.341941e+01  1.342276e+01  1.346445e+01  1.348422e+01 
##   log_lik[17]   log_lik[18]   log_lik[19]   log_lik[20]   log_lik[21] 
##  1.335915e+01  1.348345e+01  1.326469e+01  1.317618e+01  1.334588e+01 
##   log_lik[22]   log_lik[23]   log_lik[24]   log_lik[25]   log_lik[26] 
##  1.348756e+01  1.340756e+01  1.348373e+01  1.326981e+01  1.348561e+01 
##   log_lik[27]   log_lik[28]   log_lik[29]   log_lik[30]   log_lik[31] 
##  1.346174e+01  1.348731e+01  1.337548e+01  1.343007e+01  1.348366e+01 
##   log_lik[32]   log_lik[33]   log_lik[34]   log_lik[35]   log_lik[36] 
##  1.333669e+01  1.348752e+01  1.342055e+01  1.331926e+01  1.342071e+01 
##   log_lik[37]   log_lik[38]   log_lik[39]   log_lik[40]   log_lik[41] 
##  1.348268e+01  1.348710e+01  1.346313e+01  1.346517e+01  1.347789e+01 
##   log_lik[42] 
##  1.346525e+01 
## 
## $value
## [1] 595.6098
## 
## $return_code
## [1] 0
## 
## $theta_tilde
##         x0  beta[1] pro_dev[1] pro_dev[2] pro_dev[3] pro_dev[4] pro_dev[5]
## [1,] -3.46 -1.50329   0.104674 -0.1943951  0.4112205 -0.4336506  0.1943952
##      pro_dev[6] pro_dev[7] pro_dev[8] pro_dev[9] pro_dev[10] pro_dev[11]
## [1,] -0.4336507 -0.1943952 -0.4710349 -0.6205692   0.7850578   0.8074876
##      pro_dev[12] pro_dev[13] pro_dev[14] pro_dev[15] pro_dev[16] pro_dev[17]
## [1,]   -1.659836   0.4186974   0.4261744  -0.2542092  -0.1719652   0.1345815
##      pro_dev[18] pro_dev[19] pro_dev[20] pro_dev[21] pro_dev[22] pro_dev[23]
## [1,]   0.5607551  -0.2317783 -0.06729076   0.3514065  -0.8074876  0.08972066
##      pro_dev[24] pro_dev[25] pro_dev[26] pro_dev[27] pro_dev[28] pro_dev[29]
## [1,]   0.1271049  -0.1869187   0.6729066   -0.747674  -0.6729063  0.06729059
##      pro_dev[30] pro_dev[31] pro_dev[32] pro_dev[33]   pro_dev[34] pro_dev[35]
## [1,]   0.8972085 -0.05233691  -0.5308487     1.37572 -8.859995e-08   0.4261745
##      pro_dev[36] pro_dev[37] pro_dev[38] pro_dev[39] pro_dev[40] pro_dev[41]
## [1,]  -0.1046745   -0.284116 -0.04486054  -0.8897319   0.2542093  -0.5607554
##      sigma_process    sigma_obs pred[1] pred[2] pred[3] pred[4] pred[5] pred[6]
## [1,]      1.337481 5.537342e-07   -3.46   -3.32   -3.58   -3.03   -3.61   -3.35
##      pred[7] pred[8] pred[9] pred[10] pred[11] pred[12] pred[13] pred[14]
## [1,]   -3.93   -4.19   -4.82    -5.65     -4.6    -3.52    -5.74    -5.18
##      pred[15] pred[16] pred[17] pred[18] pred[19] pred[20] pred[21] pred[22]
## [1,]    -4.61    -4.95    -5.18       -5    -4.25    -4.56    -4.65    -4.18
##      pred[23] pred[24] pred[25] pred[26] pred[27] pred[28] pred[29] pred[30]
## [1,]    -5.26    -5.14    -4.97    -5.22    -4.32    -5.32    -6.22    -6.13
##      pred[31] pred[32] pred[33] pred[34] pred[35] pred[36] pred[37] pred[38]
## [1,]    -4.93       -5    -5.71    -3.87    -3.87     -3.3    -3.44    -3.82
##      pred[39] pred[40] pred[41] pred[42] intercept[1] intercept[2] intercept[3]
## [1,]    -3.88    -5.07    -4.73    -5.48        -3.46        -3.32        -3.58
##      intercept[4] intercept[5] intercept[6] intercept[7] intercept[8]
## [1,]        -3.03        -3.61        -3.35        -3.93        -4.19
##      intercept[9] intercept[10] intercept[11] intercept[12] intercept[13]
## [1,]        -4.82         -5.65          -4.6         -3.52         -5.74
##      intercept[14] intercept[15] intercept[16] intercept[17] intercept[18]
## [1,]         -5.18         -4.61         -4.95         -5.18            -5
##      intercept[19] intercept[20] intercept[21] intercept[22] intercept[23]
## [1,]         -4.25         -4.56         -4.65         -4.18         -5.26
##      intercept[24] intercept[25] intercept[26] intercept[27] intercept[28]
## [1,]         -5.14         -4.97         -5.22         -4.32         -5.32
##      intercept[29] intercept[30] intercept[31] intercept[32] intercept[33]
## [1,]         -6.22         -6.13         -4.93            -5         -5.71
##      intercept[34] intercept[35] intercept[36] intercept[37] intercept[38]
## [1,]         -3.87         -3.87          -3.3         -3.44         -3.82
##      intercept[39] intercept[40] intercept[41] intercept[42] log_lik[1]
## [1,]         -3.88         -5.07         -4.73         -5.48   13.26919
##      log_lik[2] log_lik[3] log_lik[4] log_lik[5] log_lik[6] log_lik[7]
## [1,]   13.48558   13.48143    13.4385   13.46689   13.48198   13.40442
##      log_lik[8] log_lik[9] log_lik[10] log_lik[11] log_lik[12] log_lik[13]
## [1,]   13.38355   13.41164    13.48663    13.37676    13.48383    13.41941
##      log_lik[14] log_lik[15] log_lik[16] log_lik[17] log_lik[18] log_lik[19]
## [1,]    13.42276    13.46445    13.48422    13.35915    13.48345    13.26469
##      log_lik[20] log_lik[21] log_lik[22] log_lik[23] log_lik[24] log_lik[25]
## [1,]    13.17618    13.34588    13.48756    13.40756    13.48373    13.26981
##      log_lik[26] log_lik[27] log_lik[28] log_lik[29] log_lik[30] log_lik[31]
## [1,]    13.48561    13.46174    13.48731    13.37548    13.43007    13.48366
##      log_lik[32] log_lik[33] log_lik[34] log_lik[35] log_lik[36] log_lik[37]
## [1,]    13.33669    13.48752    13.42055    13.31926    13.42071    13.48268
##      log_lik[38] log_lik[39] log_lik[40] log_lik[41] log_lik[42]
## [1,]     13.4871    13.46313    13.46517    13.47789    13.46525

MAP estimation in Stan

  • Check that model converged
map_fit$return_code
## [1] 0
  • MAP value when converged
map_fit$value
## [1] 595.6098

MAP estimation in Stan

  • grep or other string matching functions needed
grep("pred",names(map_fit$par))
##  [1] 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
## [26] 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
pred = map_fit$par[grep("pred",names(map_fit$par))]

MAP estimation in Stan

  • grep or other string matching functions needed
grep("pred",names(map_fit$par))
##  [1] 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
## [26] 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
pred = map_fit$par[grep("pred",names(map_fit$par))]

MAP estimation in Stan

  • grep or other string matching functions needed
grep("pred",names(map_fit$par))
##  [1] 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
## [26] 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
pred = map_fit$par[grep("pred",names(map_fit$par))]

MAP estimation in Stan

  • include SEs of estimates, via Hessian
map_fit = fit_stan(y = SalmonSurvCUI$logit.s, 
     model_name="dlm-intercept", 
     map_estimation=TRUE, 
     hessian=TRUE)
  • or change algorithm (“LBFGS” default) to “BFGS” or “Newton”
map_fit = fit_stan(y = SalmonSurvCUI$logit.s, 
     model_name="dlm-intercept", 
     map_estimation=TRUE, 
     hessian=TRUE, 
     algorithm="BFGS")

MAP estimation in Stan

  • Posterior means, 95% CIs

MAP estimation in Stan

  • MAP estimates = posterior

Multivariate time series models

  • Multivariate state space models
  • Stan package for fitting models similar to MARSS
devtools::install_github("atsa-es/tvvarss")

Multivariate time series models

  • Most models in tvvarss() are exactly the same as MARSS
  • Takes in data matrix y
  • family argument defaults to Gaussian, but can be many others
tvvarss::tvvarss(y = y, family = "poisson",...)

Multivariate time series models

  • Like MARSS, variances can be shared (default) or unique by species, or time series
tvvarss::tvvarss(y = y, shared_r = R, shared_q = Q,...)

Multivariate time series models

  • Optional ‘process’ argument acts as Z in MARSS and maps time series to latent state processes
tvvarss::tvvarss(y = y, shared_r = R, shared_q = Q,...)

Multivariate time series models

  • Perhaps most exciting feature is time-varying interactions (AR coefficients)
  • Time-varying vector autoregressive models
  • AR coefficients behave like a DLM
tvvarss::tvvarss(y = y, dynamicB = TRUE,...)
  • Time-varying parameters are data hungry! We’ll dive into this more later in the quarter

Bayesian DFA models

  • DFA models from atsar bundled with other DFA code we’ve developed
devtools::install_github("fate-ewi/bayesdfa")

Bayesian DFA models

  • DFA poses interesting identifiability challenges
  • Bayesian models generally involve > 1 MCMC chain
  • Direction of DFA loadings / trends don’t have meaning

Bayesian DFA models

Bayesian DFA models

  • Solution 1: priors for identifiability

  • Solution 2: post-hoc ‘chain flipping’ before convergence tests run

Bayesian DFA models

  • Fitting DFA models will give very similar answers to using MARSS
data("harborSealWA")

Bayesian DFA models

We’ll extract predictions from the best model,

fit = bayesdfa::fit_dfa(y = t(harborSealWA[,-1]), num_trends = 1)

And as an stanfit object, we can extract summaries of states x

pars = extract(fit$model)

Bayesian DFA models

bayesdfa::plot_fitted(fit) + theme_bw()

Bayesian DFA models

  • DFA extension 1: long format data (optional), replicate observations

Bayesian DFA models

  • DFA extension 2: temporal extremes

\[x_{t}=x_{t-1}+\delta_{t-1}\] \[\delta_{t-1} \sim Normal(0,q)\]

  • But do these deviations have to be normal? NO!

\[\delta_{t-1} \sim Student-t(\nu,0,q)\]

Bayesian DFA models

Bayesian DFA models

  • Extremes example: coho salmon body size in coastal fisheries

Bayesian DFA models

  • Fitting the model
coho = readRDS("coho_mean_length.rds")
# rename variables
coho=dplyr::rename(coho, obs=mean_length,time=recovery_year,ts=fishery)
fit = fit_dfa(y = coho, data_shape="long",estimate_nu = TRUE)

Bayesian DFA models

  • Let’s look at df parameter, \(\nu\)
  • Mean ~ 19.5, median ~ 16.3
pars = extract(fit$model)

Bayesian DFA models

  • \(\nu\) parameter ranges from 2 (heavy tails) to normal (30)

(source: Wikipedia)

Bayesian DFA models

  • DFA extension 2: increased flexibility in trends
  • Sampling might be ‘chunky’

Bayesian DFA models

  • AR models (conventional DFA) may not be well suited

  • Alternative approach: use Gaussian Process models to model smooth trends

  • Feddern et al. 2021 (Global Change Biology)

Bayesian DFA models

  • What’s a Gaussian Process model?

  • Instead of modeling process \[x_{t} = x_{t-1}...\]

  • We model the covariance of the \(\textbf{x}\)

  • But there are a lot of \(\textbf{x}\)!!

Bayesian DFA models

  • For model with \(T\) timesteps

  • \(\Sigma\) has \(T*(T-1)/2\) parameters

  • GP models implement covariance functions, e.g. exponential, Gaussian, Matern

  • Model covariance as function of space or distance in time

Bayesian DFA models

  • Example: Gaussian covariance (aka ‘squared exponential’)

\[\Sigma_{t,t+3} = \sigma^2exp(\frac{(t-(t+3))^2}{2\tau^2})\]

  • \(\sigma^2\) controls variability
  • \(\tau\) controls how quickly covariance decays between points

Bayesian DFA models

  • Implementation
fit_dfa(y= y, num_trends = 2, ..., trend_model = c("gp"))

Bayesian DFA models

  • DFA extension 3: alternate constraints on loadings matrix

  • Conventional constraints on loadings matrix

##             [,1]       [,2]       [,3]
## [1,] -0.56047565  0.0000000  0.0000000
## [2,] -0.23017749  0.4609162  0.0000000
## [3,]  1.55870831 -1.2650612  0.4007715
## [4,]  0.07050839 -0.6868529  0.1106827
## [5,]  0.12928774 -0.4456620 -0.5558411

Bayesian DFA models

  • Alternative: compositional DFA model
##           [,1]       [,2]       [,3]
## [1,] 0.6871757 0.24516978 0.06765452
## [2,] 0.1491745 0.06014499 0.79068047
## [3,] 0.3885419 0.59306319 0.01839489
## [4,] 0.2919682 0.30545457 0.40257721
## [5,] 0.3760087 0.33262301 0.29136825

Bayesian DFA models

  • Under compositional DFA model, time series are true mixtures of underlying trends

\[\textbf{Y} = \textbf{Zx}\]

fit_dfa(y = y, num_trends = 2, ..., z_model = "proportion")
  • When might this be useful?
  • Stable isotope data, environmental monitoring, pollutant data, etc

Writing our own Stan scripts

  • Stan scripts always start with a data block
  • Data needs to be typed
data {
  int<lower=0> N;
  int<lower=0> K;
  real y[N];
  int P;
  int y_int[N];
  matrix[N, K] x;
}

Writing our own Stan scripts

  • transformed data block is optional
transformed data {
  int zeros[N];
  for(i in 1:N) {
    zeros[i] = 0;
  }
}

Writing our own Stan scripts

  • parameters block containts any parameters
parameters {
  real x0;
  vector[K] beta0;
  vector[K] pro_dev[N-1];
  real<lower=0> sigma_process[K];
  real<lower=0> sigma_obs;
}

Writing our own Stan scripts

  • parameters block containts any parameters
parameters {
  real x0;
  vector[K] beta0;
  vector[K] pro_dev[N-1];
  real<lower=0> sigma_process[K];
  real<lower=0> sigma_obs;
}

Writing our own Stan scripts

  • transformed parameters block contains derived quantities
  • examples: predicted states for state space time series model
transformed parameters {
  vector[N] pred;
  pred[1] = x0;
  for(i in 2:N) {
    pred[i] = phi*pred[i-1] + sigma_process*pro_dev[i-1];
  }
}

Writing our own Stan scripts

  • model block contains priors on parameters
  • likelihood (or data model)
model {
  x0 ~ normal(0,10);
  phi ~ normal(0,1);
  sigma_process ~ student_t(3,0,2);
  sigma_obs ~ student_t(3,0,2);
  pro_dev ~ std_normal();
  for(i in 1:N) {
    y[i] ~ normal(pred[i], sigma_obs);
  }
}

Writing our own Stan scripts

  • generated data block (optional)
  • includes quantities useful for model selection, prediction, forecasts, etc
generated quantities {
  vector[N] log_lik;
  for (n in 1:N) log_lik[n] = normal_lpdf(y[n] | pred[n], sigma_obs);
}

Writing our own Stan scripts

  • Ok – let’s do this in practice

  • Problem # 1: let’s take our state space random walk model and modify it to better handle extreme events with the Student - t distribution

  • Mathematically

\[\delta_{t-1} \sim Normal(0,q)\] becomes

\[\delta_{t-1} \sim Student-t(\nu,0,q)\]

Writing our own Stan scripts

Writing our own Stan scripts

Writing our own Stan scripts

  • 2 things need to change

  • Change process error deviations from Normal to Student-t distribution

  • Add \(\nu\) as a parameter, with constraints (> 2)

Writing our own Stan scripts

parameters {
real<lower=2> nu;
...
}

Writing our own Stan scripts

model{
...
nu ~ student_t(3, 2, 3);
pro_dev ~ student_t(nu, 0, 1);
//pro_dev ~ std_normal();
...
}

Writing our own Stan scripts

  • Done! Now we can fit the model
y = gme$Close
N = length(y)
n_pos = length(which(!is.na(y)))
pos_indx = c(which(!is.na(y)),0,0)

fit = stan(file = "ss_ar_t.stan", 
  data = list(y = y, N = N, n_pos = n_pos, pos_indx = pos_indx))

Writing our own Stan scripts