References

Load packages

library(tidyverse)
library(ggdag)
library(biostan)
## Warning: replacing previous import 'rstan::loo' by 'rstanarm::loo' when loading 'biostan'
library(rstan)

Load data

The dataset is a case-control study on the outcome of invasive cervical cancer and the exposure of herpes simplex virus (HSV). Let \(Y\) be the binary outcome, \(X_{true}\) be the true binary covariate, and \(X_{mis}\) be the mismeasured binary covariate. The dataset contains an internal validation dataset with both \(X_{true}\) and \(X_{mis}\). The rest of the dataset only has \(X_{mis}\). Let \(R_{X}\) be the indicator to indicate that \(X_{true}\) was observed for that row.

cervical <- tribble(
    ~Y, ~X_true, ~X_mis, ~count, ~R_X,
    ## Complete data for the internal validation set
    1, 0, 0, 13, 1,
    1, 0, 1,  3, 1,
    1, 1, 0,  5, 1,
    1, 1, 1, 18, 1,
    0, 0, 0, 33, 1,
    0, 0, 1, 11, 1,
    0, 1, 0, 16, 1,
    0, 1, 1, 16, 1,
    ## Incomplete data for the rest of the dataset
    1, NA, 0, 318, 0,
    1, NA, 1, 375, 0,
    0, NA, 0, 701, 0,
    0, NA, 1, 535, 0
)
cervical
## # A tibble: 12 x 5
##        Y X_true X_mis count   R_X
##    <dbl>  <dbl> <dbl> <dbl> <dbl>
##  1     1      0     0    13     1
##  2     1      0     1     3     1
##  3     1      1     0     5     1
##  4     1      1     1    18     1
##  5     0      0     0    33     1
##  6     0      0     1    11     1
##  7     0      1     0    16     1
##  8     0      1     1    16     1
##  9     1     NA     0   318     0
## 10     1     NA     1   375     0
## 11     0     NA     0   701     0
## 12     0     NA     1   535     0
## Stan does not take NAs.
cervical <- cervical %>%
    mutate(X_true = if_else(is.na(X_true), -1, X_true))

Model

The model considered is a classical mismeasurement model represented in the following model DAG that also contains parameters. This type of DAGs with parameters in them are commonly seen in the Bayesian literature, in particular to describe hierarchical models. Here we also have the missing indicator assuming MCAR (validation subset is a random subsample within the sample).

dagify(Y ~ Xt, Xm ~ Xt, Xm ~ phi, Xt ~ psi, Y ~ beta, Rx ~ gamma) %>%
    ggdag()

Since the missing process is ignorable, we can work on the other parts only. We partition the true exposure \(X_{true}\) into the complete part \(X_{true}^{c}\) and the reduced part \(X_{true}^{r}\) (missing part). Let \(i=1,\dots,m\) index the reduced part and \(i=m+1,\dots,n\) index the complete part.

\[\begin{align*} p(X_{true}^{r}, \beta, \phi, \psi | Y, X_{mis}, X_{true}^{c}) &\propto p(\beta, \phi, \psi, Y, X_{mis}, X_{true})\\ \\ &= \underbrace{p(\beta, \phi, \psi)}_{\text{prior}} \underbrace{p(Y, X_{mis}, X_{true} | \beta, \phi, \psi)}_{\text{likelihood}}\\ \\ &~~~~\text{By iid given parameters}\\ &= p(\beta, \phi, \psi) \prod^{n}_{i=1} p(Y_{i}, X_{mis,i}, X_{true,i} | \beta, \phi, \psi)\\ \\ &~~~~\text{Partition again}\\ &= p(\beta, \phi, \psi)\\ &~~~\times \prod^{m}_{i=1} p(Y_{i}, X_{mis,i}^{r}, X_{true,i} | \beta, \phi, \psi)\\ &~~~\times \prod^{n}_{j=m+1} p(Y_{j}, X_{mis,j}^{c}, X_{true,j} | \beta, \phi, \psi)\\ \\ &~~~~\text{Factor likelihood following DAG}\\ &= p(\beta, \phi, \psi)\\ &~~~\times \prod^{m}_{i=1} \underbrace{p(Y_{i} | X_{true,i}^{r}, \beta)}_{\text{outcome model}} \underbrace{p(X_{mis,i} | X_{true,i}^{r}, \phi)}_{\text{error model}} \underbrace{p(X_{true,i}^{r} | \psi)}_{\text{covariate model}}\\ &~~~\times \prod^{n}_{j=m+1} \underbrace{p(Y_{j} | X_{true,j}^{c}, \beta)}_{\text{outcome model}} \underbrace{p(X_{mis,j} | X_{true,j}^{c}, \phi)}_{\text{error model}} \underbrace{p(X_{true,j}^{c} | \psi)}_{\text{covariate model}}\\ \end{align*} \]

Since we do not observe \(X_{true}^{r}\), we need integrate out (really sum out here) this part from the joint posterior. That is,

\[\begin{align*} p(\beta, \phi, \psi | Y, X_{mis}, X_{true}^{c}) &= \int p(X_{true}^{r}=x_{true}^{r}, \beta, \phi, \psi | Y, X_{mis}, X_{true}^{c}) d\mu(x_{true}^{r}) \end{align*} \]

The only part that depends on \(X_{true}^{r}\) in the prior \(\times\) likelihood expression is the following likelihood for the reduced part.

\[\begin{align*} \prod^{m}_{i=1} \underbrace{p(Y_{i} | X_{true,i}^{r}, \beta)}_{\text{outcome model}} \underbrace{p(X_{mis,i} | X_{true,i}^{r}, \phi)}_{\text{error model}} \underbrace{p(X_{true,i}^{r} | \psi)}_{\text{covariate model}}\\ \end{align*} \]

Conditioning on the parameters, \(X_{true,i}^{r}\) are independent from each other. Thus, we can sum out each term in the product.

\[\begin{align*} &\prod^{m}_{i=1} \sum^{1}_{x=0} p(Y_{i} | X_{true,i}^{r}=x, \beta) p(X_{mis,i} | X_{true,i}^{r}=x, \phi) p(X_{true,i}^{r}=x | \psi)\\ &= \prod^{m}_{i=1} \left[ p(Y_{i} | X_{true,i}^{r}=0, \beta) p(X_{mis,i} | X_{true,i}^{r}=0, \phi) p(X_{true,i}^{r}=0 | \psi) \right.\\ &~~~~~~~~~~~\left. + p(Y_{i} | X_{true,i}^{r}=1, \beta) p(X_{mis,i} | X_{true,i}^{r}=1, \phi) p(X_{true,i}^{r}=1 | \psi) \right]\\ \\ &= \prod^{m}_{i=1} \left[ p(Y_{i} | X_{true,i}^{r}=0, \beta) p(X_{mis,i} | X_{true,i}^{r}=0, \phi) (1-\psi) \right.\\ &~~~~~~~~~~~\left. + p(Y_{i} | X_{true,i}^{r}=1, \beta) p(X_{mis,i} | X_{true,i}^{r}=1, \phi) (\psi) \right]\\ \end{align*} \]

Therefore, the marginalized posterior and the corresponding prior \(\times\) likelihood expression are the following.

\[\begin{align*} p(\beta, \phi, \psi | Y, X_{mis}, X_{true}^{c}) &= \int p(X_{true}^{r}=x_{true}^{r}, \beta, \phi, \psi | Y, X_{mis}, X_{true}^{c}) d\mu(x_{true}^{r})\\ \\ &\propto p(\beta, \phi, \psi)\\ &~~~\times \prod^{m}_{i=1} \left[ p(Y_{i} | X_{true,i}^{r}=0, \beta) p(X_{mis,i} | X_{true,i}^{r}=0, \phi) (1-\psi) \right.\\ &~~~~~~~~~~~~~~\left. + p(Y_{i} | X_{true,i}^{r}=1, \beta) p(X_{mis,i} | X_{true,i}^{r}=1, \phi) (\psi) \right]\\ &~~~\times \prod^{n}_{j=m+1} p(Y_{j} | X_{true,j}^{c}, \beta) p(X_{mis,j} | X_{true,j}^{c}, \phi) p(X_{true,j}^{c} | \psi)\\ \end{align*} \]

The corresponding log expression is the following.

\[\begin{align*} \log p(\beta, \phi, \psi | Y, X_{mis}, X_{true}^{c}) &\propto \log p(\beta, \phi, \psi)\\ &~~~+ \sum^{m}_{i=1}\log \left[ p(Y_{i} | X_{true,i}^{r}=0, \beta) p(X_{mis,i} | X_{true,i}^{r}=0, \phi) (1-\psi) \right.\\ &~~~~~~~~~~~~~~~~~\left. + p(Y_{i} | X_{true,i}^{r}=1, \beta) p(X_{mis,i} | X_{true,i}^{r}=1, \phi) (\psi) \right]\\ &~~~+ \sum^{n}_{j=m+1} \left[ \log p(Y_{j} | X_{true,j}^{c}, \beta)+ \log p(X_{mis,j} | X_{true,j}^{c}, \phi)+ \log p(X_{true,j}^{c} | \psi) \right]\\ \end{align*} \]

Each term in the first summation can be written safely with the log_sum_exp() function in Stan.

\[\begin{align*} &\log \left[ p(Y_{i} | X_{true,i}^{r}=0, \beta) p(X_{mis,i} | X_{true,i}^{r}=0, \phi) (1-\psi) \right.\\ &~~~~~~\left. + p(Y_{i} | X_{true,i}^{r}=1, \beta) p(X_{mis,i} | X_{true,i}^{r}=1, \phi) (\psi) \right]\\ &= \text{log_sum_exp} \left( \right.\\ &~~~~~~\log p(Y_{i} | X_{true,i}^{r}=0, \beta) + \log p(X_{mis,i} | X_{true,i}^{r}=0, \phi) + \log (1-\psi),\\ &~~~~~~\log p(Y_{i} | X_{true,i}^{r}=1, \beta) + \log p(X_{mis,i} | X_{true,i}^{r}=1, \phi) + \log (\psi)\\ &~~~\left. \right)\\ \end{align*} \]

For an individual without the gold standard measurement, \(p(X_{true,i}^{r} = 1 | \psi) = \psi\) before data. Conditioning on the observed data (\(X_{mis,i}\) and \(Y_{i}\)) and parameters, the distribution of the unobserved \(X_{true,i}\) is the following.

\[\begin{align*} p(X_{true,i}^{r} = 1 | X_{mis,i}, Y_{i}, \beta, \phi, \psi) &= \frac{p(X_{true,i}^{r} = 1, X_{mis,i}, Y_{i}, \beta, \phi, \psi)}{\sum^{1}_{x=0} p(X_{true,i}^{r} = 1, X_{mis,i}, Y_{i}, \beta, \phi, \psi)}\\ &~~~~\text{Using $q$ for unnormalized density}\\ &= \frac{q(X_{true,i}^{r} = 1, X_{mis,i}, Y_{i}, \beta, \phi, \psi)}{\sum^{1}_{x=0} q(X_{true,i}^{r} = 1, X_{mis,i}, Y_{i}, \beta, \phi, \psi)}\\ \\ &~~~~\text{Also}\\ \log p(X_{true,i}^{r} = 1 | X_{mis,i}, Y_{i}, \beta, \phi, \psi) &= \log q(X_{true,i}^{r} = 1, X_{mis,i}, Y_{i}, \beta, \phi, \psi) - \log \left[ \sum^{1}_{x=0} q(X_{true,i}^{r} = 1, X_{mis,i}, Y_{i}, \beta, \phi, \psi) \right]\\ \end{align*} \]

Thus, we can reuse the log_sum_exp values for normalization to obtain \(p(X_{true,i}^{r} = 1 | X_{mis,i}, Y_{i}, \beta, \phi, \psi)\).

Stan implementation

## Load model
model1_code <- biostan::read_stan_file("./mismeasurement_bugs_book.stan")
biostan::print_stan_code(model1_code)
## data { 
##     /* N: Number of rows */ 
##     int<lower=0> N; 
##     int<lower=0,upper=1> Y[N]; 
##     int<lower=-1,upper=1> X_true[N]; 
##     int<lower=0,upper=1> X_mis[N]; 
##     int<lower=0> count[N]; 
##     int<lower=0,upper=1> R_X[N]; 
## } 
##  
## transformed data { 
##  
## } 
##  
## parameters { 
##     /* Outcome model parameters */ 
##     real beta0; 
##     real beta1; 
##     /* Error model parameter */ 
##     real<lower=0,upper=1> phi[2,2]; 
##     /* phi[1,1] = P(X_mis = 1 | X_true = 0, Y = 0) */ 
##     /* phi[1,2] = P(X_mis = 1 | X_true = 0, Y = 1) */ 
##     /* phi[2,1] = P(X_mis = 1 | X_true = 1, Y = 0) */ 
##     /* phi[2,2] = P(X_mis = 1 | X_true = 1, Y = 1) */ 
##     /* Covariate model parameter */ 
##     /* psi = P(X_true = 1) */ 
##     real<lower=0,upper=1> psi; 
## } 
##  
## transformed parameters { 
##     /* individual-level lp at X_true = 0 and X_true = 1 */ 
##     matrix[N,2] lp; 
##     /* https://mc-stan.org/docs/2_18/functions-reference/matrix-broadcast.html */ 
##     lp = rep_matrix(0, N, 2); 
##  
##     /* Loop over rows */ 
##     /* This is really a part of the model. */ 
##     /* Done here to carry over lp for generated quantities. */ 
##     /* p403 in Lambart 2018. */ 
##     for (i in 1:N) { 
##         if (R_X[i] == 1) { 
##             /* Assigned to column for observed X_true */ 
##             lp[i,X_true[i]+1] = 
##                 (/* Outcome model */ 
##                  bernoulli_lpmf(Y[i] | inv_logit(beta0 + beta1 * X_true[i])) 
##                  /* Error model */ 
##                  + bernoulli_lpmf(X_mis[i] | phi[X_true[i]+1, Y[i]+1]) 
##                  /* Covariate model */ 
##                  + bernoulli_lpmf(X_true[i] | psi)); 
##  
##         } else { 
##             /* Contribution for a row with UNOBSERVED X_true (marginalized over X_true) */ 
##  
##             /* X_true = 0 type contribution */ 
##             /* p(Yi|Xi=0,beta) p(Xi*|Xi=0,phi) p(Xi=0|psi) */ 
##             lp[i,1] = (/* Outcome model */ 
##                        bernoulli_lpmf(Y[i] | inv_logit(beta0)) 
##                        /* Error model */ 
##                        + bernoulli_lpmf(X_mis[i] | phi[1, Y[i]+1]) 
##                        /* Covariate model */ 
##                        + log(1-psi)); 
##  
##             /* X_true = 1 type contribution */ 
##             /* p(Yi|Xi=1,beta)p(Xi*|Xi=1,phi)p(Xi=1|psi) */ 
##             lp[i,2] = (/* Outcome model */ 
##                        bernoulli_lpmf(Y[i] | inv_logit(beta0 + beta1)) 
##                        /* Error model */ 
##                        + bernoulli_lpmf(X_mis[i] | phi[2, Y[i]+1]) 
##                        /* Covariate model */ 
##                        + log(psi)); 
##         } 
##     } 
## } 
##  
## model { 
##  
##     /* Priors */ 
##     /*  Outcome model parameters */ 
##     target += normal_lpdf(beta0 | 0, 100); 
##     target += normal_lpdf(beta1 | 0, 100); 
##     /*  Error model parameter */ 
##     target += uniform_lpdf(phi[1,1] | 0, 1); 
##     target += uniform_lpdf(phi[1,2] | 0, 1); 
##     target += uniform_lpdf(phi[2,1] | 0, 1); 
##     target += uniform_lpdf(phi[2,2] | 0, 1); 
##     /*  Covariate model parameter */ 
##     target += uniform_lpdf(psi | 0, 1); 
##  
##     /* Loop over rows */ 
##     for (i in 1:N) { 
##         if (R_X[i] == 1) { 
##             /* Contribution for a row with OBSERVED X */ 
##             /* count[i] to account for the sample size */ 
##             target += count[i] * lp[i,X_true[i]+1]; 
##  
##         } else { 
##             /* Sum up using log_sum_exp to marginalize over X_true */ 
##             /* count[i] to account for the sample size */ 
##             target += count[i] * log_sum_exp(lp[i]); 
##         } 
##     } 
##  
## } 
##  
## generated quantities { 
##     /* p(X_true = 1 | X_mis, Y, params) */ 
##     real p_X[N]; 
##  
##     for (i in 1:N) { 
##         if (R_X[i] == 1) { 
##             /* Observed X_true if available */ 
##             p_X[i] = X_true[i]; 
##         } else { 
##             p_X[i] = exp(lp[i,2] - log_sum_exp(lp[i])); 
##         } 
##     } 
## } 
## 
## Model fit
model1 <- stan(model_code = model1_code,
               data = c(N = nrow(cervical),
                        as.list(cervical)),
               cores = n_cores)
## Model check
pars <- c("beta0","beta1","phi","psi","lp__")
print(model1, pars = pars)
## Inference for Stan model: cdc83feb624d11817c1ebdec00a75c3b.
## 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% n_eff Rhat
## beta0       -0.91    0.00 0.20    -1.31    -1.03    -0.90    -0.78    -0.53  2491    1
## beta1        0.62    0.01 0.36    -0.08     0.38     0.61     0.86     1.31  2436    1
## phi[1,1]     0.32    0.00 0.05     0.21     0.28     0.32     0.35     0.42  2668    1
## phi[1,2]     0.22    0.00 0.08     0.08     0.16     0.21     0.27     0.39  3079    1
## phi[2,1]     0.57    0.00 0.06     0.45     0.53     0.57     0.61     0.69  2830    1
## phi[2,2]     0.76    0.00 0.06     0.64     0.72     0.76     0.81     0.89  2680    1
## psi          0.49    0.00 0.04     0.41     0.46     0.49     0.52     0.58  3406    1
## lp__     -2825.09    0.05 1.92 -2829.68 -2826.12 -2824.73 -2823.68 -2822.40  1601    1
## 
## Samples were drawn using NUTS(diag_e) at Thu Apr  4 04:47:11 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).
traceplot(model1, pars = pars)

pairs(model1, pars = pars)

The convergence looks ok. The results are consistent with the example fit in BUGS. The predicted probabilities of \(X_{true,i} = 1\) given observed data are following. Only the last four rows are meaningful.

## Data
cervical
## # A tibble: 12 x 5
##        Y X_true X_mis count   R_X
##    <dbl>  <dbl> <dbl> <dbl> <dbl>
##  1     1      0     0    13     1
##  2     1      0     1     3     1
##  3     1      1     0     5     1
##  4     1      1     1    18     1
##  5     0      0     0    33     1
##  6     0      0     1    11     1
##  7     0      1     0    16     1
##  8     0      1     1    16     1
##  9     1     -1     0   318     0
## 10     1     -1     1   375     0
## 11     0     -1     0   701     0
## 12     0     -1     1   535     0
## Prediction
print(model1, pars = "p_X")
## Inference for Stan model: cdc83feb624d11817c1ebdec00a75c3b.
## 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% n_eff Rhat
## p_X[1]  0.00     NaN 0.00 0.00 0.00 0.00 0.00  0.00   NaN  NaN
## p_X[2]  0.00     NaN 0.00 0.00 0.00 0.00 0.00  0.00   NaN  NaN
## p_X[3]  1.00     NaN 0.00 1.00 1.00 1.00 1.00  1.00   NaN  NaN
## p_X[4]  1.00     NaN 0.00 1.00 1.00 1.00 1.00  1.00   NaN  NaN
## p_X[5]  0.00     NaN 0.00 0.00 0.00 0.00 0.00  0.00   NaN  NaN
## p_X[6]  0.00     NaN 0.00 0.00 0.00 0.00 0.00  0.00   NaN  NaN
## p_X[7]  1.00     NaN 0.00 1.00 1.00 1.00 1.00  1.00   NaN  NaN
## p_X[8]  1.00     NaN 0.00 1.00 1.00 1.00 1.00  1.00   NaN  NaN
## p_X[9]  0.31       0 0.10 0.13 0.23 0.30 0.38  0.53  2583    1
## p_X[10] 0.83       0 0.08 0.66 0.79 0.84 0.89  0.95  3018    1
## p_X[11] 0.33       0 0.07 0.20 0.28 0.33 0.37  0.47  3160    1
## p_X[12] 0.58       0 0.09 0.40 0.52 0.59 0.65  0.75  2709    1
## 
## Samples were drawn using NUTS(diag_e) at Thu Apr  4 04:47:11 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).
print(sessionInfo())
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.4
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rstan_2.18.2       StanHeaders_2.18.1 biostan_0.1.0      ggdag_0.1.0        forcats_0.3.0     
##  [6] stringr_1.4.0      dplyr_0.8.0        purrr_0.3.0        readr_1.3.1        tidyr_0.8.2       
## [11] tibble_2.0.1       ggplot2_3.1.0      tidyverse_1.2.1    doRNG_1.7.1        rngtools_1.3.1    
## [16] pkgmaker_0.27      registry_0.5       doParallel_1.0.14  iterators_1.0.10   foreach_1.4.4     
## [21] knitr_1.21        
## 
## loaded via a namespace (and not attached):
##   [1] minqa_1.2.4        colorspace_1.4-0   ggridges_0.5.1     rsconnect_0.8.13   markdown_0.9      
##   [6] base64enc_0.1-3    rstudioapi_0.9.0   farver_1.1.0       ggrepel_0.8.0      DT_0.5            
##  [11] fansi_0.4.0        lubridate_1.7.4    xml2_1.2.0         codetools_0.2-16   splines_3.5.2     
##  [16] shinythemes_1.1.2  bayesplot_1.6.0    jsonlite_1.6       nloptr_1.2.1       broom_0.5.1       
##  [21] ggforce_0.1.3      shiny_1.2.0        compiler_3.5.2     httr_1.4.0         backports_1.1.3   
##  [26] assertthat_0.2.0   Matrix_1.2-15      lazyeval_0.2.1     cli_1.0.1          tweenr_1.0.1      
##  [31] later_0.8.0        htmltools_0.3.6    prettyunits_1.0.2  tools_3.5.2        igraph_1.2.4      
##  [36] gtable_0.2.0       glue_1.3.0         reshape2_1.4.3     V8_2.0             Rcpp_1.0.0        
##  [41] cellranger_1.1.0   nlme_3.1-137       crosstalk_1.0.0    ggraph_1.0.2       xfun_0.4          
##  [46] ps_1.3.0           lme4_1.1-20        rvest_0.3.2        mime_0.6           miniUI_0.1.1.1    
##  [51] gtools_3.8.1       MASS_7.3-51.1      zoo_1.8-4          scales_1.0.0       tidygraph_1.1.2   
##  [56] dagitty_0.2-2      rstanarm_2.18.2    colourpicker_1.0   hms_0.4.2          promises_1.0.1    
##  [61] inline_0.3.15      shinystan_2.5.0    curl_3.3           yaml_2.2.0         gridExtra_2.3     
##  [66] loo_2.0.0          stringi_1.3.1      dygraphs_1.1.1.6   boot_1.3-20        pkgbuild_1.0.2    
##  [71] bibtex_0.4.2       rlang_0.3.1        pkgconfig_2.0.2    matrixStats_0.54.0 evaluate_0.13     
##  [76] lattice_0.20-38    labeling_0.3       rstantools_1.5.1   htmlwidgets_1.3    processx_3.2.1    
##  [81] tidyselect_0.2.5   plyr_1.8.4         magrittr_1.5       R6_2.4.0           generics_0.0.2    
##  [86] pillar_1.3.1       haven_2.0.0        withr_2.1.2        units_0.6-2        xts_0.11-2        
##  [91] survival_2.43-3    modelr_0.1.3       crayon_1.3.4       KernSmooth_2.23-15 utf8_1.1.4        
##  [96] rmarkdown_1.11     viridis_0.5.1      grid_3.5.2         readxl_1.2.0       callr_3.1.1       
## [101] threejs_0.3.1      digest_0.6.18      xtable_1.8-3       httpuv_1.4.5.1     stats4_3.5.2      
## [106] munsell_0.5.0      viridisLite_0.3.0  shinyjs_1.0
## Record execution time and multicore use
end_time <- Sys.time()
diff_time <- difftime(end_time, start_time, units = "auto")
cat("Started  ", as.character(start_time), "\n",
    "Finished ", as.character(end_time), "\n",
    "Time difference of ", diff_time, " ", attr(diff_time, "units"), "\n",
    "Used ", foreach::getDoParWorkers(), " cores\n",
    "Used ", foreach::getDoParName(), " as backend\n",
    sep = "")
## Started  2019-04-04 04:46:27
## Finished 2019-04-04 04:47:28
## Time difference of 1.026295 mins
## Used 12 cores
## Used doParallelMC as backend