search

dmphillippo / multinma

Network meta-analysis of individual and aggregate data in Stan
https://dmphillippo.github.io/multinma
35 stars 17 forks source link

test failure with provisional StanHeaders #4

Closed bgoodri closed 10 months ago

bgoodri commented 3 years ago

One of your tests is failing with the version of StanHeaders that we are struggling to get onto CRAN. Since the test has skip_on_cran(), it is not the worst thing in the world, but I think there is a mistake somewhere when calling line 303 in test-nma.R because I now get

SAMPLING FOR MODEL 'binomial_1par' NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: vector[min] indexing: accessing element out of range. index 3 out of range; expecting index to be between 1 and 2  (in 'model_binomial_1par' at line 19)

[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                                                     
[2] "  Exception: vector[min] indexing: accessing element out of range. index 3 out of range; expecting index to be between 1 and 2  (in 'model_binomial_1par' at line 19)"
error occurred during calling the sampler; sampling not done

I am not sure if or how that was working before with the StanHeaders / rstan currently on CRAN. The list of data passed to the data block of your Stan program looks like:

List of 47
 $ ns_ipd                  : int 0
 $ ni_ipd                  : num 0
 $ ns_agd_arm              : int 1
 $ ni_agd_arm              : int 2
 $ ns_agd_contrast         : num 0
 $ ni_agd_contrast         : num 0
 $ nt                      : num 2
 $ nint                    : num 1
 $ nX                      : int 2
 $ int_thin                : num 1
 $ narm_ipd                : num 0
 $ ipd_arm                 : num(0) 
 $ ipd_trt                 : num(0) 
 $ agd_arm_trt             : Named num [1:2] 1 2
  ..- attr(*, "names")= chr [1:2] "1" "2"
 $ agd_contrast_trt        : num[0 (1d)] 
 $ agd_contrast_trt_b      : num[0 (1d)] 
 $ agd_contrast_y          : num(0) 
 $ agd_contrast_Sigma      : num [1, 1] 1
 $ RE                      : num 0
 $ RE_cor                  : num [1, 1] 1
 $ which_RE                : int(0) 
 $ QR                      : logi FALSE
 $ X                       : num [1:2, 1:2] 1 1 0 1
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "1" "2"
  .. ..$ : chr [1:2] ".studyA" ".trtb"
 $ R_inv                   : num[0 , 0 ] 
 $ has_offset              : logi FALSE
 $ offsets                 : num(0) 
 $ prior_intercept_dist    : num 1
 $ prior_intercept_location: num 0
 $ prior_intercept_scale   : num 10
 $ prior_intercept_df      : num 0
 $ prior_trt_dist          : num 1
 $ prior_trt_location      : num 0
 $ prior_trt_scale         : num 10
 $ prior_trt_df            : num 0
 $ prior_reg_dist          : num 1
 $ prior_reg_location      : num 0
 $ prior_reg_scale         : num 10
 $ prior_reg_df            : num 0
 $ prior_het_dist          : num 1
 $ prior_het_location      : num 0
 $ prior_het_scale         : num 1
 $ prior_het_df            : num 0
 $ prior_het_type          : num 1
 $ ipd_r                   : int(0) 
 $ agd_arm_r               : num [1:2] 500 500
 $ agd_arm_n               : num [1:2] 1000 1000
 $ link                    : num 1
dmphillippo commented 3 years ago

Hi @bgoodri, thanks for getting in touch and sorry for the delay - I'm away from work at the moment.

All tests (including those skipped on CRAN) have been running and passing for me, both locally and on GitHub Actions. With the CRAN versions of rstan and StanHeaders I can't recreate the issue you are seeing; both running the test in question and passing the data block you pasted directly to sampling() works fine.

The test in question runs and passes for me:

library(multinma)
library(dplyr)
library(testthat)

single_study_a <- tibble(study = "A", trt = c("a", "b"), r = 500, n = 1000)
net_ss_a <- set_agd_arm(single_study_a, study, trt, r = r, n = n)

test_that("nma() doesn't fail with a single study", {
  skip_on_cran()

  fit_ss_a <- nma(net_ss_a, prior_intercept = normal(0, 10), prior_trt = normal(0, 10))
  expect_s3_class(fit_ss_a, "stan_nma")
  d_ss_a <- relative_effects(fit_ss_a)
  expect_equivalent(d_ss_a$summary$mean, 0, tol = 0.01)
  expect_equivalent(d_ss_a$summary$sd, sqrt(4/500), tol = 0.01)

})
Output

``` SAMPLING FOR MODEL 'binomial_1par' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 0 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds. Chain 1: Adjust your expectations accordingly! Chain 1: Chain 1: Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup) Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup) Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup) Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup) Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup) Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup) Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling) Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling) Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling) Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling) Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling) Chain 1: Iteration: 2000 / 2000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 0.052 seconds (Warm-up) Chain 1: 0.046 seconds (Sampling) Chain 1: 0.098 seconds (Total) Chain 1: SAMPLING FOR MODEL 'binomial_1par' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 0 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup) Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup) Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup) Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup) Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup) Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup) Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling) Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling) Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling) Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling) Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling) Chain 2: Iteration: 2000 / 2000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.059 seconds (Warm-up) Chain 2: 0.049 seconds (Sampling) Chain 2: 0.108 seconds (Total) Chain 2: SAMPLING FOR MODEL 'binomial_1par' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 0 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup) Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup) Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup) Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup) Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup) Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup) Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling) Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling) Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling) Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling) Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling) Chain 3: Iteration: 2000 / 2000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.043 seconds (Warm-up) Chain 3: 0.063 seconds (Sampling) Chain 3: 0.106 seconds (Total) Chain 3: SAMPLING FOR MODEL 'binomial_1par' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 0 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup) Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup) Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup) Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup) Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup) Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup) Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling) Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling) Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling) Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling) Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling) Chain 4: Iteration: 2000 / 2000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.062 seconds (Warm-up) Chain 4: 0.042 seconds (Sampling) Chain 4: 0.104 seconds (Total) Chain 4: Test passed ```

Passing the reconstructed data block directly works as expected too:

library(rstan)

standat <- list(
  ns_ipd                   = 0,
  ni_ipd                   = 0,
  ns_agd_arm               = 1,
  ni_agd_arm               = 2,
  ns_agd_contrast          = 0,
  ni_agd_contrast          = 0,
  nt                       = 2,
  nint                     = 1,
  nX                       = 2,
  int_thin                 = 1,
  narm_ipd                 = 0,
  ipd_arm                  = numeric(0),
  ipd_trt                  = numeric(0),
  agd_arm_trt              = c("1" = 1, "2" = 2),
  agd_contrast_trt         = numeric(0),
  agd_contrast_trt_b       = numeric(0),
  agd_contrast_y           = numeric(0),
  agd_contrast_Sigma       = matrix(1),
  RE                       = 0,
  RE_cor                   = matrix(1),
  which_RE                 = integer(0),
  QR                       = FALSE,
  X                        = matrix(c(1, 1, 0, 1), ncol = 2, 
                                    dimnames = list(1:2, c(".studyA", ".trtb"))),
  R_inv                    = matrix(0, ncol = 0, nrow = 0),
  has_offset               = FALSE,
  offsets                  = numeric(0),
  prior_intercept_dist     = 1,
  prior_intercept_location = 0,
  prior_intercept_scale    = 10,
  prior_intercept_df       = 0,
  prior_trt_dist           = 1,
  prior_trt_location       = 0,
  prior_trt_scale          = 10,
  prior_trt_df             = 0,
  prior_reg_dist           = 1,
  prior_reg_location       = 0,
  prior_reg_scale          = 10,
  prior_reg_df             = 0,
  prior_het_dist           = 1,
  prior_het_location       = 0,
  prior_het_scale          = 1,
  prior_het_df             = 0,
  prior_het_type           = 1,
  ipd_r                    = integer(0),
  agd_arm_r                = c(500, 500),
  agd_arm_n                = c(1000, 1000),
  link                     = 1)

fit <- sampling(multinma:::stanmodels$binomial_1par, data = standat)
fit
Output

``` SAMPLING FOR MODEL 'binomial_1par' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 0 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds. Chain 1: Adjust your expectations accordingly! Chain 1: Chain 1: Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup) Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup) Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup) Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup) Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup) Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup) Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling) Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling) Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling) Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling) Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling) Chain 1: Iteration: 2000 / 2000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 0.045 seconds (Warm-up) Chain 1: 0.032 seconds (Sampling) Chain 1: 0.077 seconds (Total) Chain 1: SAMPLING FOR MODEL 'binomial_1par' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 0 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup) Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup) Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup) Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup) Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup) Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup) Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling) Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling) Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling) Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling) Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling) Chain 2: Iteration: 2000 / 2000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.037 seconds (Warm-up) Chain 2: 0.036 seconds (Sampling) Chain 2: 0.073 seconds (Total) Chain 2: SAMPLING FOR MODEL 'binomial_1par' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 0 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup) Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup) Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup) Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup) Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup) Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup) Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling) Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling) Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling) Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling) Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling) Chain 3: Iteration: 2000 / 2000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.041 seconds (Warm-up) Chain 3: 0.048 seconds (Sampling) Chain 3: 0.089 seconds (Total) Chain 3: SAMPLING FOR MODEL 'binomial_1par' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 0 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup) Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup) Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup) Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup) Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup) Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup) Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling) Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling) Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling) Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling) Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling) Chain 4: Iteration: 2000 / 2000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.046 seconds (Warm-up) Chain 4: 0.051 seconds (Sampling) Chain 4: 0.097 seconds (Total) Chain 4: > fit Inference for Stan model: binomial_1par. 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 beta_tilde[1] 0.00 0.00 0.06 -0.13 -0.04 0.00 0.04 0.13 1309 1 beta_tilde[2] 0.00 0.00 0.09 -0.18 -0.06 0.00 0.06 0.18 1230 1 theta_agd_arm_bar[1] 0.50 0.00 0.02 0.47 0.49 0.50 0.51 0.53 1309 1 theta_agd_arm_bar[2] 0.50 0.00 0.02 0.47 0.49 0.50 0.51 0.53 2956 1 allbeta[1] 0.00 0.00 0.06 -0.13 -0.04 0.00 0.04 0.13 1309 1 allbeta[2] 0.00 0.00 0.09 -0.18 -0.06 0.00 0.06 0.18 1230 1 mu[1] 0.00 0.00 0.06 -0.13 -0.04 0.00 0.04 0.13 1309 1 d[1] 0.00 0.00 0.09 -0.18 -0.06 0.00 0.06 0.18 1230 1 fitted_agd_arm[1] 499.65 0.45 16.13 467.94 488.98 499.34 510.77 531.75 1309 1 fitted_agd_arm[2] 500.09 0.29 15.85 468.91 489.26 500.48 510.75 530.81 2956 1 log_lik[1] -4.20 0.02 0.71 -6.23 -4.37 -3.92 -3.73 -3.68 1192 1 log_lik[2] -4.18 0.02 0.71 -6.14 -4.34 -3.91 -3.73 -3.68 2062 1 resdev[1] 1.04 0.04 1.43 0.00 0.11 0.48 1.39 5.09 1192 1 resdev[2] 1.01 0.03 1.42 0.00 0.10 0.46 1.32 4.91 2062 1 lp__ -1387.32 0.03 1.00 -1389.96 -1387.74 -1387.02 -1386.60 -1386.32 1335 1 Samples were drawn using NUTS(diag_e) at Tue Aug 24 16:48:41 2021. 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). ```

I am struggling to figure out where the error is occurring for you as well, since line 19 of "model_binomial_1par" seems to be the beginning of the transformed parameters block.

If you let me know how to install the development version of StanHeaders (and rstan?) you are seeing the problems with then I can try to investigate further. Should I install the versions from the https://mc-stan.org/r-packages/ repo and build the multinma package against those?