Movatterモバイル変換

Type:	Package
Title:	R Bayesian Evidence Synthesis Tools
Description:	Tool-set to support Bayesian evidence synthesis. This includes meta-analysis, (robust) prior derivation from historical data, operating characteristics and analysis (1 and 2 sample cases). Please refer to Weber et al. (2021) <doi:10.18637/jss.v100.i19> for details on applying this package while Neuenschwander et al. (2010) <doi:10.1177/1740774509356002> and Schmidli et al. (2014) <doi:10.1111/biom.12242> explain details on the methodology.
Version:	1.8-2
Date:	2025-04-25
Depends:	R (≥ 3.5.0)
Imports:	methods, Rcpp (≥ 0.12.0), RcppParallel (≥ 5.0.1), rstan (≥2.32.0), rstantools (≥ 2.4.0), assertthat, mvtnorm, Formula,checkmate, bayesplot (≥ 1.4.0), ggplot2, dplyr, stats, utils,matrixStats, abind, rlang, jsonlite, lifecycle
LinkingTo:	BH (≥ 1.72.0), Rcpp (≥ 0.12.0), RcppEigen (≥ 0.3.3.3.0),RcppParallel (≥ 5.0.1), rstan (≥ 2.32.0), StanHeaders (≥2.32.0)
License:	GPL (≥ 3)
LazyData:	true
Biarch:	true
NeedsCompilation:	yes
UseLTO:	true
URL:	https://opensource.nibr.com/RBesT/
BugReports:	https://github.com/Novartis/RBesT/issues
Suggests:	rmarkdown, knitr, testthat (≥ 2.0.0), foreach, purrr,rstanarm (≥ 2.17.2), scales, tools, broom, tidyr, parallel,brms, glue, ragg, withr
VignetteBuilder:	knitr
SystemRequirements:	GNU make, pandoc (>= 1.12.3), pngquant, C++17
Encoding:	UTF-8
RoxygenNote:	7.3.2
Config/testthat/edition:	3
Packaged:	2025-04-25 15:24:26 UTC; weberse2
Author:	Novartis Pharma AG [cph], Sebastian Weber [aut, cre], Beat Neuenschwander [ctb], Heinz Schmidli [ctb], Baldur Magnusson [ctb], Yue Li [ctb], Satrajit Roychoudhury [ctb], Lukas A. Widmer [ctb], Trustees of Columbia University [cph] (R/stanmodels.R, configure, configure.win)
Maintainer:	Sebastian Weber <sebastian.weber@novartis.com>
Repository:	CRAN
Date/Publication:	2025-04-25 18:30:02 UTC

R Bayesian Evidence Synthesis Tools

Description

The RBesT tools are designed to support in the derivation ofparametric informative priors, asses design characeristics andperform analyses. Supported endpoints include normal, binary andPoisson.

Details

For introductory material, please refer to the vignettes which include

• Introduction (binary)

• Introduction (normal)

• Customizing RBesT Plots

• Robust MAP, advanced usage

The main function of the package isgMAP(). See it'shelp page for a detailed description of the statistical model.

Global Options

Version History

SeeNEWS.md file.

Author(s)

Maintainer: Sebastian Webersebastian.weber@novartis.com

Other contributors:

• Novartis Pharma AG [copyright holder]

• Beat Neuenschwanderbeat.neuenschwander@novartis.com [contributor]

• Heinz Schmidliheinz.schmidli@novartis.com [contributor]

• Baldur Magnussonbaldur.magnusson@novartis.com [contributor]

• Yue Liyue-1.li@novartis.com [contributor]

• Satrajit Roychoudhurysatrajit.roychoudhury@novartis.com [contributor]

• Lukas A. Widmerlukas_andreas.widmer@novartis.com (ORCID) [contributor]

• Trustees of Columbia University (R/stanmodels.R, configure, configure.win) [copyright holder]

References

Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.19.3. https://mc-stan.org

See Also

Useful links:

• https://opensource.nibr.com/RBesT/

• Report bugs athttps://github.com/Novartis/RBesT/issues

Ankylosing Spondylitis.

Description

Data set containing historical information for placebo for a phaseII trial of ankylosing spondylitis patients. The primary efficacyendpoint was the percentage of patients with a 20% responseaccording to the Assessment of SpondyloArthritis internationalSociety criteria for improvement (ASAS20) at week 6.

Usage

AS

Format

A data frame with 8 rows and 3 variables:

study

study

n

study size

r

number of events

References

Baeten D. et. al,The Lancet, 2013, (382), 9906, p 1705

Examples

## Setting up dummy sampling for fast execution of example## Please use 4 chains and 20x more warmup & iter in practice.user_mc_options <- options(RBesT.MC.warmup=50, RBesT.MC.iter=100,                            RBesT.MC.chains=2, RBesT.MC.thin=1)set.seed(34563)map_AS <- gMAP(cbind(r, n - r) ~ 1 | study,  family = binomial,  data = AS,  tau.dist = "HalfNormal", tau.prior = 1,  beta.prior = 2)## Recover user set sampling defaultsoptions(.user_mc_options)

Exact Confidence interval for Binary Proportion

Description

This function calculates the exact confidendence interval for aresponse rate presented byn andr.

Usage

BinaryExactCI(r, n, alpha = 0.05, drop = TRUE)

Arguments

Details

Confidence intervals are obtained by a procedure first given inClopper and Pearson (1934). This guarantees that the confidencelevel is at least (1-\alpha).

Details can be found in the publication listed below.

Value

100 (1-\alpha)\response rate

References

Clopper, C. J. & Pearson, E. S. The use of confidence orfiducial limits illustrated in the case of the binomial. Biometrika 1934.

Examples

BinaryExactCI(3, 20, 0.05)

Functional programming utilities

Description

function from functional

Usage

Curry(FUN, ...)

Summarize Arrays

Description

The function calculates summary statistics from arbitrary arrays.

Usage

SimSum(  x,  min.max = FALSE,  n.sim = FALSE,  probs = c(0.025, 0.5, 0.975),  margin = ifelse(is.null(dim(x) | length(dim(x)) == 1), 2, length(dim(x))))

Arguments

Details

The function calculates by default the mean, standarddeviation and the specified qantiles which are by default themedian and the 95% interval.

If a mulit-dimensional array is specified asx, then thefunction will by default calculate the summaries over the margin ofthe largest dimension. For the case of a vector and a matrix, thefunction will transpose the results for better readabiliy.

Automatic Fitting of Mixtures of Conjugate Distributions to a Sample

Description

Fitting a series of mixtures of conjugate distributions to asample, using Expectation-Maximization (EM). The number ofmixture components is specified by the vectorNc. First aNc[1] component mixture is fitted, then aNc[2]component mixture, and so on. The mixture providing the best AICvalue is then selected.

Usage

automixfit(sample, Nc = seq(1, 4), k = 6, thresh = -Inf, verbose = FALSE, ...)

Arguments

Details

Thetype argument specifies the distribution ofthe mixture components, and can be a normal, beta or gammadistribution.

The penalty parameterk is 2 for the standard AICdefinition.Collet (2003) suggested to use values in therange from 2 to 6, where larger values ofk penalize morecomplex models. To favor mixtures with fewer components a value of6 is used as default.

Value

As result the best fitting mixture model is returned,i.e. the model with lowest AIC. All other models are saved in theattributemodels.

References

Collet D.Modeling Survival Data in Medical Research.2003; Chapman and Hall/CRC.

Examples

# random sample of size 1000 from a mixture of 2 beta componentsbm <- mixbeta(beta1 = c(0.4, 20, 90), beta2 = c(0.6, 35, 65))bmSamp <- rmix(bm, 1000)# fit with EM mixture models with up to 10 components and stop if# AIC increasesbmFit <- automixfit(bmSamp, Nc = 1:10, thresh = 0, type = "beta")bmFit# advanced usage: find out about all discarded modelsbmFitAll <- attr(bmFit, "models")sapply(bmFitAll, AIC, k = 6)

Scrambles the order of a mcmc array object for usage as a mcmcsample. It is advisable to set order once per mcmc run, otherwisecorrelations in the mcmc sample will be lost.

Description

Scrambles the order of a mcmc array object for usage as a mcmcsample. It is advisable to set order once per mcmc run, otherwisecorrelations in the mcmc sample will be lost.

Usage

chains2sample(chains, order, drop = TRUE)

Fast column-wise calculation of unbiased variances

Description

Fast column-wise calculation of unbiased variances

Usage

colVars(a)

Ulcerative Colitis.

Description

Data set containing historical information for placebo arm of aphase II proof-of-concept trial for the treatment of ulcerativecolitis. The primary outcome is remission at week 8 (binary).

Usage

colitis

Format

A data frame with 4 rows and 3 variables:

study

study

n

study size

r

number of events

References

Neuenschwander B, Capkun-Niggli G, Branson M,Spiegelhalter DJ. Summarizing historical information on controls inclinical trials.Clin Trials. 2010; 7(1):5-18

Crohn's disease.

Description

Data set containing historical information for placebo arm ofrelevant studies for the treatment of Crohn's disease. The primaryoutcome is change from baseline in Crohn's Disease Activity Index(CDAI) over a duration of 6 weeks. Standard deviation of changefrom baseline endpoint is approximately 88.

Usage

crohn

Format

A data frame with 4 rows and 3 variables:

study

study

n

study size

y

mean CDAI change

References

Hueber W. et. al,Gut, 2012, 61(12):1693-1700

Examples

## Setting up dummy sampling for fast execution of example## Please use 4 chains and 20x more warmup & iter in practice.user_mc_options <- options(RBesT.MC.warmup=50, RBesT.MC.iter=100,                            RBesT.MC.chains=2, RBesT.MC.thin=1)set.seed(546346)map_crohn <- gMAP(cbind(y, y.se) ~ 1 | study,  family = gaussian,  data = transform(crohn, y.se = 88 / sqrt(n)),  weights = n,  tau.dist = "HalfNormal", tau.prior = 44,  beta.prior = cbind(0, 88))## Recover user set sampling defaultsoptions(.user_mc_options)

Beta-Binomial Probabilities

Description

Beta-Binomial Probabilities

Usage

dBetaBinomial(r, n, a, b, log = FALSE)

Arguments

Details

r,n,a,b can be scalar or vectors. If vectors are used, they must be of the same length

Decision Function for 1 Sample Designs

Description

The function sets up a 1 sample one-sided decision function with anarbitrary number of conditions.

Usage

decision1S(pc = 0.975, qc = 0, lower.tail = TRUE)oc1Sdecision(pc = 0.975, qc = 0, lower.tail = TRUE)

Arguments

Details

The function creates a one-sided decision function whichtakes two arguments. The first argument is expected to be a mixture(posterior) distribution. This distribution is tested whether itfulfills all the required threshold conditions specified with thepc andqc arguments and returns 1 if all conditionsare met and 0 otherwise. Hence, forlower.tail=TRUEconditioni is equivalent to

P(\theta \leq q_{c,i}) > p_{c,i}

and the decision function is implemented as indicator function onthe basis of the heavy-side step functionH(x) which is0forx \leq 0 and1 forx > 0. As all conditionsmust be met, the final indicator function returns

\Pi_i H_i(P(\theta \leq q_{c,i}) - p_{c,i} ).

When the second argument is set toTRUE a distance metric isreturned component-wise per defined condition as

D_i = \log(P(\theta < q_{c,i})) - \log(p_{c,i}) .

These indicator functions can be used as input for 1-sampleboundary, OC or PoS calculations usingoc1S() orpos1S() .

Value

The function returns a decision function which takes twoarguments. The first argument is expected to be a mixture(posterior) distribution which is tested if the specifiedconditions are met. The logical second argument determines if thefunction acts as an indicator function or if the function returnsthe distance from the decision boundary for each condition inlog-space, i.e. the distance is 0 at the decision boundary,negative for a 0 decision and positive for a 1 decision.

Functions

• oc1Sdecision(): Deprecated old function name. Please usedecision1S instead.

References

Neuenschwander B, Rouyrre N, Hollaender H, Zuber E,Branson M. A proof of concept phase II non-inferioritycriterion.Stat. in Med.. 2011, 30:1618-1627

See Also

Other design1S:decision1S_boundary(),oc1S(),pos1S()

Examples

# see Neuenschwander et al., 2011# example is for a time-to-event trial evaluating non-inferiority# using a normal approximation for the log-hazard ratio# reference scales <- 2theta_ni <- 0.4theta_a <- 0alpha <- 0.05beta <- 0.2za <- qnorm(1 - alpha)zb <- qnorm(1 - beta)n1 <- round((s * (za + zb) / (theta_ni - theta_a))^2) # n for which design was intendednL <- 233c1 <- theta_ni - za * s / sqrt(n1)# flat priorflat_prior <- mixnorm(c(1, 0, 100), sigma = s)# standard NI designdecA <- decision1S(1 - alpha, theta_ni, lower.tail = TRUE)# for double criterion with indecision point (mean estimate must be# lower than this)theta_c <- c1# double criterion design# statistical significance (like NI design)dec1 <- decision1S(1 - alpha, theta_ni, lower.tail = TRUE)# require mean to be at least as good as theta_cdec2 <- decision1S(0.5, theta_c, lower.tail = TRUE)# combinationdecComb <- decision1S(c(1 - alpha, 0.5), c(theta_ni, theta_c), lower.tail = TRUE)theta_eval <- c(theta_a, theta_c, theta_ni)# we can display the decision function definitiondecComb# and use it to decide if a given distribution fulfills all# criterions defined# for the priordecComb(flat_prior)# or for a possible outcome of the trial# here with HR of 0.8 for 40 eventsdecComb(postmix(flat_prior, m = log(0.8), n = 40))

Decision Boundary for 1 Sample Designs

Description

Calculates the decision boundary for a 1 sample design. This is thecritical value at which the decision function will change from 0(failure) to 1 (success).

Usage

decision1S_boundary(prior, n, decision, ...)## S3 method for class 'betaMix'decision1S_boundary(prior, n, decision, ...)## S3 method for class 'normMix'decision1S_boundary(prior, n, decision, sigma, eps = 1e-06, ...)## S3 method for class 'gammaMix'decision1S_boundary(prior, n, decision, eps = 1e-06, ...)

Arguments

Details

The specification of the 1 sample design (prior, samplesize and decision function,D(y)), uniquely defines thedecision boundary

y_c = \max_y\{D(y) = 1\},

which is the maximal value ofy whenever the decisionD(y)function changes its value from 1 to 0 for a decision functionwithlower.tail=TRUE (otherwise the definition isy_c =\max_{y}\{D(y) = 0\}). The decisionfunction may change at most at a single critical value as onlyone-sided decision functions are supported. Here,y is defined for binary and Poisson endpoints as the sufficientstatisticy = \sum_{i=1}^{n} y_i and for the normalcase as the mean\bar{y} = 1/n \sum_{i=1}^n y_i.

The convention for the critical valuey_c depends on whethera left (lower.tail=TRUE) or right-sided decision function(lower.tail=FALSE) is used. Forlower.tail=TRUE thecritical valuey_c is the largest value for which thedecision is 1,D(y \leq y_c) = 1, while forlower.tail=FALSE thenD(y > y_c) = 1 holds. This isaligned with the cumulative density function definition within R(see for examplepbinom()).

Value

Returns the critical valuey_c.

Methods (by class)

• decision1S_boundary(betaMix): Applies for binomial model with a mixturebeta prior. The calculations use exact expressions.

• decision1S_boundary(normMix): Applies for the normal model with knownstandard deviation\sigma and a normal mixture prior for themean. As a consequence from the assumption of a known standarddeviation, the calculation discards sampling uncertainty of thesecond moment. The functiondecision1S_boundary has an extraargumenteps (defaults to10^{-6}). The critical valuey_c is searched in the region of probability mass1-eps fory.

• decision1S_boundary(gammaMix): Applies for the Poisson model with a gammamixture prior for the rate parameter. The functiondecision1S_boundary takes an extra argumenteps (defaults to10^{-6})which determines the region of probability mass1-eps wherethe boundary is searched fory.

See Also

Other design1S:decision1S(),oc1S(),pos1S()

Examples

# non-inferiority example using normal approximation of log-hazard# ratio, see ?decision1S for all detailss <- 2flat_prior <- mixnorm(c(1, 0, 100), sigma = s)nL <- 233theta_ni <- 0.4theta_a <- 0alpha <- 0.05beta <- 0.2za <- qnorm(1 - alpha)zb <- qnorm(1 - beta)n1 <- round((s * (za + zb) / (theta_ni - theta_a))^2)theta_c <- theta_ni - za * s / sqrt(n1)# double criterion design# statistical significance (like NI design)dec1 <- decision1S(1 - alpha, theta_ni, lower.tail = TRUE)# require mean to be at least as good as theta_cdec2 <- decision1S(0.5, theta_c, lower.tail = TRUE)# combinationdecComb <- decision1S(c(1 - alpha, 0.5), c(theta_ni, theta_c), lower.tail = TRUE)# critical value of double criterion designdecision1S_boundary(flat_prior, nL, decComb)# ... is limited by the statistical significance ...decision1S_boundary(flat_prior, nL, dec1)# ... or the indecision point (whatever is smaller)decision1S_boundary(flat_prior, nL, dec2)

Decision Function for 2 Sample Designs

Description

The function sets up a 2 sample one-sided decision function with anarbitrary number of conditions on the difference distribution.

Usage

decision2S(  pc = 0.975,  qc = 0,  lower.tail = TRUE,  link = c("identity", "logit", "log"))oc2Sdecision(  pc = 0.975,  qc = 0,  lower.tail = TRUE,  link = c("identity", "logit", "log"))

Arguments

Details

This function creates a one-sided decision function on thebasis of the difference distribution in a 2 sample situation. Tosupport double criterion designs, seeNeuenschwander et al.,2010, an arbitrary number of criterions can be given. The decisionfunction demands that the probability mass below the critical valueqc of the difference\theta_1 - \theta_2 is at leastpc. Hence, forlower.tail=TRUE conditioni isequivalent to

P(\theta_1 - \theta_2 \leq q_{c,i}) > p_{c,i}

and the decision function is implemented as indicator functionusing the heavy-side step functionH(x) which is0 forx \leq 0 and1 forx > 0. As all conditions mustbe met, the final indicator function returns

\Pi_i H_i(P(\theta_1 - \theta_2 \leq q_{c,i}) - p_{c,i} ),

which is1 if all conditions are met and0otherwise. Forlower.tail=FALSE differences must be greaterthan the given quantilesqc.

Note that whenever alink other thanidentity isrequested, then the underlying densities are first transformedusing the link function and then the probabilties for thedifferences are calculated in the transformed space. Hence, for abinary endpoint the defaultidentity link will calculaterisk differences, thelogit link will lead to decisionsbased on the differences inlogits corresponding to acriterion based on the log-odds. Thelog link will evaluateratios instead of absolute differences which could be useful for abinary endpoint or counting rates. The respective criticalquantilesqc must be given on the transformed scale.

Value

The function returns a decision function which takes threearguments. The first and second argument are expected to be mixture(posterior) distributions from which the difference distribution isformed and all conditions are tested. The third argument determinesif the function acts as an indicator function or if the functionreturns the distance from the decision boundary for each conditionin log-space. That is, the distance is 0 at the decision boundary,negative for a 0 decision and positive for a 1 decision.

Functions

• oc2Sdecision(): Deprecated old function name. Please usedecision2S instead.

References

Gsponer T, Gerber F, Bornkamp B, Ohlssen D,Vandemeulebroecke M, Schmidli H.A practical guide to Bayesian groupsequential designs.Pharm. Stat.. 2014; 13: 71-80

See Also

Other design2S:decision2S_boundary(),oc2S(),pos2S()

Examples

# see Gsponer et al., 2010priorT <- mixnorm(c(1, 0, 0.001), sigma = 88, param = "mn")priorP <- mixnorm(c(1, -49, 20), sigma = 88, param = "mn")# the success criteria is for delta which are larger than some# threshold value which is why we set lower.tail=FALSEsuccessCrit <- decision2S(c(0.95, 0.5), c(0, 50), FALSE)# the futility criterion acts in the opposite directionfutilityCrit <- decision2S(c(0.90), c(40), TRUE)print(successCrit)print(futilityCrit)# consider decision for specific outcomespostP_interim <- postmix(priorP, n = 10, m = -50)postT_interim <- postmix(priorT, n = 20, m = -80)futilityCrit(postP_interim, postT_interim)successCrit(postP_interim, postT_interim)# Binary endpoint with double criterion decision on log-odds scale# 95% certain positive difference and an odds ratio of 2 at leastdecL2 <- decision2S(c(0.95, 0.5), c(0, log(2)), lower.tail = FALSE, link = "logit")# 95% certain positive difference and an odds ratio of 3 at leastdecL3 <- decision2S(c(0.95, 0.5), c(0, log(3)), lower.tail = FALSE, link = "logit")# data scenariopost1 <- postmix(mixbeta(c(1, 1, 1)), n = 40, r = 10)post2 <- postmix(mixbeta(c(1, 1, 1)), n = 40, r = 18)# positive outcome and a median odds ratio of at least 2 ...decL2(post2, post1)# ... but not more than 3decL3(post2, post1)

Decision Boundary for 2 Sample Designs

Description

Thedecision2S_boundary function defines a 2 sample design(priors, sample sizes, decision function) for the calculation ofthe decision boundary. A function is returned which calculates thecritical value of the first sampley_{1,c} as a function ofthe outcome in the second sampley_2. At the decisionboundary, the decision function will change between 0 (failure) and1 (success) for the respective outcomes.

Usage

decision2S_boundary(prior1, prior2, n1, n2, decision, ...)## S3 method for class 'betaMix'decision2S_boundary(prior1, prior2, n1, n2, decision, eps, ...)## S3 method for class 'normMix'decision2S_boundary(  prior1,  prior2,  n1,  n2,  decision,  sigma1,  sigma2,  eps = 1e-06,  Ngrid = 10,  ...)## S3 method for class 'gammaMix'decision2S_boundary(prior1, prior2, n1, n2, decision, eps = 1e-06, ...)

Arguments

Details

For a 2 sample design the specification of the priors, thesample sizes and the decision function,D(y_1,y_2), uniquelydefines the decision boundary

D_1(y_2) = \max_{y_1}\{D(y_1,y_2) = 1\},

which is the critical value ofy_{1,c} conditional on thevalue ofy_2 whenever the decisionD(y_1,y_2) functionchanges its value from 0 to 1 for a decision function withlower.tail=TRUE (otherwise the definition isD_1(y_2) =\max_{y_1}\{D(y_1,y_2) = 0\}). The decision function may change at most at a single criticalvalue for giveny_{2} as only one-sided decision functionsare supported. Here,y_2 is defined for binary and Poissonendpoints as the sufficient statisticy_2 = \sum_{i=1}^{n_2}y_{2,i} and for the normal case as the mean\bar{y}_2 = 1/n_2\sum_{i=1}^{n_2} y_{2,i}.

Value

Returns a function with a single argument. This functioncalculates in dependence of the outcomey_2 in sample 2 thecritical valuey_{1,c} for which the defined design willchange the decision from 0 to 1 (or vice versa, depending on thedecision function).

Methods (by class)

• decision2S_boundary(betaMix): Applies for binomial model with a mixturebeta prior. The calculations use exact expressions. If theoptional argumenteps is defined, then an approximate methodis used which limits the search for the decision boundary to theregion of1-eps probability mass. This is useful for designswith large sample sizes where an exact approach is very costly tocalculate.

• decision2S_boundary(normMix): Applies for the normal model with knownstandard deviation\sigma and normal mixture priors for themeans. As a consequence from the assumption of a known standarddeviation, the calculation discards sampling uncertainty of thesecond moment. The function has two extra arguments (withdefaults):eps (10^{-6}) andNgrid (10). Thedecision boundary is searched in the region of probability mass1-eps, respectively fory_1 andy_2. Thecontinuous decision function is evaluated at a discrete grid, whichis determined by a spacing with\delta_2 =\sigma_2/\sqrt{N_{grid}}. Once the decision boundary is evaluatedat the discrete steps, a spline is used to inter-polate thedecision boundary at intermediate points.

• decision2S_boundary(gammaMix): Applies for the Poisson model with a gammamixture prior for the rate parameter. The functiondecision2S_boundary takes an extra argumenteps (defaults to10^{-6}) whichdetermines the region of probability mass1-eps where theboundary is searched fory_1 andy_2, respectively.

See Also

Other design2S:decision2S(),oc2S(),pos2S()

Examples

# see ?decision2S for details of examplepriorT <- mixnorm(c(1, 0, 0.001), sigma = 88, param = "mn")priorP <- mixnorm(c(1, -49, 20), sigma = 88, param = "mn")# the success criteria is for delta which are larger than some# threshold value which is why we set lower.tail=FALSEsuccessCrit <- decision2S(c(0.95, 0.5), c(0, 50), FALSE)# the futility criterion acts in the opposite directionfutilityCrit <- decision2S(c(0.90), c(40), TRUE)# success criterion boundarysuccessBoundary <- decision2S_boundary(priorP, priorT, 10, 20, successCrit)# futility criterion boundaryfutilityBoundary <- decision2S_boundary(priorP, priorT, 10, 20, futilityCrit)curve(successBoundary(x), -25:25 - 49, xlab = "y2", ylab = "critical y1")curve(futilityBoundary(x), lty = 2, add = TRUE)# hence, for mean in sample 2 of 10, the critical value for y1 isy1c <- futilityBoundary(-10)# around the critical value the decision for futility changesfutilityCrit(postmix(priorP, m = y1c + 1E-3, n = 10), postmix(priorT, m = -10, n = 20))futilityCrit(postmix(priorP, m = y1c - 1E-3, n = 10), postmix(priorT, m = -10, n = 20))

Transform Densities with a link function

Description

One-to-one transforms (mixture) of densities using a link function.

Usage

dlink(object) <- value

Arguments

Effective Sample Size for a Conjugate Prior

Description

Calculates the Effective Sample Size (ESS) for a mixture prior. TheESS indicates how many experimental units the prior is roughlyequivalent to.

Usage

ess(mix, method = c("elir", "moment", "morita"), ...)## S3 method for class 'betaMix'ess(mix, method = c("elir", "moment", "morita"), ..., s = 100)## S3 method for class 'gammaMix'ess(mix, method = c("elir", "moment", "morita"), ..., s = 100, eps = 1e-04)## S3 method for class 'normMix'ess(  mix,  method = c("elir", "moment", "morita"),  ...,  family = gaussian,  sigma,  s = 100)

Arguments

Details

The ESS is calculated using either the expected localinformation ratio (elir)Neuenschwander etal. (2020), the moments approach or the method byMorita et al. (2008).

The elir approach measures effective sample size in terms of theaverage curvature of the prior in relation to the Fisherinformation. Informally this corresponds to the average peakinessof the prior in relation to the information content of a singleobservation. The elir approach is the only ESS which fulfillspredictive consistency. The predictive consistency of the ESSrequires that the ESS of a prior is consistent when considering anaveraged posterior ESS of additional data distributed according tothe predictive distribution of the prior. The expectation of theposterior ESS is taken wrt to the prior predictive distribution andthe averaged posterior ESS corresponds to the sum of the prior ESSand the number of forward simulated data items. The elir approachresults in ESS estimates which are neither conservative nor liberalwhereas the moments method yields conservative and the moritamethod liberal results. See the example section for a demonstrationof predictive consistency.

For the moments method the mean and standard deviation of themixture are calculated and then approximated by the conjugatedistribution with the same mean and standard deviation. Forconjugate distributions, the ESS is well defined. See the examplesfor a step-wise calculation in the beta mixture case.

The Morita method used here evaluates the mixture prior at the modeinstead of the mean as proposed originally by Morita. The methodmay lead to very optimistic ESS values, especially if the mixturecontains many components. The calculation of the Morita approachhere follows the approach presented in Neuenschwander B. et all(2019) which avoids the need for a minimization and does notrestrict the ESS to be an integer.

The argumentssigma andfamily are specific fornormal mixture densities. These specify the sampling standarddeviation for agaussian family (the default) while alsoallowing to consider the ESS of standard one-parameter exponentialfamilies, i.e.binomial orpoisson. The functionsupports non-gaussian families with unit dispersion only.

Value

Returns the ESS of the prior as floating point number.

Methods (by class)

• ess(betaMix): ESS for beta mixtures.

• ess(gammaMix): ESS for gamma mixtures.

• ess(normMix): ESS for normal mixtures.

Supported Conjugate Prior-Likelihood Pairs

References

Morita S, Thall PF, Mueller P. Determining theeffective sample size of a parametric prior.Biometrics2008;64(2):595-602.

Neuenschwander B., Weber S., Schmidli H., O’HaganA. (2020). Predictively consistent prior effective samplesizes.Biometrics, 76(2),578–587. https://doi.org/10.1111/biom.13252

Examples

# Conjugate Beta examplea <- 5b <- 15prior <- mixbeta(c(1, a, b))ess(prior)(a + b)# Beta mixture examplebmix <- mixbeta(rob = c(0.2, 1, 1), inf = c(0.8, 10, 2))ess(bmix, "elir")ess(bmix, "moment")# moments method is equivalent to# first calculate momentsbmix_sum <- summary(bmix)# then calculate a and b of a matching betaab_matched <- ms2beta(bmix_sum["mean"], bmix_sum["sd"])# finally take the sum of a and b which are equivalent# to number of responders/non-responders respectivleyround(sum(ab_matched))ess(bmix, method = "morita")# One may also calculate the ESS on the logit scale, which# gives slightly different results due to the parameter# transformation, e.g.:prior_logit <- mixnorm(c(1, log(5 / 15), sqrt(1 / 5 + 1 / 15)))ess(prior_logit, family = binomial)bmix_logit <- mixnorm(rob = c(0.2, 0, 2), inf = c(0.8, log(10 / 2), sqrt(1 / 10 + 1 / 2)))ess(bmix_logit, family = binomial)# Predictive consistency of elirn_forward <- 1E1bmixPred <- preddist(bmix, n = n_forward)pred_samp <- rmix(bmixPred, 1E2)# use more samples here for greater accuracy, e.g.# pred_samp <- rmix(bmixPred, 1E3)pred_ess <- sapply(pred_samp, function(r) ess(postmix(bmix, r = r, n = n_forward), "elir"))ess(bmix, "elir")mean(pred_ess) - n_forward# Normal mixture examplenmix <- mixnorm(rob = c(0.5, 0, 2), inf = c(0.5, 3, 4), sigma = 10)ess(nmix, "elir")ess(nmix, "moment")# the reference scale determines the ESSsigma(nmix) <- 20ess(nmix)# we may also interpret normal mixtures as densities assigned to# parameters of a logit transformed response rate of a binomialnmix_logit <- mixnorm(c(1, logit(1 / 4), 2 / sqrt(10)))ess(nmix_logit, family = binomial)# Gamma mixture examplegmix <- mixgamma(rob = c(0.3, 20, 4), inf = c(0.7, 50, 10))ess(gmix) ## interpreted as appropriate for a Poisson likelihood (default)likelihood(gmix) <- "exp"ess(gmix) ## interpreted as appropriate for an exponential likelihood

Fill numeric objects

Description

Returns the numeric input object with the value given and respectsdimensionalty and type of input.

Usage

fill(x, value)

Arguments

Forest Plot

Description

Creates a forest plot forgMAP() analysis objects.

Usage

forest_plot(  x,  prob = 0.95,  est = c("both", "MAP", "Mean", "none"),  model = c("stratified", "both", "meta"),  point_est = c("median", "mean"),  size = 1.25,  alpha = 0.5)

Arguments

Details

The function creates a forest plot suitable forgMAP() analyses. Note that the Meta-Analytic-Predictiveprior is included by default in the plot as opposed to only showingthe estimated model mean. See the examples below to obtain standardforest plots.

Also note that the plot internally flips the x andy-axis. Therefore, if you want to manipulate the x-axis, you haveto give commands affecting the y-axis (see examples).

Value

The function returns aggplot2 plot object.

Customizingggplot2 plots

The returned plot is aggplot2 object. Please refer to the"Customizing Plots" vignette which is part ofRBesTdocumentation for an introduction. For simple modifications (changelabels, add reference lines, ...) consider the commands found inbayesplot-helpers. For more advancedcustomizations please use theggplot2 package directly. Adescription of the most common tasks can be found in theR Cookbook and a fullreference of available commands can be found at theggplot2 documentationsite.

See Also

gMAP()

Examples

# we consider the example AS MAP analysisexample(AS)# default forest plot for a gMAP analysisforest_plot(map_AS)# standard forest plot (only stratified estimate and Mean)forest_plot(map_AS, est = c("Mean"), model = "stratified")# to further customize these plots, first load bayesplot and ggplot2library(bayesplot)library(ggplot2)# to make plots with red colors, big fonts for presentations, suppress# the x axis label and add another title (with a subtitle)color_scheme_set("red")theme_set(theme_default(base_size = 16))forest_plot(map_AS, size = 2) +  yaxis_title(FALSE) +  ggtitle("Ankylosing Spondylitis Forest Plot",    subtitle = "Control Group Response Rate"  )# the defaults are set withcolor_scheme_set("blue")theme_set(theme_default(base_size = 12))

Meta-Analytic-Predictive Analysis for Generalized Linear Models

Description

Meta-Analytic-Predictive (MAP) analysis for generalized linearmodels suitable for normal, binary, or Poisson data. Modelspecification and overall syntax follows mainlystats::glm() conventions.

Usage

gMAP(  formula,  family = gaussian,  data,  weights,  offset,  tau.strata,  tau.dist = c("HalfNormal", "TruncNormal", "Uniform", "Gamma", "InvGamma", "LogNormal",    "TruncCauchy", "Exp", "Fixed"),  tau.prior,  tau.strata.pred = 1,  beta.prior,  prior_PD = FALSE,  REdist = c("normal", "t"),  t.df = 5,  contrasts = NULL,  iter = getOption("RBesT.MC.iter", 6000),  warmup = getOption("RBesT.MC.warmup", 2000),  thin = getOption("RBesT.MC.thin", 4),  init = getOption("RBesT.MC.init", 1),  chains = getOption("RBesT.MC.chains", 4),  cores = getOption("mc.cores", 1L))## S3 method for class 'gMAP'print(x, digits = 3, probs = c(0.025, 0.5, 0.975), ...)## S3 method for class 'gMAP'fitted(object, type = c("response", "link"), probs = c(0.025, 0.5, 0.975), ...)## S3 method for class 'gMAP'coef(object, probs = c(0.025, 0.5, 0.975), ...)## S3 method for class 'gMAP'as.matrix(x, ...)## S3 method for class 'gMAP'summary(  object,  type = c("response", "link"),  probs = c(0.025, 0.5, 0.975),  ...)

Arguments

Details

The meta-analytic-predictive (MAP) approach derives a prior fromhistorical data using a hierarchical model. The statistical model isformulated as a generalized linear mixed model for binary, normal(with fixed\sigma) and Poisson endpoints:

y_{ih}|\theta_{ih} \sim f(y_{ih} | \theta_{ih})

Here,i=1,\ldots,N is the index for observations, andh=1,\ldots,H is the index for the grouping (usually studies).The model assumes the linear predictor for a transformed mean as

g(\theta_{ih}; x_{ih},\beta) = x_{ih} \, \beta + \epsilon_h

withx_{ih} being the row vector ofk covariates forobservationi. The variance component is assumed by defaultnormal

\epsilon_h \sim N(0,\tau^2), \qquad h=1,\ldots,H

Lastly, the Bayesian implementation assumes independent normalpriors for thek regression coefficients and a prior for thebetween-group standard deviation\tau (seetaud.distfor available distributions).

The MAP prior will then be derived from the above model as theconditional distribution of\theta_{\star} given theavailable data and the vector of covariatesx_{\star}defining the overall intercept

\theta_{\star}| x_{\star},y .

A simple and common case arises for one observation (summarystatistic) per trial. For a normal endpoint, the model then simplifiesto the standard normal-normal hierarchical model. In the abovenotation,i=h=1,\ldots,H and

y_h|\theta_h \sim N(\theta_h,s_h^2)

\theta_h = \mu + \epsilon_h

\epsilon_h \sim N(0,\tau^2),

where the more common\mu is used for the only (intercept)parameter\beta_1. Since there are no covariates, the MAPprior is simplyPr(\theta_{\star} |y_1,\ldots,y_H).

The hierarchical model is a compromise between the two extremecases of full pooling (\tau=0, full borrowing, nodiscounting) and no pooling (\tau=\infty, no borrowing,stratification). The information content of thehistorical data grows with H (number of historical data items)indefinitely for full pooling whereas no information isgained in a stratified analysis. For a fixed\tau, the maximum effective samplesize of the MAP prior isn_\infty (H\rightarrow\infty), which for a normal endpoint with fixed\sigma is

n_\infty = \left(\frac{\tau^2}{\sigma^2}\right)^{-1},

(Neuenschwander et al., 2010). Hence, the ratio\tau/\sigma limits the amount of information a MAP prior isequivalent to. This allows for a classification of\tauvalues in relation to\sigma, which is crucial to define apriorP_\tau. The following classification is useful in aclinical trial setting:

The above formula forn_\infty assumes a known\tau. This is unrealistic as the between-trial heterogeneityparameter is often not well estimable, in particular if the numberof trials is small (H small). The above table helps to specify aprior distribution for\tau appropriate for the given contextwhich defines the crucial parameter\sigma. For binary andPoisson endpoints, normal approximations can be used to determine\sigma. See examples below for concrete cases.

The design matrixX is defined by the formula for the linearpredictor and is always of the formresponse ~ predictor | grouping, which followsstats::glm()conventions. The syntax has been extended to include aspecification of the grouping (for example study) factor of thedata with a horizontal bar,|. The bar separates theoptionally specified grouping level, i.e. in the binary endpointcasecbind(r, n-r) ~ 1 | study. By default it is assumedthat each row corresponds to an individual group (for which anindividual parameter is estimated). Specifics for the differentendpoints are:

normal

family=gaussian assumes an identity linkfunction. Theresponse should be given as matrix with twocolumns with the first column being the observed mean valuey_{ih} and the second column the standard errorse_{ih} (of the mean). Additionally, it is recommendedto specify with theweight argument the number of unitswhich contributed to the (mean) measurementy_{ih}. This information is used to estimate\sigma.

binary

family=binomial assumes a logit linkfunction. Theresponse must be given as two-column matrixwith number of respondersr (first column) and non-respondersn-r (second column).

Poisson

family=poisson assumes a log linkfunction. Theresponse is a vector of counts. The totalexposure times can be specified by anoffset, which will belinearly added to the linear predictor. Theoffset can begiven as part of the formula,y ~ 1 + offset(log(exposure))or as theoffset argument togMAP. Note that theexposure unit must be given as log-offset.

Value

The function returns a S3 object of typegMAP. Seethe methods section below for applicable functions to query theobject.

Methods (by generic)

• print(gMAP): displays a summary of the gMAP analysis.

• fitted(gMAP): returns the quantiles of the posterior shrinkageestimates for each data item used during the analysis of the givengMAP object.

• coef(gMAP): returns the quantiles of the predictivedistribution. User can choose withtype if the result is onthe response or the link scale.

• as.matrix(gMAP): extracts the posterior sample of the model.

• summary(gMAP): returns the summaries of a gMAP.analysis. Output is agMAPsummary object, which is a list containing

  tau

  posterior summary of the heterogeneity standard deviation

  beta

  posterior summary of the regression coefficients

  theta.pred

  summary of the predictive distribution (given in dependence on thetype argument either onresponse orlink scale)

  theta

  posterior summary of the mean estimate (also depends on thetype argument)

Differential Discounting

The above model assumes the same between-group standard deviation\tau, which implies that the data are equally relevant. Thisassumption can be relaxed to more than one\tau. That is,

\epsilon_h \sim N(0,\tau_{s(h)}^2)

wheres(h) assignes grouph to one ofSbetween-group heterogeneity strata.

For example, in a situation with two randomized and fourobservational studies, one may want to assume\tau_1 (fortrials 1 and 2) and\tau_2 (for trials 3-6) for thebetween-trial standard deviations of the control means. Moreheterogeneity (less relevance) for the observational studies canthen be expressed by appropriate priors for\tau_1 and\tau_2. In this case,S=2 and the strata assignments(seetau.strata argument) would bes(1)=s(2)=1,s(3)=\ldots=s(6)=2.

Prior Specification

The prior distribution for the regression coefficients\betais normal.

• If a single number is given, then this is used as the standarddeviation and the default mean of 0 is used.

• If a vector is given, it must be of the same lengthas number of covariates defined and is used as standarddeviation.

• If a matrix with a single row is given, its first row will beused as mean and the second row will be used as standard deviationfor all regression coefficients.

• Lastly, a two-column matrix (mean and standard deviation columns)with as many columns as regression coefficients can be given.

It is recommended to always specify abeta.prior. Perdefault a mean of 0 is set. The standard deviation is set to 2 forthe binary case, to 100 *sd(y) for the normal case and tosd(log(y + 0.5 + offset)) for the Poisson case.

For the between-trial heterogeniety\tau prior, a dispersionparameter must always be given for each exchangeabilitystratum. For the differenttau.prior distributions, twoparameters are needed out of which one is set to a default value ifapplicable:

For a prior distribution with a default location parameter, avector of length equal to the number of exchangability strata canbe given. Otherwise, a two-column matrix with as many rows asexchangability strata must be given, except for a single\taustratum, for which a vector of length two defines the parameters aand b.

Random seed

The MAP analysis is performed usingMarkov-Chain-Monte-Carlo (MCMC) inrstan::rstan(). MCMCis a stochastic algorithm. To obtain exactly reproducible resultsyou must use thebase::set.seed() functionbefore callinggMAP. SeeRBesT()overview page for global options on setting further MCMC simulationparameters.

References

Neuenschwander B, Capkun-Niggli G, Branson M,Spiegelhalter DJ. Summarizing historical information on controls inclinical trials.Clin Trials. 2010; 7(1):5-18

Schmidli H, Gsteiger S, Roychoudhury S, O'Hagan A, Spiegelhalter D,Neuenschwander B. Robust meta-analytic-predictive priors inclinical trials with historical control information.Biometrics 2014;70(4):1023-1032.

Weber S, Li Y, Seaman III J.W., Kakizume T, Schmidli H. ApplyingMeta-Analytic Predictive Priors with the R Bayesian evidencesynthesis tools.JSS 2021; 100(19):1-32

See Also

plot.gMAP(),forest_plot(),automixfit(),predict.gMAP()

Examples

## Setting up dummy sampling for fast execution of example## Please use 4 chains and 20x more warmup & iter in practice.user_mc_options <- options(RBesT.MC.warmup=50, RBesT.MC.iter=100,                            RBesT.MC.chains=2, RBesT.MC.thin=1)# Binary data example 1# Mean response rate is ~0.25. For binary endpoints# a conservative choice for tau is a HalfNormal(0,1) as long as# the mean response rate is in the range of 0.2 to 0.8. For# very small or large rates consider the n_infinity approach# illustrated below.# for exact reproducible results, the seed must be setset.seed(34563)map_AS <- gMAP(cbind(r, n - r) ~ 1 | study,  family = binomial,  data = AS,  tau.dist = "HalfNormal", tau.prior = 1,  beta.prior = 2)print(map_AS)# obtain numerical summariesmap_sum <- summary(map_AS)print(map_sum)names(map_sum)# [1] "tau"        "beta"       "theta.pred" "theta"map_sum$theta.pred# graphical model checks (returns list of ggplot2 plots)map_checks <- plot(map_AS)# forest plot with shrinkage estimatesmap_checks$forest_model# density of MAP prior on response scalemap_checks$densityThetaStar# density of MAP prior on link scalemap_checks$densityThetaStarLink# obtain shrinkage estimatesfitted(map_AS)# regression coefficientscoef(map_AS)# finally fit MAP prior with parametric mixturemap_mix <- mixfit(map_AS, Nc = 2)plot(map_mix)$mix# optionally select number of components automatically via AICmap_automix <- automixfit(map_AS)plot(map_automix)$mix# Normal example 2, see normal vignette# Prior considerations# The general principle to derive a prior for tau can be based on the# n_infinity concept as discussed in Neuenschwander et al., 2010.# This assumes a normal approximation which applies for the colitis# data set as:p_bar <- mean(with(colitis, r / n))s <- round(1 / sqrt(p_bar * (1 - p_bar)), 1)# s is the approximate sampling standard deviation and a# conservative prior is tau ~ HalfNormal(0,s/2)tau_prior_sd <- s / 2# Evaluate HalfNormal prior for tautau_cat <- c(  pooling = 0,  small = 0.0625,  moderate = 0.125,  substantial = 0.25,  large = 0.5,  veryLarge = 1,  stratified = Inf)# Interval probabilites (basically saying we are assuming# heterogeniety to be smaller than very large)diff(2 * pnorm(tau_cat * s, 0, tau_prior_sd))# Cumulative probabilities as 1-F1 - 2 * (pnorm(tau_cat * s, 0, tau_prior_sd) - 0.5)## Recover user set sampling defaultsoptions(.user_mc_options)

internal function used for integration of densities which appearsto be much more stable from -Inf to +Inf in the logit space whilethe density to be integrated recieves inputs from 0 to 1 such thatthe inverse distribution function must be used. The integral solvedis int_x dmix(mix,x) integrand(x) where integrand must be given aslog and we integrate over the support of mix.

Description

integrate density in logit space and split by component suchthat the quantile function of each component is used. Thisensures that the R implementation of the quantile function isalways used.

Usage

integrate_density_log(  log_integrand,  mix,  Lplower = -Inf,  Lpupper = Inf,  eps = getOption("RBesT.integrate_prob_eps", 1e-06))

Arguments

k nearest neighbor algorithm for multi-variate data

Description

k nearest neighbor algorithm for multi-variate data

Usage

knn(X, K = 2, init, Ninit = 50, verbose = FALSE, tol, Niter.max = 500)

Arguments

Read and Set Likelihood to the Corresponding Conjugate Prior

Description

Read and set the likelihood distribution corresponding to the conjugate prior distribution.

Usage

likelihood(mix)likelihood(mix) <- value

Arguments

Details

If the prior and posterior distributions are in the same family, then the prior distributionis called a conjugate prior for the likelihood function.

Supported Conjugate Prior-Likelihood Pairs

Examples

# Gamma mixturegmix <- mixgamma(c(0.3, 20, 4), c(0.7, 50, 10))# read out conjugate partnerlikelihood(gmix)ess(gmix)# set conjugate partnerlikelihood(gmix) <- "exp"# ... which changes the interpretation of the mixtureess(gmix)

Logit (log-odds) and inverse-logit function.

Description

Calculates the logit (log-odds) and inverse-logit.

Usage

logit(mu)inv_logit(eta)

Arguments

Details

Values of mu equal to 0 or 1 will return-\infty or\infty respectively.

Value

A numeric object of the same type as mu and eta containingthe logits or inverse logit of the input values. The logit andinverse transformation equates to

\text{logit}(\mu) = \log(\mu/(1-\mu))

\text{logit}^{-1}(\eta)= \exp(\eta)/(1 + \exp(\eta)).

Examples

logit(0.2)inv_logit(-1.386)

Extract log likelihood from fitted EM objects

Description

Extract log likelihood from fitted EM objects

Usage

## S3 method for class 'EM'logLik(object, ...)

Numerically stable log of the inv_logit function

Description

Numerically stable log of the inv_logit function

Usage

log_inv_logit(mat)

Mixture Distributions

Description

Density, cumulative distribution function, quantilefunction and random number generation for supported mixturedistributions. (d/p/q/r)mix are generic and work with any mixturesupported by BesT (see table below).

Usage

dmix(mix, x, log = FALSE)pmix(mix, q, lower.tail = TRUE, log.p = FALSE)qmix(mix, p, lower.tail = TRUE, log.p = FALSE)rmix(mix, n)## S3 method for class 'mix'mix[[..., rescale = FALSE]]

Arguments

Details

A mixture distribution is defined as a linearsuperposition ofK densities of the same distributionalclass. The mixture distributions supported have the form

f(x,\mathbf{w},\mathbf{a},\mathbf{b}) = \sum_{k=1}^K w_k \, f_k(x,a_k,b_k).

Thew_k are the mixing coefficients which must sum to1. Moreover, each densityf is assumed to beparametrized by two parameters such that each componentk isdefined by a triplet,(w_k,a_k,b_k).

Individual mixture components can be extracted using the[[operator, see examples below.

The supported densities are normal, beta and gamma which can beinstantiated withmixnorm(),mixbeta(), ormixgamma(), respectively. In addition, the respectivepredictive distributions are supported. These can be obtained bycallingpreddist() which returns appropriate normal,beta-binomial or Poisson-gamma mixtures.

For convenience asummary function is defined for allmixtures. It returns the mean, standard deviation and the requestedquantiles which can be specified with the argumentprobs.

Value

dmix gives the weighted sum of the densities of eachcomponent.

pmix calculates the distribution function byevaluating the weighted sum of each components distributionfunction.

qmix returns the quantile for the givenpby using that the distribution function is monotonous and hence agradient based minimization scheme can be used to find the matchingquantileq.

rmix generates a random sample of sizen by first sampling a latent component indicator in therange1..K for each draw and then the function samples fromeach component a random draw using the respective samplingfunction. Thernorm function returns the random draws asnumerical vector with an additional attributeind whichgives the sampled component indicator.

Supported Conjugate Prior-Likelihood Pairs

See Also

plot.mix()

Other mixdist:mixbeta(),mixcombine(),mixgamma(),mixjson,mixmvnorm(),mixnorm(),mixplot

Examples

## a beta mixturebm <- mixbeta(weak = c(0.2, 2, 10), inf = c(0.4, 10, 100), inf2 = c(0.4, 30, 80))## extract the two most informative componentsbm[[c(2, 3)]]## rescaling needed in order to plotplot(bm[[c(2, 3), rescale = TRUE]])summary(bm)

Beta Mixture Density

Description

The Beta mixture density and auxilary functions.

Usage

mixbeta(..., param = c("ab", "ms", "mn"))ms2beta(m, s, drop = TRUE)mn2beta(m, n, drop = TRUE)## S3 method for class 'betaMix'print(x, ...)## S3 method for class 'betaBinomialMix'print(x, ...)## S3 method for class 'betaMix'summary(object, probs = c(0.025, 0.5, 0.975), ...)## S3 method for class 'betaBinomialMix'summary(object, probs = c(0.025, 0.5, 0.975), ...)

Arguments

Details

Each entry in the... argument list is expected tobe a triplet of numbers which defines the weightw_k, firstand second parameter of the mixture componentk. A tripletcan optionally be named which will be used appropriately.

The first and second parameter can be given in differentparametrizations which is set by theparam option:

ab

Natural parametrization of Beta density (a=shape1 andb=shape2). Default.

ms

Mean and standard deviation,m=a/(a+b) ands=\sqrt{\frac{m(1-m)}{1+n}}, wheren=a+b is the number of observations. Note thats must be less than\sqrt{m(1-m)}.

mn

Mean and number of observations,n=a+b.

Value

mixbeta returns a beta mixture with the specified mixture components.ms2beta andmn2beta return the equivalent naturala andb parametrization given parametersm,s, orn.

See Also

Other mixdist:mix,mixcombine(),mixgamma(),mixjson,mixmvnorm(),mixnorm(),mixplot

Examples

## a beta mixturebm <- mixbeta(rob = c(0.2, 2, 10), inf = c(0.4, 10, 100), inf2 = c(0.4, 30, 80))# mean/standard deviation parametrizationbm2 <- mixbeta(rob = c(0.2, 0.3, 0.2), inf = c(0.8, 0.4, 0.01), param = "ms")# mean/observations parametrizationbm3 <- mixbeta(rob = c(0.2, 0.3, 5), inf = c(0.8, 0.4, 30), param = "mn")# even mixed is possiblebm4 <- mixbeta(rob = c(0.2, mn2beta(0.3, 5)), inf = c(0.8, ms2beta(0.4, 0.1)))# print methods are definedbm4print(bm4)

Combine Mixture Distributions

Description

Combining mixture distributions of the same class to form a new mixture.

Usage

mixcombine(..., weight, rescale = TRUE)

Arguments

Details

Combines mixtures of the same class of randomvariable to form a new mixture distribution.

Value

A R-object with the new mixture distribution.

See Also

robustify()

Other mixdist:mix,mixbeta(),mixgamma(),mixjson,mixmvnorm(),mixnorm(),mixplot

Examples

# beta with two informative componentsbm <- mixbeta(inf = c(0.5, 10, 100), inf2 = c(0.5, 30, 80))# robustified with mixcombine, i.e. a 10% uninformative part addedunif <- mixbeta(rob = c(1, 1, 1))mixcombine(bm, unif, weight = c(9, 1))

Difference of mixture distributions

Description

Density, cumulative distribution function, quantilefunction and random number generation for the difference of two mixturedistributions.

Usage

dmixdiff(mix1, mix2, x)pmixdiff(mix1, mix2, q, lower.tail = TRUE)qmixdiff(mix1, mix2, p, lower.tail = TRUE)rmixdiff(mix1, mix2, n)

Arguments

Details

Ifx_1 \sim f_1(x_1) andx_2 \simf_2(x_2), the density of the differenced\equiv x_1 - x_2 is given by

f_d(d) = \int f_1(u) \, f_2(u - d) \, du.

The cumulative distribution function equates to

F_d(d) = \int f_1(u) \, (1-F_2(u-d)) \, du.

Both integrals are performed over the full support of thedensities and use the numerical integration functionintegrate().

Value

Respective density, quantile, cumulative density or randomnumbers.

Examples

# 1. Difference between two beta distributions, i.e. Pr( mix1 - mix2 > 0)mix1 <- mixbeta(c(1, 11, 4))mix2 <- mixbeta(c(1, 8, 7))pmixdiff(mix1, mix2, 0, FALSE)# Interval probability, i.e. Pr( 0.3 > mix1 - mix2 > 0)pmixdiff(mix1, mix2, 0.3) - pmixdiff(mix1, mix2, 0)# 2. two distributions, one of them a mixturem1 <- mixbeta(c(1, 30, 50))m2 <- mixbeta(c(0.75, 20, 50), c(0.25, 1, 1))# random sample of differenceset.seed(23434)rM <- rmixdiff(m1, m2, 1E4)# histogram of random numbers and exact densityhist(rM, prob = TRUE, new = TRUE, nclass = 40)curve(dmixdiff(m1, m2, x), add = TRUE, n = 51)# threshold probabilities for difference, at 0 and 0.2pmixdiff(m1, m2, 0)mean(rM < 0)pmixdiff(m1, m2, 0.2)mean(rM < 0.2)# median of differencemdn <- qmixdiff(m1, m2, 0.5)mean(rM < mdn)# 95%-intervalqmixdiff(m1, m2, c(0.025, 0.975))quantile(rM, c(0.025, 0.975))

Utility function to instantiate 2 parameter mixture densities.

Description

Utility function to instantiate 2 parameter mixture densities.

Usage

mixdist3(...)

Fit of Mixture Densities to Samples

Description

Expectation-Maximization (EM) based fitting of parametric mixturedensities to numerical samples. This provides a convenient approachto approximate MCMC samples with a parametric mixture distribution.

Usage

mixfit(sample, type = c("norm", "beta", "gamma", "mvnorm"), thin, ...)## Default S3 method:mixfit(sample, type = c("norm", "beta", "gamma", "mvnorm"), thin, ...)## S3 method for class 'gMAP'mixfit(sample, type, thin, ...)## S3 method for class 'gMAPpred'mixfit(sample, type, thin, ...)## S3 method for class 'array'mixfit(sample, type, thin, ...)

Arguments

Details

Parameters of EM fitting functions

Nc

Number of mixture components. Required parameter.

mix_init

Initial mixture density. If missing (default) then a k-nearest-neighbor algorithm is used to find an initial mixture density.

Ninit

Number of data points used for initialization. Defaults to 50.

verbose

If set toTRUE the function will inform about fitting process

maxIter

Maximal number of iterations. Defaults to 500.

tol

Defines a convergence criteria as an upper bound for the change in the log-likelihood, i.e. once the derivative (with respect to iterations) of the log-likelihood falls belowtol, the function declares convergence and stops.

eps

Must be a triplet of numbers which set the desired accuracy of the inferred parameters per mixture component. See below for a description of the parameters used during EM. EM is stopped once a running mean of the absolute difference between the last successiveNeps estimates is below the giveneps for all parameters. Defaults to 5E-3 for each parameter.

Neps

Number of iterations used for the running mean of parameter estimates to test for convergence. Defaults to 5.

constrain_gt1

Logical value controlling if the Beta EM constrains all parameters a & b to be greater than 1. By default constraints are turned on (new since 1.6-0).

By default the EM convergence is declared whenthe desired accuracy of the parameters has been reached over the lastNeps estimates. Iftol andNeps is specified, thenwhatever criterion is met first will stop the EM.

The parameters per componentk used internally during fittingare for the different EM procedures:

normal

logit(w_k), \mu_k, \log(\sigma_k)

beta

logit(w_k), \log(a_k), \log(b_k)

constrained beta

logit(w_k), \log(a_k-1), \log(b_k-1)

gamma

logit(w_k), \log(\alpha_k), \log(\beta_k)

Note: Whenever nomix_init argument is given,the EM fitting routines assume that the data vector is given inrandom order. If in the unlikely event that the EM gets caught in alocal extremum, then random reordering of the data vector mayalleviate the issue.

Value

A mixture object according the requestedtype isreturned. The object has additional information attached, i.e. thelog-likelihood can be queried and diagnostic plots can begenerated. See links below.

Methods (by class)

• mixfit(default): Performs an EM fit for the givensample. Thinning is applied only if thin is specified.

• mixfit(gMAP): Fits the default predictive distribution from agMAP analysis. Automatically obtains the predictive distribution ofthe intercept only case on the response scale mixture from thegMAP() object. For the binomial case a beta mixture,for the gaussian case a normal mixture and for the Poisson case agamma mixture will be used. In the gaussian case, the resultingnormal mixture will set the reference scale to the estimatedsigma ingMAP() call.

• mixfit(gMAPpred): Fits a mixture density for each prediction fromthegMAP() prediction.

• mixfit(array): Fits a mixture density for an MCMC sample. It isrecommended to provide a thinning argument which roughly yieldsindependent draws (i.e. useacf() to identify athinning lag with small auto-correlation). The input array isexpected to have 3 dimensions which are nested as iterations,chains, and draws.

References

Dempster A.P., Laird N.M., Rubin D.B. MaximumLikelihood from Incomplete Data via the EM Algorithm.Journalof the Royal Statistical Society, Series B 1977; 39 (1): 1-38.

See Also

Other EM:plot.EM()

Examples

bmix <- mixbeta(rob = c(0.2, 1, 1), inf = c(0.8, 10, 2))bsamp <- rmix(bmix, 1000)bfit <- mixfit(bsamp, type = "beta", Nc = 2)# diagnostic plots can easily by generated from the EM fit withbfit.check <- plot(bfit)names(bfit.check)# check convergence of parametersbfit.check$mixbfit.check$mixdensbfit.check$mixecdf# obtain the log-likelihoodlogLik(bfit)# or AICAIC(bfit)

The Gamma Mixture Distribution

Description

The gamma mixture density and auxiliary functions.

Usage

mixgamma(..., param = c("ab", "ms", "mn"), likelihood = c("poisson", "exp"))ms2gamma(m, s, drop = TRUE)mn2gamma(m, n, likelihood = c("poisson", "exp"), drop = TRUE)## S3 method for class 'gammaMix'print(x, ...)## S3 method for class 'gammaPoissonMix'print(x, ...)## S3 method for class 'gammaExpMix'print(x, ...)## S3 method for class 'gammaMix'summary(object, probs = c(0.025, 0.5, 0.975), ...)## S3 method for class 'gammaPoissonMix'summary(object, probs = c(0.025, 0.5, 0.975), ...)

Arguments

Details

Each entry in the... argument list is expected tobe a triplet of numbers which defines the weightw_k, firstand second parameter of the mixture componentk. A tripletcan optionally be named which will be used appropriately.

The first and second parameter can be given in differentparametrizations which is set by theparam option:

ab

Natural parametrization of Gamma density (a=shape andb=rate). Default.

ms

Mean and standard deviation,m=a/b ands=\sqrt{a}/b.

mn

Mean and number of observations. Translation to naturalparameter depends on thelikelihood argument. For a Poissonlikelihoodn=b (anda=m \cdot n), for an Explikelihoodn=a (andb=n/m).

Value

mixgamma returns a gamma mixture with the specified mixture components.ms2gamma andmn2gamma return the equivalent naturala andb parametrization givenparametersm,s, orn.

See Also

Other mixdist:mix,mixbeta(),mixcombine(),mixjson,mixmvnorm(),mixnorm(),mixplot

Examples

# Gamma mixture with robust and informative componentgmix <- mixgamma(rob = c(0.3, 20, 4), inf = c(0.7, 50, 10))# objects can be printedgmix# or explicitlyprint(gmix)# summaries are definedsummary(gmix)# sub-components may be extracted# by component numbergmix[[2]]# or component namegmix[["inf"]]# alternative mean and standard deviation parametrizationgmsMix <- mixgamma(rob = c(0.5, 8, 0.5), inf = c(0.5, 9, 2), param = "ms")# or mean and number of observations parametrizationgmnMix <- mixgamma(rob = c(0.2, 2, 1), inf = c(0.8, 2, 5), param = "mn")# and mixed parametrizations are also possiblegfmix <- mixgamma(rob1 = c(0.15, mn2gamma(2, 1)),                  rob2 = c(0.15, ms2gamma(2, 5)),                  inf = c(0.7, 50, 10))

Write and Read a Mixture Object with JSON

Description

These functions write and read a mixture object in the JSON format.

Usage

write_mix_json(mix, con, ...)read_mix_json(con, ..., rescale = TRUE)

Arguments

Details

The mixture objects are written or read from the connectioncon, which can be a character string specifying a filepath or a connection object as detailed inbase::connections().

When writing mixture objects as JSON it is strongly recommended toexplicitly set the number of digits (argumentdigits) to beused for the numerical representation in order to control theaccuracy of the JSON representation of the mixture object. If themixture object inherits from the"EM" class (as is thecase when the mixture is created using themixfit() function), then the mixture objectwill be cast to a simple mixture object such that diagnosticsfrom the"EM" fitting procedure are dropped from theobject. For easier readability the user is encouraged to set theargumentpretty=TRUE, which is passed to thejsonlite::toJSON() function and makes the output morehuman readable.

Note that when reading in mixture objects, then these are notnecessarily equal to the mixtures passed to thewrite_mix_jsonfunction. This is a consequence of the limited precision of thetextual representation as defined by thedigits argument.

Value

Thewrite_mix_json function does not return a value whiletheread_mix_json returns the mixture object stored in theconnection specified.

See Also

Other mixdist:mix,mixbeta(),mixcombine(),mixgamma(),mixmvnorm(),mixnorm(),mixplot

Examples

## Not run: nm <- mixnorm(rob = c(0.2, 0, 2), inf = c(0.8, 2, 2), sigma = 5)write_mix_json(nm, "normal_mixture.json", pretty=TRUE, digits=1)mix <- read_mix_json("normal_mixture.json")## End(Not run)

takes x and transforms it according to the defined link function ofthe mixture

Description

takes x and transforms it according to the defined link function ofthe mixture

Usage

mixlink(mix, x)

Multivariate Normal Mixture Density

Description

The multivariate normal mixture density and auxiliaryfunctions.

Usage

mixmvnorm(..., sigma, param = c("ms", "mn", "msr"))msr2mvnorm(m, s, r, unlist = TRUE)## S3 method for class 'mvnormMix'print(x, ...)## S3 method for class 'mvnormMix'summary(object, ...)## S3 method for class 'mvnormMix'sigma(object, ...)

Arguments

Details

Each entry in the... argument list is a numericvector defining one component of the mixture multivariatenormal distribution. The first entry of the component definingvector is the weight of the mixture component followed by thevector of means in each dimension and finally a specificationof the covariance matrix, which depends on the chosenparametrization. The covariance matrix is expected to be givenas numeric vector in a column-major format, which is standardconversion applied to matrices by the vector concatenationfunctionbase::c(). Please refer to the examplessection below.

Each component defining vector can be specified in different waysas determined by theparam option:

ms

Mean vector and covariance matrixs. Default.

mn

Mean vector and number of observations.n determinesthe covariance for each component via the relation\Sigma/nwith\Sigma being the known reference covariance.

msr

Mean vector, standard deviations and correlations incolumn-major format (corresponds to order when printing multi-variatenormal mixtures).

The reference covariance\Sigma is the known covariance inthe normal-normal model (observation covariance). The functionsigma can be used to query the reference covariance and mayalso be used to assign a new reference covariance, see examplesbelow. In casesigma is not specified, the user has tosupplysigma as argument to functions which require areference covariance.

Value

Returns a multivariate normal mixture with the specifiedmixture components.

See Also

Other mixdist:mix,mixbeta(),mixcombine(),mixgamma(),mixjson,mixnorm(),mixplot

Examples

# default mean & covariance parametrizationS <- diag(c(1, 2)) %*% matrix(c(1, 0.5, 0.5, 1), 2, 2) %*% diag(c(1, 2))mvnm1 <- mixmvnorm(  rob = c(0.2, c(0, 0), diag(c(2, 2)^2)),  inf = c(0.8, c(0.5, 1), S / 4), sigma = S)print(mvnm1)summary(mvnm1)set.seed(657846)mixSamp1 <- rmix(mvnm1, 500)colMeans(mixSamp1)# alternative mean, sd and correlation parametrizationmvnm1_alt <- mixmvnorm(  rob = c(0.2, c(0, 0), c(2, 2), 0.0),  inf = c(0.8, c(0.5, 1), c(1, 2) / 2, 0.5),  sigma = msr2mvnorm(s = c(1, 2), r = 0.5, unlist = FALSE)$s,  param = "msr")print(mvnm1_alt)

Normal Mixture Density

Description

The normal mixture density and auxiliary functions.

Usage

mixnorm(..., sigma, param = c("ms", "mn"))mn2norm(m, n, sigma, drop = TRUE)## S3 method for class 'normMix'print(x, ...)## S3 method for class 'normMix'summary(object, probs = c(0.025, 0.5, 0.975), ...)## S3 method for class 'normMix'sigma(object, ...)sigma(object) <- value

Arguments

Details

Each entry in the... argument list is expected tobe a triplet of numbers which defines the weightw_k, firstand second parameter of the mixture componentk. A tripletcan optionally be named which will be used appropriately.

The first and second parameter can be given in differentparametrizations which is set by theparam option:

ms

Mean and standard deviation. Default.

mn

Mean and number of observations.n determiness via the relations=\sigma/\sqrt{n} with\sigma being the fixed reference scale.

The reference scale\sigma is the fixed standard deviation inthe one-parameter normal-normal model (observation standarddeviation). The functionsigma can be used to query thereference scale and may also be used to assign a new referencescale, see examples below. In case thesigma is notspecified, the user has to supplysigma as argument tofunctions which require a reference scale.

Value

Returns a normal mixture with the specified mixturecomponents.mn2norm returns the mean and standard deviationgiven a mean and sample size parametrization.

Functions

• sigma(object) <- value: Allows to assign a new reference scalesigma.

See Also

Other mixdist:mix,mixbeta(),mixcombine(),mixgamma(),mixjson,mixmvnorm(),mixplot

Examples

nm <- mixnorm(rob = c(0.2, 0, 2), inf = c(0.8, 2, 2), sigma = 5)print(nm)summary(nm)plot(nm)set.seed(57845)mixSamp <- rmix(nm, 500)plot(nm, samp = mixSamp)# support defined by quantilesqmix(nm, c(0.01, 0.99))# density functiondmix(nm, seq(-5, 5, by = 2))# distribution functionpmix(nm, seq(-5, 5, by = 2))# the reference scale can be changed (it determines the ESS)ess(nm)sigma(nm) <- 10ess(nm)

Plot mixture distributions

Description

Plotting for mixture distributions

Usage

## S3 method for class 'mix'plot(x, prob = 0.99, fun = dmix, log = FALSE, comp = TRUE, size = 1.25, ...)## S3 method for class 'mvnormMix'plot(x, prob = 0.99, fun = dmix, log = FALSE, comp = TRUE, size = 1.25, ...)

Arguments

Details

Plot function for mixture distribution objects. It showsthe density/quantile/cumulative distribution (corresponds tod/q/pmix function) for some specific central probabilitymass defined byprob. By default the x-axis is chosen toshow 99\% of the probability density mass.

Value

Aggplot2::ggplot() object is returned.

Customizingggplot2 plots

The returned plot is aggplot2 object. Please refer to the"Customizing Plots" vignette which is part ofRBesTdocumentation for an introduction. For simple modifications (changelabels, add reference lines, ...) consider the commands found inbayesplot-helpers. For more advancedcustomizations please use theggplot2 package directly. Adescription of the most common tasks can be found in theR Cookbook and a fullreference of available commands can be found at theggplot2 documentationsite.

See Also

Other mixdist:mix,mixbeta(),mixcombine(),mixgamma(),mixjson,mixmvnorm(),mixnorm()

Examples

# beta with two informative componentsbm <- mixbeta(inf = c(0.5, 10, 100), inf2 = c(0.5, 30, 80))plot(bm)plot(bm, fun = pmix)# for customizations of the plot we need to load ggplot2 firstlibrary(ggplot2)# show a histogram along with the densityplot(bm) + geom_histogram(  data = data.frame(x = rmix(bm, 1000)),  aes(y = ..density..), bins = 50, alpha = 0.4)# note: we can also use bayesplot for histogram plots with a density ...library(bayesplot)mh <- mcmc_hist(data.frame(x = rmix(bm, 1000)), freq = FALSE) +  overlay_function(fun = dmix, args = list(mix = bm))# ...and even add each componentfor (k in 1:ncol(bm)) {  mh <- mh + overlay_function(fun = dmix, args = list(mix = bm[[k]]), linetype = I(2))}print(mh)# normal mixturenm <- mixnorm(rob = c(0.2, 0, 2), inf = c(0.8, 6, 2), sigma = 5)plot(nm)plot(nm, fun = qmix)# obtain ggplot2 object and change titlepl <- plot(nm)pl + ggtitle("Normal 2-Component Mixture")

Mixture distributions asbrms priors

Description

Adapter function converting mixture distributions foruse withbrms::brm() models via thebrms::stanvar() facility.

Usage

mixstanvar(..., verbose = FALSE)

Arguments

Details

To declare mixture priors in abrms::brm()model requires two steps: First, the mixture densities need tobe converted by the adapter functionmixstanvar into astanvars object which is passed to thestanvarsargument of thebrms::brm() function. Doing soextends the Stan code and data generated bybrms::brm() such that the mixture densities can beused as priors within thebrms::brm() model. Thesecond step is to assign parameters of thebrms::brm() model to the mixture density as priorusing thebrms::set_prior() command ofbrms.

The adapter function translates the mixture distributions asdefined inR to the respective mixture distribution inStan. Within Stan the mixture distributions are named inaccordance to theR functions used to create the respectivemixture distributions. That is, a mixture density of normalscreated bymixnorm() is referred to asmixnorm_lpdf in Stan such that one can refer to the densityasmixnorm within thebrms::set_prior()functions (the suffix⁠_lpdf⁠ is automatically added bybrms::brm()). The arguments to these mixturedistributions depend on the specific distribution type as follows:

These arguments to the mixture densities refer to the differentdensity parameters and are automatically extracted from themixtures when converted. Important here are the argument names asthese must be used to declare the mixture prior. For each mixtureto convert as part of the... argument tomixstanvara label is determined using the name given in the list. In case noname is given, then the name of theR object itself isused. To declare a prior for a parameter the mixture distributionmust be used with it's arguments following the conventionlabel_argument. Please refer to the examples section for anillustration.

Note: Models created bybrms::brm() do use bydefault a data-dependent centering of covariates leading to a shiftof the overall intercept. This is commonly not desirable inapplications of this functionality. It is therefore stronglyrecommended to pass the optioncenter=FALSE as argument tothebrms formula created with thebrms::bf()function as demonstrated with the example below.

Value

stanvars object to be used as argument to thestanvars argument of abrms::brm() model.

Examples

## Not run: # The mixstanvar adapter requires the optional packages brms and gluestopifnot(require("brms"), require("glue"))# Assume we prefer a logistic regression MCMC analysis rather than a# beta-binomial analysis for the responder endpoint of the ankylosing# spondylitis (AS) example. Reasons to prefer a regression analysis is# to allow for baseline covariate adjustments, for example.map_AS_beta <- mixbeta(c(0.62, 19.2, 57.8), c(0.38, 3.5, 9.4))# First we need to convert the beta mixture to a respective mixture on# the log odds scale and approximate it with a normal mixture density.map_AS_samp <- rmix(map_AS_beta, 1E4)map_AS <- mixfit(logit(map_AS_samp), type = "norm", Nc = 2)# Trial results for placebo and secukinumab.trial <- data.frame(  n = c(6, 24),  r = c(1, 15),  arm = factor(c("placebo", "secukinumab")))# Define brms model such that the overall intercept corresponds to the# placebo response rate on the logit scale. NOTE: The use of# center=FALSE is required here as detailed in the note above.model <- bf(r | trials(n) ~ 1 + arm, family = binomial, center = FALSE)# to obtain detailed information on the declared model parameters use# get_prior(model, data=trial)# declare model prior with reference to mixture normal map prior...model_prior <- prior(mixnorm(map_w, map_m, map_s), coef = Intercept) +  prior(normal(0, 2), class = b)# ... which must be made available to brms using the mixstanvar adapter.# Note that the map_AS prior is labeled "map" as referred to in the# previous prior declaration.analysis <- brm(model,  data = trial, prior = model_prior,  stanvars = mixstanvar(map = map_AS),  seed = 365634, refresh = 0)# Let's compare the logistic regression estimate for the probability# of a positive treatment effect (secukinumab response rate exceeding# the response rate of placebo) to the direct beta-binomial analysis:hypothesis(analysis, "armsecukinumab > 0")post_secukinumab <- postmix(mixbeta(c(1, 0.5, 1)), r = 15, n = 24)post_placebo <- postmix(map_AS_beta, r = 1, n = 6)pmixdiff(post_secukinumab, post_placebo, 0, lower.tail = FALSE)# The posterior probability for a positive treatment effect# is very close to unity in both cases.## End(Not run)

Operating Characteristics for 1 Sample Design

Description

Theoc1S function defines a 1 sample design (prior, samplesize, decision function) for the calculation of the frequency atwhich the decision is evaluated to 1 conditional on assumingknown parameters. A function is returned which performs the actualoperating characteristics calculations.

Usage

oc1S(prior, n, decision, ...)## S3 method for class 'betaMix'oc1S(prior, n, decision, ...)## S3 method for class 'normMix'oc1S(prior, n, decision, sigma, eps = 1e-06, ...)## S3 method for class 'gammaMix'oc1S(prior, n, decision, eps = 1e-06, ...)

Arguments

Details

Theoc1S function defines a 1 sample design andreturns a function which calculates its operatingcharacteristics. This is the frequency with which the decisionfunction is evaluated to 1 under the assumption of a given truedistribution of the data defined by a known parameter\theta. The 1 sample design is defined by the prior, thesample size and the decision function,D(y). These uniquelydefine the decision boundary, seedecision1S_boundary().

When calling theoc1S function, then internally the criticalvaluey_c (usingdecision1S_boundary()) iscalculated and a function is returns which can be used tocalculated the desired frequency which is evaluated as

F(y_c|\theta).

Value

Returns a function with one argumenttheta whichcalculates the frequency at which the decision function isevaluated to 1 for the defined 1 sample design. Note that thereturned function takes vectors arguments.

Methods (by class)

• oc1S(betaMix): Applies for binomial model with a mixturebeta prior. The calculations use exact expressions.

• oc1S(normMix): Applies for the normal model with knownstandard deviation\sigma and a normal mixture prior for themean. As a consequence from the assumption of a known standarddeviation, the calculation discards sampling uncertainty of thesecond moment. The functionoc1S has an extraargumenteps (defaults to10^{-6}). The critical valuey_c is searched in the region of probability mass1-eps fory.

• oc1S(gammaMix): Applies for the Poisson model with a gammamixture prior for the rate parameter. The functionoc1S takes an extra argumenteps (defaults to10^{-6})which determines the region of probability mass1-eps wherethe boundary is searched fory.

See Also

Other design1S:decision1S(),decision1S_boundary(),pos1S()

Examples

# non-inferiority example using normal approximation of log-hazard# ratio, see ?decision1S for all detailss <- 2flat_prior <- mixnorm(c(1, 0, 100), sigma = s)nL <- 233theta_ni <- 0.4theta_a <- 0alpha <- 0.05beta <- 0.2za <- qnorm(1 - alpha)zb <- qnorm(1 - beta)n1 <- round((s * (za + zb) / (theta_ni - theta_a))^2)theta_c <- theta_ni - za * s / sqrt(n1)# standard NI designdecA <- decision1S(1 - alpha, theta_ni, lower.tail = TRUE)# double criterion design# statistical significance (like NI design)dec1 <- decision1S(1 - alpha, theta_ni, lower.tail = TRUE)# require mean to be at least as good as theta_cdec2 <- decision1S(0.5, theta_c, lower.tail = TRUE)# combinationdecComb <- decision1S(c(1 - alpha, 0.5), c(theta_ni, theta_c), lower.tail = TRUE)theta_eval <- c(theta_a, theta_c, theta_ni)# evaluate different designs at two sample sizesdesignA_n1 <- oc1S(flat_prior, n1, decA)designA_nL <- oc1S(flat_prior, nL, decA)designC_n1 <- oc1S(flat_prior, n1, decComb)designC_nL <- oc1S(flat_prior, nL, decComb)# evaluate designs at the key log-HR of positive treatment (HR<1),# the indecision point and the NI margindesignA_n1(theta_eval)designA_nL(theta_eval)designC_n1(theta_eval)designC_nL(theta_eval)# to understand further the dual criterion design it is useful to# evaluate the criterions separatley:# statistical significance criterion to warrant NI...designC1_nL <- oc1S(flat_prior, nL, dec1)# ... or the clinically determined indifference pointdesignC2_nL <- oc1S(flat_prior, nL, dec2)designC1_nL(theta_eval)designC2_nL(theta_eval)# see also ?decision1S_boundary to see which of the two criterions# will drive the decision

Operating Characteristics for 2 Sample Design

Description

Theoc2S function defines a 2 sample design (priors, samplesizes & decision function) for the calculation of operatingcharaceristics. A function is returned which calculates thecalculates the frequency at which the decision function isevaluated to 1 when assuming known parameters.

Usage

oc2S(prior1, prior2, n1, n2, decision, ...)## S3 method for class 'betaMix'oc2S(prior1, prior2, n1, n2, decision, eps, ...)## S3 method for class 'normMix'oc2S(  prior1,  prior2,  n1,  n2,  decision,  sigma1,  sigma2,  eps = 1e-06,  Ngrid = 10,  ...)## S3 method for class 'gammaMix'oc2S(prior1, prior2, n1, n2, decision, eps = 1e-06, ...)

Arguments

Details

Theoc2S function defines a 2 sample design andreturns a function which calculates its operatingcharacteristics. This is the frequency with which the decisionfunction is evaluated to 1 under the assumption of a given truedistribution of the data defined by the known parameter\theta_1 and\theta_2. The 2 sample design is definedby the priors, the sample sizes and the decision function,D(y_1,y_2). These uniquely define the decision boundary , seedecision2S_boundary().

Calling theoc2S function calculates the decision boundaryD_1(y_2) (seedecision2S_boundary()) and returnsa function which can be used to calculate the desired frequencywhich is evaluated as

\int f_2(y_2|\theta_2) F_1(D_1(y_2)|\theta_1) dy_2.

See below for examples and specifics for the supported mixturepriors.

Value

Returns a function which when called with two argumentstheta1 andtheta2 will return the frequencies atwhich the decision function is evaluated to 1 whenever the data isdistributed according to the known parameter values in eachsample. Note that the returned function takes vector arguments.

Methods (by class)

• oc2S(betaMix): Applies for binomial model with a mixturebeta prior. The calculations use exact expressions. If theoptional argumenteps is defined, then an approximate methodis used which limits the search for the decision boundary to theregion of1-eps probability mass. This is useful for designswith large sample sizes where an exact approach is very costly tocalculate.

• oc2S(normMix): Applies for the normal model with knownstandard deviation\sigma and normal mixture priors for themeans. As a consequence from the assumption of a known standarddeviation, the calculation discards sampling uncertainty of thesecond moment. The function has two extra arguments (withdefaults):eps (10^{-6}) andNgrid (10). Thedecision boundary is searched in the region of probability mass1-eps, respectively fory_1 andy_2. Thecontinuous decision function is evaluated at a discrete grid, whichis determined by a spacing with\delta_2 =\sigma_2/\sqrt{N_{grid}}. Once the decision boundary is evaluatedat the discrete steps, a spline is used to inter-polate thedecision boundary at intermediate points.

• oc2S(gammaMix): Applies for the Poisson model with a gammamixture prior for the rate parameter. The functionoc2S takes an extra argumenteps (defaults to10^{-6}) whichdetermines the region of probability mass1-eps where theboundary is searched fory_1 andy_2, respectively.

References

Schmidli H, Gsteiger S, Roychoudhury S, O'Hagan A, Spiegelhalter D, Neuenschwander B.Robust meta-analytic-predictive priors in clinical trials with historical control information.Biometrics 2014;70(4):1023-1032.

See Also

Other design2S:decision2S(),decision2S_boundary(),pos2S()

Examples

# example from Schmidli et al., 2014dec <- decision2S(0.975, 0, lower.tail = FALSE)prior_inf <- mixbeta(c(1, 4, 16))prior_rob <- robustify(prior_inf, weight = 0.2, mean = 0.5)prior_uni <- mixbeta(c(1, 1, 1))N <- 40N_ctl <- N - 20# compare designs with different priorsdesign_uni <- oc2S(prior_uni, prior_uni, N, N_ctl, dec)design_inf <- oc2S(prior_uni, prior_inf, N, N_ctl, dec)design_rob <- oc2S(prior_uni, prior_rob, N, N_ctl, dec)# type I errorcurve(design_inf(x, x), 0, 1)curve(design_uni(x, x), lty = 2, add = TRUE)curve(design_rob(x, x), lty = 3, add = TRUE)# powercurve(design_inf(0.2 + x, 0.2), 0, 0.5)curve(design_uni(0.2 + x, 0.2), lty = 2, add = TRUE)curve(design_rob(0.2 + x, 0.2), lty = 3, add = TRUE)

Diagnostic plots for EM fits

Description

Produce diagnostic plots of EM fits returned frommixfit().

Usage

## S3 method for class 'EM'plot(x, size = 1.25, link = c("identity", "logit", "log"), ...)

Arguments

Value

A list ofggplot2::ggplot() plots fordiagnostics of the EM run. Detailed EM diagnostic plots areincluded only if the global optionRBesT.verbose is set toTRUE. These include plots of the parameters of eachcomponent vs the iteration. The plot of the mixture density with ahistogram and a density of the fitted sample is always returned.

Customizingggplot2 plots

The returned plot is aggplot2 object. Please refer to the"Customizing Plots" vignette which is part ofRBesTdocumentation for an introduction. For simple modifications (changelabels, add reference lines, ...) consider the commands found inbayesplot-helpers. For more advancedcustomizations please use theggplot2 package directly. Adescription of the most common tasks can be found in theR Cookbook and a fullreference of available commands can be found at theggplot2 documentationsite.

See Also

Other EM:mixfit()

Examples

bmix <- mixbeta(rob = c(0.2, 1, 1), inf = c(0.8, 10, 2))bsamp <- rmix(bmix, 1000)bfit <- mixfit(bsamp, type = "beta", Nc = 2)pl <- plot(bfit)print(pl$mixdens)print(pl$mix)# a number of additional plots are generated in verbose mode.user_option <- options(RBesT.verbose = TRUE)pl_all <- plot(bfit)# recover previous user optionsoptions(.user_option)names(pl_all)# [1] "mixdist" "a"   "b"   "w"   "m"   "N"   "Lm"  "lN"  "Lw"  "lli" "mixdens" "mixecdf" "mix"

Diagnostic plots for gMAP analyses

Description

Diagnostic plots for gMAP analyses

Usage

## S3 method for class 'gMAP'plot(x, size = NULL, ...)

Arguments

Details

Creates MCMC diagnostics and a forest plot (includingmodel estimates) for agMAP() analysis. For acustomized forest plot, please use the dedicated functionforest_plot().

Value

The function returns a list ofggplot2::ggplot()objects.

Customizingggplot2 plots

The returned plot is aggplot2 object. Please refer to the"Customizing Plots" vignette which is part ofRBesTdocumentation for an introduction. For simple modifications (changelabels, add reference lines, ...) consider the commands found inbayesplot-helpers. For more advancedcustomizations please use theggplot2 package directly. Adescription of the most common tasks can be found in theR Cookbook and a fullreference of available commands can be found at theggplot2 documentationsite.

Probability of Success for a 1 Sample Design

Description

Thepos1S function defines a 1 sample design (prior, samplesize, decision function) for the calculation of the frequency atwhich the decision is evaluated to 1 when assuming a distributionfor the parameter. A function is returned which performs theactual operating characteristics calculations.

Usage

pos1S(prior, n, decision, ...)## S3 method for class 'betaMix'pos1S(prior, n, decision, ...)## S3 method for class 'normMix'pos1S(prior, n, decision, sigma, eps = 1e-06, ...)## S3 method for class 'gammaMix'pos1S(prior, n, decision, eps = 1e-06, ...)

Arguments

Details

Thepos1S function defines a 1 sample design andreturns a function which calculates its probability of success.The probability of success is the frequency with which the decisionfunction is evaluated to 1 under the assumption of a given truedistribution of the data implied by a distirbution of the parameter\theta.

Calling thepos1S function calculates the critical valuey_c and returns a function which can be used to evaluate thePoS for different predictive distributions and is evaluated as

\int F(y_c|\theta) p(\theta) d\theta,

whereF is the distribution function of the samplingdistribution andp(\theta) specifies the assumed truedistribution of the parameter\theta. The distributionp(\theta) is a mixture distribution and given as themix argument to the function.

Value

Returns a function that takes as single argumentmix, which is the mixture distribution of the controlparameter. Calling this function with a mixture distribution thencalculates the PoS.

Methods (by class)

• pos1S(betaMix): Applies for binomial model with a mixturebeta prior. The calculations use exact expressions.

• pos1S(normMix): Applies for the normal model with knownstandard deviation\sigma and a normal mixture prior for themean. As a consequence from the assumption of a known standarddeviation, the calculation discards sampling uncertainty of thesecond moment. The functionpos1S has an extraargumenteps (defaults to10^{-6}). The critical valuey_c is searched in the region of probability mass1-eps fory.

• pos1S(gammaMix): Applies for the Poisson model with a gammamixture prior for the rate parameter. The functionpos1S takes an extra argumenteps (defaults to10^{-6})which determines the region of probability mass1-eps wherethe boundary is searched fory.

See Also

Other design1S:decision1S(),decision1S_boundary(),oc1S()

Examples

# non-inferiority example using normal approximation of log-hazard# ratio, see ?decision1S for all detailss <- 2flat_prior <- mixnorm(c(1, 0, 100), sigma = s)nL <- 233theta_ni <- 0.4theta_a <- 0alpha <- 0.05beta <- 0.2za <- qnorm(1 - alpha)zb <- qnorm(1 - beta)n1 <- round((s * (za + zb) / (theta_ni - theta_a))^2)theta_c <- theta_ni - za * s / sqrt(n1)# assume we would like to conduct at an interim analysis# of PoS after having observed 20 events with a HR of 0.8.# We first need the posterior at the interim ...post_ia <- postmix(flat_prior, m = log(0.8), n = 20)# dual criteriondecComb <- decision1S(c(1 - alpha, 0.5), c(theta_ni, theta_c), lower.tail = TRUE)# ... and we would like to know the PoS for a successful# trial at the end when observing 10 more eventspos_ia <- pos1S(post_ia, 10, decComb)# our knowledge at the interim is just the posterior at# interim such that the PoS ispos_ia(post_ia)

Probability of Success for 2 Sample Design

Description

Thepos2S function defines a 2 sample design (priors, samplesizes & decision function) for the calculation of the probabilityof success. A function is returned which calculates the calculatesthe frequency at which the decision function is evaluated to 1 whenparameters are distributed according to the given distributions.

Usage

pos2S(prior1, prior2, n1, n2, decision, ...)## S3 method for class 'betaMix'pos2S(prior1, prior2, n1, n2, decision, eps, ...)## S3 method for class 'normMix'pos2S(  prior1,  prior2,  n1,  n2,  decision,  sigma1,  sigma2,  eps = 1e-06,  Ngrid = 10,  ...)## S3 method for class 'gammaMix'pos2S(prior1, prior2, n1, n2, decision, eps = 1e-06, ...)

Arguments

Details

Thepos2S function defines a 2 sample design andreturns a function which calculates its probability of success.The probability of success is the frequency with which the decisionfunction is evaluated to 1 under the assumption of a given truedistribution of the data implied by a distirbution of theparameters\theta_1 and\theta_2.

The calculation is analogous to the operating characeristicsoc2S() with the difference that instead of assumingknown (point-wise) true parameter values a distribution isspecified for each parameter.

Calling thepos2S function calculates the decision boundaryD_1(y_2) and returns a function which can be used to evaluate thePoS for different predictive distributions. It is evaluated as

\int\int\int f_2(y_2|\theta_2) \, p(\theta_2) \, F_1(D_1(y_2)|\theta_1) \, p(\theta_1) \, dy_2 d\theta_2 d\theta_1.

whereF is the distribution function of the samplingdistribution andp(\theta_1) andp(\theta_2) specifiesthe assumed true distribution of the parameters\theta_1 and\theta_2, respectively. Each distributionp(\theta_1)andp(\theta_2) is a mixture distribution and given as themix1 andmix2 argument to the function.

For example, in the binary case an integration of the predictivedistribution, the BetaBinomial, instead of the binomialdistribution will be performed over the data space wherever thedecision function is evaluated to 1. All other aspects of thecalculation are as for the 2-sample operating characteristics, seeoc2S().

Value

Returns a function which when called with two argumentsmix1 andmix2 will return the frequencies atwhich the decision function is evaluated to 1. Each argument isexpected to be a mixture distribution representing the assumed truedistribution of the parameter in each group.

Methods (by class)

• pos2S(betaMix): Applies for binomial model with a mixturebeta prior. The calculations use exact expressions. If theoptional argumenteps is defined, then an approximate methodis used which limits the search for the decision boundary to theregion of1-eps probability mass. This is useful for designswith large sample sizes where an exact approach is very costly tocalculate.

• pos2S(normMix): Applies for the normal model with knownstandard deviation\sigma and normal mixture priors for themeans. As a consequence from the assumption of a known standarddeviation, the calculation discards sampling uncertainty of thesecond moment. The function has two extra arguments (withdefaults):eps (10^{-6}) andNgrid (10). Thedecision boundary is searched in the region of probability mass1-eps, respectively fory_1 andy_2. Thecontinuous decision function is evaluated at a discrete grid, whichis determined by a spacing with\delta_2 =\sigma_2/\sqrt{N_{grid}}. Once the decision boundary is evaluatedat the discrete steps, a spline is used to inter-polate thedecision boundary at intermediate points.

• pos2S(gammaMix): Applies for the Poisson model with a gammamixture prior for the rate parameter. The functionpos2S takes an extra argumenteps (defaults to10^{-6}) whichdetermines the region of probability mass1-eps where theboundary is searched fory_1 andy_2, respectively.

See Also

Other design2S:decision2S(),decision2S_boundary(),oc2S()

Examples

# see ?decision2S for details of examplepriorT <- mixnorm(c(1, 0, 0.001), sigma = 88, param = "mn")priorP <- mixnorm(c(1, -49, 20), sigma = 88, param = "mn")# the success criteria is for delta which are larger than some# threshold value which is why we set lower.tail=FALSEsuccessCrit <- decision2S(c(0.95, 0.5), c(0, 50), FALSE)# example interim outcomepostP_interim <- postmix(priorP, n = 10, m = -50)postT_interim <- postmix(priorT, n = 20, m = -80)# assume that mean -50 / -80 were observed at the interim for# placebo control(n=10) / active treatment(n=20) which gives# the posteriorspostP_interimpostT_interim# then the PoS to succeed after another 20/30 patients ispos_final <- pos2S(postP_interim, postT_interim, 20, 30, successCrit)pos_final(postP_interim, postT_interim)

Conjugate Posterior Analysis

Description

Calculates the posterior distribution for datadata given a priorpriormix, where the prior is a mixture of conjugate distributions.The posterior is then also a mixture of conjugate distributions.

Usage

postmix(priormix, data, ...)## S3 method for class 'betaMix'postmix(priormix, data, n, r, ...)## S3 method for class 'normMix'postmix(priormix, data, n, m, se, ...)## S3 method for class 'gammaMix'postmix(priormix, data, n, m, ...)

Arguments

Details

A conjugate prior-likelihood pair has the convenientproperty that the posterior is in the same distributional class asthe prior. This property also applies to mixtures of conjugatepriors. Let

p(\theta;\mathbf{w},\mathbf{a},\mathbf{b})

denote a conjugate mixture prior density for data

y|\theta \sim f(y|\theta),

wheref(y|\theta) is the likelihood. Then the posterior isagain a mixture with each componentk equal to the respectiveposterior of thekth prior component and updated weightsw'_k,

p(\theta;\mathbf{w'},\mathbf{a'},\mathbf{b'}|y) = \sum_{k=1}^K w'_k \, p_k(\theta;a'_k,b'_k|y).

The weightw'_k forkth component is determined by themarginal likelihood of the new datay under thekth priordistribution which is given by the predictive distribution of thekth component,

w'_k \propto w_k \, \int p_k(\theta;a_k,b_k) \, f(y|\theta) \, d\theta \equiv w^\ast_k .

The final weightw'_k is then given by appropriatenormalization,w'_k = w^\ast_k / \sum_{k=1}^K w^\ast_k. In other words, the weight ofcomponentk is proportional to the likelihood that datay is generated from the respective component, i.e. themarginal probability; for details, see for exampleSchmidliet al., 2015.

Note: The prior weightsw_k are fixed, but theposterior weightsw'_k \neq w_k still change due to thechanging normalization.

The datay can either be given as individual data or assummary data (sufficient statistics). See below for details for theimplemented conjugate mixture prior densities.

Methods (by class)

• postmix(betaMix): Calculates the posterior beta mixturedistribution. The individual data vector is expected to be a vectorof 0 and 1, i.e. a series of Bernoulli experiments. Alternatively,the sufficient statisticsn andr can be given,i.e. number of trials and successes, respectively.

• postmix(normMix): Calculates the posterior normal mixturedistribution with the sampling likelihood being a normal with fixedstandard deviation. Either an individual data vectordatacan be given or the sufficient statistics which are the standarderrorse and sample meanm. If the sample sizen is used instead of the sample standard error, then thereference scale of the prior is used to calculate the standarderror. Should standard errorse and sample sizen begiven, then the reference scale of the prior is updated; however itis recommended to use the commandsigma() to set thereference standard deviation.

• postmix(gammaMix): Calculates the posterior gamma mixturedistribution for Poisson and exponential likelihoods. Only thePoisson case is supported in this version.

Supported Conjugate Prior-Likelihood Pairs

References

Schmidli H, Gsteiger S, Roychoudhury S, O'Hagan A, Spiegelhalter D, Neuenschwander B.Robust meta-analytic-predictive priors in clinical trials with historical control information.Biometrics 2014;70(4):1023-1032.

Examples

# binary example with individual data (1=event,0=no event), uniform priorprior.unif <- mixbeta(c(1, 1, 1))data.indiv <- c(1, 0, 1, 1, 0, 1)posterior.indiv <- postmix(prior.unif, data.indiv)print(posterior.indiv)# or with summary data (number of events and number of patients)r <- sum(data.indiv)n <- length(data.indiv)posterior.sum <- postmix(prior.unif, n = n, r = r)print(posterior.sum)# binary example with robust informative prior and conflicting dataprior.rob <- mixbeta(c(0.5, 4, 10), c(0.5, 1, 1))posterior.rob <- postmix(prior.rob, n = 20, r = 18)print(posterior.rob)# normal example with individual datasigma <- 88prior.mean <- -49prior.se <- sigma / sqrt(20)prior <- mixnorm(c(1, prior.mean, prior.se), sigma = sigma)data.indiv <- c(-46, -227, 41, -65, -103, -22, 7, -169, -69, 90)posterior.indiv <- postmix(prior, data.indiv)# or with summary data (mean and number of patients)mn <- mean(data.indiv)n <- length(data.indiv)posterior.sum <- postmix(prior, m = mn, n = n)print(posterior.sum)

Predictive Distributions for Mixture Distributions

Description

Predictive distribution for mixture of conjugate distributions(beta, normal, gamma).

Usage

preddist(mix, ...)## S3 method for class 'betaMix'preddist(mix, n = 1, ...)## S3 method for class 'normMix'preddist(mix, n = 1, sigma, ...)## S3 method for class 'gammaMix'preddist(mix, n = 1, ...)## S3 method for class 'mvnormMix'preddist(mix, ...)

Arguments

Details

Given a mixture density (either a posterior or a prior)

p(\theta,\mathbf{w},\mathbf{a},\mathbf{b})

and a data likelihood of

y|\theta \sim f(y|\theta),

the predictive distribution of a one-dimensional summaryy_nof $n$ future observations is distributed as

y_n \sim \int p(\theta,\mathbf{w},\mathbf{a},\mathbf{b}) \, f(y_n|\theta) \, d\theta .

This distribution is the marginal distribution of the data underthe mixture density. For binary and Poisson datay_n =\sum_{i=1}^n y_i is the sum over future events. For normal data,it is the mean\bar{y}_n = 1/n \sum_{i=1}^n y_i.

Value

The function returns for a normal, beta or gamma mixturethe matching predictive distribution fory_n.

Methods (by class)

• preddist(betaMix): Obtain the matching predictive distributionfor a beta distribution, the BetaBinomial.

• preddist(normMix): Obtain the matching predictive distributionfor a Normal with constant standard deviation. Note that thereference scale of the returned Normal mixture is meaninglessas the individual components are updated appropriatley.

• preddist(gammaMix): Obtain the matching predictive distributionfor a Gamma. Only Poisson likelihoods are supported.

• preddist(mvnormMix): Multivariate normal mixtures predictivedensities are not (yet) supported.

Supported Conjugate Prior-Likelihood Pairs

Examples

# Example 1: predictive distribution from uniform prior.bm <- mixbeta(c(1, 1, 1))bmPred <- preddist(bm, n = 10)# predictive proabilities and cumulative predictive probabilitiesx <- 0:10d <- dmix(bmPred, x)names(d) <- xbarplot(d)cd <- pmix(bmPred, x)names(cd) <- xbarplot(cd)# medianmdn <- qmix(bmPred, 0.5)mdn# Example 2: 2-comp Beta mixturebm <- mixbeta(inf = c(0.8, 15, 50), rob = c(0.2, 1, 1))plot(bm)bmPred <- preddist(bm, n = 10)plot(bmPred)mdn <- qmix(bmPred, 0.5)mdnd <- dmix(bmPred, x = 0:10)n.sim <- 100000r <- rmix(bmPred, n.sim)dtable(r) / n.sim# Example 3: 3-comp Normal mixturem3 <- mixnorm(c(0.50, -0.2, 0.1), c(0.25, 0, 0.2), c(0.25, 0, 0.5), sigma = 10)print(m3)summary(m3)plot(m3)predm3 <- preddist(m3, n = 2)plot(predm3)print(predm3)summary(predm3)

Predictions from gMAP analyses

Description

Produces a sample of the predictive distribution.

Usage

## S3 method for class 'gMAP'predict(  object,  newdata,  type = c("response", "link"),  probs = c(0.025, 0.5, 0.975),  na.action = na.pass,  thin,  ...)## S3 method for class 'gMAPpred'print(x, digits = 3, ...)## S3 method for class 'gMAPpred'summary(object, ...)## S3 method for class 'gMAPpred'as.matrix(x, ...)

Arguments

Details

Predictions are made using the\tau predictionstratum of the gMAP object. For details on the syntax, please refertopredict.glm() and the example below.

See Also

gMAP(),predict.glm()

Examples

## Setting up dummy sampling for fast execution of example## Please use 4 chains and 20x more warmup & iter in practice.user_mc_options <- options(RBesT.MC.warmup=50, RBesT.MC.iter=100,                            RBesT.MC.chains=2, RBesT.MC.thin=1)# create a fake data set with a covariatetrans_cov <- transform(transplant, country = cut(1:11, c(0, 5, 8, Inf), c("CH", "US", "DE")))set.seed(34246)map <- gMAP(cbind(r, n - r) ~ 1 + country | study,  data = trans_cov,  tau.dist = "HalfNormal",  tau.prior = 1,  # Note on priors: we make the overall intercept weakly-informative  # and the regression coefficients must have tighter sd as these are  # deviations in the default contrast parametrization  beta.prior = rbind(c(0, 2), c(0, 1), c(0, 1)),  family = binomial,  ## ensure fast example runtime  thin = 1, chains = 1)# posterior predictive distribution for each input data item (shrinkage estimates)pred_cov <- predict(map)pred_cov# extract sample as matrixsamp <- as.matrix(pred_cov)# predictive distribution for each input data item (if the input studies were new ones)pred_cov_pred <- predict(map, trans_cov)pred_cov_pred# a summary function returns the results as matrixsummary(pred_cov)# obtain a prediction for new data with specific covariatespred_new <- predict(map, data.frame(country = "CH", study = 12))pred_new## Recover user set sampling defaultsoptions(.user_mc_options)

Robustify Mixture Priors

Description

Add a non-informative component to a mixture prior.

Usage

robustify(priormix, weight, mean, n = 1, ...)## S3 method for class 'betaMix'robustify(priormix, weight, mean, n = 1, ...)## S3 method for class 'gammaMix'robustify(priormix, weight, mean, n = 1, ...)## S3 method for class 'normMix'robustify(priormix, weight, mean, n = 1, ..., sigma)

Arguments

Details

It is recommended to robustify informative priors derivedwithgMAP() using unit-information priors . Thisprotects against prior-data conflict, see for exampleSchmidli et al., 2015.

The procedure can be used with beta, gamma and normal mixturepriors. A unit-information prior (seeKass and Wasserman,1995) corresponds to a prior which represents the observation ofn=1 at the null hypothesis. As the null is problem dependent westrongly recommend to make use of themean argumentaccordingly. See below for the definition of the default mean.

The weights of the mixture priors are rescaled to(1-weight)while the non-informative prior is assigned theweightgiven.

Value

New mixture with an extra non-informative component namedrobust.

Methods (by class)

• robustify(betaMix): The defaultmean is set to 1/2 whichrepresents no difference between the occurrence rates for one of thetwo outcomes. As the uniformBeta(1,1) is more appropriate inpractical applications,RBesT usesn+1 as the samplesize such that the default robust prior is the uniform instead oftheBeta(1/2,1/2) which strictly defined would be the unitinformation prior in this case.

• robustify(gammaMix): The defaultmean is set to the mean of theprior mixture. It is strongly recommended to explicitly set themean to the location of the null hypothesis.

• robustify(normMix): The defaultmean is set to the meanof the prior mixture. It is strongly recommended to explicitly setthe mean to the location of the null hypothesis, which is veryoften equal to 0. It is also recommended to explicitly set thesampling standard deviation using thesigma argument.

References

Schmidli H, Gsteiger S, Roychoudhury S, O'Hagan A,Spiegelhalter D, Neuenschwander B. Robust meta-analytic-predictivepriors in clinical trials with historical control information.Biometrics 2014;70(4):1023-1032.

Kass RE, Wasserman L A Reference Bayesian Test for NestedHypotheses and its Relationship to the Schwarz CriterionJAmer Statist Assoc 1995; 90(431):928-934.

See Also

mixcombine()

Examples

bmix <- mixbeta(inf1 = c(0.2, 8, 3), inf2 = c(0.8, 10, 2))plot(bmix)rbmix <- robustify(bmix, weight = 0.1, mean = 0.5)rbmixplot(rbmix)gmnMix <- mixgamma(inf1 = c(0.2, 2, 3), inf2 = c(0.8, 2, 5), param = "mn")plot(gmnMix)rgmnMix <- robustify(gmnMix, weight = 0.1, mean = 2)rgmnMixplot(rgmnMix)nm <- mixnorm(inf1 = c(0.2, 0.5, 0.7), inf2 = c(0.8, 2, 1), sigma = 2)plot(nm)rnMix <- robustify(nm, weight = 0.1, mean = 0, sigma = 2)rnMixplot(rnMix)

Support of Distributions

Description

Returns the support of a distribution.

Usage

support(mix)

Arguments

Transplant.

Description

Data set containing historical information for standard treatmentfor a phase IV trial in de novo transplant patients. The primaryoutcome is treament failure (binary).

Usage

transplant

Format

A data frame with 4 rows and 3 variables:

study

study

n

study size

r

number of events

References

Neuenschwander B, Capkun-Niggli G, Branson M,Spiegelhalter DJ. Summarizing historical information on controls inclinical trials.Clin Trials. 2010; 7(1):5-18

Find root of univariate function of integers

Description

Uses a bisectioning algorithm to search the give interval for achange of sign and returns the integer which is closest to 0.

Usage

uniroot_int(  f,  interval,  ...,  f.lower = f(interval[1], ...),  f.upper = f(interval[2], ...),  maxIter = 1000)