Supplemental Materials

This website serves as the supplemental materials for the paper titled “Hierarchical Point Process Models for Recurring Safety Critical Events Involving Commercial Truck Drivers: A Reliability Framework for Human Performance Modeling”, published on Journal of Quality Technology.

Two models are proposed in this paper:

1. A Bayesian hierarchical non-homogeneous Poisson process (NHPP) with the power law process (PLP) intensity function,
2. A Bayesian hierarchical jump power law process (JPLP)

The first section shows the R code to simulate data following the NHPP and JPLP data generating process. Only 10 drivers are assumed to minimize the simulation and Bayesian estimation time. Then, Stan code to perform Bayesian hierarchical PLP and JPLP estimation is presented in the next section. In the third section, R code for performing Hamiltonian Monte Carlo for the simulated data are demonstrated. In the last section, we present the R code to get summary statistics (posterior mean, 95% credible interval, Gelman-Rubin statistics \(\hat{R}\), and effective sample size ESS) from the posterior distribution are shown.

Simulate data

Simulate NHPP data with PLP intensity function

pacman::p_load(dplyr, rstan, broom)
source('Functions/NHPP_functions.R')
set.seed(123)

df_NHPP = sim_hier_nhpp(D = 10, beta = 1.2)
str(df_NHPP$hier_dat)

List of 9
 $ N           : int 255
 $ K           : num 3
 $ S           : int 95
 $ D           : int 10
 $ id          : int [1:95] 1 1 1 1 1 1 1 1 1 1 ...
 $ tau         : num [1:95] 11.09 9.62 10.42 8.22 10.38 ...
 $ event_time  : num [1:255] 1.04 3.29 4.15 5.16 6.49 ...
 $ group_size  : int [1:95] 8 7 1 0 6 16 5 2 3 8 ...
 $ X_predictors:'data.frame':   95 obs. of  3 variables:
  ..$ x1: num [1:95] 0.6651 -0.0857 0.9146 2.0706 0.8546 ...
  ..$ x2: num [1:95] 1.9446 0.0545 1.0107 2.6988 1.0172 ...
  ..$ x3: int [1:95] 0 3 3 1 1 1 2 4 2 2 ...

NHPP is estimated at shift-level. Here are the definition of the simulated NHPP data passed to Stan:

• N: the total number of events,
• K: the number of predictor variables for \(\theta\),
• S: the total number of shifts \(S\),
• D: the total number of drivers \(D\),
• id: driver ID for each shift \(d = \{1, 2, \ldots, D\}\),
• tau: shift end time \(\tau\),
• event_time: time of the events (SCEs) \(t_{d, s, i}\),
• group_size: the number of events for a driver in a shift \(n_{d, s}\),
• X_predictors: the predictor matrix.

Simulate JPLP data

source('Functions/JPLP_functions.R')
set.seed(123)
df_JPLP = sim_hier_JPLP(D = 10, beta = 1.2)
str(df_JPLP$stan_dt)

List of 11
 $ N           : int 162
 $ K           : num 3
 $ S           : int 331
 $ D           : num 10
 $ id          : int [1:331] 1 1 1 1 1 1 1 1 1 1 ...
 $ r_trip      : int [1:331] 1 2 3 4 5 1 2 3 4 1 ...
 $ t_trip_start: num [1:331] 0 1.66 4.71 6.44 9.16 0 2.94 5.29 6.79 0 ...
 $ t_trip_end  : num [1:331] 1.66 4.71 6.44 9.16 11.09 ...
 $ event_time  : num [1:162] 2.71 4.57 5.7 6.27 6.66 ...
 $ group_size  : int [1:331] 0 2 2 1 0 1 1 1 0 1 ...
 $ X_predictors: num [1:331, 1:3] 0.665 0.665 0.665 0.665 0.665 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:331] "1" "1.1" "1.2" "1.3" ...
  .. ..$ : chr [1:3] "x1" "x2" "x3"

In contrast to the NHPP, JPLP is estimated at trip level. Here are the definition of the simulated JPLP data passed to Stan:

• N: the total number of events,
• K: the number of predictor variables for \(\theta\),
• S: the total number of trips \(S\),
• D: the total number of drivers \(D\),
• id: driver ID for each shift \(d = \{1, 2, \ldots, D\}\),
• r_trip: the index of the trip within a shift \(r\) (the first one, second one, …)
• t_trip_start: trip start time \(a_{d, s, r-1}\),
• t_trip_end: trip end time \(a_{d, s, r}\),
• event_time: time of the events (SCEs) \(t_{d, s, i}\),
• group_size: the number of events for a driver in a shift \(n_{d, s}\),
• X_predictors: the predictor matrix.

Note that the trip start and end time are counted from the beginning of the shift, and the rest time is excluded from calculation.

Stan code

The Stan codes for both models are written using self-define likelihood functions. These likelihood functions have been derived in the manuscript.

Stan code to estimate a hierarchical PLP

functions{
  real nhpp_log(vector t, real beta, real theta, real tau){
    vector[num_elements(t)] loglik_part;
    real loglikelihood;
    for (i in 1:num_elements(t)){
      loglik_part[i] = log(beta) - beta*log(theta) + (beta - 1)*log(t[i]);
    }
    loglikelihood = sum(loglik_part) - (tau/theta)^beta;
    return loglikelihood;
  }
  real nhppnoevent_lp(real tau, real beta, real theta){
    real loglikelihood = - (tau/theta)^beta;
    return(loglikelihood);
  }
}
data {
  int<lower=1> N;                // total # of failures
  int<lower=1> K;                // number of predictors
  int<lower=1> S;                // total # of shifts
  int<lower=1> D;                // total # of drivers
  int<lower=1> id[S];            // driver index, must be an array
  vector<lower=0>[S] tau;        // truncated time
  vector<lower=0>[N] event_time; // failure time
  int group_size[S];             // group sizes
  matrix[S, K] X_predictors;     // predictor variable matrix
}
transformed data{
  matrix[S, K] X_centered;
  vector[K] X_means;
  for(k0 in 1:K){
    X_means[k0] = mean(X_predictors[, k0]);
    X_centered[,k0] = X_predictors[, k0] - X_means[k0];
  }
}
parameters{
  real mu0;             // hyperparameter: mean
  real<lower=0> sigma0; // hyperparameter: s.e.
  real<lower=0> beta;   // shape parameter
  vector[K] R1_K;       // fixed parameters each of K predictors
  vector[D] R0;         // random intercept for each of D drivers
}
model{
  int position = 1;
  vector[S] theta_temp;

  for (s0 in 1:S){
    theta_temp[s0] = exp(R0[id[s0]] + X_centered[s0,]*R1_K);
  }

  for (s1 in 1:S){
    if(group_size[s1] == 0) {
      target += nhppnoevent_lp(tau[s1], beta, theta_temp[s1]);
    }else{
      segment(event_time, position, group_size[s1]) ~ nhpp(beta, theta_temp[s1], tau[s1]);
      position += group_size[s1];
    }
  }
  beta ~ gamma(1, 1);
  R0 ~ normal(mu0, sigma0);
  R1_K  ~ normal(0, 10);
  mu0 ~ normal(0, 10);
  sigma0 ~ gamma(1, 1);
}
generated quantities{
  real mu0_true = mu0 - dot_product(X_means, R1_K);
  vector[D] R0_true = R0 - dot_product(X_means, R1_K);
  //real theta_correct = theta_temp - dot_product(X_centered, R1_K);
}

Stan code to estimate a hierarchical JPLP

functions{
  // LogLikelihood function for shifts with events (N_{event} > 0)
  real jplp_log(vector t_event, // time of SCEs
                real trip_start,
                real trip_end,
                int r,// trip index
                real beta,
                real theta,
                real kappa)
  {
    vector[num_elements(t_event)] loglik;
    real loglikelihood;
    for (i in 1:num_elements(t_event))
    {
      loglik[i] = (r - 1)*log(kappa) + log(beta) - beta*log(theta) + (beta - 1)*log(t_event[i]);
    }
    loglikelihood = sum(loglik) - kappa^(r - 1)*theta^(-beta)*(trip_end^beta - trip_start^beta);
    return loglikelihood;
  }
  // LogLikelihood function for shifts with no event (N_{event} = 0)
  real jplpoevent_lp(real trip_start,
                     real trip_end,
                     int r,
                     real beta,
                     real theta,
                     real kappa)
  {
    real loglikelihood = - kappa^(r - 1)*theta^(-beta)*(trip_end^beta - trip_start^beta);
    return(loglikelihood);
  }
}
data {
  int<lower=0> N; //total # of events
  int<lower=1> D; //total # of drivers
  int<lower=1> K; //number of predictors
  int<lower=0> S; //total # of trips, not shifts!!
  int<lower=1> id[S];//driver index, must be an array
  int r_trip[S];//index of trip $r$
  vector<lower=0>[S] t_trip_start;//trip start time
  vector<lower=0>[S] t_trip_end;//trip end time
  vector<lower=0>[N] event_time; //failure time
  int group_size[S]; //group sizes
  matrix[S, K] X_predictors;//predictor variable matrix
}
transformed data{
  matrix[S, K] X_centered;
  vector[K] X_means;
  for(k0 in 1:K){
    X_means[k0] = mean(X_predictors[, k0]);
    X_centered[,k0] = X_predictors[, k0] - X_means[k0];
  }
}
parameters{
  real mu0; // hyperparameter
  real<lower=0> sigma0;// hyperparameter
  real<lower=0> beta;
  real<lower=0, upper=1> kappa;
  vector[K] R1_K; // fixed parameters for K predictors
  vector[D] R0; // random intercept for D drivers
}
model{
  int position = 1;
  vector[S] theta_temp;

  for (s0 in 1:S){
    theta_temp[s0] = exp(R0[id[s0]] + X_centered[s0,]*R1_K);
  }

  for (s1 in 1:S){ // Likelihood estimation for JPLP based on trips, not shifts
    if(group_size[s1] == 0){
      target += jplpoevent_lp(t_trip_start[s1], t_trip_end[s1], r_trip[s1], beta, theta_temp[s1], kappa);
      }else{
      segment(event_time, position, group_size[s1]) ~ jplp_log(t_trip_start[s1], t_trip_end[s1], r_trip[s1], beta, theta_temp[s1], kappa);
      position += group_size[s1];
    }
  }
//PRIORS
  beta ~ gamma(1, 1);
  kappa ~ uniform(0, 1);
  R0 ~ normal(mu0, sigma0);
  R1_K  ~ normal(0, 10);
  mu0 ~ normal(0, 10);
  sigma0 ~ gamma(1, 1);
}
generated quantities{
  real mu0_true = mu0 - dot_product(X_means, R1_K);
  vector[D] R0_true = R0 - dot_product(X_means, R1_K);
  //real theta_correct = theta_temp - dot_product(X_centered, R1_K);
}

Bayesian estimation for simulated data

NHPP with PLP intensity function

fit_NHPP = stan("Stan/nhpp_plp_hierarchical.stan",
                chains = 4, iter = 5000, warmup = 1000, data = df_NHPP$hier_dat)

SAMPLING FOR MODEL 'nhpp_plp_hierarchical' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 1: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 1: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 1: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.905 seconds (Warm-up)
Chain 1:                7.388 seconds (Sampling)
Chain 1:                9.293 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'nhpp_plp_hierarchical' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 2: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 2: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 2: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.754 seconds (Warm-up)
Chain 2:                6.767 seconds (Sampling)
Chain 2:                8.521 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'nhpp_plp_hierarchical' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 3: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 3: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 3: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 2.019 seconds (Warm-up)
Chain 3:                6.751 seconds (Sampling)
Chain 3:                8.77 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'nhpp_plp_hierarchical' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 4: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 4: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 4: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.938 seconds (Warm-up)
Chain 4:                7.214 seconds (Sampling)
Chain 4:                9.152 seconds (Total)
Chain 4: 

JPLP

fit_JPLP = stan("Stan/jplp_hierarchical.stan",
                chains = 4, iter = 5000, warmup = 1000, data = df_JPLP$stan_dt)

SAMPLING FOR MODEL 'jplp_hierarchical' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 1: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 1: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 1: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 8.542 seconds (Warm-up)
Chain 1:                31.131 seconds (Sampling)
Chain 1:                39.673 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'jplp_hierarchical' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 2: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 2: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 2: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 7.733 seconds (Warm-up)
Chain 2:                29.298 seconds (Sampling)
Chain 2:                37.031 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'jplp_hierarchical' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.001 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 3: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 3: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 3: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 8.16 seconds (Warm-up)
Chain 3:                28.924 seconds (Sampling)
Chain 3:                37.084 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'jplp_hierarchical' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 4: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 4: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 4: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 8.453 seconds (Warm-up)
Chain 4:                27.719 seconds (Sampling)
Chain 4:                36.172 seconds (Total)
Chain 4: 

Posterior summary and diagnostic statistics

NHPP with PLP intensity function

Posterior mean, 95% credible interval, ESS, and \(\hat{R}\)

est_NHPP = broom.mixed::tidy(fit_NHPP, conf.int = T, rhat = T, ess = T)
est_NHPP

# A tibble: 27 x 7
   term    estimate std.error conf.low conf.high  rhat   ess
   <chr>      <dbl>     <dbl>    <dbl>     <dbl> <dbl> <int>
 1 mu0        2.01     0.207     1.60      2.43   1.00 12477
 2 sigma0     0.527    0.183     0.299     1.01   1.00  9779
 3 beta       1.26     0.0785    1.12      1.42   1.00  9832
 4 R1_K[1]    0.933    0.0929    0.764     1.13   1.00  9549
 5 R1_K[2]    0.271    0.0706    0.139     0.415  1.00 16499
 6 R1_K[3]    0.211    0.0504    0.116     0.313  1.00 14545
 7 R0[1]      1.64     0.128     1.38      1.88   1.00 12865
 8 R0[2]      1.69     0.150     1.40      1.99   1.00 14187
 9 R0[3]      2.87     0.301     2.37      3.55   1.00 13207
10 R0[4]      1.94     0.181     1.60      2.31   1.00 20213
# ... with 17 more rows

Trace plots for selected parameters

t_NHPP = stan_trace(fit_NHPP, pars = c('mu0', 'sigma0', 'beta'), ncol = 1)
t_NHPP

JPLP

Posterior mean, 95% credible interval, ESS, and \(\hat{R}\)

est_JPLP = broom.mixed::tidy(fit_JPLP, conf.int = T, rhat = T, ess = T)
est_JPLP

# A tibble: 28 x 7
   term    estimate std.error conf.low conf.high  rhat   ess
   <chr>      <dbl>     <dbl>    <dbl>     <dbl> <dbl> <int>
 1 mu0        2.11     0.222    1.68       2.57   1.00  6575
 2 sigma0     0.506    0.183    0.266      0.987  1.00  6908
 3 beta       1.18     0.119    0.969      1.43   1.00  5024
 4 kappa      0.800    0.0787   0.651      0.955  1.00  4704
 5 R1_K[1]    0.995    0.136    0.762      1.29   1.00  4210
 6 R1_K[2]    0.216    0.0915   0.0507     0.411  1.00  9781
 7 R1_K[3]    0.235    0.0693   0.110      0.381  1.00  7657
 8 R0[1]      1.77     0.187    1.42       2.16   1.00  6752
 9 R0[2]      1.99     0.220    1.60       2.46   1.00  6687
10 R0[3]      2.55     0.299    2.04       3.22   1.00  6678
# ... with 18 more rows

Trace plots for selected parameters

t_JPLP = stan_trace(fit_JPLP, pars = c('mu0', 'sigma0', 'beta', 'kappa'), ncol = 1)
t_JPLP

Supplementary trace plots for selected parameters in the manuscript

Trace plots below are generated from the 496 large commercial truck drivers, which is demonstrated as a case study in the main manuscript.

knitr::include_graphics("Figures/Aim3_trace_plot.jpeg")

Session Information

sessionInfo()

R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=Chinese (Simplified)_China.936 
[2] LC_CTYPE=Chinese (Simplified)_China.936   
[3] LC_MONETARY=Chinese (Simplified)_China.936
[4] LC_NUMERIC=C                              
[5] LC_TIME=Chinese (Simplified)_China.936    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
[1] broom_0.7.6          rstan_2.21.2         ggplot2_3.3.3       
[4] StanHeaders_2.21.0-7 dplyr_1.0.6         

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.6         lattice_0.20-44    tidyr_1.1.3       
 [4] prettyunits_1.1.1  ps_1.6.0           assertthat_0.2.1  
 [7] digest_0.6.27      utf8_1.2.1         V8_3.4.2          
[10] R6_2.5.0           plyr_1.8.6         backports_1.2.1   
[13] stats4_4.1.0       evaluate_0.14      coda_0.19-4       
[16] highr_0.9          pillar_1.6.1       rlang_0.4.11      
[19] curl_4.3.1         rstudioapi_0.13    callr_3.7.0       
[22] jquerylib_0.1.4    Matrix_1.3-3       rmarkdown_2.8     
[25] labeling_0.4.2     stringr_1.4.0      TMB_1.7.20        
[28] loo_2.4.1          munsell_0.5.0      compiler_4.1.0    
[31] xfun_0.23          pkgconfig_2.0.3    pkgbuild_1.2.0    
[34] htmltools_0.5.1.1  broom.mixed_0.2.6  downlit_0.2.1     
[37] tidyselect_1.1.1   tibble_3.1.2       gridExtra_2.3     
[40] codetools_0.2-18   matrixStats_0.59.0 fansi_0.5.0       
[43] crayon_1.4.1       withr_2.4.2        grid_4.1.0        
[46] nlme_3.1-152       jsonlite_1.7.2     gtable_0.3.0      
[49] lifecycle_1.0.0    DBI_1.1.1          magrittr_2.0.1    
[52] scales_1.1.1       RcppParallel_5.1.4 cli_2.5.0         
[55] stringi_1.6.2      farver_2.1.0       reshape2_1.4.4    
[58] bslib_0.2.5.1      ellipsis_0.3.2     generics_0.1.0    
[61] vctrs_0.3.8        distill_1.2        tools_4.1.0       
[64] glue_1.4.2         purrr_0.3.4        processx_3.5.2    
[67] parallel_4.1.0     yaml_2.2.1         inline_0.3.19     
[70] colorspace_2.0-1   knitr_1.33         sass_0.4.0