Skip to content

Commit

Permalink
Merge branch 'feature/fix_multinomial_int_thin' into develop
Browse files Browse the repository at this point in the history
  • Loading branch information
@dmphillippo
dmphillippo committed Jan 10, 2024
2 parents 2ef5a50 + 5592c63 commit 5e12870
Show file tree
Hide file tree
Showing 39 changed files with 4,448 additions and 3,737 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: multinma
Title: Bayesian Network Meta-Analysis of Individual and Aggregate Data
Version: 0.5.1.9006
Version: 0.5.1.9007
Authors@R:
person(given = c("David", "M."),
family = "Phillippo",
Expand Down
2 changes: 2 additions & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,8 @@ a data frame of layout coordinates) now works as expected when `nudge > 0`.
* Fix: Documentation corrections (PR #24).
* Fix: Added missing `as.tibble.stan_nma()` and `as_tibble.stan_nma()` methods,
to complement the existing `as.data.frame.stan_nma()`.
* Fix: Bug in ordered multinomial models where data in studies with missing
categories could be assigned the wrong category (#28).

# multinma 0.5.1

Expand Down
7 changes: 6 additions & 1 deletion R/nma.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1170,7 +1170,12 @@ nma <- function(network,
agd_arm_cumint_labels <- rep(agd_arm_data_labels, each = n_int_thin)
agd_arm_cumint_labels <- paste0(agd_arm_cumint_labels, ", ", rep_len(1:n_int_thin * int_thin, length.out = length(agd_arm_cumint_labels)))

fnames_oi[grepl("^theta_bar_cum_agd_arm\\[[0-9]+\\]$", fnames_oi)] <- paste0("theta_bar_cum_agd_arm[", agd_arm_cumint_labels, "]")
if (likelihood == "ordered") {
l_cat <- colnames(y_agd_arm$.r)
agd_arm_cumint_labels <- paste0(rep(agd_arm_cumint_labels, times = length(l_cat)), ", ", rep(l_cat, each = length(agd_arm_cumint_labels)))
}

fnames_oi[grepl("^theta_bar_cum_agd_arm\\[", fnames_oi)] <- paste0("theta_bar_cum_agd_arm[", agd_arm_cumint_labels, "]")
if (likelihood %in% c("bernoulli2", "binomial2"))
fnames_oi[grepl("^theta2_bar_cum\\[[0-9]+\\]$", fnames_oi)] <- paste0("theta2_bar_cum[", agd_arm_cumint_labels, "]")
}
Expand Down
19 changes: 13 additions & 6 deletions R/stan_nma-class.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -566,6 +566,8 @@ plot_integration_error <- function(x, ...,

# Get cumulative integration points
twoparbin <- x$likelihood %in% c("binomial2", "bernoulli2")
multi <- x$likelihood == "ordered"

ipars <- c()
if (has_agd_arm(x$network)) {
ipars <- c(ipars, "theta_bar_cum_agd_arm")
Expand All @@ -588,23 +590,24 @@ plot_integration_error <- function(x, ...,

n_int <- x$network$n_int

rx <- "^(theta2?)\\[(.+): (.+), ([0-9]+)\\]$"
rx <- if (multi) "^(theta2?)\\[(.+): (.+), ([0-9]+), (.+)\\]$" else "^(theta2?)\\[(.+): (.+), ([0-9]+)\\]$"

int_dat <- tidyr::pivot_longer(int_dat, cols = -dplyr::one_of(".draw"),
names_pattern = rx,
names_to = c("parameter", "study", "treatment", "n_int"),
names_to = if (multi) c("parameter", "study", "treatment", "n_int", "category") else c("parameter", "study", "treatment", "n_int"),
names_transform = list(n_int = as.integer),
values_to = "value")

int_dat$study <- factor(int_dat$study, levels = levels(x$network$studies))
int_dat$treatment <- factor(int_dat$treatment, levels = levels(x$network$treatments))
if (multi) int_dat$category <- forcats::fct_inorder(factor(int_dat$category))

# Estimate integration error by subtracting final value
int_dat <- dplyr::left_join(dplyr::filter(int_dat, .data$n_int != max(.data$n_int)),
dplyr::filter(int_dat, .data$n_int == max(.data$n_int)) %>%
dplyr::rename(final_value = "value") %>%
dplyr::select(-"n_int"),
by = c("parameter", "study", "treatment", ".draw")) %>%
by = if (multi) c("parameter", "study", "treatment", "category", ".draw") else c("parameter", "study", "treatment", ".draw")) %>%
dplyr::mutate(diff = .data$value - .data$final_value)

int_thin <- min(int_dat$n_int)
Expand Down Expand Up @@ -675,9 +678,13 @@ plot_integration_error <- function(x, ...,
args = rlang::dots_list(orientation = orientation, ..., !!! v_args, .homonyms = "first"))
}

p <- p +
ggplot2::facet_wrap(~ study + treatment) +
theme_multinma()
if (multi) {
p <- p + ggplot2::facet_grid(study + treatment ~ category)
} else {
p <- p + ggplot2::facet_wrap(~ study + treatment)
}

p <- p + theme_multinma()

return(p)
}
Expand Down
72 changes: 37 additions & 35 deletions inst/stan/ordered_multinomial.stan
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,26 +59,20 @@ transformed parameters {
// Is this only necessary if link > 2? Since ordered_(logistic|probit) are available
for (i in 1:ni_ipd) {
vector[ipd_ncat[i] - 1] q_temp;
for (k in 1:(ipd_ncat[i] - 1)) {
if (link == 1) // logit link
q_temp[k] = inv_logit(eta_ipd[i] - cc[ipd_cat[i, k+1]-1]);
else if (link == 2) // probit link
q_temp[k] = Phi(eta_ipd[i] - cc[ipd_cat[i, k+1]-1]);
else if (link == 3) // cloglog link
q_temp[k] = inv_cloglog(eta_ipd[i] - cc[ipd_cat[i, k+1]-1]);
}

// Category 1
if (link == 1) // logit link
q_temp[1] = inv_logit(eta_ipd[i] - cc[1]);
else if (link == 2) // probit link
q_temp[1] = Phi(eta_ipd[i] - cc[1]);
else if (link == 3) // cloglog link
q_temp[1] = inv_cloglog(eta_ipd[i] - cc[1]);

theta_ipd[i, 1] = 1 - q_temp[1];

// Categories 2:(ipd_ncat - 1)
for (k in 2:(ipd_ncat[i] - 1)) {
if (link == 1) // logit link
q_temp[k] = inv_logit(eta_ipd[i] - cc[ipd_cat[i, k]]);
else if (link == 2) // probit link
q_temp[k] = Phi(eta_ipd[i] - cc[ipd_cat[i, k]]);
else if (link == 3) // cloglog link
q_temp[k] = inv_cloglog(eta_ipd[i] - cc[ipd_cat[i, k]]);

// Store predictor in actual category column, rather than left-aligned
theta_ipd[i, ipd_cat[i, k]] = q_temp[k - 1] - q_temp[k];
}
Expand Down Expand Up @@ -107,7 +101,7 @@ transformed parameters {
eta_agd_arm_RE = eta_agd_arm_noRE[(1 + (i-1)*nint_max):((i-1)*nint_max + nint)];

for (k in 1:(agd_arm_ncat[i] - 1)) {
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), agd_arm_cat[i, k]] = inv_logit(eta_agd_arm_RE - cc[agd_arm_cat[i, k]]);
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), k] = inv_logit(eta_agd_arm_RE - cc[agd_arm_cat[i, k+1]-1]);
}
}
} else if (link == 2) { // probit link
Expand All @@ -118,7 +112,7 @@ transformed parameters {
eta_agd_arm_RE = eta_agd_arm_noRE[(1 + (i-1)*nint_max):((i-1)*nint_max + nint)];

for (k in 1:(agd_arm_ncat[i] - 1)) {
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), agd_arm_cat[i, k]] = Phi(eta_agd_arm_RE - cc[agd_arm_cat[i, k]]);
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), k] = Phi(eta_agd_arm_RE - cc[agd_arm_cat[i, k+1]-1]);
}
}
} else if (link == 3) { // cloglog link
Expand All @@ -129,7 +123,7 @@ transformed parameters {
eta_agd_arm_RE = eta_agd_arm_noRE[(1 + (i-1)*nint_max):((i-1)*nint_max + nint)];

for (k in 1:(agd_arm_ncat[i] - 1)) {
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), agd_arm_cat[i, k]] = inv_cloglog(eta_agd_arm_RE - cc[agd_arm_cat[i, k]]);
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), k] = inv_cloglog(eta_agd_arm_RE - cc[agd_arm_cat[i, k+1]-1]);
}
}
}
Expand All @@ -139,19 +133,19 @@ transformed parameters {
if (link == 1) { // logit link
for (i in 1:ni_agd_arm) {
for (k in 1:(agd_arm_ncat[i] - 1)) {
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), agd_arm_cat[i, k]] = inv_logit(eta_agd_arm_noRE[(1 + (i-1)*nint_max):((i-1)*nint_max + nint)] - cc[agd_arm_cat[i, k]]);
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), k] = inv_logit(eta_agd_arm_noRE[(1 + (i-1)*nint_max):((i-1)*nint_max + nint)] - cc[agd_arm_cat[i, k+1]-1]);
}
}
} else if (link == 2) { // probit link
for (i in 1:ni_agd_arm) {
for (k in 1:(agd_arm_ncat[i] - 1)) {
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), agd_arm_cat[i, k]] = Phi(eta_agd_arm_noRE[(1 + (i-1)*nint_max):((i-1)*nint_max + nint)] - cc[agd_arm_cat[i, k]]);
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), k] = Phi(eta_agd_arm_noRE[(1 + (i-1)*nint_max):((i-1)*nint_max + nint)] - cc[agd_arm_cat[i, k+1]-1]);
}
}
} else if (link == 3) { // cloglog link
for (i in 1:ni_agd_arm) {
for (k in 1:(agd_arm_ncat[i] - 1)) {
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), agd_arm_cat[i, k]] = inv_cloglog(eta_agd_arm_noRE[(1 + (i-1)*nint_max):((i-1)*nint_max + nint)] - cc[agd_arm_cat[i, k]]);
theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), k] = inv_cloglog(eta_agd_arm_noRE[(1 + (i-1)*nint_max):((i-1)*nint_max + nint)] - cc[agd_arm_cat[i, k+1]-1]);
}
}
}
Expand All @@ -160,7 +154,7 @@ transformed parameters {

for (i in 1:ni_agd_arm) {
for (k in 1:(agd_arm_ncat[i] - 1)) {
q_agd_arm_bar[i, agd_arm_cat[i, k]] = mean(theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), agd_arm_cat[i, k]]);
q_agd_arm_bar[i, k] = mean(theta_agd_arm_ii[(1 + (i-1)*nint):(i*nint), k]);
}
}

Expand All @@ -178,7 +172,7 @@ transformed parameters {
eta_agd_arm_RE = eta_agd_arm_noRE[i];

for (k in 1:(agd_arm_ncat[i] - 1)) {
q_agd_arm_bar[i, agd_arm_cat[i, k]] = inv_logit(eta_agd_arm_RE - cc[agd_arm_cat[i, k]]);
q_agd_arm_bar[i, k] = inv_logit(eta_agd_arm_RE - cc[agd_arm_cat[i, k+1]-1]);
}
}
} else if (link == 2) { // probit link
Expand All @@ -189,7 +183,7 @@ transformed parameters {
eta_agd_arm_RE = eta_agd_arm_noRE[i];

for (k in 1:(agd_arm_ncat[i] - 1)) {
q_agd_arm_bar[i, agd_arm_cat[i, k]] = Phi(eta_agd_arm_RE - cc[agd_arm_cat[i, k]]);
q_agd_arm_bar[i, k] = Phi(eta_agd_arm_RE - cc[agd_arm_cat[i, k+1]-1]);
}
}
} else if (link == 3) { // cloglog link
Expand All @@ -200,7 +194,7 @@ transformed parameters {
eta_agd_arm_RE = eta_agd_arm_noRE[i];

for (k in 1:(agd_arm_ncat[i] - 1)) {
q_agd_arm_bar[i, agd_arm_cat[i, k]] = inv_cloglog(eta_agd_arm_RE - cc[agd_arm_cat[i, k]]);
q_agd_arm_bar[i, k] = inv_cloglog(eta_agd_arm_RE - cc[agd_arm_cat[i, k+1]-1]);
}
}
}
Expand All @@ -210,19 +204,19 @@ transformed parameters {
if (link == 1) { // logit link
for (i in 1:ni_agd_arm) {
for (k in 1:(agd_arm_ncat[i] - 1)) {
q_agd_arm_bar[i, agd_arm_cat[i, k]] = inv_logit(eta_agd_arm_noRE[i] - cc[agd_arm_cat[i, k]]);
q_agd_arm_bar[i, k] = inv_logit(eta_agd_arm_noRE[i] - cc[agd_arm_cat[i, k+1]-1]);
}
}
} else if (link == 2) { // probit link
for (i in 1:ni_agd_arm) {
for (k in 1:(agd_arm_ncat[i] - 1)) {
q_agd_arm_bar[i, agd_arm_cat[i, k]] = Phi(eta_agd_arm_noRE[i] - cc[agd_arm_cat[i, k]]);
q_agd_arm_bar[i, k] = Phi(eta_agd_arm_noRE[i] - cc[agd_arm_cat[i, k+1]-1]);
}
}
} else if (link == 3) { // cloglog link
for (i in 1:ni_agd_arm) {
for (k in 1:(agd_arm_ncat[i] - 1)) {
q_agd_arm_bar[i, agd_arm_cat[i, k]] = inv_cloglog(eta_agd_arm_noRE[i] - cc[agd_arm_cat[i, k]]);
q_agd_arm_bar[i, k] = inv_cloglog(eta_agd_arm_noRE[i] - cc[agd_arm_cat[i, k+1]-1]);
}
}
}
Expand All @@ -237,10 +231,10 @@ transformed parameters {

// Categories 2:(agd_arm_ncat - 1)
for (k in 2:(agd_arm_ncat[i] - 1))
theta_agd_arm_bar[i, agd_arm_cat[i, k]] = q_agd_arm_bar[i, agd_arm_cat[i, k - 1]] - q_agd_arm_bar[i, agd_arm_cat[i, k]];
theta_agd_arm_bar[i, agd_arm_cat[i, k]] = q_agd_arm_bar[i, k - 1] - q_agd_arm_bar[i, k];

// Category agd_arm_ncat
theta_agd_arm_bar[i, agd_arm_cat[i, agd_arm_ncat[i]]] = q_agd_arm_bar[i, agd_arm_cat[i, agd_arm_ncat[i] - 1]];
theta_agd_arm_bar[i, agd_arm_cat[i, agd_arm_ncat[i]]] = q_agd_arm_bar[i, agd_arm_ncat[i] - 1];
}

}
Expand Down Expand Up @@ -299,13 +293,21 @@ generated quantities {
}
}

// Cumulative integration - note this is for the intermediate q
// Cumulative integration
for (i in 1:ni_agd_arm) {
for (k in 1:(ncat - 1)) {
for (j in 1:n_int_thin) {
theta_bar_cum_agd_arm[(i - 1)*n_int_thin + j, agd_arm_cat[i, k]] = mean(theta_agd_arm_ii[(1 + (i - 1)*nint):((i - 1)*nint + j*int_thin), agd_arm_cat[i, k]]);
}
for (j in 1:n_int_thin) {
vector[ncat - 1] q_agd_arm_bar_cum;
for (k in 1:(agd_arm_ncat[i] - 1)) q_agd_arm_bar_cum[k] = mean(theta_agd_arm_ii[(1 + (i - 1)*nint):((i - 1)*nint + j*int_thin), k]);

// Category 1
theta_bar_cum_agd_arm[(i - 1)*n_int_thin + j, 1] = 1 - q_agd_arm_bar_cum[1];

// Categories 2:(agd_arm_ncat - 1)
for (k in 2:(agd_arm_ncat[i] - 1))
theta_bar_cum_agd_arm[(i - 1)*n_int_thin + j, agd_arm_cat[i, k]] = q_agd_arm_bar_cum[k - 1] - q_agd_arm_bar_cum[k];

// Category agd_arm_ncat
theta_bar_cum_agd_arm[(i - 1)*n_int_thin + j, agd_arm_cat[i, agd_arm_ncat[i]]] = q_agd_arm_bar_cum[agd_arm_ncat[i] - 1];
}
}

}
Loading

0 comments on commit 5e12870

Please sign in to comment.