Example 1: single vs multi-chains

There is some overhead in parallel processing when using the future package, so the payoff when using parallel processing may vary by the problem. Here it is about a 2-4x speedup, but you can see more benefit at higher iterations. Note that this may not yield as many usable iterations as a single large chain if a substantial burnin period is needed, but it will enable useful convergence diagnostics. Note that the future package can implement sequential processing, which effectively turns the kmbayes_parallel into a loop, but still has all other advantages of multiple chains.

# enable parallel processing (up to 4 simultaneous processes here)
future::plan(strategy = future::multisession)

# single run of 4000 observations from bkmr package
set.seed(111)
system.time(kmfit <- suppressMessages(kmbayes(y = y, Z = Z, X = X, iter = 4000, verbose = FALSE, varsel = FALSE)))

##    user  system elapsed 
##  42.521   0.025  11.118

# 4 runs of 1000 observations from bkmrhat package
set.seed(111)
system.time(kmfit5 <- suppressMessages(kmbayes_parallel(nchains=4, y = y, Z = Z, X = X, iter = 1000, verbose = FALSE, varsel = FALSE)))

## Chain 1 
## Chain 2 
## Chain 3 
## Chain 4

##    user  system elapsed 
##  15.307   0.123 340.016

Example 2: Diagnostics

The diagnostics from the rstan package come from the monitor function (see the help files for that function in the rstan pacakge)

# Using rstan functions (set burnin/warmup to zero for comparability with coda numbers given later
#  posterior summaries should be performed after excluding warmup/burnin)
singlediag = kmbayes_diagnose(kmfit, warmup=0, digits_summary=2)

## Single chain

## Inference for the input samples (1 chains: each with iter = 4000; warmup = 0):
## 
##            Q5  Q50  Q95 Mean  SD  Rhat Bulk_ESS Tail_ESS
## beta1     1.9  2.0  2.0  2.0 0.0  1.00     2665     2969
## beta2     0.0  0.1  0.3  0.1 0.1  1.00     3041     3693
## lambda    4.3 11.1 27.6 13.1 8.2  1.01      259      196
## r1        0.0  0.0  0.1  0.0 0.1  1.00      170      162
## r2        0.0  0.0  0.1  0.0 0.1  1.00      310      244
## r3        0.0  0.0  0.0  0.0 0.0  1.00      155      128
## r4        0.0  0.0  0.1  0.0 0.1  1.00      170      155
## r5        0.0  0.0  0.0  0.0 0.1  1.01       90      131
## sigsq.eps 0.2  0.3  0.5  0.4 0.1  1.00     1326     1631
## 
## For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
## effective sample size for bulk and tail quantities respectively (an ESS > 100 
## per chain is considered good), and Rhat is the potential scale reduction 
## factor on rank normalized split chains (at convergence, Rhat <= 1.05).

# Using rstan functions (multiple chains enable R-hat)
multidiag = kmbayes_diagnose(kmfit5, warmup=0, digits_summary=2)

## Parallel chains

## Inference for the input samples (4 chains: each with iter = 1000; warmup = 0):
## 
##            Q5  Q50  Q95 Mean  SD  Rhat Bulk_ESS Tail_ESS
## beta1     1.9  2.0  2.0  2.0 0.0  1.00     2565     1385
## beta2     0.0  0.1  0.3  0.1 0.1  1.00     2789     2854
## lambda    4.3 10.7 25.0 12.3 6.6  1.01      323      265
## r1        0.0  0.0  0.1  0.1 0.2  1.03       99       61
## r2        0.0  0.0  0.1  0.1 0.1  1.02      195      232
## r3        0.0  0.0  0.1  0.0 0.1  1.04       91       82
## r4        0.0  0.0  0.1  0.0 0.1  1.03      116       96
## r5        0.0  0.0  0.2  0.1 0.3  1.05       67       22
## sigsq.eps 0.2  0.3  0.5  0.4 0.1  1.01     1062      883
## 
## For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
## effective sample size for bulk and tail quantities respectively (an ESS > 100 
## per chain is considered good), and Rhat is the potential scale reduction 
## factor on rank normalized split chains (at convergence, Rhat <= 1.05).

# using coda functions, not using any burnin (for demonstration only)
kmfitcoda = as.mcmc(kmfit, iterstart = 1)
kmfit5coda = as.mcmc.list(kmfit5, iterstart = 1)

# single chain trace plot
traceplot(kmfitcoda)

The trace plots look typical, and fine, but trace plots don’t give a full picture of convergence. Note that there is apparent quick convergence for a couple of parameters demonstrated by movement away from the starting value and concentration of the rest of the samples within a narrow band.

Seeing visual evidence that different chains are sampling from the same marginal distributions is reassuring about the stability of the results.

# multiple chain trace plot
traceplot(kmfit5coda)

Now examine “cross correlation”, which can help identify highly correlated parameters in the posterior, which can be problematic for MCMC sampling. Here there is a block {r3,r4,r5} which appear to be highly correlated. All other things equal, having highly correlated parameters in the posterior means that more samples are needed than would be needed with uncorrelated parameters.

# multiple cross-correlation plot (combines all samples)
crosscorr(kmfit5coda)

##                 beta1       beta2      lambda         r1          r2
## beta1      1.00000000  0.04252082 -0.01597349  0.2003976  0.16893771
## beta2      0.04252082  1.00000000 -0.06910184 -0.1130723 -0.11967535
## lambda    -0.01597349 -0.06910184  1.00000000 -0.0609043 -0.06292417
## r1         0.20039759 -0.11307231 -0.06090430  1.0000000  0.69761056
## r2         0.16893771 -0.11967535 -0.06292417  0.6976106  1.00000000
## r3         0.15362529 -0.09215259 -0.02952802  0.7775382  0.82344238
## r4         0.22913128 -0.10932048 -0.04174805  0.9045248  0.78294291
## r5         0.09391064 -0.05213624  0.01700295  0.2544064  0.30114999
## sigsq.eps -0.06070043  0.07624883 -0.34190260 -0.1174797 -0.14333533
##                    r3          r4          r5   sigsq.eps
## beta1      0.15362529  0.22913128  0.09391064 -0.06070043
## beta2     -0.09215259 -0.10932048 -0.05213624  0.07624883
## lambda    -0.02952802 -0.04174805  0.01700295 -0.34190260
## r1         0.77753815  0.90452482  0.25440642 -0.11747969
## r2         0.82344238  0.78294291  0.30114999 -0.14333533
## r3         1.00000000  0.75420334  0.29350443 -0.12028459
## r4         0.75420334  1.00000000  0.33543998 -0.12125084
## r5         0.29350443  0.33543998  1.00000000 -0.03852572
## sigsq.eps -0.12028459 -0.12125084 -0.03852572  1.00000000

crosscorr.plot(kmfit5coda)

Now examine “autocorrelation” to identify parameters that have high correlation between subsequent iterations of the MCMC sampler, which can lead to inefficient MCMC sampling. All other things equal, having highly autocorrelated parameters in the posterior means that more samples are needed than would be needed with low-autocorrelation parameters.

# multiple chain trace plot
#autocorr(kmfit5coda) # lots of output
autocorr.plot(kmfit5coda)

Graphical tools can be limited, and are sometimes difficult to use effectively with scale parameters (of which bkmr has many). Additionally, no single diagnostic is perfect, leading many authors to advocate the use of multiple, complementary diagnostics. Thus, more formal diagnostics are helpful.

Gelman’s r-hat diagnostic gives an interpretable diagnostic: the expected reduction in the standard error of the posterior means if you could run the chains to an infinite size. These give some idea about when is a fine idea to stop sampling. There are rules of thumb about using r-hat to stop sampling that are available from several authors (for example you can consult the help files for rstan and coda).

Effective sample size is also useful - it estimates the amount of information in your chain, expressed in terms of the number of independent posterior samples it would take to match that information (e.g. if we could just sample from the posterior directly).

# Gelman's r-hat using coda estimator (will differ from rstan implementation)
gelman.diag(kmfit5coda)

## Potential scale reduction factors:
## 
##           Point est. Upper C.I.
## beta1           1.00       1.00
## beta2           1.00       1.01
## lambda          1.00       1.00
## r1              1.05       1.09
## r2              1.12       1.30
## r3              1.06       1.16
## r4              1.05       1.07
## r5              1.03       1.03
## sigsq.eps       1.01       1.02
## 
## Multivariate psrf
## 
## 1.07

# effective sample size
effectiveSize(kmfitcoda)

##     beta1     beta2    lambda        r1        r2        r3        r4        r5 
## 2892.7570 3198.3261  202.6782  118.1804  284.6426  295.7300  146.9502  115.5348 
## sigsq.eps 
## 1242.4771

effectiveSize(kmfit5coda)

##     beta1     beta2    lambda        r1        r2        r3        r4        r5 
## 2592.0433 3341.5839  424.5452  109.5655  175.1005  115.6253  159.1675  181.0398 
## sigsq.eps 
## 1333.3042

Example 3: Posterior summaries

Posterior kernel marginal densities, 1 chain

# posterior kernel marginal densities using `mcmc` and `mcmc` objects
densplot(kmfitcoda)

Posterior kernel marginal densities, multiple chains combined. Look for multiple modes that may indicate non-convergence of some chains

# posterior kernel marginal densities using `mcmc` and `mcmc` objects
densplot(kmfit5coda)

Other diagnostics from the coda package are available here.

Finally, the chains from the original kmbayes_parallel fit can be combined into a single chain (see the help files for how to deal with burn-in, the default in bkmr is to use the first half of the chain, which is respected here). The kmbayes_combine function smartly first combines the burn-in iterations and then combines the iterations after burnin, such that the burn-in rules of subsequent functions within the bkmr package are respected. Note that unlike the as.mcmc.list function, this function combines all iterations into a single chain, so trace plots will not be good diagnotistics in this combined object, and it should be used once one is assured that all chains have converged and the burn-in is acceptable.

With this combined set of samples, you can follow any of the post-processing functions from the bkmr functions, which are described here: https://jenfb.github.io/bkmr/overview.html. For example, see below the estimation of the posterior mean difference along a series of quantiles of all exposures in Z.

# posterior summaries using `mcmc` and `mcmc` objects
summary(kmfitcoda)

## 
## Iterations = 1:4000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 4000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##               Mean      SD  Naive SE Time-series SE
## beta1      1.97357 0.04598 0.0007270      0.0008549
## beta2      0.13518 0.08710 0.0013771      0.0015401
## lambda    13.08699 8.23403 0.1301915      0.5783746
## r1         0.02829 0.06333 0.0010013      0.0058256
## r2         0.03552 0.05200 0.0008222      0.0030823
## r3         0.02140 0.04096 0.0006476      0.0023818
## r4         0.02931 0.06623 0.0010472      0.0054636
## r5         0.03055 0.10618 0.0016788      0.0098783
## sigsq.eps  0.35970 0.08495 0.0013432      0.0024101
## 
## 2. Quantiles for each variable:
## 
##               2.5%     25%      50%      75%    97.5%
## beta1      1.88651 1.94226  1.97310  2.00380  2.06407
## beta2     -0.03357 0.07694  0.13345  0.19279  0.30592
## lambda     3.95746 7.68835 11.09318 16.25233 33.16143
## r1         0.01031 0.01229  0.01676  0.02568  0.08740
## r2         0.01056 0.01474  0.02290  0.04175  0.10979
## r3         0.01025 0.01194  0.01437  0.02179  0.06415
## r4         0.01021 0.01286  0.01763  0.02860  0.07571
## r5         0.01014 0.01198  0.01451  0.02031  0.06762
## sigsq.eps  0.23236 0.29752  0.34779  0.40793  0.55917

summary(kmfit5coda)

## 
## Iterations = 1:1000
## Thinning interval = 1 
## Number of chains = 4 
## Sample size per chain = 1000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##               Mean      SD  Naive SE Time-series SE
## beta1      1.97157 0.04777 0.0007553       0.001162
## beta2      0.13206 0.09025 0.0014270       0.001600
## lambda    12.27249 6.63455 0.1049014       0.347197
## r1         0.05839 0.17650 0.0027907       0.022236
## r2         0.05182 0.10180 0.0016095       0.007928
## r3         0.04333 0.13064 0.0020657       0.012293
## r4         0.04463 0.13129 0.0020759       0.012878
## r5         0.07322 0.25238 0.0039905       0.048491
## sigsq.eps  0.35697 0.08628 0.0013642       0.002368
## 
## 2. Quantiles for each variable:
## 
##               2.5%     25%      50%      75%   97.5%
## beta1      1.88082 1.93959  1.97017  2.00159  2.0664
## beta2     -0.04983 0.07316  0.13350  0.18970  0.3091
## lambda     3.71673 7.50762 10.71569 15.42226 28.9153
## r1         0.01025 0.01238  0.01690  0.02882  0.7633
## r2         0.01060 0.01657  0.02839  0.05083  0.2150
## r3         0.01035 0.01245  0.01658  0.02387  0.4567
## r4         0.01021 0.01224  0.01563  0.02416  0.5177
## r5         0.01021 0.01229  0.01585  0.02261  1.1716
## sigsq.eps  0.21938 0.29655  0.34822  0.40572  0.5591

# highest posterior density intervals using `mcmc` and `mcmc` objects
HPDinterval(kmfitcoda)

##                 lower       upper
## beta1      1.88272469  2.05913505
## beta2     -0.03379651  0.30462773
## lambda     2.62734881 27.66498240
## r1         0.01003506  0.05828581
## r2         0.01000265  0.08925073
## r3         0.01009384  0.04319179
## r4         0.01001438  0.06122670
## r5         0.01003380  0.03968790
## sigsq.eps  0.21094543  0.52732628
## attr(,"Probability")
## [1] 0.95

HPDinterval(kmfit5coda)

## [[1]]
##                 lower       upper
## beta1      1.86951656  2.06341096
## beta2     -0.03839616  0.30647850
## lambda     3.39259225 28.34283839
## r1         0.01001457  0.16930580
## r2         0.01028361  0.16154920
## r3         0.01045287  0.06977346
## r4         0.01014374  0.13509464
## r5         0.01020661  0.06500149
## sigsq.eps  0.21505633  0.53032414
## attr(,"Probability")
## [1] 0.95
## 
## [[2]]
##                 lower       upper
## beta1      1.89169671  2.06839274
## beta2     -0.05307858  0.31160505
## lambda     3.25617836 22.66751016
## r1         0.01053872  0.76527137
## r2         0.01008093  0.10932136
## r3         0.01003452  0.07355298
## r4         0.01008777  0.16628667
## r5         0.01035072  0.05746120
## sigsq.eps  0.20985426  0.51534072
## attr(,"Probability")
## [1] 0.95
## 
## [[3]]
##                 lower       upper
## beta1      1.87657975  2.05685234
## beta2     -0.03508781  0.31833285
## lambda     3.04977671 23.93857170
## r1         0.01003505  0.08707388
## r2         0.01022525  0.11415988
## r3         0.01017528  0.14258896
## r4         0.01020691  0.07132056
## r5         0.01007435  0.05801491
## sigsq.eps  0.18913969  0.50397884
## attr(,"Probability")
## [1] 0.95
## 
## [[4]]
##                 lower       upper
## beta1      1.87858665  2.06076279
## beta2     -0.04643783  0.31027316
## lambda     3.08683817 27.73022418
## r1         0.01002781  0.08317944
## r2         0.01018996  0.08783009
## r3         0.01040574  0.08914759
## r4         0.01004444  0.07922687
## r5         0.01000356  1.47345520
## sigsq.eps  0.19615666  0.55578747
## attr(,"Probability")
## [1] 0.95

# combine multiple chains into a single chain
fitkmccomb = kmbayes_combine(kmfit5)


# For example:
summary(fitkmccomb)

## Fitted object of class 'bkmrfit'
## Iterations: 4000 
## Outcome family: gaussian  
## Model fit on: 2026-10-20 08:04:31.235701 
## Running time:  5.62377 mins 
## 
## Acceptance rates for Metropolis-Hastings algorithm:
##    param      rate
## 1 lambda 0.4361090
## 2     r1 0.1827957
## 3     r2 0.3060765
## 4     r3 0.1587897
## 5     r4 0.1737934
## 6     r5 0.1592898
## 
## Parameter estimates (based on iterations 2001-4000):
##       param     mean      sd    q_2.5   q_97.5
## 1     beta1  1.96893 0.04537  1.87980  2.05907
## 2     beta2  0.13814 0.08995 -0.03586  0.31167
## 3 sigsq.eps  0.36261 0.08354  0.23126  0.56016
## 4        r1  0.02125 0.01529  0.01025  0.06188
## 5        r2  0.03967 0.03487  0.01062  0.13518
## 6        r3  0.01930 0.01231  0.01022  0.05548
## 7        r4  0.02001 0.01601  0.01009  0.06208
## 8        r5  0.01870 0.01261  0.01021  0.05746
## 9    lambda 11.67207 6.01689  3.61709 26.68135

## NULL

mean.difference <- suppressWarnings(OverallRiskSummaries(fit = fitkmccomb, y = y, Z = Z, X = X, 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact"))
mean.difference

##    quantile         est         sd
## 1      0.25 -0.43174980 0.09054293
## 2      0.30 -0.37574135 0.06685229
## 3      0.35 -0.17904107 0.04131393
## 4      0.40 -0.14035080 0.03453813
## 5      0.45 -0.06161091 0.03071210
## 6      0.50  0.00000000 0.00000000
## 7      0.55  0.16914871 0.05057507
## 8      0.60  0.31162209 0.07039285
## 9      0.65  0.57650940 0.11070336
## 10     0.70  0.70282460 0.13214568
## 11     0.75  0.84098199 0.17352702

with(mean.difference, {
  plot(quantile, est, pch=19, ylim=c(min(est - 1.96*sd), max(est + 1.96*sd)), 
       axes=FALSE, ylab= "Mean difference", xlab = "Joint quantile")
  segments(x0=quantile, x1=quantile, y0 = est - 1.96*sd, y1 = est + 1.96*sd)
  abline(h=0)
  axis(1)
  axis(2)
  box(bty='l')
})

Example 5: Parallel posterior summaries as diagnostics

bkmrhat also has ported versions of the native posterior summarization functions to compare how these summaries vary across parallel chains. Note that these should serve as diagnostics, and final posterior inference should be done on the combined chain. The easiest of these functions to demonstrate is the OverallRiskSummaries_parallel function, which simply runs OverallRiskSummaries (from the bkmr package) on each chain and combines the results. Notably, this function fixes the y-axis at zero for the median, so it under-represents overall predictive variation across chains, but captures variation in effect estimates across the chains. Ideally, that variation is negligible - e.g. if you see differences between chains that would result in different interpretations, you should re-fit the model with more iterations. In this example, the results are reasonably consistent across chains, but one might want to run more iterations if, say, the differences seen across the upper error bounds are of such a magnitude as to be practically meaningful.

set.seed(111)
system.time(kmfitbma.list <- suppressWarnings(kmbayes_parallel(nchains=4, y = y, Z = Z, X = X, iter = 1000, verbose = FALSE, varsel = TRUE)))

## Chain 1

## Iteration: 100 (10% completed; 10.5647 secs elapsed)

## Iteration: 200 (20% completed; 21.55969 secs elapsed)

## Iteration: 300 (30% completed; 33.6677 secs elapsed)

## Iteration: 400 (40% completed; 44.91266 secs elapsed)

## Iteration: 500 (50% completed; 54.43071 secs elapsed)

## Iteration: 600 (60% completed; 1.06014 mins elapsed)

## Iteration: 700 (70% completed; 1.23671 mins elapsed)

## Iteration: 800 (80% completed; 1.43028 mins elapsed)

## Iteration: 900 (90% completed; 1.63798 mins elapsed)

## Iteration: 1000 (100% completed; 1.81484 mins elapsed)

## Chain 2

## Iteration: 100 (10% completed; 12.52805 secs elapsed)

## Iteration: 200 (20% completed; 24.95635 secs elapsed)

## Iteration: 300 (30% completed; 38.06532 secs elapsed)

## Iteration: 400 (40% completed; 47.85933 secs elapsed)

## Iteration: 500 (50% completed; 56.17085 secs elapsed)

## Iteration: 600 (60% completed; 1.1039 mins elapsed)

## Iteration: 700 (70% completed; 1.30731 mins elapsed)

## Iteration: 800 (80% completed; 1.49777 mins elapsed)

## Iteration: 900 (90% completed; 1.69636 mins elapsed)

## Iteration: 1000 (100% completed; 1.84151 mins elapsed)

## Chain 3

## Iteration: 100 (10% completed; 10.50215 secs elapsed)

## Iteration: 200 (20% completed; 22.61535 secs elapsed)

## Iteration: 300 (30% completed; 35.22633 secs elapsed)

## Iteration: 400 (40% completed; 45.7622 secs elapsed)

## Iteration: 500 (50% completed; 55.8743 secs elapsed)

## Iteration: 600 (60% completed; 1.12271 mins elapsed)

## Iteration: 700 (70% completed; 1.30781 mins elapsed)

## Iteration: 800 (80% completed; 1.50354 mins elapsed)

## Iteration: 900 (90% completed; 1.70104 mins elapsed)

## Iteration: 1000 (100% completed; 1.84019 mins elapsed)

## Chain 4

## Iteration: 100 (10% completed; 10.67531 secs elapsed)

## Iteration: 200 (20% completed; 20.76254 secs elapsed)

## Iteration: 300 (30% completed; 32.06135 secs elapsed)

## Iteration: 400 (40% completed; 44.02886 secs elapsed)

## Iteration: 500 (50% completed; 50.81131 secs elapsed)

## Iteration: 600 (60% completed; 1.0452 mins elapsed)

## Iteration: 700 (70% completed; 1.24567 mins elapsed)

## Iteration: 800 (80% completed; 1.43899 mins elapsed)

## Iteration: 900 (90% completed; 1.64897 mins elapsed)

## Iteration: 1000 (100% completed; 1.80384 mins elapsed)

##    user  system elapsed 
##   4.476   0.031 110.623

meandifference_par = OverallRiskSummaries_parallel(kmfitbma.list, y = y, Z = Z, X = X ,qs = seq(0.25, 0.75, by = 0.05), q.fixed = 0.5, method = "exact")

## Chain 1 
## Chain 2 
## Chain 3 
## Chain 4

head(meandifference_par)

##   quantile        est         sd chain
## 1     0.25 -0.4237405 0.09912107     1
## 2     0.30 -0.3601165 0.07425543     1
## 3     0.35 -0.1734402 0.04549098     1
## 4     0.40 -0.1390188 0.03842276     1
## 5     0.45 -0.0538018 0.03317052     1
## 6     0.50  0.0000000 0.00000000     1

nchains = length(unique(meandifference_par$chain))

with(meandifference_par, {
  plot.new()
  plot.window(ylim=c(min(est - 1.96*sd), max(est + 1.96*sd)), 
              xlim=c(min(quantile), max(quantile)),
       ylab= "Mean difference", xlab = "Joint quantile")
  for(cch in seq_len(nchains)){
    width = diff(quantile)[1]
    jit = runif(1, -width/5, width/5)
   points(jit+quantile[chain==cch], est[chain==cch], pch=19, col=cch) 
   segments(x0=jit+quantile[chain==cch], x1=jit+quantile[chain==cch], y0 = est[chain==cch] - 1.96*sd[chain==cch], y1 = est[chain==cch] + 1.96*sd[chain==cch], col=cch)
  }
  abline(h=0)
  axis(1)
  axis(2)
  box(bty='l')
  legend("bottom", col=1:nchains, pch=19, lty=1, legend=paste("chain", 1:nchains), bty="n")
})

regfuns_par = PredictorResponseUnivar_parallel(kmfitbma.list, y = y, Z = Z, X = X ,qs = seq(0.25, 0.75, by = 0.05), q.fixed = 0.5, method = "exact")

## Chain 1 
## Chain 2 
## Chain 3 
## Chain 4

head(regfuns_par)

##   variable         z        est        se chain
## 1       z1 -2.159373 -0.7372148 0.6261796     1
## 2       z1 -2.048986 -0.7033279 0.5981104     1
## 3       z1 -1.938600 -0.6691741 0.5700151     1
## 4       z1 -1.828214 -0.6347894 0.5419344     1
## 5       z1 -1.717827 -0.6002116 0.5139157     1
## 6       z1 -1.607441 -0.5654801 0.4860135     1

nchains = length(unique(meandifference_par$chain))

# single variable
with(regfuns_par[regfuns_par$variable=="z1",], {
  plot.new()
  plot.window(ylim=c(min(est - 1.96*se), max(est + 1.96*se)), 
              xlim=c(min(z), max(z)),
       ylab= "Predicted Y", xlab = "Z")
  pc = c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#999999")
  pc2 = c("#0000001A", "#E69F001A", "#56B4E91A", "#009E731A", "#F0E4421A", "#0072B21A", "#D55E001A", "#CC79A71A", "#9999991A")
  for(cch in seq_len(nchains)){
   ribbonX = c(z[chain==cch], rev(z[chain==cch]))
   ribbonY = c(est[chain==cch] + 1.96*se[chain==cch], rev(est[chain==cch] - 1.96*se[chain==cch]))
   polygon(x=ribbonX, y = ribbonY, col=pc2[cch], border=NA)
   lines(z[chain==cch], est[chain==cch], pch=19, col=pc[cch]) 
  }
  axis(1)
  axis(2)
  box(bty='l')
  legend("bottom", col=1:nchains, pch=19, lty=1, legend=paste("chain", 1:nchains), bty="n")
})

Example 6: Continuing a fit

Sometimes you just need to run more samples in an existing chain. For example, you run a bkmr fit for 3 days, only to find you don’t have enough samples. A “continued” fit just means that you can start off at the last iteration you were at and just keep building on an existing set of results by lengthening the Markov chain. Unfortunately, due to how the kmbayes function accepts starting values (for the official install version), you can’t quite do this exactly in many cases (The function will relay a message and possible solutions, if any. bkmr package authors are aware of this issue). The kmbayes_continue function continues a bkmr fit as well as the bkmr package will allow. The r parameters from the fit must all be initialized at the same value, so kmbayes_continue starts a new MCMC fit at the final values of all parameters from the prior bkmr fit, but sets all of the r parameters to the mean at the last iteration from the previous fit. Additionally, if h.hat parameters are estimated, these are fixed to be above zero to meet similar constraints, either by fixing them at their posterior mean or setting to a small positive value. One should inspect trace plots to see whether this will cause issues (e.g. if the traceplots demonstrate different patterns in the samples before and after the continuation). Here’s an example with a quick check of diagnostics of the first part of the chain, and the combined chain (which could be used for inference or extended again, if necessary). We caution users that this function creates 2 distinct, if very similar Markov chains, and to use appropriate caution if traceplots differ before and after each continuation. Nonetheless, in many cases one can act as though all samples are from the same Markov chain.

Note that if you install the developmental version of the bkmr package you can continue fits from exactly where they left off, so you get a true, single Markov chain. You can install that via the commented code below

# install dev version of bkmr to allow true continued fits.
#install.packages("devtools")
#devtools::install_github("jenfb/bkmr")

set.seed(111)
# run 100 initial iterations for a model with only 2 exposures
Z2 = Z[,1:2]
kmfitbma.start <- suppressWarnings(kmbayes(y = y, Z = Z2, X = X, iter = 500, verbose = FALSE, varsel = FALSE))

## Iteration: 50 (10% completed; 0.06584 secs elapsed)

## Iteration: 100 (20% completed; 0.12186 secs elapsed)

## Iteration: 150 (30% completed; 0.18961 secs elapsed)

## Iteration: 200 (40% completed; 0.24566 secs elapsed)

## Iteration: 250 (50% completed; 0.31253 secs elapsed)

## Iteration: 300 (60% completed; 0.3785 secs elapsed)

## Iteration: 350 (70% completed; 0.43527 secs elapsed)

## Iteration: 400 (80% completed; 0.50197 secs elapsed)

## Iteration: 450 (90% completed; 0.5579 secs elapsed)

## Iteration: 500 (100% completed; 0.6241 secs elapsed)

kmbayes_diag(kmfitbma.start)

## Single chain

## Inference for the input samples (1 chains: each with iter = 500; warmup = 250):
## 
##             Q5  Q50  Q95 Mean  SD  Rhat Bulk_ESS Tail_ESS
## beta1      1.9  1.9  2.0  1.9 0.0  1.00      279      217
## beta2     -0.1  0.1  0.3  0.1 0.1  1.00      238      206
## lambda     4.6 13.0 31.3 14.9 7.8  1.08       13       16
## r1         0.0  0.0  0.0  0.0 0.0  1.17        5        7
## r2         0.0  0.0  0.1  0.1 0.0  1.05       10       16
## sigsq.eps  0.3  0.4  0.6  0.4 0.1  1.00       87      219
## 
## For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
## effective sample size for bulk and tail quantities respectively (an ESS > 100 
## per chain is considered good), and Rhat is the potential scale reduction 
## factor on rank normalized split chains (at convergence, Rhat <= 1.05).

##                  mean     se_mean         sd        2.5%        25%         50%
## beta1      1.94302414 0.002744458 0.04643461  1.84782583 1.91330142  1.94655854
## beta2      0.10483405 0.006191034 0.09620251 -0.07303313 0.04304370  0.10467199
## lambda    14.90060019 1.869438290 7.81070475  4.60230691 8.65519030 13.04053754
## r1         0.01476940 0.002347690 0.00449794  0.01078574 0.01235942  0.01320025
## r2         0.05291092 0.009623519 0.04332210  0.01010751 0.01924294  0.03815578
## sigsq.eps  0.41328122 0.010095801 0.09214694  0.27259916 0.35051490  0.39841682
##                   75%       97.5% n_eff      Rhat valid          Q5         Q50
## beta1      1.97193285  2.02693739   282 1.0049308     1  1.86400517  1.94655854
## beta2      0.16804699  0.27700828   236 0.9977248     1 -0.05155040  0.10467199
## lambda    18.65686378 34.19469565    18 1.0843443     1  4.60230691 13.04053754
## r1         0.01703200  0.02877183     9 1.1678074     1  0.01078574  0.01320025
## r2         0.07782864  0.15639005    21 1.0523022     1  0.01010751  0.03815578
## sigsq.eps  0.46665521  0.62932750    94 1.0023168     1  0.28064363  0.39841682
##                   Q95    MCSE_Q2.5     MCSE_Q25    MCSE_Q50    MCSE_Q75
## beta1      2.01504706 0.0098237593 0.0048432127 0.003083228 0.004034983
## beta2      0.26110337 0.0166696373 0.0093358891 0.007387626 0.010829577
## lambda    31.32637623 1.1971720562 2.0194134613 2.532679950 1.535664875
## r1         0.02153674 0.0002641772 0.0009430751 0.001429907 0.004168245
## r2         0.14102786 0.0028040064 0.0071858659 0.011789241 0.013548261
## sigsq.eps  0.60258149 0.0068027724 0.0052486067 0.005123756 0.008923788
##            MCSE_Q97.5     MCSE_SD Bulk_ESS Tail_ESS
## beta1     0.006818301 0.001946731      279      217
## beta2     0.008917280 0.005057796      238      206
## lambda    1.781705411 1.344630129       13       16
## r1        0.004563786 0.001813399        5        7
## r2        0.010363171 0.006905270       10       16
## sigsq.eps 0.013662390 0.007248003       87      219

# run 2000 additional iterations
moreiterations = kmbayes_continue(kmfitbma.start, iter=2000)

## Validating control.params...

## Validating starting.values...

## Iteration: 201 (10% completed; 0.28091 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param  rate
## 1 lambda 0.425
## 2     r1 0.140
## 3     r2 0.475

## Iteration: 401 (20% completed; 0.55099 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param   rate
## 1 lambda 0.4500
## 2     r1 0.1475
## 3     r2 0.3975

## Iteration: 601 (30% completed; 0.84436 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param      rate
## 1 lambda 0.4750000
## 2     r1 0.1400000
## 3     r2 0.4083333

## Iteration: 801 (40% completed; 1.12426 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param    rate
## 1 lambda 0.48250
## 2     r1 0.13875
## 3     r2 0.40125

## Iteration: 1001 (50% completed; 1.39722 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param  rate
## 1 lambda 0.457
## 2     r1 0.137
## 3     r2 0.412

## Iteration: 1201 (60% completed; 1.68068 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param      rate
## 1 lambda 0.4683333
## 2     r1 0.1483333
## 3     r2 0.3975000

## Iteration: 1401 (70% completed; 1.96359 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param      rate
## 1 lambda 0.4885714
## 2     r1 0.1500000
## 3     r2 0.4028571

## Iteration: 1601 (80% completed; 2.27651 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param     rate
## 1 lambda 0.484375
## 2     r1 0.155625
## 3     r2 0.408750

## Iteration: 1801 (90% completed; 2.55114 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param      rate
## 1 lambda 0.4911111
## 2     r1 0.1627778
## 3     r2 0.4061111

## Iteration: 2001 (100% completed; 2.83461 secs elapsed)

## Acceptance rates for Metropolis-Hastings algorithm:
##    param   rate
## 1 lambda 0.4935
## 2     r1 0.1590
## 3     r2 0.4055

kmbayes_diag(moreiterations)

## Single chain

## Inference for the input samples (1 chains: each with iter = 2500; warmup = 1250):
## 
##            Q5  Q50  Q95 Mean  SD  Rhat Bulk_ESS Tail_ESS
## beta1     1.9  1.9  2.0  1.9 0.0  1.00     1195     1239
## beta2     0.0  0.1  0.3  0.1 0.1  1.00      861     1219
## lambda    4.9 12.1 30.5 14.3 8.5  1.00      126       93
## r1        0.0  0.0  0.1  0.0 0.0  1.03       77       78
## r2        0.0  0.0  0.2  0.1 0.0  1.00      112      100
## sigsq.eps 0.3  0.4  0.6  0.4 0.1  1.00      514      890
## 
## For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
## effective sample size for bulk and tail quantities respectively (an ESS > 100 
## per chain is considered good), and Rhat is the potential scale reduction 
## factor on rank normalized split chains (at convergence, Rhat <= 1.05).

##                  mean     se_mean         sd        2.5%        25%         50%
## beta1      1.94451464 0.001354599 0.04679152  1.84957428 1.91337229  1.94425420
## beta2      0.10591205 0.003098008 0.09190734 -0.07244492 0.04307369  0.10417774
## lambda    14.34371537 0.741761549 8.52169987  4.05070504 7.75788463 12.05628895
## r1         0.02298157 0.001565890 0.01454607  0.01069964 0.01263807  0.01684080
## r2         0.05795423 0.004171398 0.04711284  0.01155283 0.02140563  0.04287896
## sigsq.eps  0.40958211 0.004090066 0.09282983  0.26471559 0.34187272  0.39924364
##                   75%       97.5% n_eff      Rhat valid          Q5         Q50
## beta1      1.97704630  2.03445739  1186 0.9998051     1  1.86691892  1.94425420
## beta2      0.16721410  0.29665914   874 0.9996921     1 -0.04406332  0.10417774
## lambda    18.43684053 35.00798030   130 1.0031465     1  4.88809530 12.05628895
## r1         0.02691796  0.06024277    87 1.0285881     1  0.01105215  0.01684080
## r2         0.08141999  0.18286772   128 1.0017863     1  0.01242185  0.04287896
## sigsq.eps  0.46162498  0.62034160   515 0.9992725     1  0.28188526  0.39924364
##                   Q95    MCSE_Q2.5     MCSE_Q25     MCSE_Q50    MCSE_Q75
## beta1      2.02154219 0.0053480030 0.0018254915 0.0017327069 0.001809085
## beta2      0.25197027 0.0060693244 0.0032074359 0.0037835712 0.004042631
## lambda    30.54061347 0.6391431563 0.5063700931 0.7650849805 1.072408138
## r1         0.05704602 0.0002311747 0.0003811785 0.0007469754 0.002080573
## r2         0.16096251 0.0004350924 0.0020817839 0.0040577461 0.005804157
## sigsq.eps  0.58415925 0.0046659948 0.0037588012 0.0037328009 0.005981460
##            MCSE_Q97.5      MCSE_SD Bulk_ESS Tail_ESS
## beta1     0.003404986 0.0009580805     1195     1239
## beta2     0.010762886 0.0021913487      861     1219
## lambda    2.300751507 0.5314745343      126       93
## r1        0.001167341 0.0011110153       77       78
## r2        0.005557337 0.0029563949      112      100
## sigsq.eps 0.004102877 0.0028937522      514      890

TracePlot(moreiterations, par="beta")
TracePlot(moreiterations, par="r")


# move to sequential processing (end of vignette)
future::plan(strategy = future::sequential)

Acknowledgments

Thanks to Haotian “Howie” Wu for invaluable feedback on early versions of the package.