/
Functions.R
2170 lines (2003 loc) · 109 KB
/
Functions.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#DESCRIPTION: self-explanatory
expit = function(x) { 1/(1+exp(-x));}
logit = function(x) { log(x/(1-x));}
#DESCRIPTION: Error-handling function
#Borrowed from the R package simsalapar
#https://www.rdocumentation.org/packages/simsalapar/versions/1.0-9/topics/tryCatch.W.E
tryCatch.W.E <- function(expr)
{
W <- NULL
w.handler <- function(w){ # warning handler
W <<- c(W,w)
invokeRestart("muffleWarning")
}
list(value = withCallingHandlers(tryCatch(expr, error = function(e) e),
warning = w.handler),
warning = W)
}
#DESCRIPTION: Simulator function for drawing binary outcomes (historical, current, and new [for validation]),
#the probabilities of which are logistic-linear functions of normal and/or bernoulli distributed
#predictors.
#
#
#ARGUMENTS:
#
#n_hist (pos. integer) size of historical data; n_hist in the paper
#
#n_curr (pos. integer) size of current data; n_curr in the paper
#
#n_new (pos. integer) size of testing data (for prediction)
#
#true_mu_hist (real) true intercept for generating historical model. mu_hist in Boonstra and Barbaro
#
#true_mu_curr (real) true intercept for generating current / new model. mu in Boonstra and Barbaro
#
#true_betas_orig (vector) true regression coefficients corresponding to original covariates. Beta^o in Boonstra and Barbaro
#
#true_betas_aug (vector) true regression coefficients corresponding to augmented covariates. Beta^a in Boonstra and Barbaro
#
#covariate_args (list) the named arguments are simple ways to govern the distribution of the predictors.
#x_correlation is the normal correlation between any pair of predictors; x_orig_binom is the
#integer indices (any subset of 1, ..., length(true_betas_orig)) indicating which of the original
#covariates should be transformed to binary values based upon being less than or greater than zero;
#x_aug_binom is the analogous set of indices for the augmented predictors (any subset of 1, ..., length(true_betas_aug))
draw_data = function(n_hist = 150,
n_curr = 50,
n_new = 100,
true_mu_hist = 0,
true_mu_curr = 0,
true_betas_orig = 0,
true_betas_aug = 0,
covariate_args = list(x_correlation = 0,
x_orig_binom = NULL,
x_aug_binom = NULL)
) {
p = length(true_betas_orig);
q = length(true_betas_aug);
stopifnot(length(true_mu_hist) == 1 && length(true_mu_curr) == 1);
x_all = matrix(rnorm((n_hist + n_curr + n_new)*(p+q)),nrow=n_hist + n_curr + n_new)%*%chol(diag(1 - covariate_args$x_correlation,p+q) + covariate_args$x_correlation);#original covariates are N(0,1)
x_all_orig = x_all[,1:p,drop = F];
#Binary covariates will be -1 or 1 with equal prevalence (assuming the latent normal is mean zero),
#which will result in a random variable with mean zero and variance 1
if(length(covariate_args$x_orig_binom)) {
x_all_orig[,covariate_args$x_orig_binom] = 2 * (x_all_orig[,covariate_args$x_orig_binom,drop = F] > 0) - 1;
}
x_all_aug = x_all[,(p + 1):(p + q),drop = F];
if(length(covariate_args$x_aug_binom)) {
x_all_aug[,covariate_args$x_aug_binom] = 2 * (x_all_aug[,covariate_args$x_aug_binom,drop = F] > 0) - 1;
}
#Linear predictors (two for each observation: contribution from original covariates and augmented covariates)
lin_pred_x_orig = drop(x_all_orig%*%true_betas_orig);
lin_pred_x_aug = drop(x_all_aug%*%true_betas_aug);
#Note the difference in intercepts between historical data and current data
risk_all = 1/(1+exp(-c(rep(true_mu_hist,n_hist),rep(true_mu_curr,n_curr+n_new)) - lin_pred_x_orig - lin_pred_x_aug));
y_all = rbinom(n_hist + n_curr + n_new, 1, risk_all);
x_hist_orig = x_all_orig[1:n_hist,,drop=F];
x_curr_orig = x_all_orig[(n_hist + 1):(n_hist + n_curr),,drop=F];
x_new_orig = x_all_orig[(n_hist + n_curr + 1):(n_hist + n_curr + n_new),,drop=F];
#
x_hist_aug = x_all_aug[1:n_hist,,drop=F];
x_curr_aug = x_all_aug[(n_hist + 1):(n_hist + n_curr),,drop=F];
x_new_aug = x_all_aug[(n_hist + n_curr + 1):(n_hist + n_curr + n_new),,drop=F];
#
y_hist = y_all[1:n_hist];
y_curr = y_all[(n_hist+1):(n_hist + n_curr)];
y_new = y_all[(n_hist + n_curr + 1):(n_hist + n_curr + n_new)];
risk_new = risk_all[(n_hist + n_curr + 1):(n_hist + n_curr + n_new)];
list(x_hist_orig = x_hist_orig,
x_hist_aug = x_hist_aug,
y_hist = y_hist,
x_curr_orig = x_curr_orig,
x_curr_aug = x_curr_aug,
y_curr = y_curr,
x_new_orig = x_new_orig,
x_new_aug = x_new_aug,
y_new = y_new,
lin_pred_x_orig = lin_pred_x_orig,
lin_pred_x_aug = lin_pred_x_aug,
risk_new = risk_new)
}
#DESCRIPTION: Program for fitting a GLM equipped with the 'standard' prior evaluated
#in Boonstra and Barbaro, which is the regularized horseshoe. It has two intended uses:
#compile stan scripts or analyze data. First, if the user provides nothing but a valid
#'stan_path', then the stan script is compiled. Second, the user provides both a compiled
#stanfit object as well asvalues for y, x_standardized, #, q, and any other desired
#arguments to actually fit a regression.
#
#
#ARGUMENTS:
#
#stan_fit: an R object of class stanfit, which allows the function to run without recompiling the stan code.
#
#stan_path: (character) a path pointing to a .stan file, which indicates the stan code to compile and run. If
#both stan_fit and stan_path are provided, stan_fit takes precedence.
#
#y (vector) outcomes corresponding to the type of glm desired. This should match whatever datatype is expected
#by the stan program.
#
#x_standardized (matrix) matrix of numeric values with number of rows equal to the length of y and number of columns
#equal to p+q. It is assumed without verification that each column is standardized to whatever scale the prior
#expects - in Boonstra and Barbaro, all predictors are marginally generated to have mean zero and unit variance, so no
#standardization is conducted. In practice, all data should be standardized to have a common scale before model fitting.
#If regression coefficients on the natural scale are desired, they be easily obtained through unstandardizing.
#
#p, q (nonneg. integers) numbers, the sum of which add up to the number of columns in x_standardized. For the standard
#prior, this distinction is only needed if a different constant scale parameter (beta_orig_scale, beta_aug_scale), which is
#the constant 'c' in the notation of Boonstra and Barbaro, is used.
#
#beta_orig_scale, beta_aug_scale (pos. real) constants indicating the prior scale of the horseshoe. Both values correspond
#to 'c' in the notation of Boonstra and Barbaro, because that paper never considers beta_orig_scale!=beta_aug_scale
#
#local_dof, global_dof (pos. integer) numbers indicating the degrees of freedom for lambda_j and tau, respectively. Boonstra,
#et al. never considered local_dof != 1 or global_dof != 1.
#
#slab_precision (pos. real) the slab-part of the regularized horseshoe, this is equivalent to (1/d)^2 in the notation of
#Boonstra and Barbaro
#
#intercept_offset (vector) vector of 0's and 1's equal having the same length as y. Those observations with a value of 1
#have an additional constant offset in their linear predictor, effectively a different intercept. This is useful to jointly
#regress two datasets in which it is believed that the regression coefficients are the same but not the intercepts and could be
#useful (but was not used) in the simulation study to compare to a benchmark, namely if both the historical and current datasets
#were available but there is a desire to adjust for potentially different baseline prevalences.
#
#only_prior (logical) should all data be ignored, sampling only from the prior?
#
#ntries (pos. integer) the stan function will run up to this many times, stopping either when the number of
#*divergent transitions* is zero or when ntries has been reached. The reported fit will be that with the fewest number of
#divergent iterations.
glm_standard = function(stan_fit = NA,
stan_path,
y = c(0,1),
x_standardized = matrix(0,length(y),6),
p = 3,
q = 3,
beta_orig_scale = 1,
beta_aug_scale = 1,
local_dof = 1,
global_dof = 1,
slab_precision = (1/15)^2,
intercept_offset = NULL,
only_prior = F,
mc_warmup = 50,
mc_iter_after_warmup = 50,
mc_chains = 1,
mc_thin = 1,
mc_stepsize = 0.1,
mc_adapt_delta = 0.9,
mc_max_treedepth = 15,
ntries = 1) {
stopifnot(ncol(x_standardized) == (p+q));
max_divergences = -Inf;
accepted_divergences = Inf;
curr_try = 1;
if(is.null(intercept_offset)) {intercept_offset = numeric(length(y));}
while(curr_try <= ntries) {
assign("curr_fit",tryCatch.W.E(stan(file = stan_path,
fit = stan_fit,
data = list(n_stan = length(y),
p_stan = p,
q_stan = q,
y_stan = y,
x_standardized_stan = x_standardized,
local_dof_stan = local_dof,
global_dof_stan = global_dof,
beta_orig_scale_stan = beta_orig_scale,
beta_aug_scale_stan = beta_aug_scale,
slab_precision_stan = slab_precision,
intercept_offset_stan = intercept_offset,
only_prior = as.integer(only_prior)),
warmup = mc_warmup,
iter = mc_iter_after_warmup + mc_warmup,
chains = mc_chains,
thin = mc_thin,
control = list(stepsize = mc_stepsize,
adapt_delta = mc_adapt_delta,
max_treedepth = mc_max_treedepth))));
if("simpleError"%in%class(curr_fit$value) || "error"%in%class(curr_fit$value)) {
stop(curr_fit$value);
}
if(!"stanfit"%in%class(stan_fit)) {
break;
}
divergent_check = unlist(lapply(curr_fit$warning,grep,pattern="divergent transitions",value=T));
rhat_check = max(summary(curr_fit$value)$summary[,"Rhat"],na.rm=T);
#Originally, the break conditions were baesd upon having both no divergent transitions as well as a max Rhat (i.e. gelman-rubin
#diagnostic) sufficiently close to 1. I subsequently changed the conditions to be based only upon the first, which is reflected
#by setting rhat = T immediately below.
break_conditions = c(divergence = F, rhat = T);
if(length(divergent_check) == 0) {#corresponds to zero divergent transitions
curr_divergences = 0;
max_divergences = max(max_divergences,curr_divergences,na.rm=T);
break_conditions["divergence"] = T;
} else {#corresponds to > zero divergent transitions
curr_divergences <- max(as.numeric(strsplit(divergent_check," ")$message),na.rm=T);
max_divergences = max(max_divergences,curr_divergences,na.rm=T);
curr_try = curr_try + 1;
}
#update if fewer divergent transitions were found
if(curr_divergences < accepted_divergences) {
accepted_divergences = curr_divergences;
max_rhat = rhat_check;
foo = rstan::extract(curr_fit$value);
hist_beta0 = as.numeric(foo$mu);
curr_beta0 = as.numeric(foo$mu) + as.numeric(foo$mu_offset);
curr_beta = foo$beta;
theta_orig = foo$theta_orig;
theta_aug = foo$theta_aug;
}
if(all(break_conditions)) {
break;
}
}
if(!"stanfit"%in%class(stan_fit)) {
curr_fit$value;
} else {
list(accepted_divergences = accepted_divergences,
max_divergences = max_divergences,
max_rhat = max_rhat,
hist_beta0 = hist_beta0,
curr_beta0 = curr_beta0,
curr_beta = curr_beta,
theta_orig = theta_orig,
theta_aug = theta_aug);
}
}
#DESCRIPTION: Program for fitting a GLM equipped with the 'naive adaptive bayes' prior evaluated in the manuscript
#
#
#ARGUMENTS: (only those distinct from glm_standard are discussed)
#
#alpha_prior_mean (vector) p-length vector giving the mean of alpha from the historical analysis,
#corresponds to m_alpha in Boonstra and Barbaro
#
#alpha_prior_cov (matrix) pxp positive definite matrix giving the variance of alpha from the historical
#analysis, corresponds to S_alpha in Boonstra and Barbaro
#
#phi_mean (real) mean of phi corresponding to a normal distribution, support is truncated to [0,1]
#
#phi_sd (pos. real) sd of phi corresponding to a normal distribution, support is truncated to [0,1]
#
#beta_aug_scale_tilde (pos. real) constant indicating the prior scale of the horseshoe for the augmented
#covariates when phi = 1, i.e. when the historical analysis is fully used. This corresponds to tilde_c in
#Boonstra and Barbaro
glm_nab = function(stan_fit = NA,
stan_path,
y = c(0,1),
x_standardized = matrix(0,length(y),6),
alpha_prior_mean = rep(0, 3),
alpha_prior_cov = diag(1, 3),
phi_mean = 0.5,
phi_sd = 2.5,
beta_orig_scale = 1,
beta_aug_scale = 1,
beta_aug_scale_tilde = 1,
local_dof = 1,
global_dof = 1,
slab_precision = (1/15)^2,
only_prior = F,
mc_warmup = 50,
mc_iter_after_warmup = 50,
mc_chains = 1,
mc_thin = 1,
mc_stepsize = 0.1,
mc_adapt_delta = 0.9,
mc_max_treedepth = 15,
ntries = 1,
eigendecomp_hist_var = NULL,
scale_to_variance225 = NULL
) {
if(is.null(eigendecomp_hist_var)) {
eigendecomp_hist_var = eigen(alpha_prior_cov);
}
eigenvec_hist_var = t(eigendecomp_hist_var$vectors);
sqrt_eigenval_hist_var = sqrt(eigendecomp_hist_var$values);
if(is.null(scale_to_variance225)) {
scale_to_variance225 = diag(alpha_prior_cov) / 225;
}
p = length(alpha_prior_mean);
q = ncol(x_standardized) - p;
if(p == 1) {
alpha_prior_mean = array(alpha_prior_mean,dim=1);
sqrt_eigenval_hist_var = array(sqrt_eigenval_hist_var,dim=1);
scale_to_variance225 = array(scale_to_variance225,dim=1);
}
max_divergences = -Inf;
accepted_divergences = Inf;
curr_try = 1;
while(curr_try <= ntries) {
assign("curr_fit",tryCatch.W.E(stan(file = stan_path,
fit = stan_fit,
data = list(n_stan = length(y),
p_stan = p,
q_stan = q,
y_stan = y,
x_standardized_stan = x_standardized,
alpha_prior_mean_stan = alpha_prior_mean,
alpha_prior_cov_stan = alpha_prior_cov,
sqrt_eigenval_hist_var_stan = sqrt_eigenval_hist_var,
eigenvec_hist_var_stan = eigenvec_hist_var,
local_dof_stan = local_dof,
global_dof_stan = global_dof,
beta_orig_scale_stan = beta_orig_scale,
beta_aug_scale_stan = beta_aug_scale,
beta_aug_scale_tilde_stan = beta_aug_scale_tilde,
slab_precision_stan = slab_precision,
scale_to_variance225 = scale_to_variance225,
phi_mean_stan = phi_mean,
phi_sd_stan = phi_sd,
only_prior = as.integer(only_prior)),
warmup = mc_warmup,
iter = mc_iter_after_warmup + mc_warmup,
chains = mc_chains,
thin = mc_thin,
control = list(stepsize = mc_stepsize,
adapt_delta = mc_adapt_delta,
max_treedepth = mc_max_treedepth))));
if("simpleError"%in%class(curr_fit$value) || "error"%in%class(curr_fit$value)) {
stop(curr_fit$value);
}
if(!"stanfit"%in%class(stan_fit)) {
break;
}
divergent_check = unlist(lapply(curr_fit$warning,grep,pattern="divergent transitions",value=T));
rhat_check = max(summary(curr_fit$value)$summary[,"Rhat"],na.rm=T);
#Originally, the break conditions were baesd upon having both no divergent transitions as well as a max Rhat (i.e. gelman-rubin
#diagnostic) sufficiently close to 1. I subsequently changed the conditions to be based only upon the first, which is reflected
#by setting rhat = T immediately below.
break_conditions = c(divergence = F, rhat = T);
if(length(divergent_check) == 0) {#corresponds to zero divergent transitions
curr_divergences = 0;
max_divergences = max(max_divergences,curr_divergences,na.rm=T);
break_conditions["divergence"] = T;
} else {#corresponds to > zero divergent transitions
curr_divergences <- max(as.numeric(strsplit(divergent_check," ")$message),na.rm=T);
max_divergences = max(max_divergences,curr_divergences,na.rm=T);
curr_try = curr_try + 1;
}
#update if fewer divergent transitions were found
if(curr_divergences < accepted_divergences) {
accepted_divergences = curr_divergences;
max_rhat = rhat_check;
foo = rstan::extract(curr_fit$value);
curr_beta0 = as.numeric(foo$mu);
curr_beta = foo$beta;
theta_orig = foo$theta_orig;
theta_aug = foo$theta_aug;
phi = foo$phi_copy;
eta = foo$eta;
}
if(all(break_conditions)) {
break;
}
}
if(!"stanfit"%in%class(stan_fit)) {
curr_fit$value;
} else {
list(accepted_divergences = accepted_divergences,
max_divergences = max_divergences,
max_rhat = max_rhat,
curr_beta0 = curr_beta0,
curr_beta = curr_beta,
theta_orig = theta_orig,
theta_aug = theta_aug,
phi = phi,
eta = eta);
}
}
#DESCRIPTION: Program for fitting a GLM equipped with the 'sensible adaptive bayes' prior evaluated in the manuscript
#
#
#ARGUMENTS: (only those distinct from glm_standard and glm_nab are discussed)
#
#aug_projection: (matrix) pxq matrix that approximately projects the regression coefficients of
#the augmented predictors onto the space of the regression coefficients for the original predictors.
#This is the matrix P in the notation of Boonstra and Barbaro. It can be calculated using the function
#'create_projection'.
#
#eigendecomp_hist_var: R object of class 'eigen' containing a pxp matrix of eigenvectors in each row
#(equivalent to v_0 in Boonstra and Barbaro) and a p-length vector of eigenvalues. This is by default
#equal to eigen(alpha_prior_cov)
#
#scale_to_variance225: a vector assumed to be such that, when multiplied by the diagonal elements of
#alpha_prior_cov, the result is a vector of elements each equal to 225. This is explicitly calculated
#if it is not provided
glm_sab = function(stan_fit = NA,
stan_path,
y = c(0,1),
x_standardized = matrix(0,length(y),6),
alpha_prior_mean = rep(0, 3),
alpha_prior_cov = diag(1, 3),
aug_projection = diag(1, 3),
phi_mean = 0.5,
phi_sd = 2.5,
beta_orig_scale = 1,
beta_aug_scale = 1,
local_dof = 1,
global_dof = 1,
slab_precision = (1/15)^2,
only_prior = F,
mc_warmup = 50,
mc_iter_after_warmup = 50,
mc_chains = 1,
mc_thin = 1,
mc_stepsize = 0.1,
mc_adapt_delta = 0.9,
mc_max_treedepth = 15,
ntries = 1,
eigendecomp_hist_var = NULL,
scale_to_variance225 = NULL
) {
if(is.null(eigendecomp_hist_var)) {
eigendecomp_hist_var = eigen(alpha_prior_cov);
}
eigenvec_hist_var = t(eigendecomp_hist_var$vectors);
sqrt_eigenval_hist_var = sqrt(eigendecomp_hist_var$values);
if(is.null(scale_to_variance225)) {
scale_to_variance225 = diag(alpha_prior_cov) / 225;
}
p = length(alpha_prior_mean);
q = ncol(x_standardized) - p;
if(p == 1) {
alpha_prior_mean = array(alpha_prior_mean,dim=1);
sqrt_eigenval_hist_var = array(sqrt_eigenval_hist_var,dim=1);
scale_to_variance225 = array(scale_to_variance225,dim=1);
}
max_divergences = -Inf;
accepted_divergences = Inf;
curr_try = 1;
while(curr_try <= ntries) {
assign("curr_fit",tryCatch.W.E(stan(file = stan_path,
fit = stan_fit,
data = list(n_stan = length(y),
p_stan = p,
q_stan = q,
y_stan = y,
x_standardized_stan = x_standardized,
aug_projection_stan = aug_projection,
alpha_prior_mean_stan = alpha_prior_mean,
alpha_prior_cov_stan = alpha_prior_cov,
sqrt_eigenval_hist_var_stan = sqrt_eigenval_hist_var,
eigenvec_hist_var_stan = eigenvec_hist_var,
local_dof_stan = local_dof,
global_dof_stan = global_dof,
beta_orig_scale_stan = beta_orig_scale,
beta_aug_scale_stan = beta_aug_scale,
slab_precision_stan = slab_precision,
scale_to_variance225 = scale_to_variance225,
phi_mean_stan = phi_mean,
phi_sd_stan = phi_sd,
only_prior = as.integer(only_prior)),
warmup = mc_warmup,
iter = mc_iter_after_warmup + mc_warmup,
chains = mc_chains,
thin = mc_thin,
control = list(stepsize = mc_stepsize,
adapt_delta = mc_adapt_delta,
max_treedepth = mc_max_treedepth))));
if("simpleError"%in%class(curr_fit$value) || "error"%in%class(curr_fit$value)) {
stop(curr_fit$value);
}
if(!"stanfit"%in%class(stan_fit)) {
break;
}
divergent_check = unlist(lapply(curr_fit$warning,grep,pattern="divergent transitions",value=T));
rhat_check = max(summary(curr_fit$value)$summary[,"Rhat"],na.rm=T);
#Originally, the break conditions were baesd upon having both no divergent transitions as well as a max Rhat (i.e. gelman-rubin
#diagnostic) sufficiently close to 1. I subsequently changed the conditions to be based only upon the first, which is reflected
#by setting rhat = T immediately below.
break_conditions = c(divergence = F, rhat = T);
if(length(divergent_check) == 0) {#corresponds to zero divergent transitions
curr_divergences = 0;
max_divergences = max(max_divergences,curr_divergences,na.rm=T);
break_conditions["divergence"] = T;
} else {#corresponds to > zero divergent transitions
curr_divergences <- max(as.numeric(strsplit(divergent_check," ")$message),na.rm=T);
max_divergences = max(max_divergences,curr_divergences,na.rm=T);
curr_try = curr_try + 1;
}
#update if fewer divergent transitions were found
if(curr_divergences < accepted_divergences) {
accepted_divergences = curr_divergences;
max_rhat = rhat_check;
foo = rstan::extract(curr_fit$value);
curr_beta0 = as.numeric(foo$mu);
curr_beta = foo$beta;
theta_orig = foo$theta_orig;
theta_aug = foo$theta_aug;
phi = foo$phi_copy;
eta = foo$eta;
}
if(all(break_conditions)) {
break;
}
}
if(!"stanfit"%in%class(stan_fit)) {
curr_fit$value;
} else {
list(accepted_divergences = accepted_divergences,
max_divergences = max_divergences,
max_rhat = max_rhat,
curr_beta0 = curr_beta0,
curr_beta = curr_beta,
theta_orig = theta_orig,
theta_aug = theta_aug,
phi = phi,
eta = eta);
}
}
#DESCRIPTION: Helper function to make projection matrix (or list of projection matrices), which is P in the notation
#of Boonstra and Barbaro
#
#VALUE: A list as long as the the length of imputes_list, with each element containing a different projection matrix
#using the indices of the imputations specified in the corresponding element of imputes_list.
#
#ARGUMENTS:
#
#x_curr_orig, x_curr_aug (matrices) matrices with equal numbers of rows and p & q columns,
#respectively. These are used to estimate the joint association between the original
#and augmented covariates, which is needed for the imputation model
#
#eigenvec_hist_var (matrix) pxp matrix with each row corresponding to an eigenvector.
#This is v_0 in Boonstra and Barbaro.
#
#imputes_list (list) list of length-2 vectors, with each vector containing the lower and upper indices of the imputations to use
#for a projection matrix in the SAB method. This is best explained with an example: if imputes_list = list(c(1,15),c(16,100),c(1,100)),
#then three projection matrices will be returned. One will be based upon the first 15 imputations from a run of MICE, the second based upon
#the last 85 imputations from that same run (i.e. the 16th-100th imputations), and the third will be based upon all 100 imputations from this
#same run. This is coded as such to allow for flexible exploration of the impact of number of imputations or variability due to imputations.
#
#seed_start (pos. integer) random seed to start each imputation
#
#predictorMatrix: (matrix) (p+q)x(p+q) matrix equivalent to argument of the same name in
#in the 'mice()' function (type '?mice'). It is specially calculated based upon a monotone
#missingness pattern (x^o is fully observed, x^a is not) Thus it samples from
#[X^a_{1...q}|X^o] = [X^a_1|X^o]*[X^a_2|X^o,X^a_1]*...
create_projection = function(x_curr_orig,
x_curr_aug,
eigenvec_hist_var,
imputes_list = list(c(1,15)),
seed_start = sample(.Machine$integer.max,1),
predictorMatrix = NULL) {
require(mice);
require(magrittr);
p = ncol(x_curr_orig);
q = ncol(x_curr_aug);
orig_covariates = colnames(x_curr_orig);
aug_covariates = colnames(x_curr_aug);
x_all = cbind(x_curr_orig,x_curr_aug);
stopifnot(class(imputes_list) == "list");
n_imputes = max(unlist(lapply(imputes_list,max)));
if(is.null(predictorMatrix)) {
predictorMatrix11 = matrix(0,nrow = p,ncol = p);
predictorMatrix12 = matrix(0,nrow = p,ncol = q);
predictorMatrix21 = matrix(1,nrow = q,ncol = p);
predictorMatrix22 = matrix(1,nrow = q,ncol = q);
predictorMatrix22[upper.tri(predictorMatrix22,diag = T)] = 0;
predictorMatrix = rbind(cbind(predictorMatrix11,predictorMatrix12),
cbind(predictorMatrix21,predictorMatrix22));
colnames(predictorMatrix) = rownames(predictorMatrix) = c(orig_covariates,aug_covariates);
rm(predictorMatrix11,predictorMatrix12,predictorMatrix21,predictorMatrix22);
} else {
stopifnot(dim(predictorMatrix) == c(p+q,p+q));
}
#Create data to be marginalized
#Column 1 identifies each eigenvector, including an additional row for the intercept offset;
#Columns 2:(p+1) are original covariates;
#Columns (p+2):(p+q+1) are the augmented covariates to be imputed and averaged;
dat_for_marg = cbind(c(1:p,nrow(x_all)+1),
rbind(eigenvec_hist_var,0),
matrix(NA,p+1,q));
colnames(dat_for_marg) = c("ID",orig_covariates,aug_covariates);
dat_for_marg = data.matrix(dat_for_marg);
#Store one copy for each of the different imputations to use
dat_for_marg = lapply(1:length(imputes_list), function(x) dat_for_marg);
#Now impute the augmented covariates based upon the artificially created original covariates, i.e. the eigenvectors, using
#the empiric associations in current data.
curr_impute = mice(rbind(x_all,
#original covariates are constant across the different imputations, so we can just use the first
dat_for_marg[[1]][,c(orig_covariates,aug_covariates)]),
printFlag = F,
predictorMatrix = predictorMatrix,#user-provided, or if missing, created above
m = n_imputes,#
maxit = 1,#only one iteration is needed, due to the monotone missingness pattern
method = "pmm",
seed = seed_start);
curr_impute = data.matrix(mice::complete(curr_impute,"long"));
for(k in 1:(p+1)) {
for(j in 1:length(imputes_list)) {
#Fill in the data_for marg
dat_for_marg[[j]][dat_for_marg[[j]][,"ID"] == dat_for_marg[[j]][k,"ID"], aug_covariates] =
colMeans(curr_impute[which(curr_impute[,".id"] == k + nrow(x_all))[imputes_list[[j]][1]:imputes_list[[j]][2]],aug_covariates,drop=F]);
}
}
#If there is any perfect correlation between predictors, the imputer will return NA. In this case, we fill in the missing data
#with whatever variable was perfectly correlated with it
for(j in 1:length(imputes_list)) {
while(any(is.na(dat_for_marg[[j]][,aug_covariates]))) {
correlated_column = which(colSums(is.na(dat_for_marg[[j]][,aug_covariates]))>0)[1];
dat_for_marg[[j]][,names(correlated_column)] = dat_for_marg[[j]][,setdiff(colnames(x_all)[abs(cor(x_all,x_curr_aug[,correlated_column]) - 1) < sqrt(.Machine$double.eps)],names(correlated_column))[1]];
}
}
#v0_inv is constant across the different imputation sets, so we can just use the first set
v0_inv = solve(dat_for_marg[[1]][1:p,orig_covariates,drop=F]);
projections = vector("list",length(imputes_list));
for(j in 1:length(imputes_list)) {
projections[[j]] = v0_inv %*% (dat_for_marg[[j]][1:p,aug_covariates,drop=F] - dat_for_marg[[j]][rep(p+1,p),aug_covariates,drop=F]);
}
projections;
}
#DESCRIPTION: This is the parent function in the simulation study. For a given data-generating mechanism (characterized by choices
#of n_hist, n_curr, true_mu_hist, true_mu_curr, true_betas_orig, true_betas_aug, and covariate_args) and modeling choice (characterized
#by choices of local_dof, global_dof, slab_precision, nab_augmented_scale, power_prop_nonzero_prior, and sab_imputes_list, phi_params),
#all of the methods in Boonstra and Barbaro are run against an arbitrary number of simulated datasets. The user can modify various
#characteristics of the underlying HMC chain. A number of operating characteristics are returned, based both on estimation and prediction.
#
#
#VALUE: A list of various results.
#
#
#ARGUMENTS:
#sim_number (arbitrary) This is a label to help the user keep track of multiple different scenario settings
#It is simply returned at the end of the function
#
#array_id (pos. integer) This is intended to be the slurm array id
#
#n_sim (pos. integer) number of simulated datasets to construct
#
#n_hist (pos. integer) size of historical data; n_hist in the paper
#
#n_curr (pos. integer) size of current data; n_curr in the paper
#
#n_new (pos. integer) size of testing data (for prediction)
#
#true_mu_hist (real) true intercept for generating historical model. mu_hist in Boonstra and Barbaro.
#
#true_mu_curr (real) true intercept for generating current / new model. mu in Boonstra and Barbaro.
#
#true_betas_orig (vector) true regression coefficients corresponding to original covariates. Beta^o in Boonstra and Barbaro
#
#true_betas_aug (vector) true regression coefficients corresponding to augmented covariates. Beta^a in Boonstra and Barbaro
#
#covariate_args (list) the named arguments are simple ways to govern the distribution of the predictors. Beta^a in Boonstra and Barbaro.
#
#beta_label (arbitrary) This is a label to help the user keep track of multiple different true regression
#coefficients. It is simply returned at the end of the function.
#
#covariate_args (list) the named arguments are simple ways to govern the distribution of the predictors.
#x_correlation is the normal correlation between any pair of predictors; x_orig_binom is the
#integer indices (any subset of 1, ..., length(true_betas_orig)) indicating which of the original
#covariates should be transformed to binary values based upon being less than or greater than zero;
#x_aug_binom is the analogous set of indices for the augmented predictors (any subset of 1, ..., length(true_betas_aug))
#
#covariate_label (arbitrary) This is a label to help the user keep track of multiple different true covariate
#distributions It is simply returned at the end of the function
#
#local_dof, global_dof (pos. integer) numbers indicating the degrees of freedom for lambda_j and tau, respectively. Boonstra,
#et al. never considered local_dof != 1 or global_dof != 1.
#
#slab_precision (pos. real) the slab-part of the regularized horseshoe, this is equivalent to (1/d)^2 in the notation of
#Boonstra and Barbaro
#
#nab_augmented_scale (pos. real) the scale parameter to accompany the tilde_lambda shrinkage of the augmented covariates
#when phi = 1 under NAB. This is equivalent to tilde_c in the notation of Boonstra and Barbaro.
#
#power_prop_nonzero_prior (pos. real in [0,1]) exponent for number covariates that are assumed to be non-zero minus 1/2,
#e.g. p^(1/3) - 0.5
#
#sab_imputes_list (list) list of length-2 vectors, with each vector containing the lower and upper indices of the imputations to use
#for a projection matrix in the SAB method. This is best explained with an example: if sab_imputes_list = list(c(1,15),c(16,100),c(1,100)),
#then three projection matrices will be constructed. One will be based upon the first 15 imputations from a run of MICE, the second based upon
#the last 85 imputations from that same run (i.e. the 16th-100th imputations), and the third will be based upon all 100 imputations from this
#same run. This is coded as such to allow for flexible exploration of the impact of number of imputations or variability due to imputations.
#Note that each of the projection matrices is crossed with each of the hyperpriors on phi (determined by the length of phi_params), i.e.
#i.e. if phi_params has length 3 and sab_imputes_list has length 2, then there will be six versions of SAB run and reported. Because NAB
#doesn't depend upon any imputation, there would be just 3 versions of NAB run and reported in this example.
#
#stan_file_path (character) local path to directory containing stan files
#
#standard_stan_filename (character) file name for standard prior
#
#sab_stan_filename (character) file name for sab prior
#
#sab_dev_stan_filename (character) file name for development version of sab prior (for testing)
#
#nab_stan_filename (character) file name for nab prior
#
#nab_dev_stan_filename (character) file name for development version of nab prior (for testing)
#
#phi_params (list) list of named lists, each with named components 'mean' and 'sd', corresponding to the mean and standard deviation, respectively,
#of a normal hyperprior on phi truncated to the [0,1] interval. Its length is the number of different hyperpriors on phi to be explored
#for each adaptive bayesian update. In Boonstra and Barbaro, 'Agnostic' corresponds to list(mean = 0.5, sd = 2.5), which is approximately uniform,
#and 'Optimistic' corresponds to list(mean = 1, sd = 0.25). Note that each of the hyperpriors on phi is crossed with each of the constructed
#projection matrices (determined by the length of sab_imputes_list), i.e. if phi_params has length 3 and sab_imputes_list has length 2,
#then there will be six versions of SAB run and reported. Because NAB doesn't depend upon any imputation, there would be just 3 versions of
#NAB run and reported in this example.
#
#mc_warmup (pos. integer) equivalent to warmup in 'stan()' function (type '?stan')
#
#mc_iter_after_warmup (pos. integer) equivalent to iter - warmup in 'stan()' function (type '?stan')
#
#mc_chains (pos. integer) equivalent to chains in 'stan()' function (type '?stan')
#
#mc_thin (pos. integer) equivalent to thin in 'stan()' function (type '?stan')
#
#mc_stepsize (pos. real) equivalent to stepsize in 'stan()' function (type '?stan')
#
#mc_adapt_delta_relaxed, mc_adapt_delta,strict (pos. real in [0,1]) Two alternatives to use as values for adapt_delta
#in 'stan()' function (type '?stan'). The relaxed version is presumably smaller and will be used for the
#standard prior, which has fewer numerical issues. The strict version is for the nab and sab priors.
#
#mc_max_treedepth (pos. integer) equivalent to max_treedepth in 'stan()' function (type '?stan')
#
#ntries_per_iter (pos. integer) each method will run up to this many times, stopping either when the number of *divergent transitions*
#is zero or when ntries has been reached. The reported fit will be that with the fewest number of divergent iterations.
#
#random_seed (pos. integer) where to initialize the simulator
#
#fit_marginal (logical) Set to TRUE to include an asymptotic estimate of the misspecified alpha_orig
#
#fit_methods (logical) Set to FALSE to do a dry run of the simulator
#
#skip_methods (vector) any methods falling in the intersection of skip_methods and
#c("Benchmark","Historical","Standard","NAB","NAB_dev","SAB","SAB_dev") will be skipped
#
#dynamic_run (logical) set to TRUE if stepping through this function and store more results
run.sim <- function(sim_number,
array_id,
n_sim,
n_hist,
n_curr,
n_new,
true_mu_hist,
true_mu_curr,
true_betas_orig,
true_betas_aug,
beta_label,
covariate_args,
covariate_label,
local_dof = 1,
global_dof = 1,
slab_precision = (1.0/15.0)^2,
nab_augmented_scale = 0.05,
power_prop_nonzero_prior = 1/3,
sab_imputes_list = list(c(1,100)),
stan_file_path = "",
standard_stan_filename = "RegHS_stable.stan",
sab_stan_filename = "SAB_stable.stan",
sab_dev_stan_filename = "SAB_dev.stan",
nab_stan_filename = "NAB_stable.stan",
nab_dev_stan_filename = "NAB_dev.stan",
phi_params = list("Agnostic" = c(mean = 0.5, sd = 2.5),
"Optimist" = c(mean = 1, sd = 0.25)),
mc_warmup = 1e3,
mc_iter_after_warmup = 5e3,
mc_chains = 2,
mc_thin = 1,
mc_stepsize = 0.01,
mc_adapt_delta_relaxed = 0.999,
mc_adapt_delta_strict = 0.9999999,
mc_max_treedepth = 20,
ntries_per_iter = 4,
random_seed = sample(.Machine$integer.max,1),
fit_marginal = T,#Set to TRUE to use all data to come up with an asymptotic estimate of the misspecified alpha_orig
fit_methods = T,#Set to FALSE to do a dry run of the simulator
skip_methods = c("NAB_dev","SAB_dev"),#methods falling in the intersection of skip_methods with c("Benchmark","Historical","Standard","NAB","NAB_dev","SAB","SAB_dev") will be skipped
dynamic_run = T) {
begin_all = Sys.time();
set.seed(random_seed);
data_seeds = sample(.Machine$integer.max,n_sim);
informational_messages = list();
stopifnot(class(sab_imputes_list) == "list");
sab_num_imputes_each = unlist(lapply(sab_imputes_list,diff)) + 1;
max_sab_index = max(unlist(lapply(sab_imputes_list,max)));
min_sab_index = min(unlist(lapply(sab_imputes_list,min)));
if(min_sab_index != 1) {
stop("The argument 'sab_imputes_list' requires that the smallest index overall be 1");
}
base_meth_names = c("Benchmark",
"Historical",
"Standard",
"NAB",
"NAB_dev",
"SAB",
"SAB_dev");
sab_suffix = expand.grid(paste0(".phi",names(phi_params)),
paste0(".imp",1:length(sab_imputes_list)));
expanded_meth_names = c("Benchmark",
"Historical",
"Standard",
paste0("NAB",names(phi_params)),
paste0("NAB_dev",names(phi_params)),
paste0("SAB",paste0(sab_suffix[,1],sab_suffix[,2])),
paste0("SAB_dev",paste0(sab_suffix[,1],sab_suffix[,2])));
rm(sab_suffix);
if("Historical" %in% skip_methods) {
skip_methods = union(skip_methods, c("NAB","NAB_dev","SAB","SAB_dev"));
informational_messages = c(informational_messages, paste0("Skipping all adaptive Bayesian methods because because Historical was skipped"));
}
if(("NAB" %in% skip_methods && !"NAB_dev" %in% skip_methods) || ("SAB" %in% skip_methods && !"SAB_dev" %in% skip_methods)) {
stop("Code assumes that NAB and SAB are to be fit if their corresponding development versions are to be fit");
}
mse_beta = #mean squared error for all betas
mse_beta_orig = #mean squared error for original betas
mse_beta_aug = #mean squared error for augmented betas
rmspe = #root mean squared prediction error
mape = #mean absolute predict error
mins_per_method = #minutes devoted to each method
max_rhat = #
matrix(NA,nrow=n_sim,ncol=length(expanded_meth_names),dimnames = list(NULL, expanded_meth_names));
store_impute_time = rep(NA, n_sim);
posterior_eff_orig = #posterior effective number of original parameters = mean(rowSums(1-kappa[orig]))
posterior_eff_aug = #posterior effective number of augmented parameters = mean(rowSums(1-kappa[aug]))
matrix(NA,nrow=n_sim,ncol=length(expanded_meth_names),dimnames = list(NULL,expanded_meth_names));
#Store characteristics of true risk in the population
population_params =
matrix(NA,
nrow=n_sim,
ncol=12,
dimnames = list(NULL,c("mean_hist",#historical prevalence
"mean_curr",#current prevalence
"signal_comp",#(n_hist+n_curr) * var(lin_pred_all[historical data+current data])
"signal_miss",#(n_hist) * var(lin_pred_aug[historical data])
"historical_signal_obs",#(n_hist) * var(lin_pred_orig[historical data])
"current_signal_obs",#(n_curr) * var(lin_pred_all[current data])
"mean_pop_risk",#mean of true population risk
"sd_pop_risk",#standard deviation of true population risk
"lowerq_pop_risk",#lower and upper quartiles of risk
"upperq_pop_risk",
"log_det_xprod_x_comp",
"log_det_xprod_x_miss")));
num_orig = length(true_betas_orig);
num_aug = length(true_betas_aug);
orig_covariates = paste0("orig",1:num_orig);
aug_covariates = paste0("aug",1:num_aug);
#####For evaluating overall MSE
matrix_true_beta = matrix(c(true_betas_orig,true_betas_aug),
nrow = mc_iter_after_warmup * mc_chains,
ncol = num_orig + num_aug,
byrow = T);
#####For evaluating MSE of original betas only
matrix_true_beta_orig = matrix(true_betas_orig,
nrow = mc_iter_after_warmup * mc_chains,
ncol = num_orig,
byrow = T);
#####For evaluating MSE of augmented betas only
matrix_true_beta_aug = matrix(true_betas_aug,
nrow = mc_iter_after_warmup * mc_chains,
ncol = num_aug,
byrow = T);
store_mean_beta = store_sd_beta = vector("list",length(expanded_meth_names));
names(store_mean_beta) = names(store_sd_beta) = expanded_meth_names;
for(k in 1:length(store_mean_beta)) {
store_mean_beta[[k]] = matrix(NA,nrow = n_sim,ncol=num_orig + num_aug,dimnames = list(NULL,c(orig_covariates,aug_covariates)));
store_sd_beta[[k]] = matrix(NA,nrow = n_sim,ncol=num_orig + num_aug,dimnames = list(NULL,c(orig_covariates,aug_covariates)));
}
rm(k);
#store posterior means and sds of phi (shrinkage weight) and eta (scale factor on S_alpha)
store_phi_mean =
store_phi_sd =
store_eta_mean =
store_eta_sd = matrix(NA, nrow = n_sim,
ncol = length(grep("NAB",expanded_meth_names)) +
length(grep("SAB",expanded_meth_names)),
dimnames = list(NULL, c(grep("NAB",expanded_meth_names,value=T),
grep("SAB",expanded_meth_names,value=T))));
#Calculate scale parameters so that E[1-sum kappa_j] = (# of parameters to fit)^(power_prop_nonzero_prior) - 0.5;
store_hierarchical_scales =
prior_eff = #prior effective number of original parameters = mean(rowSums(1-kappa[orig]))
vector("list",length(base_meth_names));
names(store_hierarchical_scales) = names(prior_eff) = base_meth_names;
if(fit_methods) {
#Benchmark: full access to data (n_hist + n_curr), but agnostic about covariates, so use skeptical hierarchical shrinkage
#with most assumed to be zero.
foo = solve_for_hiershrink_scale(target_mean1 = -0.5 + (num_orig + num_aug) ^ power_prop_nonzero_prior,
target_mean2 = NA,
npar1 = num_orig + num_aug,
npar2 = 0,
local_dof = local_dof,
regional_dof = -Inf,
global_dof = global_dof,
slab_precision = slab_precision,
n = n_hist + n_curr,
sigma = 2,
n_sim = round(2e6/(num_orig + num_aug)));
store_hierarchical_scales$Benchmark = foo$scale1;
prior_eff$Benchmark = foo$prior_num1;
rm(foo);
#Historical: access to historical data alone, with a skeptical prior. This is both a method in and of itself as well as
#the 'prior' analysis that will be provided to the SAB methods.
foo = solve_for_hiershrink_scale(target_mean1 = -0.5 + num_orig ^ power_prop_nonzero_prior,
target_mean2 = NA,
npar1 = num_orig,
npar2 = 0,
local_dof = local_dof,
regional_dof = -Inf,
global_dof = global_dof,
slab_precision = slab_precision,
n = n_hist,
sigma = 2,
n_sim = round(2e6/(num_orig + num_aug)));
store_hierarchical_scales$Historical = foo$scale1;
prior_eff$Historical = foo$prior_num1;
rm(foo);
#Standard: access to current data alone, with a skeptical prior.
#This also corresponds to the standard scales for SAB and NAB, where standard means that
#no historical information is used
foo = solve_for_hiershrink_scale(target_mean1 = -0.5 + (num_orig + num_aug) ^ power_prop_nonzero_prior,
target_mean2 = NA,
npar1 = num_orig + num_aug,
npar2 = 0,
local_dof = local_dof,
regional_dof = -Inf,
global_dof = global_dof,
slab_precision = slab_precision,
n = n_curr,
sigma = 2,
n_sim = round(2e6/(num_orig + num_aug)));
store_hierarchical_scales$Standard =
store_hierarchical_scales$NAB =
store_hierarchical_scales$SAB =
foo$scale1;
prior_eff$Standard =
prior_eff$NAB =
prior_eff$SAB =
foo$prior_num1;
rm(foo);
#
store_hierarchical_scales$NAB_aug_tilde = nab_augmented_scale;
}