Bayesian RE fractional polynomial NMA for grouped survival data

Tested on BEE R 3.6.1

Introduction

This vignette provides a short example of a Bayesian random effects fractional polynomial NMA model for grouped survival data. A random effect is put on the intercept term only (i.e. on the proportional hazards part of the polynomial) but not on the time-varying terms of the hazard ratio.

Prepare the environment

library(dplyr)
library(gemtc)        
library(gemtcPlus)     
library(ggmcmc)

Load in the data

data("grouped_TTE")

First order random intercept model (log-normal distribution for RE variance)

Plan the model

model_plan <-  plan_fp(model.pars = list(exponents = 0, t.eval = "midpoint"),
                       bth.model = "REINT", ref.std = "STUDY2", nma.ref.trt = "B",
                       bth.prior = list(type = "var", distr = "dlnorm", 
                                        meanlog = -2.71, 
                                        sdlog = 1.74),
                       n.chains = 2,
                       n.iter = 6000,
                       n.burnin = 1000,
                       n.thin = 1,
                       descr = "First order random intercept model (LN prior for RE Var)",
                       descr_s = "FP (1o, REint, LN)"
                       )

Ready the data

# Returns list list contaiing a jags list ready for input to `nma_fit` and a network object
model_input <- nma_pre_proc(grouped_TTE, model_plan)

Figure Network graph

plot(model_input$network, displaylabels = TRUE)

Fit the model

model  <- nma_fit(model_input = model_input)
## module glm loaded
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 675
##    Unobserved stochastic nodes: 20
##    Total graph size: 9191
## 
## Initializing model

Post processing

# Prepare plot data
nodes <- colnames(as.mcmc(model)[[1]])
sel <- grep("d[2,", nodes, fixed = TRUE)
plot_data <- ggs(as.mcmc(model)[, sel])

Produce diagnostic plots to further assess convergence. Here: select the contrasts trt 2 vs trt 1 for visibility.

Figure Traceplot

ggs_traceplot(plot_data)

Figure Densityplot

ggs_density(plot_data)

Figure Auto-correlation plot

ggs_autocorrelation(plot_data)

Figure Running means

ggs_running(plot_data)

The diagnostic plots show that the chains are much (!) too short in this example (for real use, run longer chains and consider thinning).

Save the results for later use.

res_re1 <- model
rm(model)

Second order random intercept model (log-normal distribution for RE variance)

Plan the model

model_plan <-  plan_fp(model.pars = list(exponents = c(0, 1), 
                                         t.eval = "midpoint"),
                       bth.model = "REINT", ref.std = "STUDY2", nma.ref.trt = "B",
                       bth.prior = list(type = "var", distr = "dlnorm", 
                                        meanlog = -2.71, 
                                        sdlog = 1.74),
                       n.chains = 2,
                       n.iter = 6000,
                       n.burnin = 1000,
                       n.thin = 1,
                       descr = "Second order random intercept model (LN prior for RE Var)",
                       descr_s = "FP (2o, REint, LN)"
                       )

Ready the data

# Returns list list contaiing a jags list ready for input to `nma_fit` and a network object
model_input <- nma_pre_proc(grouped_TTE, model_plan)

Fit the model

model  <- nma_fit(model_input = model_input)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 675
##    Unobserved stochastic nodes: 20
##    Total graph size: 10482
## 
## Initializing model

Post processing

# Prepare plot data
nodes <- colnames(as.mcmc(model)[[1]])
sel <- grep("d[2,", nodes, fixed = TRUE)
plot_data <- ggs(as.mcmc(model)[, sel])

Figure Traceplot

ggs_traceplot(plot_data)

The chains are obviously MUCH too short!

Save the results for later use.

res_re2 <- model
rm(model)

First order random intercept model (uniform distribution for RE standard deviation)

Plan the model

model_plan <-  plan_fp(model.pars = list(exponents = 0, t.eval = "midpoint"),
                       bth.model = "REINT", ref.std = "STUDY2", nma.ref.trt = "B",
                       bth.prior = list(type = "sd", distr = "unif", min = 0, max = 2),
                       n.chains = 2,
                       n.iter = 6000,
                       n.burnin = 1000,
                       n.thin = 1,
                       descr = "First order random intercept model (unif[0,2] prior)",
                       descr_s = "FP (1o, REint, U[0,2])"
                       )

Ready the data

# Returns list list contaiing a jags list ready for input to `nma_fit` and a network object
model_input <- nma_pre_proc(grouped_TTE, model_plan)

Fit the model

model  <- nma_fit(model_input = model_input)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 675
##    Unobserved stochastic nodes: 20
##    Total graph size: 9187
## 
## Initializing model

Post processing

Do convergence assessment as above (not shown).

The chains are obviously MUCH too short!

Save the results for later use.

res_re3 <- model
rm(model)

Second order random intercept model (uniform distribution for RE standard deviation)

Plan the model

model_plan <-  plan_fp(model.pars = list(exponents = c(0, 1), t.eval = "midpoint"),
                       bth.model = "REINT", ref.std = "STUDY2", nma.ref.trt = "B",
                       bth.prior = list(type = "sd", distr = "unif", min = 0, max = 2),
                       n.chains = 2,
                       n.iter = 6000,
                       n.burnin = 1000,
                       n.thin = 1,
                       descr = "Second order random intercept model (unif[0,2] prior)",
                       descr_s = "FP (2o, REint, U[0,2])"
                       )

Ready the data

# Returns list list contaiing a jags list ready for input to `nma_fit` and a network object
model_input <- nma_pre_proc(grouped_TTE, model_plan)

Fit the model

model  <- nma_fit(model_input = model_input)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 675
##    Unobserved stochastic nodes: 20
##    Total graph size: 10478
## 
## Initializing model

Post processing

Do convergence assessment as above (not shown).

The chains are obviously MUCH too short!

Save the results for later use.

res_re4 <- model
rm(model)

Produce outputs of interest

Start with an object collecting all fits.

all_res <- list(res_re1, res_re2, res_re3, res_re4)

Model comparison

dcompare <- get_fp_comparison(all_res)
cat("__Table__ Model comparison")

Table Model comparison

pander::pandoc.table(dcompare, row.names = FALSE, split.tables = Inf)
Model Order Exponents REINT RE DIC pD meanDev
FP (1o, REint, LN) 1 0 TRUE FALSE 2706.6 23.9 2683
FP (2o, REint, LN) 2 0, 1 TRUE FALSE 2643.6 32.2 2631.2
FP (1o, REint, U[0,2]) 1 0 TRUE FALSE 2712.5 27.2 2684.6
FP (2o, REint, U[0,2]) 2 0, 1 TRUE FALSE 2649.3 39.7 2624.7

Hazard ratio estimates

# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")

  ## Tables: Hazard ratio estimates for each segment

  HR_rev <- get_fp_HR(x = seq(0.1, 24, 0.1), 
                      fit = res_i, 
                      trt.nos = 1:res_i$data.jg$Ntrt, 
                      ref.no = 2, # use this treatment in the list of treatments as reference for the HRs 
                      revert = TRUE # revert to get HRs ref vs other treatments
                      ) 
  
  fig1 <- plot_fp_HR(HR_rev, xlab = "Month", breaks = c(0.25, 0.5, 1, 2))
  fig2 <- plot_fp_HR(HR_rev, xlab = "Month", breaks = c(0.25, 0.5, 1, 2), facet = TRUE) # !! HERE: ADD CIs !!
  
  
  dHRtab <- HR_rev %>%
    mutate(Month = round(x, 1)) %>%        # trick, otherwise equality testing fails to pick out all timepoints in filter step
    filter(Month %in% seq(3, 24, 3)) %>%
    mutate(Comparison = lab, Month, HR = round(median, 3), lCI = round(lCI, 3), uCI = round(uCI, 3)) %>%
    select(Comparison, Month, HR, lCI, uCI)
  
  
  cat("__Figure__ Hazard ratios treatment A vs other treatments\n")
  plot(fig1)
  cat("\n\n")
  cat("__Figure__ Hazard ratios treatment A vs other treatments (multi-panel)\n")
  plot(fig2)
  cat("\n\n")
  cat("__Table__ Hazard ratios treatment A vs other treatments\n")
  pander::pandoc.table(dHRtab, row.names = FALSE, split.tables = Inf)
  cat("\n\n")
  cat("\n\n")
  
  rm(HR_rev)
  rm(dHRtab)
  rm(fig1)
  rm(fig2)
  rm(res_i)
}

First order random intercept model (LN prior for RE Var)

Figure Hazard ratios treatment A vs other treatments

Figure Hazard ratios treatment A vs other treatments (multi-panel)

Table Hazard ratios treatment A vs other treatments

Comparison Month HR lCI uCI
A vs B 3 0.821 0.596 1.097
A vs B 6 0.859 0.632 1.123
A vs B 9 0.884 0.641 1.159
A vs B 12 0.904 0.653 1.182
A vs B 15 0.914 0.66 1.199
A vs B 18 0.923 0.662 1.212
A vs B 21 0.932 0.666 1.226
A vs B 24 0.942 0.67 1.247
A vs C 3 0.912 0.621 1.37
A vs C 6 0.888 0.644 1.25
A vs C 9 0.878 0.647 1.23
A vs C 12 0.87 0.641 1.23
A vs C 15 0.867 0.63 1.236
A vs C 18 0.861 0.619 1.258
A vs C 21 0.859 0.606 1.27
A vs C 24 0.857 0.595 1.284
A vs D 3 0.55 0.356 0.843
A vs D 6 0.658 0.438 0.972
A vs D 9 0.733 0.485 1.083
A vs D 12 0.789 0.519 1.174
A vs D 15 0.829 0.542 1.262
A vs D 18 0.869 0.566 1.325
A vs D 21 0.903 0.581 1.394
A vs D 24 0.932 0.594 1.467
A vs E 3 0.707 0.404 1.14
A vs E 6 0.799 0.493 1.247
A vs E 9 0.856 0.541 1.324
A vs E 12 0.899 0.572 1.412
A vs E 15 0.933 0.591 1.485
A vs E 18 0.961 0.605 1.565
A vs E 21 0.988 0.62 1.656
A vs E 24 1.011 0.628 1.733
A vs F 3 0.681 0.35 1.423
A vs F 6 0.732 0.412 1.351
A vs F 9 0.762 0.446 1.346
A vs F 12 0.78 0.465 1.392
A vs F 15 0.798 0.475 1.438
A vs F 18 0.808 0.479 1.489
A vs F 21 0.822 0.481 1.539
A vs F 24 0.833 0.482 1.578

Second order random intercept model (LN prior for RE Var)

Figure Hazard ratios treatment A vs other treatments

Figure Hazard ratios treatment A vs other treatments (multi-panel)

Table Hazard ratios treatment A vs other treatments

Comparison Month HR lCI uCI
A vs B 3 2.693 1.254 4.393
A vs B 6 2.079 0.947 4.232
A vs B 9 1.889 0.844 4.166
A vs B 12 1.831 0.782 4.214
A vs B 15 1.818 0.764 4.291
A vs B 18 1.845 0.778 4.399
A vs B 21 1.916 0.81 4.519
A vs B 24 2.011 0.831 4.596
A vs C 3 1.519 0.805 2.938
A vs C 6 1.205 0.458 2.588
A vs C 9 1.105 0.374 2.47
A vs C 12 1.078 0.347 2.456
A vs C 15 1.09 0.347 2.436
A vs C 18 1.114 0.367 2.458
A vs C 21 1.154 0.403 2.525
A vs C 24 1.211 0.446 2.637
A vs D 3 1.628 0.545 2.98
A vs D 6 1.34 0.308 2.948
A vs D 9 1.324 0.243 3.278
A vs D 12 1.416 0.221 3.823
A vs D 15 1.586 0.217 4.646
A vs D 18 1.784 0.224 5.644
A vs D 21 2.028 0.238 6.912
A vs D 24 2.346 0.261 8.817
A vs E 3 1.987 0.952 8.13
A vs E 6 1.407 0.442 7.077
A vs E 9 1.277 0.301 7.522
A vs E 12 1.296 0.253 8.782
A vs E 15 1.412 0.237 10.71
A vs E 18 1.592 0.242 13.39
A vs E 21 1.847 0.258 17.11
A vs E 24 2.151 0.291 22.25
A vs F 3 0.763 0.081 1.319
A vs F 6 0.456 0.026 0.963
A vs F 9 0.392 0.018 0.964
A vs F 12 0.394 0.016 1.038
A vs F 15 0.429 0.018 1.182
A vs F 18 0.49 0.023 1.406
A vs F 21 0.574 0.031 1.699
A vs F 24 0.697 0.043 2.214

First order random intercept model (unif[0,2] prior)

Figure Hazard ratios treatment A vs other treatments

Figure Hazard ratios treatment A vs other treatments (multi-panel)

Table Hazard ratios treatment A vs other treatments

Comparison Month HR lCI uCI
A vs B 3 0.934 0.668 6.592
A vs B 6 0.954 0.674 7.2
A vs B 9 0.961 0.659 7.634
A vs B 12 0.966 0.657 7.915
A vs B 15 0.968 0.651 8.088
A vs B 18 0.97 0.652 8.236
A vs B 21 0.976 0.649 8.361
A vs B 24 0.981 0.639 8.489
A vs C 3 0.93 0.549 9.612
A vs C 6 0.919 0.544 9.832
A vs C 9 0.906 0.537 10.17
A vs C 12 0.903 0.534 10.46
A vs C 15 0.9 0.526 10.66
A vs C 18 0.898 0.523 10.75
A vs C 21 0.894 0.516 11.05
A vs C 24 0.891 0.508 11.01
A vs D 3 0.626 0.339 6.098
A vs D 6 0.713 0.423 7.995
A vs D 9 0.78 0.475 9.47
A vs D 12 0.829 0.506 10.75
A vs D 15 0.877 0.517 11.89
A vs D 18 0.916 0.531 12.91
A vs D 21 0.952 0.536 13.75
A vs D 24 0.982 0.538 14.63
A vs E 3 0.807 0.399 3.19
A vs E 6 0.866 0.474 3.821
A vs E 9 0.902 0.488 4.37
A vs E 12 0.94 0.501 4.903
A vs E 15 0.972 0.509 5.39
A vs E 18 1.005 0.51 5.785
A vs E 21 1.029 0.512 6.161
A vs E 24 1.054 0.509 6.423
A vs F 3 0.679 0.254 74.03
A vs F 6 0.74 0.311 83.27
A vs F 9 0.779 0.338 88.36
A vs F 12 0.823 0.354 93.63
A vs F 15 0.847 0.357 97.82
A vs F 18 0.871 0.365 100
A vs F 21 0.886 0.374 103.5
A vs F 24 0.903 0.38 104.9

Second order random intercept model (unif[0,2] prior)

Figure Hazard ratios treatment A vs other treatments

Figure Hazard ratios treatment A vs other treatments (multi-panel)

Table Hazard ratios treatment A vs other treatments

Comparison Month HR lCI uCI
A vs B 3 0.903 0.619 1.238
A vs B 6 0.859 0.479 1.586
A vs B 9 0.833 0.428 1.77
A vs B 12 0.838 0.407 1.864
A vs B 15 0.84 0.402 1.916
A vs B 18 0.857 0.403 1.959
A vs B 21 0.869 0.413 1.957
A vs B 24 0.89 0.424 1.924
A vs C 3 0.553 0.272 0.94
A vs C 6 0.461 0.187 0.929
A vs C 9 0.445 0.158 0.946
A vs C 12 0.449 0.149 1.004
A vs C 15 0.461 0.147 1.065
A vs C 18 0.489 0.151 1.125
A vs C 21 0.523 0.157 1.197
A vs C 24 0.563 0.165 1.315
A vs D 3 0.325 0.097 1.232
A vs D 6 0.371 0.072 2.081
A vs D 9 0.404 0.068 2.82
A vs D 12 0.43 0.07 3.453
A vs D 15 0.455 0.076 4.041
A vs D 18 0.48 0.085 4.542
A vs D 21 0.506 0.097 4.99
A vs D 24 0.532 0.112 5.4
A vs E 3 0.77 0.372 1.666
A vs E 6 0.715 0.224 2.301
A vs E 9 0.732 0.19 2.73
A vs E 12 0.754 0.185 3.024
A vs E 15 0.772 0.194 3.271
A vs E 18 0.822 0.212 3.583
A vs E 21 0.851 0.238 3.881
A vs E 24 0.92 0.278 4.191
A vs F 3 0.067 0.007 2.485
A vs F 6 0.034 0.003 1.467
A vs F 9 0.029 0.002 1.287
A vs F 12 0.03 0.002 1.299
A vs F 15 0.035 0.002 1.415
A vs F 18 0.042 0.003 1.609
A vs F 21 0.054 0.003 1.887
A vs F 24 0.073 0.004 2.333

Survivor function estimates

The NMA baseline estimate from the ref_trt arm from ref_std is used. These are combined with the time-varying hazard-ratio functions from the NMA to obtain the survivor functions for the other interventions.

ref_trt <- "B"
ref_std <- "STUDY2"
hor <- 60

# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")
  
  ## Plots of survivor functions over time ("NMA result"), ref study/arm and timehorizons specified in settings function
  sel_ref <- which(attr(res_i$data.jg, "d_arms")$study == ref_std & attr(res_i$data.jg, "d_arms")$treatment == ref_trt)
  id_ref_std <- attr(res_i$data.jg, "d_arms")$studyn[sel_ref]
  id_ref_arm <- attr(res_i$data.jg, "d_arms")$arm[sel_ref]
    
  S_extrap <- get_fp_S(fit = res_i, 
                       ref.std = id_ref_std, 
                       ref.arm = id_ref_arm, 
                       time = seq(0.1, hor, 0.1))

  fig <- ggplot(data = S_extrap) +        
    geom_line(aes(x = time, y = S, col = treatment, linetype = treatment)) +
    ylim(0, 1) +
    xlab("Month") + ylab("Survival probability") +
    theme_bw() +
    theme(legend.title = element_blank())
  cat("__Figure__ Survivor function estimates (time horizon:", hor, "months) \n")
  plot(fig)
  cat("\n\n")
  
  fig <- ggplot(data = S_extrap) + 
    facet_wrap(~treatment) +
    geom_ribbon(aes(x = time, ymin = lCrI, ymax = uCrI), fill = "lightblue", alpha = 0.8) +
    geom_line(aes(x = time, y = S)) +
    ylim(0, 1) +
    xlab("Month") + ylab("Survival probability") +
    theme_bw()
  cat("__Figure__ Survivor function estimates by treatment (time horizon:", hor, "months) \n")
  plot(fig)
  cat("\n\n")
  
  rm(list = c("S_extrap", "fig"))
  rm(res_i)
}

First order random intercept model (LN prior for RE Var)

Figure Survivor function estimates (time horizon: 60 months)

Figure Survivor function estimates by treatment (time horizon: 60 months)

Second order random intercept model (LN prior for RE Var)

Figure Survivor function estimates (time horizon: 60 months)

Figure Survivor function estimates by treatment (time horizon: 60 months)

First order random intercept model (unif[0,2] prior)

Figure Survivor function estimates (time horizon: 60 months)

Figure Survivor function estimates by treatment (time horizon: 60 months)

Second order random intercept model (unif[0,2] prior)

Figure Survivor function estimates (time horizon: 60 months)

Figure Survivor function estimates by treatment (time horizon: 60 months)

Model fit: observed KM data vs estimated S(t)

For every arm in every study, the study baseline hazard estimate is combined with the corresponding contrast estimate (both from the NMA) to obtain the estimated survivor functions.

hor <- 36
# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")
  
  gof <- get_fp_GoF(fit = res_i, time = seq(0.1, hor, 0.1))

  fig <- ggplot() + 
    geom_line(data = gof %>% filter(type == "nma"), aes(x = time, y = S, col = treatment)) +
    geom_line(data = gof %>% filter(type == "obs"), aes(x = time, y = S, col = treatment), linetype = "dashed") +
    facet_wrap(~study, ncol = 2) +
    ylim(0, 1) + xlim(0, 36) +
    xlab("Month") + ylab("Survival probability") +
    theme_bw() +
    theme(legend.position = "top", legend.title = element_blank())

  cat("__Figure__ Goodness-of-fit: estimated (solid lines) and observed (dashed) survivor functions for each study\n")
  plot(fig)
  cat("\n\n")
  
  rm(list = c("gof", "fig"))
  rm(res_i)
}

First order random intercept model (LN prior for RE Var)

Figure Goodness-of-fit: estimated (solid lines) and observed (dashed) survivor functions for each study

Second order random intercept model (LN prior for RE Var)

Figure Goodness-of-fit: estimated (solid lines) and observed (dashed) survivor functions for each study

First order random intercept model (unif[0,2] prior)

Figure Goodness-of-fit: estimated (solid lines) and observed (dashed) survivor functions for each study

Second order random intercept model (unif[0,2] prior)

Figure Goodness-of-fit: estimated (solid lines) and observed (dashed) survivor functions for each study

Parameter estimates (baseline and contrasts)

# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")

  cest <- get_fp_contrasts(res_i)

  cat("\n\n")
  cat("__Table__ Contrast estimates in fractional polynomial vs network reference\n")
  pander::pandoc.table(cest, row.names = FALSE, split.tables = Inf)
  cat("\n\n")

  rm(cest)
  rm(res_i)
}

First order random intercept model (LN prior for RE Var)

Table Contrast estimates in fractional polynomial vs network reference

Treatment Int.med Int.lCI Int.uCI Slope1.med Slope1.lCI Slope1.uCI
A -0.272 -0.646 0.097 0.063 -0.06 0.185
C -0.209 -0.78 0.36 0.086 -0.103 0.292
D 0.614 0.117 1.097 -0.197 -0.359 -0.015
E 0.27 -0.272 0.855 -0.122 -0.311 0.079
F 0.267 -0.772 1.156 -0.04 -0.367 0.278

Second order random intercept model (LN prior for RE Var)

Table Contrast estimates in fractional polynomial vs network reference

Treatment Int.med Int.lCI Int.uCI Slope1.med Slope1.lCI Slope1.uCI Slope2.med Slope2.lCI Slope2.uCI
A 1.46 0.666 1.642 -0.425 -0.912 -0.052 0.03 0.002 0.063
C 0.521 0.074 0.855 0.106 -0.245 0.391 -0.007 -0.044 0.024
D 0.421 0.057 0.983 0.306 0.017 0.644 -0.062 -0.102 -0.019
E -0.256 -1.05 0.195 0.58 0.239 1.087 -0.079 -0.125 -0.036
F 0.642 -0.08 1.908 1.1 -0.11 1.88 -0.11 -0.195 -0.024

First order random intercept model (unif[0,2] prior)

Table Contrast estimates in fractional polynomial vs network reference

Treatment Int.med Int.lCI Int.uCI Slope1.med Slope1.lCI Slope1.uCI
A -0.094 -0.45 1.756 0.028 -0.142 0.166
C -0.028 -0.803 0.616 0.041 -0.087 0.283
D 0.579 -0.167 1.216 -0.189 -0.462 0.006
E 0.325 -0.301 1.469 -0.107 -0.459 0.158
F 0.45 -2.571 1.622 -0.115 -0.448 0.272

Second order random intercept model (unif[0,2] prior)

Table Contrast estimates in fractional polynomial vs network reference

Treatment Int.med Int.lCI Int.uCI Slope1.med Slope1.lCI Slope1.uCI Slope2.med Slope2.lCI Slope2.uCI
A -0.09 -0.483 0.38 -0.214 -0.6 0.477 0.013 -0.031 0.047
C 0.256 -0.403 1.039 0.349 -0.106 0.717 -0.037 -0.076 0.003
D 0.994 0.298 1.896 -0.082 -0.476 0.358 -0.03 -0.073 0.042
E 0.006 -0.266 0.202 0.18 -0.269 0.744 -0.041 -0.093 0.033
F 1.27 -2.171 3.325 1.615 1.121 2.129 -0.163 -0.212 -0.093

Posterior correlation plots

# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")

  corrs <- get_fp_corrs(res_i)

  for(j in 1:res_i$data.jg$Ntrt){
    cat("\n\n")
    cat("__Table__ Posterior correlations of (multivariate) contrasts for", dimnames(corrs)$treatment[j],"vs reference\n")
    pander::pandoc.table(corrs[j,,], row.names = dimnames(corrs[j,,])[[1]], split.tables = Inf)
    cat("\n\n")
  }
  
  rm(corrs)
  rm(res_i)
}

First order random intercept model (LN prior for RE Var)

Table Posterior correlations of (multivariate) contrasts for B vs reference

  Int Slope1
Int 1 NA
Slope1 NA 1

Table Posterior correlations of (multivariate) contrasts for A vs reference

  Int Slope1
Int 1 -0.69
Slope1 -0.69 1

Table Posterior correlations of (multivariate) contrasts for C vs reference

  Int Slope1
Int 1 -0.847
Slope1 -0.847 1

Table Posterior correlations of (multivariate) contrasts for D vs reference

  Int Slope1
Int 1 -0.752
Slope1 -0.752 1

Table Posterior correlations of (multivariate) contrasts for E vs reference

  Int Slope1
Int 1 -0.753
Slope1 -0.753 1

Table Posterior correlations of (multivariate) contrasts for F vs reference

  Int Slope1
Int 1 -0.856
Slope1 -0.856 1

Second order random intercept model (LN prior for RE Var)

Table Posterior correlations of (multivariate) contrasts for B vs reference

  Int Slope1 Slope2
Int 1 NA NA
Slope1 NA 1 NA
Slope2 NA NA 1

Table Posterior correlations of (multivariate) contrasts for A vs reference

  Int Slope1 Slope2
Int 1 0.05 -0.073
Slope1 0.05 1 -0.941
Slope2 -0.073 -0.941 1

Table Posterior correlations of (multivariate) contrasts for C vs reference

  Int Slope1 Slope2
Int 1 -0.014 -0.087
Slope1 -0.014 1 -0.822
Slope2 -0.087 -0.822 1

Table Posterior correlations of (multivariate) contrasts for D vs reference

  Int Slope1 Slope2
Int 1 0.843 -0.533
Slope1 0.843 1 -0.733
Slope2 -0.533 -0.733 1

Table Posterior correlations of (multivariate) contrasts for E vs reference

  Int Slope1 Slope2
Int 1 -0.446 0.617
Slope1 -0.446 1 -0.716
Slope2 0.617 -0.716 1

Table Posterior correlations of (multivariate) contrasts for F vs reference

  Int Slope1 Slope2
Int 1 0.16 -0.217
Slope1 0.16 1 -0.925
Slope2 -0.217 -0.925 1

First order random intercept model (unif[0,2] prior)

Table Posterior correlations of (multivariate) contrasts for B vs reference

  Int Slope1
Int 1 NA
Slope1 NA 1

Table Posterior correlations of (multivariate) contrasts for A vs reference

  Int Slope1
Int 1 -0.015
Slope1 -0.015 1

Table Posterior correlations of (multivariate) contrasts for C vs reference

  Int Slope1
Int 1 -0.202
Slope1 -0.202 1

Table Posterior correlations of (multivariate) contrasts for D vs reference

  Int Slope1
Int 1 -0.644
Slope1 -0.644 1

Table Posterior correlations of (multivariate) contrasts for E vs reference

  Int Slope1
Int 1 -0.659
Slope1 -0.659 1

Table Posterior correlations of (multivariate) contrasts for F vs reference

  Int Slope1
Int 1 -0.483
Slope1 -0.483 1

Second order random intercept model (unif[0,2] prior)

Table Posterior correlations of (multivariate) contrasts for B vs reference

  Int Slope1 Slope2
Int 1 NA NA
Slope1 NA 1 NA
Slope2 NA NA 1

Table Posterior correlations of (multivariate) contrasts for A vs reference

  Int Slope1 Slope2
Int 1 -0.853 0.82
Slope1 -0.853 1 -0.979
Slope2 0.82 -0.979 1

Table Posterior correlations of (multivariate) contrasts for C vs reference

  Int Slope1 Slope2
Int 1 -0.768 0.581
Slope1 -0.768 1 -0.838
Slope2 0.581 -0.838 1

Table Posterior correlations of (multivariate) contrasts for D vs reference

  Int Slope1 Slope2
Int 1 0.958 -0.908
Slope1 0.958 1 -0.94
Slope2 -0.908 -0.94 1

Table Posterior correlations of (multivariate) contrasts for E vs reference

  Int Slope1 Slope2
Int 1 0.538 -0.515
Slope1 0.538 1 -0.936
Slope2 -0.515 -0.936 1

Table Posterior correlations of (multivariate) contrasts for F vs reference

  Int Slope1 Slope2
Int 1 -0.38 -0.001
Slope1 -0.38 1 -0.764
Slope2 -0.001 -0.764 1

Appendix

# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("## ", title, "  \n\n")

  jginfo <- get_jags_info(res_i, include.comments = TRUE)
  cat("```\n", jginfo, "\n```\n\n")
  
  rm(jginfo)
  rm(out)
}

First order random intercept model (LN prior for RE Var)

 ##############################################
# DATA                                       #
##############################################
list(
 Na  =  c(2, 2, 2, 2, 2, 2, 2) ,
 Nobs  =  675 ,
 Ns  =  7 ,
 Ntrt  =  6 ,
 P1  =  0 ,
 a  =  c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2) ,
 dt  =  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1) ,
 feprior_mean  =  c(0, 0) ,
 feprior_prec  = structure(.Data =  c(1e-04, 0, 0, 1e-04) ,.Dim =  c(2, 2) ) ,
 n  =  c(182, 177, 172, 164, 155, 148, 143, 138, 132, 128, 125, 120, 117, 110, 108, 106, 101, 98, 96, 90, 87, 83, 79, 75, 70, 66, 59, 56, 53, 50, 47, 45, 39, 37, 34, 32, 31, 26, 25, 24, 23, 21, 19, 16, 14, 13, 11, 10, 8, 6, 5, 3, 1, 183, 179, 172, 159, 152, 147, 139, 133, 131, 124, 119, 115, 110, 104, 101, 95, 93, 92, 88, 83, 80, 76, 73, 68, 62, 59, 53, 50, 43, 40, 37, 35, 33, 31, 30, 29, 26, 24, 22, 19, 18, 18, 18, 18, 17, 17, 15, 14, 12, 10, 8, 4, 2, 322, 319, 309, 301, 288, 277, 262, 255, 248, 238, 229, 
226, 216, 206, 200, 195, 189, 183, 176, 171, 163, 160, 150, 148, 142, 136, 130, 127, 121, 116, 113, 112, 108, 105, 101, 101, 96, 95, 91, 89, 86, 84, 327, 324, 318, 307, 297, 283, 275, 263, 256, 248, 241, 235, 225, 217, 212, 201, 194, 184, 178, 169, 161, 156, 149, 142, 134, 131, 125, 117, 115, 107, 101, 96, 93, 89, 83, 79, 76, 74, 70, 68, 65, 63, 375, 374, 368, 360, 352, 339, 326, 319, 315, 308, 295, 289, 283, 271, 262, 259, 247, 238, 229, 225, 213, 203, 197, 190, 180, 161, 135, 121, 98, 78, 61, 55, 
52, 45, 39, 35, 375, 372, 356, 338, 319, 308, 295, 282, 276, 269, 260, 251, 242, 228, 221, 210, 203, 200, 187, 180, 172, 168, 161, 155, 149, 132, 115, 104, 82, 66, 53, 48, 43, 39, 34, 30, 181, 181, 180, 176, 172, 169, 167, 158, 153, 145, 143, 141, 137, 133, 129, 126, 123, 121, 118, 113, 112, 108, 105, 101, 99, 99, 98, 95, 92, 90, 88, 87, 86, 83, 82, 82, 78, 78, 78, 78, 181, 179, 176, 175, 173, 172, 169, 166, 165, 164, 159, 155, 150, 144, 140, 135, 132, 132, 132, 129, 127, 126, 123, 121, 116, 109, 
107, 105, 105, 104, 101, 97, 94, 93, 93, 90, 557, 552, 541, 532, 521, 505, 487, 473, 462, 452, 436, 423, 409, 388, 377, 368, 357, 348, 338, 328, 320, 312, 305, 295, 280, 273, 262, 255, 243, 235, 230, 226, 219, 215, 208, 204, 197, 184, 171, 156, 147, 147, 114, 114, 99, 99, 80, 77, 57, 48, 38, 27, 21, 20, 17, 14, 8, 6, 6, 1, 1, 553, 544, 530, 518, 501, 480, 468, 448, 435, 420, 406, 387, 377, 366, 354, 348, 339, 329, 321, 315, 304, 295, 283, 280, 272, 263, 254, 247, 237, 230, 224, 218, 212, 203, 192, 
187, 180, 174, 162, 149, 138, 131, 120, 111, 100, 92, 80, 63, 56, 48, 38, 31, 20, 17, 14, 11, 9, 7, 5, 5, 1, 189, 184, 178, 166, 156, 144, 138, 131, 123, 117, 111, 106, 103, 98, 93, 88, 85, 79, 74, 68, 64, 59, 55, 50, 46, 43, 40, 37, 35, 31, 28, 26, 25, 20, 19, 15, 15, 12, 10, 9, 9, 3, 3, 3, 3, 188, 179, 175, 167, 159, 151, 141, 136, 131, 127, 123, 117, 112, 106, 101, 96, 94, 88, 84, 77, 70, 64, 57, 52, 44, 39, 39, 30, 30, 26, 23, 21, 20, 16, 16, 12, 11, 11, 10, 8, 7, 4, 4, 2, 2, 2, 363, 360, 344, 
332, 314, 305, 286, 269, 261, 248, 236, 226, 220, 214, 199, 197, 188, 180, 177, 168, 162, 155, 151, 150, 148, 144, 140, 136, 128, 122, 118, 116, 113, 113, 109, 105, 98, 94, 87, 83, 75, 69, 64, 57, 53, 49, 46, 40, 37, 31, 25, 22, 18, 14, 10, 9, 7, 6, 3, 2, 1, 369, 366, 358, 347, 333, 322, 314, 302, 294, 277, 267, 254, 240, 236, 225, 218, 210, 197, 190, 179, 179, 173, 166, 163, 160, 154, 152, 148, 145, 143, 139, 136, 127, 123, 119, 118, 116, 112, 109, 106, 102, 99, 94, 85, 74, 66, 58, 50, 42, 35, 31, 
28, 25, 20, 17, 15, 13, 11, 6, 3, 2, 1, 1) ,
 r  =  c(2, 2, 4, 6, 4, 5, 4, 4, 4, 2, 5, 1, 7, 1, 1, 3, 1, 0, 3, 2, 1, 3, 2, 2, 2, 3, 1, 0, 0, 0, 1, 4, 1, 1, 1, 0, 4, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 4, 10, 4, 2, 6, 4, 0, 5, 3, 3, 3, 5, 1, 5, 1, 0, 3, 3, 3, 0, 0, 2, 2, 0, 3, 1, 3, 0, 1, 2, 0, 0, 0, 0, 2, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 10, 7, 13, 10, 15, 7, 7, 9, 9, 3, 9, 10, 5, 5, 6, 5, 7, 5, 7, 3, 10, 1, 6, 6, 5, 3, 5, 5, 3, 0, 4, 3, 3, 0, 4, 1, 4, 1, 3, 2, 3, 2, 4, 9, 8, 12, 7, 11, 5, 7, 5, 4, 8, 7, 2, 10, 
6, 8, 4, 8, 6, 3, 5, 5, 6, 2, 5, 5, 1, 6, 5, 3, 1, 2, 5, 2, 1, 1, 2, 1, 1, 0, 0, 0, 5, 6, 7, 11, 12, 7, 3, 6, 12, 5, 6, 11, 8, 1, 11, 8, 8, 3, 10, 9, 4, 5, 9, 2, 6, 3, 2, 1, 2, 1, 0, 3, 0, 0, 0, 0, 13, 14, 15, 7, 10, 12, 5, 6, 7, 8, 8, 14, 7, 10, 7, 3, 13, 5, 5, 2, 4, 4, 3, 3, 1, 2, 5, 1, 1, 0, 2, 0, 0, 0, 0, 0, 1, 4, 4, 3, 2, 9, 5, 8, 2, 2, 4, 4, 4, 3, 3, 2, 3, 5, 1, 4, 3, 4, 2, 0, 1, 3, 3, 2, 2, 1, 1, 3, 1, 0, 4, 0, 0, 0, 0, 2, 3, 1, 2, 1, 3, 3, 1, 1, 5, 4, 5, 6, 4, 5, 3, 0, 0, 3, 2, 1, 3, 2, 5, 
7, 2, 2, 0, 1, 3, 4, 3, 1, 0, 3, 2, 3, 8, 6, 9, 14, 15, 10, 9, 8, 13, 11, 12, 21, 9, 8, 11, 9, 10, 9, 8, 6, 4, 6, 12, 4, 8, 4, 10, 6, 3, 1, 4, 2, 3, 2, 5, 5, 4, 6, 1, 0, 0, 0, 15, 0, 0, 3, 4, 1, 0, 0, 1, 0, 2, 0, 4, 0, 0, 0, 0, 0, 4, 8, 6, 12, 19, 9, 17, 11, 15, 12, 18, 10, 11, 11, 5, 9, 9, 6, 4, 10, 7, 9, 0, 5, 6, 6, 5, 7, 5, 5, 4, 3, 8, 9, 2, 6, 0, 5, 7, 3, 0, 3, 3, 6, 0, 3, 5, 1, 0, 0, 0, 2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 10, 7, 8, 3, 6, 4, 5, 4, 4, 0, 3, 3, 1, 1, 4, 3, 2, 2, 2, 2, 1, 1, 2, 
2, 0, 1, 3, 2, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 6, 2, 5, 5, 3, 8, 2, 4, 1, 2, 4, 3, 4, 4, 0, 2, 0, 3, 4, 1, 1, 0, 2, 3, 3, 0, 0, 0, 2, 0, 1, 1, 0, 0, 2, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 3, 15, 11, 18, 8, 19, 17, 8, 13, 11, 10, 6, 6, 15, 2, 8, 8, 3, 9, 6, 7, 4, 1, 2, 4, 4, 3, 8, 6, 4, 2, 3, 0, 4, 4, 7, 2, 4, 2, 4, 4, 3, 4, 1, 0, 0, 3, 0, 3, 2, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 3, 7, 11, 14, 10, 8, 12, 8, 17, 10, 13, 14, 4, 11, 7, 8, 13, 7, 11, 0, 6, 7, 3, 3, 6, 2, 3, 3, 2, 4, 3, 9, 4, 3, 
1, 2, 3, 1, 2, 2, 2, 4, 3, 5, 1, 2, 2, 2, 4, 0, 0, 0, 1, 0, 1, 1, 0, 3, 2, 0, 0, 0, 0) ,
 reprior_meanlog  =  -2.71 ,
 reprior_sdlog  =  1.74 ,
 s  =  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2) ,
 t  = structure(.Data =  c(1, 3, 4, 5, 4, 5, 1, 4, 3, 6, 2, 1, 2, 3) ,.Dim =  c(7, 2) ) ,
 time  =  c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 
34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 
26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 
35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 
12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 
37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 
18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 
58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 61.5, 62.5) 
)

##############################################
# MODEL                                      #
##############################################
# Fractional polynomial, 1st order, random effects model: a RE is put only on the scale parameter 
#                                                          i.e. on the intercept in the frac poly
#                                                         a log-normal prior is used for the RE Var
#                                                           similar to the informative Turner priors
#                                                           using the standard parameterization via meanlog, sdlog
# -------------------------------------------------------------------------------
# Data: grouped survival data, binomial likelihood, linear predictor on log-hazard
#         Nobs     number of observations
#         n[i]     patients at risk in interval i
#         r[i]     events during interval i
#         time[i]  mid-point of interval i
#         dt[i]    length of interval i
#         Ns       number of studies
#         Na[j]    number of arms in study j
#         Ntrt     number of treatments
#         s[i]     study number for obs i
#         a[i]     arm number (within study) for obs i
#         t[i,j]   treatment in study i arm j
#         P1       exponent of the time varying term in the fractional polynomial
#         feprior_mean[1:2]         prior mean (for contrasts d and baselines mu)
#         feprior_prec[1:2, 1:2]    prior precision (for d and mu)
#         reprior_meanlog           Y~LN(meanlog, sdlog) then log(Y)~N(meanlog, sdlog^2)
#         reprior_sdlog           
# -------------------------------------------------------------------------------

model{

## Sampling model (likelihood)
for (i in 1:Nobs){
  time1[i] <- (equals(P1, 0) * log(time[i]) + (1 - equals(P1, 0)) * pow(time[i], P1)) 

  # likelihood: digitized KM curves
  r[i] ~ dbin(p[i], n[i])
  p[i] <- 1 - exp(-h[i] * dt[i])  # cumulative hazard over interval [t,t+dt] expressed as deaths per person-month

  # fractional polynomial
  log(h[i]) <- Beta[s[i], a[i], 1] + Beta[s[i], a[i], 2] * time1[i]
  }


## Arm level parameters = study effect + trt effect (RE model, consistency eq for pop pars)
for (i in 1:Ns){
  w[i, 1] <- 0
  delta[i, 1] <- 0

  for (j in 1:Na[i]){
    Beta[i, j, 1] <- mu[i, 1] + delta[i, j]
    Beta[i, j, 2] <- mu[i, 2] + d[t[i, j], 2] - d[t[i, 1], 2]
    }

  for (j in 2:Na[i]){
    delta[i, j] ~ dnorm(md[i, j], taud[i, j])
    md[i, j] <- d[t[i, j], 1] - d[t[i, 1], 1] + sw[i, j]
    w[i, j] <- (delta[i, j] - d[t[i, j], 1] + d[t[i, 1], 1])
    sw[i, j] <- sum(w[i, 1:(j - 1)]) / (j - 1) 
    taud[i, j] <- tau * 2 * (j - 1) / j 
    }
    
  }


## Priors
for (j in 1:Ns){
  mu[j, 1:2] ~ dmnorm(feprior_mean[1:2], feprior_prec[,]) 
  }

d[1, 1] <- 0
d[1, 2] <- 0
for (i in 2:Ntrt){
  d[i, 1:2] ~ dmnorm(feprior_mean[1:2], feprior_prec[,]) 
  }

reprior_preclog <- 1 / (reprior_sdlog * reprior_sdlog)
logVar ~ dnorm(reprior_meanlog, reprior_preclog) 
sd <- sqrt(exp(logVar))
tau <- 1 / (sd * sd)

} # end of model
##############################################

Second order random intercept model (LN prior for RE Var)

 ##############################################
# DATA                                       #
##############################################
list(
 Na  =  c(2, 2, 2, 2, 2, 2, 2) ,
 Nobs  =  675 ,
 Ns  =  7 ,
 Ntrt  =  6 ,
 P1  =  0 ,
 P2  =  1 ,
 a  =  c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2) ,
 dt  =  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1) ,
 feprior_mean  =  c(0, 0, 0) ,
 feprior_prec  = structure(.Data =  c(1e-04, 0, 0, 0, 1e-04, 0, 0, 0, 1e-04) ,.Dim =  c(3, 3) ) ,
 n  =  c(182, 177, 172, 164, 155, 148, 143, 138, 132, 128, 125, 120, 117, 110, 108, 106, 101, 98, 96, 90, 87, 83, 79, 75, 70, 66, 59, 56, 53, 50, 47, 45, 39, 37, 34, 32, 31, 26, 25, 24, 23, 21, 19, 16, 14, 13, 11, 10, 8, 6, 5, 3, 1, 183, 179, 172, 159, 152, 147, 139, 133, 131, 124, 119, 115, 110, 104, 101, 95, 93, 92, 88, 83, 80, 76, 73, 68, 62, 59, 53, 50, 43, 40, 37, 35, 33, 31, 30, 29, 26, 24, 22, 19, 18, 18, 18, 18, 17, 17, 15, 14, 12, 10, 8, 4, 2, 322, 319, 309, 301, 288, 277, 262, 255, 248, 238, 229, 
226, 216, 206, 200, 195, 189, 183, 176, 171, 163, 160, 150, 148, 142, 136, 130, 127, 121, 116, 113, 112, 108, 105, 101, 101, 96, 95, 91, 89, 86, 84, 327, 324, 318, 307, 297, 283, 275, 263, 256, 248, 241, 235, 225, 217, 212, 201, 194, 184, 178, 169, 161, 156, 149, 142, 134, 131, 125, 117, 115, 107, 101, 96, 93, 89, 83, 79, 76, 74, 70, 68, 65, 63, 375, 374, 368, 360, 352, 339, 326, 319, 315, 308, 295, 289, 283, 271, 262, 259, 247, 238, 229, 225, 213, 203, 197, 190, 180, 161, 135, 121, 98, 78, 61, 55, 
52, 45, 39, 35, 375, 372, 356, 338, 319, 308, 295, 282, 276, 269, 260, 251, 242, 228, 221, 210, 203, 200, 187, 180, 172, 168, 161, 155, 149, 132, 115, 104, 82, 66, 53, 48, 43, 39, 34, 30, 181, 181, 180, 176, 172, 169, 167, 158, 153, 145, 143, 141, 137, 133, 129, 126, 123, 121, 118, 113, 112, 108, 105, 101, 99, 99, 98, 95, 92, 90, 88, 87, 86, 83, 82, 82, 78, 78, 78, 78, 181, 179, 176, 175, 173, 172, 169, 166, 165, 164, 159, 155, 150, 144, 140, 135, 132, 132, 132, 129, 127, 126, 123, 121, 116, 109, 
107, 105, 105, 104, 101, 97, 94, 93, 93, 90, 557, 552, 541, 532, 521, 505, 487, 473, 462, 452, 436, 423, 409, 388, 377, 368, 357, 348, 338, 328, 320, 312, 305, 295, 280, 273, 262, 255, 243, 235, 230, 226, 219, 215, 208, 204, 197, 184, 171, 156, 147, 147, 114, 114, 99, 99, 80, 77, 57, 48, 38, 27, 21, 20, 17, 14, 8, 6, 6, 1, 1, 553, 544, 530, 518, 501, 480, 468, 448, 435, 420, 406, 387, 377, 366, 354, 348, 339, 329, 321, 315, 304, 295, 283, 280, 272, 263, 254, 247, 237, 230, 224, 218, 212, 203, 192, 
187, 180, 174, 162, 149, 138, 131, 120, 111, 100, 92, 80, 63, 56, 48, 38, 31, 20, 17, 14, 11, 9, 7, 5, 5, 1, 189, 184, 178, 166, 156, 144, 138, 131, 123, 117, 111, 106, 103, 98, 93, 88, 85, 79, 74, 68, 64, 59, 55, 50, 46, 43, 40, 37, 35, 31, 28, 26, 25, 20, 19, 15, 15, 12, 10, 9, 9, 3, 3, 3, 3, 188, 179, 175, 167, 159, 151, 141, 136, 131, 127, 123, 117, 112, 106, 101, 96, 94, 88, 84, 77, 70, 64, 57, 52, 44, 39, 39, 30, 30, 26, 23, 21, 20, 16, 16, 12, 11, 11, 10, 8, 7, 4, 4, 2, 2, 2, 363, 360, 344, 
332, 314, 305, 286, 269, 261, 248, 236, 226, 220, 214, 199, 197, 188, 180, 177, 168, 162, 155, 151, 150, 148, 144, 140, 136, 128, 122, 118, 116, 113, 113, 109, 105, 98, 94, 87, 83, 75, 69, 64, 57, 53, 49, 46, 40, 37, 31, 25, 22, 18, 14, 10, 9, 7, 6, 3, 2, 1, 369, 366, 358, 347, 333, 322, 314, 302, 294, 277, 267, 254, 240, 236, 225, 218, 210, 197, 190, 179, 179, 173, 166, 163, 160, 154, 152, 148, 145, 143, 139, 136, 127, 123, 119, 118, 116, 112, 109, 106, 102, 99, 94, 85, 74, 66, 58, 50, 42, 35, 31, 
28, 25, 20, 17, 15, 13, 11, 6, 3, 2, 1, 1) ,
 r  =  c(2, 2, 4, 6, 4, 5, 4, 4, 4, 2, 5, 1, 7, 1, 1, 3, 1, 0, 3, 2, 1, 3, 2, 2, 2, 3, 1, 0, 0, 0, 1, 4, 1, 1, 1, 0, 4, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 4, 10, 4, 2, 6, 4, 0, 5, 3, 3, 3, 5, 1, 5, 1, 0, 3, 3, 3, 0, 0, 2, 2, 0, 3, 1, 3, 0, 1, 2, 0, 0, 0, 0, 2, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 10, 7, 13, 10, 15, 7, 7, 9, 9, 3, 9, 10, 5, 5, 6, 5, 7, 5, 7, 3, 10, 1, 6, 6, 5, 3, 5, 5, 3, 0, 4, 3, 3, 0, 4, 1, 4, 1, 3, 2, 3, 2, 4, 9, 8, 12, 7, 11, 5, 7, 5, 4, 8, 7, 2, 10, 
6, 8, 4, 8, 6, 3, 5, 5, 6, 2, 5, 5, 1, 6, 5, 3, 1, 2, 5, 2, 1, 1, 2, 1, 1, 0, 0, 0, 5, 6, 7, 11, 12, 7, 3, 6, 12, 5, 6, 11, 8, 1, 11, 8, 8, 3, 10, 9, 4, 5, 9, 2, 6, 3, 2, 1, 2, 1, 0, 3, 0, 0, 0, 0, 13, 14, 15, 7, 10, 12, 5, 6, 7, 8, 8, 14, 7, 10, 7, 3, 13, 5, 5, 2, 4, 4, 3, 3, 1, 2, 5, 1, 1, 0, 2, 0, 0, 0, 0, 0, 1, 4, 4, 3, 2, 9, 5, 8, 2, 2, 4, 4, 4, 3, 3, 2, 3, 5, 1, 4, 3, 4, 2, 0, 1, 3, 3, 2, 2, 1, 1, 3, 1, 0, 4, 0, 0, 0, 0, 2, 3, 1, 2, 1, 3, 3, 1, 1, 5, 4, 5, 6, 4, 5, 3, 0, 0, 3, 2, 1, 3, 2, 5, 
7, 2, 2, 0, 1, 3, 4, 3, 1, 0, 3, 2, 3, 8, 6, 9, 14, 15, 10, 9, 8, 13, 11, 12, 21, 9, 8, 11, 9, 10, 9, 8, 6, 4, 6, 12, 4, 8, 4, 10, 6, 3, 1, 4, 2, 3, 2, 5, 5, 4, 6, 1, 0, 0, 0, 15, 0, 0, 3, 4, 1, 0, 0, 1, 0, 2, 0, 4, 0, 0, 0, 0, 0, 4, 8, 6, 12, 19, 9, 17, 11, 15, 12, 18, 10, 11, 11, 5, 9, 9, 6, 4, 10, 7, 9, 0, 5, 6, 6, 5, 7, 5, 5, 4, 3, 8, 9, 2, 6, 0, 5, 7, 3, 0, 3, 3, 6, 0, 3, 5, 1, 0, 0, 0, 2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 10, 7, 8, 3, 6, 4, 5, 4, 4, 0, 3, 3, 1, 1, 4, 3, 2, 2, 2, 2, 1, 1, 2, 
2, 0, 1, 3, 2, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 6, 2, 5, 5, 3, 8, 2, 4, 1, 2, 4, 3, 4, 4, 0, 2, 0, 3, 4, 1, 1, 0, 2, 3, 3, 0, 0, 0, 2, 0, 1, 1, 0, 0, 2, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 3, 15, 11, 18, 8, 19, 17, 8, 13, 11, 10, 6, 6, 15, 2, 8, 8, 3, 9, 6, 7, 4, 1, 2, 4, 4, 3, 8, 6, 4, 2, 3, 0, 4, 4, 7, 2, 4, 2, 4, 4, 3, 4, 1, 0, 0, 3, 0, 3, 2, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 3, 7, 11, 14, 10, 8, 12, 8, 17, 10, 13, 14, 4, 11, 7, 8, 13, 7, 11, 0, 6, 7, 3, 3, 6, 2, 3, 3, 2, 4, 3, 9, 4, 3, 
1, 2, 3, 1, 2, 2, 2, 4, 3, 5, 1, 2, 2, 2, 4, 0, 0, 0, 1, 0, 1, 1, 0, 3, 2, 0, 0, 0, 0) ,
 reprior_meanlog  =  -2.71 ,
 reprior_sdlog  =  1.74 ,
 s  =  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2) ,
 t  = structure(.Data =  c(1, 3, 4, 5, 4, 5, 1, 4, 3, 6, 2, 1, 2, 3) ,.Dim =  c(7, 2) ) ,
 time  =  c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 
34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 
26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 
35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 
12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 
37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 
18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 
58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 61.5, 62.5) 
)

##############################################
# MODEL                                      #
##############################################
# Fractional polynomial, 2nd order, random effects model: a RE is put only on the scale parameter 
#                                                          i.e. on the intercept in the frac poly
#                                                         a log-normal prior is used for the RE Var
#                                                           similar to the informative Turner priors
#                                                           using the standard parameterization via meanlog, sdlog
# -------------------------------------------------------------------------------
# Data: grouped survival data, binomial likelihood, linear predictor on log-hazard
#         Nobs     number of observations
#         n[i]     patients at risk in interval i
#         r[i]     events during interval i
#         time[i]  mid-point of interval i
#         dt[i]    length of interval i
#         Ns       number of studies
#         Na[j]    number of arms in study j
#         Ntrt     number of treatments
#         s[i]     study number for obs i
#         a[i]     arm number (within study) for obs i
#         t[i,j]   treatment in study i arm j
#         P1       exponent of the time varying term in the fractional polynomial
#         feprior_mean[1:3]         prior mean (for contrasts d and baselines mu)
#         feprior_prec[,]           prior precision (for d and mu)
#         reprior_meanlog           Y~LN(meanlog, sdlog) then log(Y)~N(meanlog, sdlog^2)
#         reprior_sdlog           
# -------------------------------------------------------------------------------

model{

## Sampling model (likelihood)
for (i in 1:Nobs){
  time1[i] <- (equals(P1, 0) * log(time[i]) + (1 - equals(P1, 0)) * pow(time[i], P1)) 
  time2[i] <- ((1 - equals(P2, P1)) * (equals(P2, 0) * log(time[i]) + (1 - equals(P2, 0)) * pow(time[i], P2)) + 
                  equals(P2, P1) * (equals(P2, 0) * log(time[i]) * log(time[i]) + (1 - equals(P2, 0)) * pow(time[i], P2) * log(time[i])))


  # likelihood: digitized KM curves
  r[i] ~ dbin(p[i], n[i])
  p[i] <- 1 - exp(-h[i] * dt[i])  # cumulative hazard over interval [t,t+dt] expressed as deaths per person-month

  # fractional polynomial
  log(h[i]) <- Beta[s[i], a[i], 1] + Beta[s[i], a[i], 2] * time1[i] + Beta[s[i], a[i], 3] * time2[i]
  }


  ## Arm level parameters = study effect + trt effect (RE model, consistency eq for pop pars)
  for (i in 1:Ns){
    w[i, 1] <- 0
    delta[i, 1] <- 0

    for (j in 1:Na[i]){
      Beta[i, j, 1] <- mu[i, 1] + delta[i, j]
      Beta[i, j, 2] <- mu[i, 2] + d[t[i, j], 2] - d[t[i, 1], 2]
      Beta[i, j, 3] <- mu[i, 3] + d[t[i, j], 3] - d[t[i, 1], 3]
      }

    for (j in 2:Na[i]){
      delta[i, j] ~ dnorm(md[i, j], taud[i, j])
      md[i, j] <- d[t[i, j], 1] - d[t[i, 1], 1] + sw[i, j]
      w[i, j] <- (delta[i, j] - d[t[i, j], 1] + d[t[i, 1], 1])
      sw[i, j] <- sum(w[i, 1:(j - 1)]) / (j - 1) 
      taud[i, j] <- tau * 2 * (j - 1) / j 
      }
    
    }


     
## Priors
for (j in 1:Ns){
  mu[j, 1:3] ~ dmnorm(feprior_mean[1:3], feprior_prec[,]) 
  }

d[1, 1] <- 0
d[1, 2] <- 0
d[1, 3] <- 0
for (i in 2:Ntrt){
  d[i, 1:3] ~ dmnorm(feprior_mean[1:3], feprior_prec[,]) 
  }


reprior_preclog <- 1 / (reprior_sdlog * reprior_sdlog)
logVar ~ dnorm(reprior_meanlog, reprior_preclog) 
sd <- sqrt(exp(logVar))
tau <- 1 / (sd * sd)

} # end of model
##############################################

First order random intercept model (unif[0,2] prior)

 ##############################################
# DATA                                       #
##############################################
list(
 Na  =  c(2, 2, 2, 2, 2, 2, 2) ,
 Nobs  =  675 ,
 Ns  =  7 ,
 Ntrt  =  6 ,
 P1  =  0 ,
 a  =  c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2) ,
 dt  =  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1) ,
 feprior_mean  =  c(0, 0) ,
 feprior_prec  = structure(.Data =  c(1e-04, 0, 0, 1e-04) ,.Dim =  c(2, 2) ) ,
 n  =  c(182, 177, 172, 164, 155, 148, 143, 138, 132, 128, 125, 120, 117, 110, 108, 106, 101, 98, 96, 90, 87, 83, 79, 75, 70, 66, 59, 56, 53, 50, 47, 45, 39, 37, 34, 32, 31, 26, 25, 24, 23, 21, 19, 16, 14, 13, 11, 10, 8, 6, 5, 3, 1, 183, 179, 172, 159, 152, 147, 139, 133, 131, 124, 119, 115, 110, 104, 101, 95, 93, 92, 88, 83, 80, 76, 73, 68, 62, 59, 53, 50, 43, 40, 37, 35, 33, 31, 30, 29, 26, 24, 22, 19, 18, 18, 18, 18, 17, 17, 15, 14, 12, 10, 8, 4, 2, 322, 319, 309, 301, 288, 277, 262, 255, 248, 238, 229, 
226, 216, 206, 200, 195, 189, 183, 176, 171, 163, 160, 150, 148, 142, 136, 130, 127, 121, 116, 113, 112, 108, 105, 101, 101, 96, 95, 91, 89, 86, 84, 327, 324, 318, 307, 297, 283, 275, 263, 256, 248, 241, 235, 225, 217, 212, 201, 194, 184, 178, 169, 161, 156, 149, 142, 134, 131, 125, 117, 115, 107, 101, 96, 93, 89, 83, 79, 76, 74, 70, 68, 65, 63, 375, 374, 368, 360, 352, 339, 326, 319, 315, 308, 295, 289, 283, 271, 262, 259, 247, 238, 229, 225, 213, 203, 197, 190, 180, 161, 135, 121, 98, 78, 61, 55, 
52, 45, 39, 35, 375, 372, 356, 338, 319, 308, 295, 282, 276, 269, 260, 251, 242, 228, 221, 210, 203, 200, 187, 180, 172, 168, 161, 155, 149, 132, 115, 104, 82, 66, 53, 48, 43, 39, 34, 30, 181, 181, 180, 176, 172, 169, 167, 158, 153, 145, 143, 141, 137, 133, 129, 126, 123, 121, 118, 113, 112, 108, 105, 101, 99, 99, 98, 95, 92, 90, 88, 87, 86, 83, 82, 82, 78, 78, 78, 78, 181, 179, 176, 175, 173, 172, 169, 166, 165, 164, 159, 155, 150, 144, 140, 135, 132, 132, 132, 129, 127, 126, 123, 121, 116, 109, 
107, 105, 105, 104, 101, 97, 94, 93, 93, 90, 557, 552, 541, 532, 521, 505, 487, 473, 462, 452, 436, 423, 409, 388, 377, 368, 357, 348, 338, 328, 320, 312, 305, 295, 280, 273, 262, 255, 243, 235, 230, 226, 219, 215, 208, 204, 197, 184, 171, 156, 147, 147, 114, 114, 99, 99, 80, 77, 57, 48, 38, 27, 21, 20, 17, 14, 8, 6, 6, 1, 1, 553, 544, 530, 518, 501, 480, 468, 448, 435, 420, 406, 387, 377, 366, 354, 348, 339, 329, 321, 315, 304, 295, 283, 280, 272, 263, 254, 247, 237, 230, 224, 218, 212, 203, 192, 
187, 180, 174, 162, 149, 138, 131, 120, 111, 100, 92, 80, 63, 56, 48, 38, 31, 20, 17, 14, 11, 9, 7, 5, 5, 1, 189, 184, 178, 166, 156, 144, 138, 131, 123, 117, 111, 106, 103, 98, 93, 88, 85, 79, 74, 68, 64, 59, 55, 50, 46, 43, 40, 37, 35, 31, 28, 26, 25, 20, 19, 15, 15, 12, 10, 9, 9, 3, 3, 3, 3, 188, 179, 175, 167, 159, 151, 141, 136, 131, 127, 123, 117, 112, 106, 101, 96, 94, 88, 84, 77, 70, 64, 57, 52, 44, 39, 39, 30, 30, 26, 23, 21, 20, 16, 16, 12, 11, 11, 10, 8, 7, 4, 4, 2, 2, 2, 363, 360, 344, 
332, 314, 305, 286, 269, 261, 248, 236, 226, 220, 214, 199, 197, 188, 180, 177, 168, 162, 155, 151, 150, 148, 144, 140, 136, 128, 122, 118, 116, 113, 113, 109, 105, 98, 94, 87, 83, 75, 69, 64, 57, 53, 49, 46, 40, 37, 31, 25, 22, 18, 14, 10, 9, 7, 6, 3, 2, 1, 369, 366, 358, 347, 333, 322, 314, 302, 294, 277, 267, 254, 240, 236, 225, 218, 210, 197, 190, 179, 179, 173, 166, 163, 160, 154, 152, 148, 145, 143, 139, 136, 127, 123, 119, 118, 116, 112, 109, 106, 102, 99, 94, 85, 74, 66, 58, 50, 42, 35, 31, 
28, 25, 20, 17, 15, 13, 11, 6, 3, 2, 1, 1) ,
 r  =  c(2, 2, 4, 6, 4, 5, 4, 4, 4, 2, 5, 1, 7, 1, 1, 3, 1, 0, 3, 2, 1, 3, 2, 2, 2, 3, 1, 0, 0, 0, 1, 4, 1, 1, 1, 0, 4, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 4, 10, 4, 2, 6, 4, 0, 5, 3, 3, 3, 5, 1, 5, 1, 0, 3, 3, 3, 0, 0, 2, 2, 0, 3, 1, 3, 0, 1, 2, 0, 0, 0, 0, 2, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 10, 7, 13, 10, 15, 7, 7, 9, 9, 3, 9, 10, 5, 5, 6, 5, 7, 5, 7, 3, 10, 1, 6, 6, 5, 3, 5, 5, 3, 0, 4, 3, 3, 0, 4, 1, 4, 1, 3, 2, 3, 2, 4, 9, 8, 12, 7, 11, 5, 7, 5, 4, 8, 7, 2, 10, 
6, 8, 4, 8, 6, 3, 5, 5, 6, 2, 5, 5, 1, 6, 5, 3, 1, 2, 5, 2, 1, 1, 2, 1, 1, 0, 0, 0, 5, 6, 7, 11, 12, 7, 3, 6, 12, 5, 6, 11, 8, 1, 11, 8, 8, 3, 10, 9, 4, 5, 9, 2, 6, 3, 2, 1, 2, 1, 0, 3, 0, 0, 0, 0, 13, 14, 15, 7, 10, 12, 5, 6, 7, 8, 8, 14, 7, 10, 7, 3, 13, 5, 5, 2, 4, 4, 3, 3, 1, 2, 5, 1, 1, 0, 2, 0, 0, 0, 0, 0, 1, 4, 4, 3, 2, 9, 5, 8, 2, 2, 4, 4, 4, 3, 3, 2, 3, 5, 1, 4, 3, 4, 2, 0, 1, 3, 3, 2, 2, 1, 1, 3, 1, 0, 4, 0, 0, 0, 0, 2, 3, 1, 2, 1, 3, 3, 1, 1, 5, 4, 5, 6, 4, 5, 3, 0, 0, 3, 2, 1, 3, 2, 5, 
7, 2, 2, 0, 1, 3, 4, 3, 1, 0, 3, 2, 3, 8, 6, 9, 14, 15, 10, 9, 8, 13, 11, 12, 21, 9, 8, 11, 9, 10, 9, 8, 6, 4, 6, 12, 4, 8, 4, 10, 6, 3, 1, 4, 2, 3, 2, 5, 5, 4, 6, 1, 0, 0, 0, 15, 0, 0, 3, 4, 1, 0, 0, 1, 0, 2, 0, 4, 0, 0, 0, 0, 0, 4, 8, 6, 12, 19, 9, 17, 11, 15, 12, 18, 10, 11, 11, 5, 9, 9, 6, 4, 10, 7, 9, 0, 5, 6, 6, 5, 7, 5, 5, 4, 3, 8, 9, 2, 6, 0, 5, 7, 3, 0, 3, 3, 6, 0, 3, 5, 1, 0, 0, 0, 2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 10, 7, 8, 3, 6, 4, 5, 4, 4, 0, 3, 3, 1, 1, 4, 3, 2, 2, 2, 2, 1, 1, 2, 
2, 0, 1, 3, 2, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 6, 2, 5, 5, 3, 8, 2, 4, 1, 2, 4, 3, 4, 4, 0, 2, 0, 3, 4, 1, 1, 0, 2, 3, 3, 0, 0, 0, 2, 0, 1, 1, 0, 0, 2, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 3, 15, 11, 18, 8, 19, 17, 8, 13, 11, 10, 6, 6, 15, 2, 8, 8, 3, 9, 6, 7, 4, 1, 2, 4, 4, 3, 8, 6, 4, 2, 3, 0, 4, 4, 7, 2, 4, 2, 4, 4, 3, 4, 1, 0, 0, 3, 0, 3, 2, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 3, 7, 11, 14, 10, 8, 12, 8, 17, 10, 13, 14, 4, 11, 7, 8, 13, 7, 11, 0, 6, 7, 3, 3, 6, 2, 3, 3, 2, 4, 3, 9, 4, 3, 
1, 2, 3, 1, 2, 2, 2, 4, 3, 5, 1, 2, 2, 2, 4, 0, 0, 0, 1, 0, 1, 1, 0, 3, 2, 0, 0, 0, 0) ,
 reprior_max  =  2 ,
 reprior_min  =  0 ,
 s  =  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2) ,
 t  = structure(.Data =  c(1, 3, 4, 5, 4, 5, 1, 4, 3, 6, 2, 1, 2, 3) ,.Dim =  c(7, 2) ) ,
 time  =  c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 
34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 
26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 
35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 
12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 
37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 
18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 
58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 61.5, 62.5) 
)

##############################################
# MODEL                                      #
##############################################
# Fractional polynomial, 1st order, random effects model: a RE is put only on the scale parameter 
#                                                          i.e. on the intercept in the frac poly
#                                                         a uniform prior is used for the RE standard deviation
# -------------------------------------------------------------------------------
# Data: grouped survival data, binomial likelihood, linear predictor on log-hazard
#         Nobs     number of observations
#         n[i]     patients at risk in interval i
#         r[i]     events during interval i
#         time[i]  mid-point of interval i
#         dt[i]    length of interval i
#         Ns       number of studies
#         Na[j]    number of arms in study j
#         Ntrt     number of treatments
#         s[i]     study number for obs i
#         a[i]     arm number (within study) for obs i
#         t[i,j]   treatment in study i arm j
#         P1       exponent of the time varying term in the fractional polynomial
#         feprior_mean[1:2]         prior mean (for contrasts d and baselines mu)
#         feprior_prec[1:2, 1:2]    prior precision (for d and mu)
#         reprior_min               lower bound of the unif prior for the RE SD
#         reprior_max               upper bound of the unif prior for the RE SD
# -------------------------------------------------------------------------------

model{

## Sampling model (likelihood)
for (i in 1:Nobs){
  time1[i] <- (equals(P1, 0) * log(time[i]) + (1 - equals(P1, 0)) * pow(time[i], P1)) 

  # likelihood: digitized KM curves
  r[i] ~ dbin(p[i], n[i])
  p[i] <- 1 - exp(-h[i] * dt[i])  # cumulative hazard over interval [t,t+dt] expressed as deaths per person-month

  # fractional polynomial
  log(h[i]) <- Beta[s[i], a[i], 1] + Beta[s[i], a[i], 2] * time1[i]
  }


## Arm level parameters = study effect + trt effect (RE model, consistency eq for pop pars)
for (i in 1:Ns){
  w[i, 1] <- 0
  delta[i, 1] <- 0

  for (j in 1:Na[i]){
    Beta[i, j, 1] <- mu[i, 1] + delta[i, j]
    Beta[i, j, 2] <- mu[i, 2] + d[t[i, j], 2] - d[t[i, 1], 2]
    }

  for (j in 2:Na[i]){
    delta[i, j] ~ dnorm(md[i, j], taud[i, j])
    md[i, j] <- d[t[i, j], 1] - d[t[i, 1], 1] + sw[i, j]
    w[i, j] <- (delta[i, j] - d[t[i, j], 1] + d[t[i, 1], 1])
    sw[i, j] <- sum(w[i, 1:(j - 1)]) / (j - 1) 
    taud[i, j] <- tau * 2 * (j - 1) / j 
    }
    
  }


## Priors
for (j in 1:Ns){
  mu[j, 1:2] ~ dmnorm(feprior_mean[1:2], feprior_prec[,]) 
  }

d[1, 1] <- 0
d[1, 2] <- 0
for (i in 2:Ntrt){
  d[i, 1:2] ~ dmnorm(feprior_mean[1:2], feprior_prec[,]) 
  }

sd ~ dunif(reprior_min, reprior_max)
tau <- 1 / (sd * sd)

} # end of model
##############################################

Second order random intercept model (unif[0,2] prior)

 ##############################################
# DATA                                       #
##############################################
list(
 Na  =  c(2, 2, 2, 2, 2, 2, 2) ,
 Nobs  =  675 ,
 Ns  =  7 ,
 Ntrt  =  6 ,
 P1  =  0 ,
 P2  =  1 ,
 a  =  c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2) ,
 dt  =  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1) ,
 feprior_mean  =  c(0, 0, 0) ,
 feprior_prec  = structure(.Data =  c(1e-04, 0, 0, 0, 1e-04, 0, 0, 0, 1e-04) ,.Dim =  c(3, 3) ) ,
 n  =  c(182, 177, 172, 164, 155, 148, 143, 138, 132, 128, 125, 120, 117, 110, 108, 106, 101, 98, 96, 90, 87, 83, 79, 75, 70, 66, 59, 56, 53, 50, 47, 45, 39, 37, 34, 32, 31, 26, 25, 24, 23, 21, 19, 16, 14, 13, 11, 10, 8, 6, 5, 3, 1, 183, 179, 172, 159, 152, 147, 139, 133, 131, 124, 119, 115, 110, 104, 101, 95, 93, 92, 88, 83, 80, 76, 73, 68, 62, 59, 53, 50, 43, 40, 37, 35, 33, 31, 30, 29, 26, 24, 22, 19, 18, 18, 18, 18, 17, 17, 15, 14, 12, 10, 8, 4, 2, 322, 319, 309, 301, 288, 277, 262, 255, 248, 238, 229, 
226, 216, 206, 200, 195, 189, 183, 176, 171, 163, 160, 150, 148, 142, 136, 130, 127, 121, 116, 113, 112, 108, 105, 101, 101, 96, 95, 91, 89, 86, 84, 327, 324, 318, 307, 297, 283, 275, 263, 256, 248, 241, 235, 225, 217, 212, 201, 194, 184, 178, 169, 161, 156, 149, 142, 134, 131, 125, 117, 115, 107, 101, 96, 93, 89, 83, 79, 76, 74, 70, 68, 65, 63, 375, 374, 368, 360, 352, 339, 326, 319, 315, 308, 295, 289, 283, 271, 262, 259, 247, 238, 229, 225, 213, 203, 197, 190, 180, 161, 135, 121, 98, 78, 61, 55, 
52, 45, 39, 35, 375, 372, 356, 338, 319, 308, 295, 282, 276, 269, 260, 251, 242, 228, 221, 210, 203, 200, 187, 180, 172, 168, 161, 155, 149, 132, 115, 104, 82, 66, 53, 48, 43, 39, 34, 30, 181, 181, 180, 176, 172, 169, 167, 158, 153, 145, 143, 141, 137, 133, 129, 126, 123, 121, 118, 113, 112, 108, 105, 101, 99, 99, 98, 95, 92, 90, 88, 87, 86, 83, 82, 82, 78, 78, 78, 78, 181, 179, 176, 175, 173, 172, 169, 166, 165, 164, 159, 155, 150, 144, 140, 135, 132, 132, 132, 129, 127, 126, 123, 121, 116, 109, 
107, 105, 105, 104, 101, 97, 94, 93, 93, 90, 557, 552, 541, 532, 521, 505, 487, 473, 462, 452, 436, 423, 409, 388, 377, 368, 357, 348, 338, 328, 320, 312, 305, 295, 280, 273, 262, 255, 243, 235, 230, 226, 219, 215, 208, 204, 197, 184, 171, 156, 147, 147, 114, 114, 99, 99, 80, 77, 57, 48, 38, 27, 21, 20, 17, 14, 8, 6, 6, 1, 1, 553, 544, 530, 518, 501, 480, 468, 448, 435, 420, 406, 387, 377, 366, 354, 348, 339, 329, 321, 315, 304, 295, 283, 280, 272, 263, 254, 247, 237, 230, 224, 218, 212, 203, 192, 
187, 180, 174, 162, 149, 138, 131, 120, 111, 100, 92, 80, 63, 56, 48, 38, 31, 20, 17, 14, 11, 9, 7, 5, 5, 1, 189, 184, 178, 166, 156, 144, 138, 131, 123, 117, 111, 106, 103, 98, 93, 88, 85, 79, 74, 68, 64, 59, 55, 50, 46, 43, 40, 37, 35, 31, 28, 26, 25, 20, 19, 15, 15, 12, 10, 9, 9, 3, 3, 3, 3, 188, 179, 175, 167, 159, 151, 141, 136, 131, 127, 123, 117, 112, 106, 101, 96, 94, 88, 84, 77, 70, 64, 57, 52, 44, 39, 39, 30, 30, 26, 23, 21, 20, 16, 16, 12, 11, 11, 10, 8, 7, 4, 4, 2, 2, 2, 363, 360, 344, 
332, 314, 305, 286, 269, 261, 248, 236, 226, 220, 214, 199, 197, 188, 180, 177, 168, 162, 155, 151, 150, 148, 144, 140, 136, 128, 122, 118, 116, 113, 113, 109, 105, 98, 94, 87, 83, 75, 69, 64, 57, 53, 49, 46, 40, 37, 31, 25, 22, 18, 14, 10, 9, 7, 6, 3, 2, 1, 369, 366, 358, 347, 333, 322, 314, 302, 294, 277, 267, 254, 240, 236, 225, 218, 210, 197, 190, 179, 179, 173, 166, 163, 160, 154, 152, 148, 145, 143, 139, 136, 127, 123, 119, 118, 116, 112, 109, 106, 102, 99, 94, 85, 74, 66, 58, 50, 42, 35, 31, 
28, 25, 20, 17, 15, 13, 11, 6, 3, 2, 1, 1) ,
 r  =  c(2, 2, 4, 6, 4, 5, 4, 4, 4, 2, 5, 1, 7, 1, 1, 3, 1, 0, 3, 2, 1, 3, 2, 2, 2, 3, 1, 0, 0, 0, 1, 4, 1, 1, 1, 0, 4, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 4, 10, 4, 2, 6, 4, 0, 5, 3, 3, 3, 5, 1, 5, 1, 0, 3, 3, 3, 0, 0, 2, 2, 0, 3, 1, 3, 0, 1, 2, 0, 0, 0, 0, 2, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 10, 7, 13, 10, 15, 7, 7, 9, 9, 3, 9, 10, 5, 5, 6, 5, 7, 5, 7, 3, 10, 1, 6, 6, 5, 3, 5, 5, 3, 0, 4, 3, 3, 0, 4, 1, 4, 1, 3, 2, 3, 2, 4, 9, 8, 12, 7, 11, 5, 7, 5, 4, 8, 7, 2, 10, 
6, 8, 4, 8, 6, 3, 5, 5, 6, 2, 5, 5, 1, 6, 5, 3, 1, 2, 5, 2, 1, 1, 2, 1, 1, 0, 0, 0, 5, 6, 7, 11, 12, 7, 3, 6, 12, 5, 6, 11, 8, 1, 11, 8, 8, 3, 10, 9, 4, 5, 9, 2, 6, 3, 2, 1, 2, 1, 0, 3, 0, 0, 0, 0, 13, 14, 15, 7, 10, 12, 5, 6, 7, 8, 8, 14, 7, 10, 7, 3, 13, 5, 5, 2, 4, 4, 3, 3, 1, 2, 5, 1, 1, 0, 2, 0, 0, 0, 0, 0, 1, 4, 4, 3, 2, 9, 5, 8, 2, 2, 4, 4, 4, 3, 3, 2, 3, 5, 1, 4, 3, 4, 2, 0, 1, 3, 3, 2, 2, 1, 1, 3, 1, 0, 4, 0, 0, 0, 0, 2, 3, 1, 2, 1, 3, 3, 1, 1, 5, 4, 5, 6, 4, 5, 3, 0, 0, 3, 2, 1, 3, 2, 5, 
7, 2, 2, 0, 1, 3, 4, 3, 1, 0, 3, 2, 3, 8, 6, 9, 14, 15, 10, 9, 8, 13, 11, 12, 21, 9, 8, 11, 9, 10, 9, 8, 6, 4, 6, 12, 4, 8, 4, 10, 6, 3, 1, 4, 2, 3, 2, 5, 5, 4, 6, 1, 0, 0, 0, 15, 0, 0, 3, 4, 1, 0, 0, 1, 0, 2, 0, 4, 0, 0, 0, 0, 0, 4, 8, 6, 12, 19, 9, 17, 11, 15, 12, 18, 10, 11, 11, 5, 9, 9, 6, 4, 10, 7, 9, 0, 5, 6, 6, 5, 7, 5, 5, 4, 3, 8, 9, 2, 6, 0, 5, 7, 3, 0, 3, 3, 6, 0, 3, 5, 1, 0, 0, 0, 2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 10, 7, 8, 3, 6, 4, 5, 4, 4, 0, 3, 3, 1, 1, 4, 3, 2, 2, 2, 2, 1, 1, 2, 
2, 0, 1, 3, 2, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 6, 2, 5, 5, 3, 8, 2, 4, 1, 2, 4, 3, 4, 4, 0, 2, 0, 3, 4, 1, 1, 0, 2, 3, 3, 0, 0, 0, 2, 0, 1, 1, 0, 0, 2, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 3, 15, 11, 18, 8, 19, 17, 8, 13, 11, 10, 6, 6, 15, 2, 8, 8, 3, 9, 6, 7, 4, 1, 2, 4, 4, 3, 8, 6, 4, 2, 3, 0, 4, 4, 7, 2, 4, 2, 4, 4, 3, 4, 1, 0, 0, 3, 0, 3, 2, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 3, 7, 11, 14, 10, 8, 12, 8, 17, 10, 13, 14, 4, 11, 7, 8, 13, 7, 11, 0, 6, 7, 3, 3, 6, 2, 3, 3, 2, 4, 3, 9, 4, 3, 
1, 2, 3, 1, 2, 2, 2, 4, 3, 5, 1, 2, 2, 2, 4, 0, 0, 0, 1, 0, 1, 1, 0, 3, 2, 0, 0, 0, 0) ,
 reprior_max  =  2 ,
 reprior_min  =  0 ,
 s  =  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2) ,
 t  = structure(.Data =  c(1, 3, 4, 5, 4, 5, 1, 4, 3, 6, 2, 1, 2, 3) ,.Dim =  c(7, 2) ) ,
 time  =  c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 
34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 
26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 
35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 
12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 
37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 
18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 
58.5, 59.5, 60.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 61.5, 62.5) 
)

##############################################
# MODEL                                      #
##############################################
# Fractional polynomial, 2nd order, random effects model: a RE is put only on the scale parameter 
#                                                          i.e. on the intercept in the frac poly
#                                                         a uniform prior is used for the RE standard deviation
# -------------------------------------------------------------------------------
# Data: grouped survival data, binomial likelihood, linear predictor on log-hazard
#         Nobs     number of observations
#         n[i]     patients at risk in interval i
#         r[i]     events during interval i
#         time[i]  mid-point of interval i
#         dt[i]    length of interval i
#         Ns       number of studies
#         Na[j]    number of arms in study j
#         Ntrt     number of treatments
#         s[i]     study number for obs i
#         a[i]     arm number (within study) for obs i
#         t[i,j]   treatment in study i arm j
#         P1       exponent of the time varying term in the fractional polynomial
#         feprior_mean[1:3]         prior mean (for contrasts d and baselines mu)
#         feprior_prec[,]           prior precision (for d and mu)
#         reprior_min               lower bound of the unif prior for the RE SD
#         reprior_max               upper bound of the unif prior for the RE SD
# -------------------------------------------------------------------------------

model{

## Sampling model (likelihood)
for (i in 1:Nobs){
  time1[i] <- (equals(P1, 0) * log(time[i]) + (1 - equals(P1, 0)) * pow(time[i], P1)) 
  time2[i] <- ((1 - equals(P2, P1)) * (equals(P2, 0) * log(time[i]) + (1 - equals(P2, 0)) * pow(time[i], P2)) + 
                  equals(P2, P1) * (equals(P2, 0) * log(time[i]) * log(time[i]) + (1 - equals(P2, 0)) * pow(time[i], P2) * log(time[i])))


  # likelihood: digitized KM curves
  r[i] ~ dbin(p[i], n[i])
  p[i] <- 1 - exp(-h[i] * dt[i])  # cumulative hazard over interval [t,t+dt] expressed as deaths per person-month

  # fractional polynomial
  log(h[i]) <- Beta[s[i], a[i], 1] + Beta[s[i], a[i], 2] * time1[i] + Beta[s[i], a[i], 3] * time2[i]
  }


  ## Arm level parameters = study effect + trt effect (RE model, consistency eq for pop pars)
  for (i in 1:Ns){
    w[i, 1] <- 0
    delta[i, 1] <- 0

    for (j in 1:Na[i]){
      Beta[i, j, 1] <- mu[i, 1] + delta[i, j]
      Beta[i, j, 2] <- mu[i, 2] + d[t[i, j], 2] - d[t[i, 1], 2]
      Beta[i, j, 3] <- mu[i, 3] + d[t[i, j], 3] - d[t[i, 1], 3]
      }

    for (j in 2:Na[i]){
      delta[i, j] ~ dnorm(md[i, j], taud[i, j])
      md[i, j] <- d[t[i, j], 1] - d[t[i, 1], 1] + sw[i, j]
      w[i, j] <- (delta[i, j] - d[t[i, j], 1] + d[t[i, 1], 1])
      sw[i, j] <- sum(w[i, 1:(j - 1)]) / (j - 1) 
      taud[i, j] <- tau * 2 * (j - 1) / j 
      }
    
    }


     
## Priors
for (j in 1:Ns){
  mu[j, 1:3] ~ dmnorm(feprior_mean[1:3], feprior_prec[,]) 
  }

d[1, 1] <- 0
d[1, 2] <- 0
d[1, 3] <- 0
for (i in 2:Ntrt){
  d[i, 1:3] ~ dmnorm(feprior_mean[1:3], feprior_prec[,]) 
  }


sd ~ dunif(reprior_min, reprior_max)
tau <- 1 / (sd * sd)

} # end of model
##############################################

Session info

date()
## [1] "Fri Oct 11 04:18:36 2024"
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggmcmc_1.5.1.1  ggplot2_3.5.1   tidyr_1.3.1     gemtcPlus_1.0.0
##  [5] R2jags_0.8-5    rjags_4-16      gemtc_1.0-2     coda_0.19-4.1  
##  [9] dplyr_1.1.4     rmarkdown_2.28 
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.5          xfun_0.48             bslib_0.8.0          
##  [4] GGally_2.2.1          CompQuadForm_1.4.3    lattice_0.22-6       
##  [7] tzdb_0.4.0            numDeriv_2016.8-1.1   mathjaxr_1.6-0       
## [10] vctrs_0.6.5           tools_4.4.1           generics_0.1.3       
## [13] parallel_4.4.1        tibble_3.2.1          fansi_1.0.6          
## [16] highr_0.11            pkgconfig_2.0.3       Matrix_1.7-0         
## [19] RColorBrewer_1.1-3    truncnorm_1.0-9       lifecycle_1.0.4      
## [22] farver_2.1.2          compiler_4.4.1        stringr_1.5.1        
## [25] munsell_0.5.1         htmltools_0.5.8.1     sys_3.4.3            
## [28] buildtools_1.0.0      sass_0.4.9            yaml_2.3.10          
## [31] crayon_1.5.3          pillar_1.9.0          nloptr_2.1.1         
## [34] jquerylib_0.1.4       MASS_7.3-61           cachem_1.1.0         
## [37] meta_7.0-0            metadat_1.2-0         boot_1.3-31          
## [40] abind_1.4-8           nlme_3.1-166          ggstats_0.7.0        
## [43] network_1.18.2        tidyselect_1.2.1      digest_0.6.37        
## [46] slam_0.1-53           stringi_1.8.4         reshape2_1.4.4       
## [49] pander_0.6.5          purrr_1.0.2           labeling_0.4.3       
## [52] maketools_1.3.1       forcats_1.0.0         splines_4.4.1        
## [55] fastmap_1.2.0         grid_4.4.1            colorspace_2.1-1     
## [58] cli_3.6.3             metafor_4.6-0         magrittr_2.0.3       
## [61] utf8_1.2.4            readr_2.1.5           withr_3.0.1          
## [64] scales_1.3.0          igraph_2.0.3          lme4_1.1-35.5        
## [67] hms_1.1.3             evaluate_1.0.1        knitr_1.48           
## [70] Rglpk_0.6-5.1         rlang_1.1.4           Rcpp_1.0.13          
## [73] glue_1.8.0            xml2_1.3.6            minqa_1.2.8          
## [76] jsonlite_1.8.9        R6_2.5.1              plyr_1.8.9           
## [79] statnet.common_4.10.0 R2WinBUGS_2.1-22.1