text:
#### Vol. 2, No. 4, 2009

 Recent Issues
 The Journal About the Journal Editorial Board Editors’ Interests Subscriptions Submission Guidelines Submission Form Policies for Authors Ethics Statement ISSN: 1944-4184 (e-only) ISSN: 1944-4176 (print) Author Index Coming Soon Other MSP Journals
Automatic growth series for right-angled Coxeter groups

### Rebecca Glover and Richard Scott

Vol. 2 (2009), No. 4, 371–385
##### Abstract

Right-angled Coxeter groups have a natural automatic structure induced by their action on a CAT($0$) cube complex. The normal form for this structure is defined with respect to the generating set consisting of all cliques in the defining graph for the group. In this paper we study the growth series for right-angled Coxeter groups with respect to this automatic generating set. In particular, we show that there exist nonisomorphic Coxeter groups with the same automatic growth series, and give a comparison with the usual growth series defined with respect to the standard generating set.

##### Keywords
Coxeter groups, growth series
##### Mathematical Subject Classification 2000
Primary: 05A15, 20F10, 20F55

format: plain


text:
378 articles – 175 references  [version française]
 HAL: hal-00080397, version 2
 Available versions: v1 (2006-06-28) v2 (2007-02-16)
 Homogenization of first order equations with $u/\epsilon$-periodic Hamiltonians. Part II: application to dislocations dynamics
 Cyril Imbert 1, Régis Monneau 2
 For the ACI JC 1025 collaboration(s)
 (2006-06-16)
 This paper is concerned with a result of homogenization of a non-local first order Hamilton-Jacobi equations describing the dislocations dynamics. Our model for the interaction between dislocations involve both an integro-differential operator and a (local) Hamiltonian depending periodicly on $u/\eps$. The first two authors studied in a previous work homogenization problems involving such local Hamiltonians. Two main ideas of this previous work are used: on the one hand, we prove an ergodicity property of this equation by constructing approximate correctors which are necessarily non periodic in space in general; on the other hand, the proof of the convergence of the solution uses here a twisted perturbed test function for a higher dimensional problem. The limit equation is a nonlinear diffusion equation involving a first order Lévy operator; the nonlinearity keeps memory of the short range interaction, while the Lévy operator keeps memory of long ones. The homogenized equation is a kind of effective plastic law for densities of dislocations moving in a single slip plane.
 1: Institut de Mathématiques et de Modélisation de Montpellier (I3M) CNRS : UMR5149 – Université Montpellier II - Sciences et techniques 2: Centre d'Enseignement et de Recherche en Mathématiques, Informatique et Calcul Scientifique (CERMICS) INRIA – Ecole des Ponts ParisTech 3: École Centrale de Lyon (ECL) Ministère de l'Enseignement Supérieur et de la Recherche Scientifique
 Subject : Mathematics/Analysis of PDEs
 Keyword(s): periodic homogenization – Hamilton-Jacobi equations – integro-differential operators – dislocations dynamics – non-periodic approximate correctors
Attached file list to this document:
 PDF
 homog-dislo.pdf(336.4 KB)
 hal-00080397, version 2 http://hal.archives-ouvertes.fr/hal-00080397 oai:hal.archives-ouvertes.fr:hal-00080397 From: Cyril Imbert <> Submitted on: Friday, 16 February 2007 14:13:11 Updated on: Friday, 16 February 2007 15:04:06


text:
# Generalization of a theorem of Mahler to higher dimensions

A seminal theorem of Kurt Mahler, in his papers Zur Approximation algebraischer Zahlen. I-III., is the following:

Let $F(x,y) \in \mathbb{Z}[x,y]$ be a binary form of degree $d \geq 3$ and non-zero discriminant $\Delta(F)$. Put $A_F$ for the area of the region

$\displaystyle \{(x,y) \in \mathbb{R}^2 : |F(x,y)| \leq 1\}$

and put $N_F(Z)$ for the number of points $(x,y) \in \mathbb{Z}^2$ such that $|F(x,y)| \leq Z$ and $F(x,y) \ne 0$. Then he proved that

$\displaystyle N_F(Z) = A_F Z^{\frac{2}{d}} + O_F \left(Z^{\frac{1}{d-1}}\right).$

Has this statement been generalized to higher dimensions? The most conservative generalization would be to decomposable forms, meaning homogeneous polynomials in $x_1, \cdots, x_n$ with integer coefficients which factor completely into linear factors over $\mathbb{C}$, say. In particular, the set-up is as follows: Let

$F(x_1, \cdots, x_n) = L_1(\mathbf{x}) \cdots L_d(\mathbf{x})$

be a decomposable form with integer coefficients, where $L_j(\mathbf{x})$ are linear forms with coefficients in $\overline{\mathbb{Q}}$, and $d > n$. Put $V_F$ for the Lebesgue measure of the set

$\displaystyle \{(x_1, \cdots, x_n) \in \mathbb{R}^n : |F(x_1, \cdots, x_n)| \leq 1\}$

and $N_F(Z)$ for the number of points $(x_1, \cdots, x_n) \in \mathbb{Z}^n$ such that $|F(x_1, \cdots, x_n)| \leq Z, F(x_1, \cdots, x_n) \ne 0$. Then has it been proven that

$\displaystyle N_F(Z) = V_F Z^{\frac{n}{d}} + o_F\left(Z^{\frac{n}{d}}\right)$

in general?

I did a search of papers which cites Mahler's Zur Approximation algebraischer Zahlen. I and did not find any of them that deals with this problem (at least from looking at their titles and abstracts).

A more optimistic and aggressive generalization would be to consider general homogeneous polynomials in $n$ variables. I am not sure what additional complications could arise.

Any help would be appreciated.

• Where Mahler's proof seems to fail in proving your generalization? Also, a small remark: You want to assume $d>n$ (to avoid, for instance, the case where $F$ is the norm form of some number field, which implies $N_F$ is infinite - think of units). Dec 30, 2015 at 7:49
• @OfirGorodetsky Mahler's papers are in German, a language which I am not fluent in, so it has been very hard to read his papers and I am sure my understanding is not quite right in any case. Further, you are right, I meant to phrase my question around incomplete decomposable forms Dec 30, 2015 at 14:40

However, as far as I know, the situation for general homogeneous $F$ in many variables has not been addressed.


text:
## Evaluating with Probabilistic Truth: Log Loss vs. O/1 Loss

I’ve been advocating an approach to classifier evaluation in which the “truth” is itself a probability distribution over possible labels and system responses take the same form.

### Is the Truth Out There?

The main motivation is that we don’t actually know what the truth is for some items being classified. Not knowing the truth can come about because our labelers are biased or inconsistent, because our coding standard is underspecified, or because the items being classified are truly marginal. That is, the truth may not actually be out there, at least in any practically useful sense.

For example, we can squint all we want at the use of a pronoun in text and still not be able to come to any consensus about its antecedent, no matter how many experts in anaphora we assemble. Similarly, we might look at a phrase like fragile X in a paper on autism and truly not know if it refers to the gene or the disease. We may not know a searcher’s intent when we try to do spelling correction on a query like Montona — maybe they want Montana the state, maybe they want the pop star, or maybe they really want the west central region of the island of São Vicente.

The truth’s not always out there. At least for language data collected from the wild.

### Partial Credit for Uncertainty

The secondary motivation is that many classifiers, like most of LingPipe’s, are able to return probabilistic answers. The basic motivation of log loss is that we want to evaluate the probability assigned by the system to the correct answers. The less uncertain the system is about the true answer, the better.

If we want to be able to trade precision (or specificity) for recall (i.e., sensitivity), it’s nice to have probabilistically normalized outputs over which to threshold.

### Probabilistic Reference Truth

Let’s just consider binary (two category) classification problems, where a positive outcome is coded as $1$ and a negative one as $0$.. Suppose there are $N$ items being classified, $x_1, \ldots, x_N$. A probabilistic corpus then consists of a sequence of probabilistic categorizations $y_1, \ldots, y_N$, where $y_n \in [0,1]$ for $1 \leq n \leq N$.

A non-probabilistic corpus is just an edge case, with $y_n \in \{ 0, 1 \}$ for $1 \leq n \leq N$.

### Probabilistic System Responses

Next, suppose that the systems we’re evaluating are themselves allowed to return conditional probabilities of positive responses, which we will write as $\hat{y}_n \in [0,1]$ for $1 \leq n \leq N$.

A non-probabilistic response is again an edge case, with $\hat{y}_n \in \{ 0, 1 \}$.

### Log Loss

The log loss for response $\hat{y}$ for reference $y$ is the negative expected log probability of the response given the reference, namely

${\cal L}(y,\hat{y}) = - \sum_{n=1}^N y_n \times \log \hat{y}_n + (1-y_n) \times \log (1 - \hat{y}_n)$.

If we think of $y$ and $\hat{y}$ as parameters to Bernoulli distributions, it’s just KL divergence.

If the reference $y$ is restricted to being $0$ or $1$, this reduces to the usual definition of log loss for a probabilistic classifier (that is, $\hat{y}$ may be any value in $[0,1]$.

### Plotting a Single Response

Because we just sum over all the responses, we can see what the result is for a single item $n$. Suppose $y_n = 0.7$, which means the reference says there’s a 70% chance of item $n$ being positive, and a 30% chance of it being negative. The plot of the loss for response $\hat{y}_n$ given reference $y_n = 0.7$ is

It’s easy to see that the response $\hat{y} = 1$ is better than $\hat{y} = 0$ for the case where the reference is $y = 0.7$.

The red vertical dashed line is at the point $\hat{y} = 0.7$, which is where $y = \hat{y}$. The blue horizontal line is at ${\cal L}(y,\hat{y}) = 0.61$, which is the minimum loss for any $\hat{y}$ when $y = 0.7$.

### Classifier Uncertainty Matches Reference Uncertainty

It’s easy to see from the graph (and prove via differentiation) that the log loss function is convex and the minimum loss occurs when $y = \hat{y} = 0.7$. That is, best performance occurs when the system response has the same uncertainty as the gold standard. In other words, we want to train the system to be uncertain in the same places as the reference (training) data is uncertain.

### Log Loss Degenerates at the Edges

If the reference is $y = 0$ and the response $\hat{y} = 1$, the log loss is infinite. Similarly for $y = 1, \hat{y} = 0$. As long as the response is not 0 or 1, the results will be finite.

If the response is 0 or 1, the loss will only be finite if the reference category is the same value, $y = \hat{y}$. In fact, having the reference equal to the response, $y = \hat{y}$, with both being 0 or 1, is the only way to get zero loss. In all other situations, loss will be non-zero. Loss is never negative because we’re negating logs of values between 0 and 1.

This is why you don’t see log loss used for bakeoffs. If you’re using a system that doesn’t have probabilistic output and the truth is binary 0/1, then every system will have infinite error.

### Relation to Regularization via Prior Data

This matching of uncertainty can be thought of as a kind of regularization. One common way to implement priors if you only have maximum likelihood estimators is to add some artificial data. If I add a new data point $x$ which has all features turned on, and give it one positive instance and one negative instance, and I’m doing logistic regression, it acts as a regularizer pushing all the coefficients closer to zero (the further they are away from zero, the worse they do on the combination of the one negative and one positive example). If I add more of this kind of data, the regularization effect is stronger. Here, we can think of adding 3 negative examples and 7 positive examples, or 0.3 and 0.7 if we allow fractional training data.

### Comparison to 0/1 Loss

Consider what happens if instead of log loss, we consider the standard 0/1 loss as expressed in a confusion matrix (true positive, true negative, false positive, and false negative counts).

With a probabilistic reference category $y$ and probabilistic system response $\hat{y}$, the expected confusion matrix counts are

$\mathbb{E}[\mbox{TP}] = y \times \hat{y}$,

$\mathbb{E}[\mbox{TN}] = (1-y) \times (1-\hat{y})$,

$\mathbb{E}[\mbox{FP}] = (1-y) \times \hat{y}$, and

$\mathbb{E}[\mbox{FN}] = y \times (1 - \hat{y})$.

With this, we can compute expected accuracy, precision, recall (= sensitivity), specificity, and F-measure as

$\mathbb{E}[\mbox{acc}] = \mathbb{E}[(\mbox{TP}+\mbox{TN})/(\mbox{TP}+\mbox{TN}+\mbox{FP}+\mbox{FN})]$
${} = y \hat{y} + (1 - y)(1-\hat{y})$.

$\mathbb{E}[\mbox{recall}] = \mathbb{E}[\mbox{TP}/(\mbox{TP}+\mbox{FN})]$
${} = y \hat{y}/(y\hat{y} + y(1 - \hat{y})) = \hat{y}$,

$\mathbb{E}[\mbox{precision}] = \mathbb{E}[\mbox{TP}/(\mbox{TP}+\mbox{FP})]$
${} = y \hat{y}/(y\hat{y} + (1-y)(\hat{y})) = y$,

$\mathbb{E}[\mbox{specificity}] = \mathbb{E}[\mbox{TN}/(\mbox{TN}+\mbox{FP})]$
${} = (1-y)(1-\hat{y})/((1-y)(1-\hat{y}) + (1-y)\hat{y}) = 1-\hat{y}$, and

$\mathbb{E}[\mbox{F-measure}] = \mathbb{E}[2 \times \mbox{prec} \times \mbox{recall} / (\mbox{prec} + \mbox{recall})]$
${} = 2 y \hat{y} / (y + \hat{y})$.

All of these statistical measures are degenerate for probabilistic evaluation. For instance, if $y > 0.5$, then accuracy is maximized with $\hat{y} = 1$, if $y < 0.5$, accuracy is maximized with $\hat{y} = 0$, and if $y = 0$, accuracy is the same for all $\hat{y} \in [0,1]$! The rest of the measures are no better behaved.

### R Source

Plotting functions is easy in R with the curve() function. Here’s how I generated the one above.

p <- 0.7;
f <- function(q) {
- (p * log(q) + (1-p) * log(1-q));
}
curve(f, xlim=c(0,1), ylim=c(0,4),
xlab="yHat", ylab="LogLoss(y, yHat)",
main="Log Loss  ( y=0.7 )",
lwd=2, cex.lab=1.5, cex.axis=1.5, cex.main=1.5,
axes=FALSE, frame=FALSE)

axis(1,at=c(0,0.5,0.7,1))
axis(2,at=c(0,f(0.7),2,4),labels=c("0", ".61", "2", "4"))
abline(h=f(0.7),col="blue",lty=3)
abline(v=0.7,col="red",lty=2)

Just be careful with that dynamic lexical scope if you want to redefine the reference result p ($y$ in the definitions above).

### 4 Responses to “Evaluating with Probabilistic Truth: Log Loss vs. O/1 Loss”

1. Vikas Raykar Says:

You may find this paper relevant

Beta Regression for Modelling Rates and Proportions
Silvia Ferraria and Francisco Cribari-Netob
http://www.ime.usp.br/~sferrari/beta.pdf

They even have a R package called betareg.

2. lingpipe Says:

Ferraria and Cribari-Netob use a generalized linear model with a Beta-type link. Translating what they do into notation I can more easily follow, they’re assuming responses $y_n \in [0,1]$ with predictors $X_n \in \mathbb{R}^K$. They then model the responses as beta distributed,

$y_n \sim \mbox{\sf Beta}(\mu_n \tau, (1-\mu_n) \tau)$,

where the mean parameter is modeled using logistic regression

$\mu_n = \mbox{logit}^{-1}(\beta X_n^T)$.

This’d be straightforward to program in BUGS or JAGS — you’d just need a prior on the regression coefficients $\beta_k$ and on the precision parameter $\tau$.

Now that I’ve translated it into terms I understand, I’ll have to get my head around the implications. My first reaction is that the deterministic part of the prediction will be the same as if you just took $y_n = \mu_n$, because the mean of $\mbox{\sf Beta}(\mu_n \tau, (1-\mu_n) \tau)$ is $\mu_n$.

3. lingpipe Says:

Phil Resnik pointed me to this paper, section 5 of which discusses training logistic regression using probabilistic “truth”:

Limin Yao, David Mimno and Andrew McCallum. 2009. Efficient Methods for Topic Model Inference on Streaming Document Collections. KDD.

This made me realize that, of course, EM uses exactly the same kind of “probabilistic training” in its E step. It’s not much more complicated to plug in logistic regression in place of the more standard naive Bayes.

4. Competitive Predictive Modeling – How Useful is it? at So much to do, so little time Says:

[…] The contest was open to anybody and who ever got the best classification model (as measured by log loss) was selected as the winner. You can read more about the details of the competition and also on […]


text:
# 1 Scope of this document

This document describes the statistical methods implemented in the rbmi R package for standard and reference-based multiple imputation of continuous longitudinal outcomes. The package implements three classes of multiple imputation (MI) approaches:

1. Conventional MI methods based on Bayesian (or approximate Bayesian) posterior draws of model parameters combined with Rubin’s rules to make inferences as described in Carpenter, Roger, and Kenward (2013) and Cro et al. (2020).

2. Conditional mean imputation methods combined with re-sampling techniques as described in Wolbers et al. (2022).

3. Bootstrapped MI methods as described in von Hippel and Bartlett (2021).

The document is structured as follows: we first provide an informal introduction to estimands and corresponding treatment effect estimation based on MI (section 2). The core of this document consists of section 3 which describes the statistical methodology in detail and also contains a comparison of the implemented approaches (section 3.10). The link between theory and the functions included in package rbmi is described in section 4. We conclude with a comparison of our package to some alternative software implementations of reference-based imputation methods (section 5).

# 2 Introduction to estimands and estimation methods

## 2.1 Estimands

The ICH E9(R1) addendum on estimands and sensitivity analyses describes a systematic approach to ensure alignment among clinical trial objectives, trial execution/conduct, statistical analyses, and interpretation of results (ICH E9 working group (2019)). As per the addendum, an estimand is a precise description of the treatment effect reflecting the clinical question posed by the trial objective which summarizes at a population-level what the outcomes would be in the same patients under different treatment conditions being compared. One important attribute of an estimand is a list of possible intercurrent events (ICEs), i.e. of events occurring after treatment initiation that affect either the interpretation or the existence of the measurements associated with the clinical question of interest, and the definition of appropriate strategies to deal with ICEs. The three most relevant strategies for the purpose of this document are the hypothetical strategy, the treatment policy strategy, and the composite strategy. For the hypothetical strategy, a scenario is envisaged in which the ICE would not occur. Under this scenario, endpoint values after the ICE are not directly observable and treated using models for missing data. For the treatment policy strategy, the treatment effect in the presence of the ICEs is targeted and analyses are based on the observed outcomes regardless whether the subject had an ICE or not. For the composite strategy, the ICE itself is included as a component of the endpoint.

## 2.2 Alignment between the estimand and the estimation method

The ICH E9(R1) addendum distinguishes between ICEs and missing data (ICH E9 working group (2019)). Whereas ICEs such as treatment discontinuations reflect clinical practice, the amount of missing data can be minimized in the conduct of a clinical trial. However, there are many connections between missing data and ICEs. For example, it is often difficult to retain subjects in a clinical trial after treatment discontinuation and a subject’s dropout from the trial leads to missing data. As another example, outcome values after ICEs addressed using a hypothetical strateg are not directly observable under the hypothetical scenario. Consequently, any observed outcome values after such ICEs are typically discarded and treated as missing data.

The addendum proposes that estimation methods to address the problem presented by missing data should be selected to align with the estimand. A recent overview of methods to align the estimator with the estimand is Mallinckrodt et al. (2020). A short introduction on estimation methods for studies with longitudinal endpoints can also be found in Wolbers et al. (2022). One prominent statistical method for this purpose is multiple imputation (MI), which is the target of the rbmi package.

### 2.2.1 Missing data prior to ICEs

Missing data may occur in subjects without an ICE or prior to the occurrence of an ICE. As such missing outcomes are not associated with an ICE, it is often plausible to impute them under a missing-at-random (MAR) assumption using a standard MMRM imputation model of the longitudinal outcomes. Informally, MAR occurs if the missing data can be fully accounted for by the baseline variables included in the model and the observed longitudinal outcomes, and if the model is correctly specified.

### 2.2.2 Implementation of the hypothetical strategy

The MAR imputation model described above is often also a good starting point for imputing data after an ICE handled using a hypothetical strategy (Mallinckrodt et al. (2020)). Informally, this assumes that unobserved values after the ICE would have been similar to the observed data from subjects who did not have the ICE and remained under follow-up. However, in some situations, it may be more reasonable to assume that missingness is “informative” and indicates a systematically better or worse outcome than in observed subjects. In such situations, MNAR imputation with a $$\delta$$-adjustment could be explored as a sensitivity analysis. $$\delta$$-adjustments add a fixed or random quantity to the imputations in order to make the imputed outcomes systematically worse or better than those observed as described in Cro et al. (2020). In rbmi only fixed $$\delta$$-adjustments are implemented.

### 2.2.3 Implementation of the treatment policy strategy

Ideally, data collection continues after an ICE handled with a treatment policy strategy and no missing data arises. Indeed, such post-ICE data are increasingly systematically collected in RCTs. However, despite best efforts, missing data after an ICE such as study treatment discontinuation may still occur because the subject drops out from the study after discontinuation. It is difficult to give definite recommendations regarding the implementation of the treatment policy strategy in the presence of missing data at this stage because the optimal method is highly context dependent and a topic of ongoing statistical research.

For ICEs which are thought to have a negligible effect on efficacy outcomes, standard MAR-based imputation may be appropriate. In contrast, an ICE such as treatment discontinuation may be expected to have a more substantial impact on efficacy outcomes. In such settings, the MAR assumption may still be plausible after conditioning on the subject’s time-varying treatment status (Guizzaro et al. (2021)). In this case, one option is to impute missing post-discontinuation data based on subjects who also discontinued treatment but continued to be followed up (Polverejan and Dragalin (2020)). Another option which may require somewhat less post-discontinuation data is to include all subjects in the imputation procedure but to model post-discontinuation data by using a time-varying treatment status indicators (e.g. time-varying indicators of treatment compliance, discontinuation, or initiation of rescue treatment) (Guizzaro et al. (2021)). In this approach, post-ICE outcomes are included in every step of the analysis, including in the fitting of the imputation model. It assumes that ICEs may impact post-ICE outcomes but that otherwise missingness is non-informative. The approach also assumes that the time-varying covariates do not contain missing values, deviations in outcomes after the ICE are correctly modeled by these time-varying covariates, and that sufficient post-ICE data are available to inform the regression coefficients of the time-varying covariates. These proposals are relatively recent and there remain open questions regarding the appropriate trade-off between model complexity (e.g. should the model account for a potentially differential effect on post-ICE outcomes depending on the timing of the ICE?) and the variance in the resulting treatment effect estimate. More generally, it is not yet established how much post-discontinuation data is required to implement such methods robustly and without the risk of substantial inflation of variance.

In some trial settings, only few subjects discontinue the randomized treatment. In other settings, treatment discontinuation rates are higher but it is difficult to retain subjects in the trial after treatment discontinuation leading to sparse data collection after treatment discontinuation. In both settings, the amount of available data after treatment discontinuation may be insufficient to inform an imputation model which explicitly models post-discontinuation data. Depending on the disease area and the anticipated mechanism of action of the intervention, it may be plausible to assume that subjects in the intervention group behave similarly to subjects in the control group after the ICE treatment discontinuation. In this case, reference-based imputation methods are an option (Mallinckrodt et al. (2020)). Reference-based imputation methods formalize the idea to impute missing data in the intervention group based on data from a control or reference group. For a general description and review of reference-based imputation methods, we refer to Carpenter, Roger, and Kenward (2013), Cro et al. (2020), I. White, Royes, and Best (2020) and Wolbers et al. (2022). For a technical description of the implemented statistical methodology for reference-based imputation, we refer to section 3 (in particular section 3.4).

### 2.2.4 Implementation of the composite strategy

The composite strategy is typically applied to binary or time-to-event outcomes but it can also be used for continuous outcomes by ascribing a suitably unfavorable value to patients who experience ICEs for which a composite strategy has been defined. One possibility to implement this is to use MI with a $$\delta$$-adjustment for post-ICE data as described in Darken et al. (2020).

# 3 Statistical methodology

## 3.1 Overview of the imputation procedure

Analyses of datasets with missing data always rely on missing data assumptions. The methods described here can be used to produce valid imputations under a MAR assumption or under reference-based imputation assumptions. MNAR imputation based on fixed $$\delta$$-adjustments as typically used in sensitivity analyses such as tipping-point analyses are also supported.

Three general imputation approaches are implemented in rbmi:

1. Conventional MI based on Bayesian (or approximate Bayesian) posterior draws from the imputation model combined with Rubin’s rules for inference as described in Carpenter, Roger, and Kenward (2013) and Cro et al. (2020).

2. Conditional mean imputation based on the REML estimate of the imputation model combined with resampling techniques (the jackknife or the bootstrap) for inference as described in Wolbers et al. (2022).

3. Bootstrapped MI methods based on REML estimates of the imputation model as described in von Hippel and Bartlett (2021).

### 3.1.1 Conventional MI

Conventional MI approaches include the following steps:

1. Base imputation model fitting step (Section 3.3)
• Fit a Bayesian multivariate normal mixed model for repeated measures (MMRM) to the observed longitudinal outcomes after exclusion of data after ICEs for which reference-based missing data imputation is desired (Section 3.3.3). Draw $$M$$ posterior samples of the estimated parameters (regression coefficients and covariance matrices) from this model.

• Alternatively, $$M$$ approximate posterior draws from the posterior distribution can be sampled by repeatedly applying conventional restricted maximum-likelihood (REML) parameter estimation of the MMRM model to nonparametric bootstrap samples from the original dataset (Section 3.3.4).

1. Imputation step (Section 3.4)
• Take a single sample $$m$$ ($$m\in 1,\ldots, M)$$ from the posterior distribution of the imputation model parameters.

• For each subject, use the sampled parameters and the defined imputation strategy to determine the mean and covariance matrix describing the subject’s marginal outcome distribution for all longitudinal outcome assessments (i.e. observed and missing outcomes).

• For each subjects, construct the conditional multivariate normal distribution of their missing outcomes given their observed outcomes (including observed outcomes after ICEs for which a reference-based assumption is desired).

• For each subject, draw a single sample from this conditional distribution to impute their missing outcomes leading to a complete imputed dataset.

• For sensitivity analyses, a pre-defined $$\delta$$-adjustment may be applied to the imputed data prior to the analysis step. (Section 3.5).

1. Analysis step (Section 3.6)
• Analyze the imputed dataset using an analysis model (e.g. ANCOVA) resulting in a point estimate and a standard error (with corresponding degrees of freedom) of the treatment effect.
1. Pooling step for inference (Section 3.7)
• Repeat steps 2. and 3. for each posterior sample $$m$$, resulting in $$M$$ complete datasets, $$M$$ point estimates of the treatment effect, and $$M$$ standard errors (with corresponding degrees of freedom). Pool the $$M$$ treatment effect estimates, standard errors, and degrees of freedom using the rules by Barnard and Rubin to obtain the final pooled treatment effect estimator, standard error, and degrees of freedom.

### 3.1.2 Conditional mean imputation

The conditional mean imputation approach includes the following steps:

1. Base imputation model fitting step (Section 3.3)
• Fit a conventional multivariate normal/MMRM model using restricted maximum likelihood (REML) to the observed longitudinal outcomes after exclusion of data after ICEs for which reference-based missing data imputation is desired (Section 3.3.2).
1. Imputation step (Section 3.4)
• For each subject, use the fitted parameters from step 1. to construct the conditional distribution of missing outcomes given observed outcomes (including observed outcomes after ICEs for which reference-based missing data imputation is desired) as described above.

• For each subject, impute their missing data deterministically by the mean of this conditional distribution leading to a complete imputed dataset.

• For sensitivity analyses, a pre-defined $$\delta$$-adjustment may be applied to the imputed data prior to the analysis step. (Section 3.5).

1. Analysis step (Section 3.6)
• Apply an analysis model (e.g. ANCOVA) to the completed dataset resulting in a point estimate of the treatment effect.
1. Jackknife or bootstrap inference step (Section 3.8)
• Inference for the treatment effect estimate from 3. is based on re-sampling techniques. Both the jackknife and the bootstrap are supported. Importantly, these methods require repeating all steps of the imputation procedure (i.e. imputation, conditional mean imputation, and analysis steps) on each of the resampled datasets.

### 3.1.3 Bootstrapped MI

The bootstrapped MI approach includes the following steps:

1. Base imputation model fitting step (Section 3.3)
• Apply conventional restricted maximum-likelihood (REML) parameter estimation of the MMRM model to $$B$$ nonparametric bootstrap samples from the original dataset using the observed longitudinal outcomes after exclusion of data after ICEs for which reference-based missing data imputation is desired.
1. Imputation step (Section 3.4)
• Take a bootstrapped dataset $$b$$ ($$b\in 1,\ldots, B)$$ and its corresponding imputation model parameter estimates.

• For each subject (from the bootstrapped dataset), use the parameter estimates and the defined strategy for dealing with their ICEs to determine the mean and covariance matrix describing the subject’s marginal outcome distribution for all longitudinal outcome assessments (i.e. observed and missing outcomes).

• For each subjects (from the bootstrapped dataset), construct the conditional multivariate normal distribution of their missing outcomes given their observed outcomes (including observed outcomes after ICEs for which reference-based missing data imputation is desired).

• For each subject (from the bootstrapped dataset), draw $$D$$ samples from this conditional distributions to impute their missing outcomes leading to $$D$$ complete imputed dataset for bootstrap sample $$b$$.

• For sensitivity analyses, a pre-defined $$\delta$$-adjustment may be applied to the imputed data prior to the analysis step. (Section 3.5).

1. Analysis step (Section 3.6)
• Analyze each of the $$B\times D$$ imputed datasets using an analysis model (e.g. ANCOVA) resulting in $$B\times D$$ point estimates of the treatment effect.
1. Pooling step for inference (Section 3.9)
• Pool the $$B\times D$$ treatment effect estimates as described in von Hippel and Bartlett (2021) to obtain the final pooled treatment effect estimate, standard error, and degrees of freedom.

## 3.2 Setting, notation, and missing data assumptions

Assume that the data are from a study with $$n$$ subjects in total and that each subject $$i$$ ($$i=1,\ldots,n$$) has $$J$$ scheduled follow-up visits at which the outcome of interest is assessed. In most applications, the data will be from a randomized trial of an intervention vs a control group and the treatment effect of interest is a comparison in outcomes at a specific visit between these randomized groups. However, single-arm trials or multi-arm trials are in principle also supported by the rbmi implementation.

Denote the observed outcome vector of length $$J$$ for subject $$i$$ by $$Y_i$$ (with missing assessments coded as NA (not available)) and its non-missing and missing components by $$Y_{i!}$$ and $$Y_{i?}$$, respectively. By default, imputation of missing outcomes in $$Y_{i}$$ is performed under a MAR assumption in rbmi. Therefore, if missing data following an ICE are to be handled using MAR imputation, this is compatible with the default assumption. As discussed in Section 2, the MAR assumption is often a good starting point for implementing a hypothetical strategy. But also note that observed outcome data after an ICE handled using a hypothetical strategy is not compatible with this strategy. Therefore, we assume that all post-ICE data after ICEs handled using a hypothetical strategy are already set to NA in $$Y_i$$ prior calling any rbmi functions. However, any observed outcomes after ICEs handled using a treatment policy strategy should be included in $$Y_i$$ as they are compatible with this strategy.

Subjects may also experience up to one ICE after which missing data imputation according to a reference-based imputation method is foreseen. For a subject $$i$$ with such an ICE, denote their first visit which is affected by the ICE by $$\tilde{t}_i \in \{1,\ldots,J\}$$. For all other subjects, set $$\tilde{t}_i=\infty$$. A subject’s outcome vector after setting observed outcomes from visit $$\tilde{t}_i$$ onwards to missing (i.e. NA) is denoted as $$Y'_i$$ and the corresponding data vector after removal of NA elements as $$Y'_{i!}$$.

MNAR $$\delta$$-adjustments are added to the imputed datasets after the formal imputation steps. This is covered in a separate section (Section 3.5).

## 3.3 The base imputation model

### 3.3.1 Included data and model specification

The purpose of the imputation model is to estimate (covariate-dependent) mean trajectories and covariance matrices for each group in the absence of ICEs handled using reference-based imputation methods. Conventionally, publications on reference-based imputation methods have implicitly assumed that the corresponding post-ICE data is missing for all subjects (Carpenter, Roger, and Kenward (2013)). We also allow the situation where post-ICE data is available for some subjects but needs to be imputed using reference-based methods for others. However, any observed data after ICEs for which reference-based imputation methods are specified is not compatible with the imputation model described below and they are therefore removed and considered as missing for the purpose of estimating the imputation model, and for this purpose only. For example, if a patient has an ICE addressed with a reference-based method but outcomes after the ICE are collected, these post-ICE outcomes will be excluded when fitting the base imputation model (but they will be included again in the following steps). That is, the base imputation model is fitted to $$Y'_{i!}$$ and not to $$Y_{i!}$$. If we did not exclude these data, then the imputation model would mistakenly estimate mean trajectories based on a mixture of observed pre- and post-ICE data which are not relevant for reference-based imputations.

Observed post-ICE outcomes in the control or reference group are also excluded from the base imputation model if the user specifies a reference-based imputation strategy for such ICEs. This ensures that an ICE has the same impact on the data included in the imputation model regardless whether the ICE occurred in the control or the intervention group. On the other hand, imputation in the reference group is based on a MAR assumption even for reference-based imputation methods and it may be preferable in some settings to include post-ICE data from the control group in the base imputation model. This can be implemented by specifying a MAR strategy for the ICE in the control group and a reference-based strategy for the same ICE in the intervention group.

The base imputation model of the longitudinal outcomes $$Y'_i$$ assumes that the mean structure is a linear function of covariates. Full flexibility for the specification of the linear predictor of the model is supported. At a minimum the covariates should include the treatment group, the (categorical) visit, and treatment-by-visit interactions. Typically, other covariates including the baseline outcome are also included. External time-varying covariates (e.g. calendar time of the visit) as well as internal time-varying (e.g. time-varying indicators of treatment discontinuation or initiation of rescue treatment) may in principle also be included if indicated (Guizzaro et al. (2021)). Missing covariate values are not allowed. This means that the values of time-varying covariates must be non-missing at every visit regardless of whether the outcome is measured or missing.

Denote the $$J\times p$$ design matrix for subject $$i$$ corresponding to the mean structure model by $$X_i$$ and the same matrix after removal of rows corresponding to missing outcomes in $$Y'_{i!}$$ by $$X'_{i!}$$. Here $$p$$ is the number of parameters in the mean structure of the model for the elements of $$Y'_{i!}$$. The base imputation model for the observed outcomes is defined as: $Y'_{i!} = X'_{i!}\beta + \epsilon_{i!} \mbox{ with } \epsilon_{i!}\sim N(0,\Sigma_{i!!})$ where $$\beta$$ is the vector of regression coefficients and $$\Sigma_{i!!}$$ is a covariance matrix which is obtained from the complete-data $$J\times J$$-covariance matrix $$\Sigma$$ by omitting rows and columns corresponding to missing outcome assessments for subject $$i$$.

Typically, a common unstructured covariance matrix for all subjects is assumed for $$\Sigma$$ but separate covariate matrices per treatment group are also supported. Indeed, the implementation also supports the specification of separate covariate matrices according to an arbitrarily defined categorical variable which groups the subjects into disjoint subset. For example, this could be useful if different covariance matrices are suspected in different subject strata. Finally, for all imputation methods described below that do not rely on Bayesian model fitting through MCMC, there is further flexibility in the choice of the covariance structure, i.e. unstructured (default), heterogeneous Toeplitz, heterogeneous compound symmetry, and AR(1) covariance structures are supported.

### 3.3.2 Restricted maximum likelihood estimation (REML)

Frequentist parameter estimation for the base imputation is based on REML. The use of REML as an improved alternative to maximum likelihood (ML) for covariance parameter estimation was originally proposed by Patterson and Thompson (1971). Since then, it has become the default method for parameter estimation in linear mixed effects models. rbmi allows to choose between ML and REML methods to estimate the model parameters, with REML being the default option.

### 3.3.3 Bayesian model fitting

The Bayesian imputation model is fitted with the R package rstan (Stan Development Team (2020)). rstan is the R interface of Stan. Stan is a powerful and flexible statistical software developed by a dedicated team and implements Bayesian inference with state-of-the-art MCMC sampling procedures. The multivariate normal model with missing data specified in section 3.3.1 can be considered a generalization of the models described in the Stan user’s guide (see Stan Development Team (2020, sec. 3.5)).

The same prior distributions as in the SAS implementation of the “five macros” are used (Roger (2021)), i.e. an improper flat priors for the regression coefficients and a weakly informative inverse Wishart prior for the covariance matrix (or matrices). Specifically, let $$S \in \mathbb{R}^{J \times J}$$ be a symmetric positive definite matrix and $$\nu \in (J-1, \infty)$$. Then the symmetric positive definite matrix $$x \in \mathbb{R}^{J \times J}$$ has density: $\text{InvWish}(x \vert \nu, S) = \frac{1}{2^{\nu J/2}} \frac{1}{\Gamma_J(\frac{\nu}{2})} \vert S \vert^{\nu/2} \vert x \vert ^{-(\nu + J + 1)/2} \text{exp}(-\frac{1}{2} \text{tr}(Sx^{-1})).$ For $$\nu > J+1$$ the mean is given by: $E[x] = \frac{S}{\nu - J - 1}.$ We choose $$S$$ equal to the estimated covariance matrix from the frequentist REML fit and $$\nu = J+2$$ as these are the lowest degrees of freedom that guarantee a finite mean. Setting the degrees of freedom with such a low $$\nu$$ ensures that the prior has little impact on the posterior. Moreover, this choice allows to interpret the parameter $$S$$ as the mean of the prior distribution.

As in the “five macros”, the MCMC algorithm is initialized at the parameters from a frequentist REML fit (see section 3.3.2). As described above, we are using only weakly informative priors for the parameters. Therefore, the Markov chain is essentially starting from the targeted stationary posterior distribution and only a minimal amount of burn-in of the chain is required.

### 3.3.4 Approximate Bayesian posterior draws via the bootstrap

Several authors have suggested that a stabler way to get Bayesian posterior draws from the imputation model is to bootstrap the incomplete data and to calculate REML estimates for each bootstrap sample (Little and Rubin (2002), Efron (1994), Honaker and King (2010), von Hippel and Bartlett (2021)). This method is proper in that the REML estimates from the bootstrap samples are asymptotically equivalent to a sample from the posterior distribution and may provide additional robustness to model misspecification (Little and Rubin (2002, sec. 10.2.3, part 6), Honaker and King (2010)). In order to retain balance between treatment groups and stratification factors across bootstrap samples, the user is able to provide stratification variables for the bootstrap in the rbmi implementation.

## 3.4 Imputation step

### 3.4.1 Marginal imputation distribution for a subject - MAR case

For each subject $$i$$, the marginal distribution of the complete $$J$$-dimensional outcome vector from all assessment visits according to the imputation model is a multivariate normal distribution. Its mean $$\tilde{\mu}_i$$ is given by the predicted mean from the imputation model conditional on the subject’s baseline characteristics, group, and, optionally, time-varying covariates. Its covariance matrix $$\tilde{\Sigma}_i$$ is given by the overall estimated covariance matrix or, if different covariance matrices are assumed for different groups, the covariance matrix corresponding to subject $$i$$’s group.

### 3.4.2 Marginal imputation distribution for a subject - reference-based imputation methods

For each subject $$i$$, we calculate the mean and covariance matrix of the complete $$J$$-dimensional outcome vector from all assessment visits as for the MAR case and denote them by $$\mu_i$$ and $$\Sigma_i$$. For reference-based imputation methods, a corresponding reference group is also required for each group. Typically, the reference group for the intervention group will be the control group. The reference mean $$\mu_{ref,i}$$ is defined as the predicted mean from the imputation model conditional on the reference group (rather than the actual group subject $$i$$ belongs to) and the subject’s baseline characteristics. The reference covariance matrix $$\Sigma_{ref,i}$$ is the overall estimated covariance matrix or, if different covariance matrices are assumed for different groups, the estimated covariance matrix corresponding to the reference group. In principle, time-varying covariates could also be included in reference-based imputation methods. However, this is only sensible for external time-varying covariates (e.g. calendar time of the visit) and not for internal time-varying covariates (e.g. treatment discontinuation) because the latter likely depend on the actual treatment group and it is typically not sensible to assume the same trajectory of the time-varying covariate for the reference group.

Based on these means and covariance matrices, the subject’s marginal imputation distribution for the reference-based imputation methods is then calculated as detailed in Carpenter, Roger, and Kenward (2013, sec. 4.3). Denote the mean and covariance matrix of this marginal imputation distribution by $$\tilde{\mu}_i$$ and $$\tilde{\Sigma}_i$$. Recall that the subject’s first visit which is affected by the ICE is denoted by $$\tilde{t}_i \in \{1,\ldots,J\}$$ (and visit $$\tilde{t}_i-1$$ is the last visit unaffected by the ICE). The marginal distribution for the patient $$i$$ is then built according to the specific assumption for the data up to and post the ICE as follows:

1. Jump to reference (JR): the patient’s outcome distribution is normally distributed with the following mean: $\tilde{\mu}_i = (\mu_i[1], \dots, \mu_i[\tilde{t}_i-1], \mu_{ref,i}[\tilde{t}_i], \dots, \mu_{ref,i}[J])^T.$ The covariance matrix is constructed as follows. First, we partition the covariance matrices $$\Sigma_i$$ and $$\Sigma_{ref,i}$$ in blocks according to the time of the ICE $$\tilde{t}_i$$: $\Sigma_{i} = \begin{bmatrix} \Sigma_{i, 11} & \Sigma_{i, 12} \\ \Sigma_{i, 21} & \Sigma_{i,22} \\ \end{bmatrix}$ $\Sigma_{ref,i} = \begin{bmatrix} \Sigma_{ref, i, 11} & \Sigma_{ref, i, 12} \\ \Sigma_{ref, i, 21} & \Sigma_{ref, i,22} \\ \end{bmatrix}.$ We want the covariance matrix $$\tilde{\Sigma}_i$$ to match $$\Sigma_i$$ for the pre-deviation measurements, and $$\Sigma_{ref,i}$$ for the conditional components for the post-deviation given the pre-deviation measurements. The solution is derived in Carpenter, Roger, and Kenward (2013, sec. 4.3) and is given by: $\begin{matrix} \tilde{\Sigma}_{i,11} = \Sigma_{i, 11} \\ \tilde{\Sigma}_{i, 21} = \Sigma_{ref,i, 21} \Sigma^{-1}_{ref,i, 11} \Sigma_{i, 11} \\ \tilde{\Sigma}_{i, 22} = \Sigma_{ref, i, 22} - \Sigma_{ref,i, 21} \Sigma^{-1}_{ref,i, 11} (\Sigma_{ref,i, 11} - \Sigma_{i,11}) \Sigma^{-1}_{ref,i, 11} \Sigma_{ref,i, 12}. \end{matrix}$

2. Copy increments in reference (CIR): the patient’s outcome distribution is normally distributed with the following mean: $\begin{split} \tilde{\mu}_i =& (\mu_i[1], \dots, \mu_i[\tilde{t}_i-1], \mu_i[\tilde{t}_i-1] + (\mu_{ref,i}[\tilde{t}_i] - \mu_{ref,i}[\tilde{t}_i-1]), \dots,\\ & \mu_i[\tilde{t}_i-1]+(\mu_{ref,i}[J] - \mu_{ref,i}[\tilde{t}_i-1]))^T. \end{split}$ The covariance matrix is derived as for the JR method.

3. Copy reference (CR): the patient’s outcome distribution is normally distributed with mean and covariance matrix taken from the reference group: $\tilde{\mu}_i = \mu_{ref,i}$ $\tilde{\Sigma}_i = \Sigma_{ref,i}.$

4. Last mean carried forward (LMCF): the patient’s outcome distribution is normally distributed with the following mean: $\tilde{\mu}_i = (\mu_i[1], \dots, \mu_i[\tilde{t}_i-1], \mu_i[\tilde{t}_i-1], \dots, \mu_i[\tilde{t}_i-1])'$ and covariance matrix: $\tilde{\Sigma}_i = \Sigma_i.$

### 3.4.3 Imputation of missing outcome data

The joint marginal multivariate normal imputation distribution of subject $$i$$’s observed and missing outcome data has mean $$\tilde{\mu}_i$$ and covariance matrix $$\tilde{\Sigma}_i$$ as defined above. The actual imputation of the missing outcome data is obtained by conditioning this marginal distribution on the subject’s observed outcome data. Of note, this approach is valid regardless whether the subject has intermittent or terminal missing data.

The conditional distribution used for the imputation is again a multivariate normal distribution and explicit formulas for the conditional mean and covariance are readily available. For completeness, we report them here with the notation and terminology of our setting. The marginal distribution for the outcome of patient $$i$$ is $$Y_i \sim N(\tilde{\mu}_i, \tilde{\Sigma}_i)$$ and the outcome $$Y_i$$ can be decomposed in the observed ($$Y_{i,!}$$) and the unobserved ($$Y_{i,?}$$) components. Analogously the mean $$\tilde{\mu}_i$$ can be decomposed as $$(\tilde{\mu}_{i,!},\tilde{\mu}_{i,?})$$ and the covariance $$\tilde{\Sigma}_i$$ as: $\tilde{\Sigma}_i = \begin{bmatrix} \tilde{\Sigma}_{i, !!} & \tilde{\Sigma}_{i,!?} \\ \tilde{\Sigma}_{i, ?!} & \tilde{\Sigma}_{i, ??} \end{bmatrix}.$ The conditional distribution of $$Y_{i,?}$$ conditional on $$Y_{i,!}$$ is then a multivariate normal distribution with expectation $E(Y_{i,?} \vert Y_{i,!})= \tilde{\mu}_{i,?} + \tilde{\Sigma}_{i, ?!} \tilde{\Sigma}_{i,!!}^{-1} (Y_{i,!} - \tilde{\mu}_{i,!})$ and covariance matrix $Cov(Y_{i,?} \vert Y_{i,!}) = \tilde{\Sigma}_{i,??} - \tilde{\Sigma}_{i,?!} \tilde{\Sigma}_{i,!!}^{-1} \tilde{\Sigma}_{i,!?}.$

Conventional random imputation consists in sampling from this conditional multivariate normal distribution. Conditional mean imputation imputes missing values with the deterministic conditional expectation $$E(Y_{i,?} \vert Y_{i,!})$$.

## 3.5$$\delta$$-adjustment

A marginal $$\delta$$-adjustment approach similar to the “five macros” in SAS is implemented (Roger (2021)), i.e. fixed non-stochastic values are added after the multivariate normal imputation step and prior to the analysis. This is relevant for sensitivity analyses in order to make imputed data systematically worse or better, respectively, than observed data. In addition, some authors have suggested $$\delta$$-type adjustments to implement a composite strategy for continuous outcomes (Darken et al. (2020)).

The implementation provides full flexibility regarding the specific implementation of the $$\delta$$-adjustment, i.e. the value that is added may depend on the randomized treatment group, the timing of the subject’s ICE, and other factors. For suggestions and case studies regarding this topic, we refer to Cro et al. (2020).

## 3.6 Analysis step

After data imputation, a standard analysis model can be applied to the completed data resulting in a treatment effect estimate. As the imputed data no longer contains missing values, the analysis model is often simple. For example, it can be an analysis of covariance (ANCOVA) model with the outcome (or the change in the outcome from baseline) at a specific visit j as the dependent variable, the randomized treatment group as the primary covariate and, typically, adjustment for the same baseline covariates as for the imputation model.

## 3.7 Pooling step for inference of (approximate) Bayesian MI and Rubin’s rules

Assume that the analysis model has been applied to $$M$$ multiple imputed random datasets which resulted in $$m$$ treatment effect estimates $$\hat{\theta}_m$$ ($$m=1,\ldots,M$$) with corresponding standard error $$SE_m$$ and (if available) degrees of freedom $$\nu_{com}$$. If degrees of freedom are not available for an analysis model, set $$\nu_{com}=\infty$$ for inference based on the normal distribution.

Rubin’s rules are used for pooling the treatment effect estimates and corresponding variances estimates from the analysis steps across the $$M$$ multiple imputed datasets. According to Rubin’s rules, the final estimate of the treatment effect is calculated as the sample mean over the $$M$$ treatment effect estimates: $\hat{\theta} = \frac{1}{M} \sum_{m = 1}^M \hat{\theta}_m.$ The pooled variance is based on two components that reflect the within and the between variance of the treatment effects across the multiple imputed datasets: $V(\hat{\theta}) = V_W(\hat{\theta}) + (1 + \frac{1}{M}) V_B(\hat{\theta})$ where $$V_W(\hat{\theta}) = \frac{1}{M}\sum_{m = 1}^M SE^2_m$$ is the within-variance and $$V_B(\hat{\theta}) = \frac{1}{M-1} \sum_{m = 1}^M (\hat{\theta}_m - \hat{\theta})^2$$ is the between-variance.

Confidence intervals and tests of the null hypothesis $$H_0: \theta=\theta_0$$ are based on the $$t$$-statistics $$T$$:

$T= (\hat{\theta}-\theta_0)/\sqrt{V(\hat{\theta})}.$ Under the null hypothesis, $$T$$ has an approximate $$t$$-distribution with $$\nu$$ degrees of freedom. $$\nu$$ is calculated according to the Barnard and Rubin approximation, see Barnard and Rubin (1999) (formula 3) or Little and Rubin (2002) (formula (5.24), page 87):

$\nu = \frac{\nu_{old}* \nu_{obs}}{\nu_{old} + \nu_{obs}}$ with $\nu_{old} = \frac{M-1}{\lambda^2} \quad\mbox{and}\quad \nu_{obs} = \frac{\nu_{com} + 1}{\nu_{com} + 3} \nu_{com} (1 - \lambda)$ where $$\lambda = \frac{(1 + \frac{1}{M})V_B(\hat{\theta})}{V(\hat{\theta})}$$ is the fraction of missing information.

## 3.8 Bootstrap and jackknife inference for conditional mean imputation

### 3.8.1 Point estimate of the treatment effect

The point estimator is obtained by applying the analysis model (Section 3.6) to a single conditional mean imputation of the missing data (see Section 3.4.3) based on the REML estimator of the parameters of the imputation model (see Section 3.3.2). We denote this treatment effect estimator by $$\hat{\theta}$$.

As demonstrated in Wolbers et al. (2022) (Section 2.4), this treatment effect estimator is valid if the analysis model is an ANCOVA model or, more generally, if the treatment effect estimator is a linear function of the imputed outcome vector. Indeed, if this is the case, then the estimator is identical to the pooled treatment effect across multiple random REML imputation with an infinite number of imputations and corresponds to a computationally efficient implementation of a proposal by von Hippel and Bartlett (2021). We expect that the conditional mean imputation method is also applicable to some other analysis models (e.g. for general MMRM analysis models) but this has not been formally justified.

### 3.8.2 Jackknife standard errors, confidence intervals (CI) and tests for the treatment effect

For a dataset containing $$n$$ subjects, the jackknife standard error depends on treatment effect estimates $$\hat{\theta}_{(-b)}$$ ($$b=1,\ldots,n$$) from samples of the original dataset which leave out the observation from subject $$b$$. As described previously, to obtain treatment effect estimates for leave-one-subject-out datasets, all steps of the imputation procedure (i.e. imputation, conditional mean imputation, and analysis steps) need to be repeated on this new dataset.

Then, the jackknife standard error is defined as $\hat{se}_{jack}=[\frac{(n-1)}{n}\cdot\sum_{b=1}^{n} (\hat{\theta}_{(-b)}-\bar{\theta}_{(.)})^2]^{1/2}$ where $$\bar{\theta}_{(.)}$$ denotes the mean of all jackknife estimates (Efron and Tibshirani (1994), chapter 10). The corresponding two-sided normal approximation $$1-\alpha$$ CI is defined as $$\hat{\theta}\pm z^{1-\alpha/2}\cdot \hat{se}_{jack}$$ where $$\hat{\theta}$$ is the treatment effect estimate from the original dataset. Tests of the null hypothesis $$H_0: \theta=\theta_0$$ are then based on the $$Z$$-score $$Z=(\hat{\theta}-\theta_0)/\hat{se}_{jack}$$ using a standard normal approximation.

A simulation study reported in Wolbers et al. (2022) demonstrated exact protection of the type I error for jackknife-based inference with a relatively low sample size (n = 100 per group) and a substantial amount of missing data (>25% of subjects with an ICE).

### 3.8.3 Bootstrap standard errors, confidence intervals (CI) and tests for the treatment effect

As an alternative to the jackknife, the bootstrap has also been implemented in rbmi (Efron and Tibshirani (1994), Davison and Hinkley (1997)).

Two different bootstrap methods are implemented in rbmi: Methods based on the bootstrap standard error and the normal approximation and percentile bootstrap methods. Denote the treatment effect estimates from $$B$$ bootstrap samples by $$\hat{\theta}^*_b$$ ($$b=1,\ldots,B$$). The bootstrap standard error $$\hat{se}_{boot}$$ is defined as the empirical standard deviation of the bootstrapped treatment effect estimates. Confidence intervals and tests based on the bootstrap standard error can then be constructed in the same way as for the jackknife. Confidence intervals using the percentile bootstrap are based on empirical quantiles of the bootstrap distribution and corresponding statistical tests are implemented in rbmi via inversion of the confidence interval. Explicit formulas for bootstrap inference as implemented in the rbmi package and some considerations regarding the required number of bootstrap samples are included in the Appendix of Wolbers et al. (2022).

A simulation study reported in Wolbers et al. (2022) demonstrated a small inflation of the type I error rate for inference based on the bootstrap standard error (up to $$5.3\%$$ for a nominal type I error rate of $$5\%$$) for a sample size of n = 100 per group and a substantial amount of missing data (>25% of subjects with an ICE). Based on this simulations, we recommend the jackknife over the bootstrap for inference because it performed better in our simulation study and is typically much faster to compute than the bootstrap.

## 3.9 Pooling step for inference of the bootstrapped MI methods

Assume that the analysis model has been applied to $$B\times D$$ multiple imputed random datasets which resulted in $$B\times D$$ treatment effect estimates $$\hat{\theta}_{bd}$$ ($$b=1,\ldots,B$$; $$d=1,\ldots,D$$).

The final estimate of the treatment effect is calculated as the sample mean over the $$B*D$$ treatment effect estimates: $\hat{\theta} = \frac{1}{BD} \sum_{b = 1}^B \sum_{d = 1}^D \hat{\theta}_{bd}.$ The pooled variance is based on two components that reflect the variability within and between imputed bootstrap samples (von Hippel and Bartlett (2021), formula 8.4): $V(\hat{\theta}) = (1 + \frac{1}{B})\frac{MSB - MSW}{D} + \frac{MSW}{BD}$

where $$MSB$$ is the mean square between the bootstrapped datasets, and $$MSW$$ is the mean square within the bootstrapped datasets and between the imputed datasets:

\begin{align*} MSB &= \frac{D}{B-1} \sum_{b = 1}^B (\bar{\theta_{b}} - \hat{\theta})^2 \\ MSW &= \frac{1}{B(D-1)} \sum_{b = 1}^B \sum_{d = 1}^D (\theta_{bd} - \bar{\theta_b})^2 \end{align*} where $$\bar{\theta_{b}}$$ is the mean across the $$D$$ estimates obtained from random imputation of the $$b$$-th bootstrap sample.

The degrees of freedom are estimated with the following formula (von Hippel and Bartlett (2021), formula 8.6):

$\nu = \frac{(MSB\cdot (B+1) - MSW\cdot B)^2}{\frac{MSB^2\cdot (B+1)^2}{B-1} + \frac{MSW^2\cdot B}{D-1}}$

Confidence intervals and tests of the null hypothesis $$H_0: \theta=\theta_0$$ are based on the $$t$$-statistics $$T$$:

$T= (\hat{\theta}-\theta_0)/\sqrt{V(\hat{\theta})}.$ Under the null hypothesis, $$T$$ has an approximate $$t$$-distribution with $$\nu$$ degrees of freedom.

## 3.10 Comparison between the implemented approaches

### 3.10.1 Treatment effect estimation

All approaches provide consistent treatment effect estimates for standard and reference-based imputation methods in case the analysis model of the completed datasets is a general linear model such as ANCOVA. Methods other than conditional mean imputation should also be valid for other analysis models. The validity of conditional mean imputation has only been formally demonstrated for analyses using the general linear model (Wolbers et al. (2022, sec. 2.4)) though it may also be applicable more widely (e.g. for general MMRM analysis models).

Treatment effects based on conditional mean imputation are deterministic. All other methods are affected by Monte Carlo sampling error and the precision of estimates depends on the number of imputations or bootstrap samples, respectively.

### 3.10.2 Standard errors of the treatment effect

All approaches provide frequentist consistent estimates of the standard error for imputation under a MAR assumption. For reference-based imputation methods, methods based on conditional mean imputation or bootstrapped MI provide frequentist consistent estimates of the standard error whereas Rubin’s rules applied to conventional MI methods provides so-called information anchored inference (Bartlett (2021), Cro, Carpenter, and Kenward (2019), von Hippel and Bartlett (2021), Wolbers et al. (2022)). Frequentist consistent estimates of the standard error lead to confidence intervals and tests which have (asymptotically) correct coverage and type I error control under the assumption that the reference-based assumption reflects the true data-generating mechanism. For finite samples, simulations for a sample size of $$n=100$$ per group reported in Wolbers et al. (2022) demonstrated that conditional mean imputation combined with the jackknife provided exact protection of the type one error rate whereas the bootstrap was associated with a small type I error inflation (between 5.1% to 5.3% for a nominal level of 5%).

It is well known that Rubin’s rules do not provide frequentist consistent estimates of the standard error for reference-based imputation methods (Seaman, White, and Leacy (2014), Liu and Pang (2016), Tang (2017), Cro, Carpenter, and Kenward (2019), Bartlett (2021)). Standard errors from Rubin’s rule are typically larger than frequentist standard error estimates leading to conservative inference and a corresponding loss of statistical power, see e.g. the simulations reported in Wolbers et al. (2022). Intuitively, this occurs because reference-based imputation methods borrow information from the reference group for imputations in the intervention group leading to a reduction in the frequentist variance of the resulting treatment effect contrast which is not captured by Rubin’s variance estimator. Formally, this occurs because the imputation and analysis models are uncongenial for reference-based imputation methods (Meng (1994), Bartlett (2021)). Cro, Carpenter, and Kenward (2019) argued that Rubin’s rule is nevertheless valid for reference-based imputation methods because it is approximately information-anchored, i.e. that the proportion of information lost due to missing data under MAR is approximately preserved in reference-based analyses. In contrast, frequentist standard errors for reference based imputation are not information anchored for reference-based imputation and standard errors under reference-based assumptions are typically smaller than those for MAR imputation.

Information anchoring is a sensible concept for sensitivity analyses, whereas for a primary analyses, it may be more important to adhere to the principles of frequentist inference. Analyses of data with missing observations generally rely on unverifiable missing data assumptions and the assumptions for reference-based imputation methods are relatively strong. Therefore, these assumptions need to be clinically justified as appropriate or at least conservative for the considered disease area and the anticipated mechanism of action of the intervention.

Conditional mean imputation combined with the jackknife is the only method which leads to deterministic standard error estimates and, consequently, confidence intervals and $$p$$-values are also deterministic. This is particularly important in a regulatory setting where it is important to ascertain whether a calculated $$p$$-value which is close to the critical boundary of 5% is truly below or above that threshold rather than being uncertain about this because of Monte Carlo error.

### 3.10.3 Computational complexity

Bayesian MI methods rely on the specification of prior distributions and the usage of Markov chain Monte Carlo (MCMC) methods. All other methods based on multiple imputation or bootstrapping require no other tuning parameters than the specification of the number of imputations $$M$$ or bootstrap samples $$B$$ and rely on numerical optimization for fitting the MMRM imputation models via REML. Conditional mean imputation combined with the jackknife has no tuning parameters.

In our rbmi implementation, the fitting of the MMRM imputation model via REML is computationally most expensive. MCMC sampling using rstan (Stan Development Team (2020)) is typically relatively fast in our setting and requires only a small burn-in and burn-between of the chains. In addition, the number of random imputations for reliable inference using Rubin’s rules is often smaller than the number of resamples required for the jackknife or the bootstrap (see e.g. the discussions in I. R. White, Royston, and Wood (2011, sec. 7) for Bayesian MI and the Appendix of Wolbers et al. (2022) for the bootstrap). Thus, for many applications, we expect that conventional MI based on Bayesian posterior draws will be fastest, followed by conventional MI using approximate Bayesian posterior draws and conditional mean imputation combined with the jackknife. Conditional mean imputation combined with the bootstrap and bootstrapped MI methods will typically be most computationally demanding. Of note, all implemented methods are conceptually straightforward to parallelize and some parallelization support is provided by rbmi.

# 4 Mapping of statistical methods to rbmi functions

For a full documentation of the rbmi package functionality we refer to the help pages of all functions and to the other package vignettes. Here we only give a brief overview of how the different steps of the imputation procedure are mapped to rbmi functions:

• The base imputation model fitting step is implemented in the function draws(). The chosen MI approach can be set using the argument method and should be one of the following:
• Bayesian posterior parameter draws from the imputation model are obtained via the argument method = method_bayes().
• Approximate Bayesian posterior parameter draws from the imputation model are obtained via argument method = method_approxbayes().
• ML or REML parameter estimates of the imputation model parameters for the original dataset and all leave-one-subject-out datasets (as required for the jackknife) are obtained via argument method = method_condmean(type = "jackknife").
• ML or REML parameter estimates of the imputation model parameters for the original dataset and bootstrapped datasets are obtained via argument method = method_condmean(type = "bootstrap").
• Bootstrapped MI methods are obtained via argument method = method_bmlmi(B=B, D=D) where $$B$$ refers to the number of bootstrap samples and $$D$$ to the number of random imputations for each bootstrap sample.
• The imputation step using random imputation or deterministic conditional mean imputation, respectively, is implemented in function impute(). Imputation can be performed assuming the already implemented imputation strategies as presented in section 3.4. Additionally, user-defined imputation strategies are also supported.
• The analysis step is implemented in function analyse() and applies the analysis model to all imputed datasets. By default, the analysis model (argument fun) is the ancova() function but alternative analysis functions can also be provided by the user. The analyse() function also allows for $$\delta$$-adjustments to the imputed datasets prior to the analysis via argument delta.
• The inference step is implemented in function pool() which pools the results across imputed datasets. The Rubin and Bernard rule is applied in case of (approximate) Bayesian MI. For conditional mean imputation, jackknife and bootstrap (normal approximation or percentile) inference is supported. For BMLMI, the pooling and inference steps are performed via pool() which in this case implements the method described in Section 3.9.

# 5 Comparison to other software implementations

An established software implementation of reference-based imputation in SAS are the so-called “five macros” by James Roger (Roger (2021)). An alternative R implementation which is also currently under development is the R package RefBasedMI (McGrath and White (2021)).

rbmi has several features which are not supported by the other implementations:

1. In addition to the Bayesian MI approach implemented also in the other packages, our implementation provides three alternative MI approaches: approximate Bayesian MI, conditional mean imputation combined with resampling, and bootstrapped MI.

2. rbmi allows for the usage of data collected after an ICE. For example, suppose that we want to adopt a treatment policy strategy for the ICE “treatment discontinuation”. A possible implementation of this strategy is to use the observed outcome data for subjects who remain in the study after the ICE and to use reference-based imputation in case the subject drops out. In our implementation, this is implemented by excluding observed post ICE data from the imputation model which assumes MAR missingness but including them in the analysis model. To our knowledge, this is not directly supported by the other implementations.

3. RefBasedMI fits the imputation model to data from each treatment group separately which implies covariate-treatment group interactions for all covariates for the pooled data from both treatment groups. In contrast, Roger’s five macros assume a joint model including data from all the randomized groups and covariate-treatment interactions covariates are not allowed. We also chose to implement a joint model but use a flexible model for the linear predictor which may or may not include an interaction term between any covariate and the treatment group. In addition, our imputation model also allows for the inclusion of time-varying covariates.

4. In our implementation, the grouping of the subjects for the purpose of the imputation model (and the definition of the reference group) does not need to correspond to the assigned treatment groups. This provides additional flexibility for the imputation procedure. It is not clear to us whether this feature is supported by Roger’s five macros or RefBasedMI.

5. We believe that our R-based implementation is more modular than RefBasedMI which should facilitate further package enhancements.

In contrast, the more general causal model introduced by I. White, Royes, and Best (2020) is available in the other implementations but is currently not supported by ours.

# References

Barnard, John, and Donald B Rubin. 1999. “Miscellanea. Small-Sample Degrees of Freedom with Multiple Imputation.” Biometrika 86 (4): 948–55.
Bartlett, Jonathan W. 2021. “Reference-Based Multiple Imputation - What Is the Right Variance and How to Estimate It.” Statistics in Biopharmaceutical Research. https://doi.org/10.1080/19466315.2021.1983455.
Carpenter, James R, James H Roger, and Michael G Kenward. 2013. “Analysis of Longitudinal Trials with Protocol Deviation: A Framework for Relevant, Accessible Assumptions, and Inference via Multiple Imputation.” Journal of Biopharmaceutical Statistics 23 (6): 1352–71.
Cro, Suzie, James R Carpenter, and Michael G Kenward. 2019. “Information-Anchored Sensitivity Analysis: Theory and Application.” Journal of the Royal Statistical Society: Series A (Statistics in Society) 182 (2): 623–45.
Cro, Suzie, Tim P Morris, Michael G Kenward, and James R Carpenter. 2020. “Sensitivity Analysis for Clinical Trials with Missing Continuous Outcome Data Using Controlled Multiple Imputation: A Practical Guide.” Statistics in Medicine 39 (21): 2815–42.
Darken, Patrick, Jack Nyberg, Shaila Ballal, and David Wright. 2020. “The Attributable Estimand: A New Approach to Account for Intercurrent Events.” Pharmaceutical Statistics 19 (5): 626–35.
Davison, Anthony C, and David V Hinkley. 1997. Bootstrap Methods and Their Application. Cambridge University Press.
Efron, Bradley. 1994. “Missing Data, Imputation, and the Bootstrap.” Journal of the American Statistical Association 89 (426): 463–75.
Efron, Bradley, and Robert J Tibshirani. 1994. An Introduction to the Bootstrap. CRC press.
Guizzaro, Lorenzo, Frank Pétavy, Robin Ristl, and Ciro Gallo. 2021. “The Use of a Variable Representing Compliance Improves Accuracy of Estimation of the Effect of Treatment Allocation Regardless of Discontinuation in Trials with Incomplete Follow-up.” Statistics in Biopharmaceutical Research 13 (1): 119–27.
Honaker, James, and Gary King. 2010. “What to Do about Missing Values in Time-Series Cross-Section Data.” American Journal of Political Science 54 (2): 561–81.
ICH E9 working group. 2019. ICH E9 (R1): Addendum on estimands and sensitivity analysis in clinical trials to the guideline on statistical principles for clinical trials. International Council for Harmonisation of Technical Requirements for Pharmaceuticals for Human Use. 2019. https://database.ich.org/sites/default/files/E9-R1_Step4_Guideline_2019_1203.pdf.
Little, Roderick JA, and Donald B Rubin. 2002. Statistical Analysis with Missing Data, Second Edition. John Wiley & Sons.
Liu, G Frank, and Lei Pang. 2016. “On Analysis of Longitudinal Clinical Trials with Missing Data Using Reference-Based Imputation.” Journal of Biopharmaceutical Statistics 26 (5): 924–36.
Mallinckrodt, CH, J Bell, G Liu, B Ratitch, M O’Kelly, I Lipkovich, P Singh, L Xu, and G Molenberghs. 2020. “Aligning Estimators with Estimands in Clinical Trials: Putting the ICH E9 (R1) Guidelines into Practice.” Therapeutic Innovation & Regulatory Science 54 (2): 353–64.
McGrath, Kevin, and Ian White. 2021. “RefBasedMI: Reference-Based Imputation for Longitudinal Clinical Trials with Protocol Deviation.” https://github.com/UCL/RefbasedMI.
Meng, Xiao-Li. 1994. “Multiple-Imputation Inferences with Uncongenial Sources of Input.” Statistical Science 9 (4): 538–58.
Patterson, H Desmond, and Robin Thompson. 1971. “Recovery of Inter-Block Information When Block Sizes Are Unequal.” Biometrika 58 (3): 545–54.
Polverejan, Elena, and Vladimir Dragalin. 2020. “Aligning Treatment Policy Estimands and Estimators—a Simulation Study in Alzheimer’s Disease.” Statistics in Biopharmaceutical Research 12 (2): 142–54.
Roger, James. 2021. “Reference-Based MI via Multivariate Normal RM (the ‘Five Macros’ and MIWithD).” https://www.lshtm.ac.uk/research/centres-projects-groups/missing-data#dia-missing-data.
Seaman, Shaun R, Ian R White, and Finbarr P Leacy. 2014. “Comment on Analysis of Longitudinal Trials with Protocol Deviations: A Framework for Relevant, Accessible Assumptions, and Inference via Multiple Imputation,’ by Carpenter, Roger, and Kenward.” Journal of Biopharmaceutical Statistics 24 (6): 1358–62.
Stan Development Team. 2020. RStan: The R Interface to Stan.” https://mc-stan.org/.
Tang, Yongqiang. 2017. “On the Multiple Imputation Variance Estimator for Control-Based and Delta-Adjusted Pattern Mixture Models.” Biometrics 73 (4): 1379–87.
von Hippel, Paul T, and Jonathan W Bartlett. 2021. “Maximum Likelihood Multiple Imputation: Faster Imputations and Consistent Standard Errors Without Posterior Draws.” Statistical Science 36 (3): 400–420.
White, Ian R, Patrick Royston, and Angela M Wood. 2011. “Multiple Imputation Using Chained Equations: Issues and Guidance for Practice.” Statistics in Medicine 30 (4): 377–99.
White, Ian, Joseph Royes, and Nicky Best. 2020. “A Causal Modelling Framework for Reference-Based Imputation and Tipping Point Analysis in Clinical Trials with Quantitative Outcome.” Journal of Biopharmaceutical Statistics 30 (2): 334–50.
Wolbers, Marcel, Alessandro Noci, Paul Delmar, Craig Gower-Page, Sean Yiu, and Jonathan W Bartlett. 2022. “Standard and Reference-Based Conditional Mean Imputation.” https://arxiv.org/abs/2109.11162.


text:
# Finding Power of Factorial Divisor

You are given two numbers $n$ and $k$. Find the largest power of $k$ $x$ such that $n!$ is divisible by $k^x$.

## Prime $k$

Let's first consider the case of prime $k$. The explicit expression for factorial

$$n! = 1 \cdot 2 \cdot 3 \ldots (n-1) \cdot n$$

Note that every $k$-th element of the product is divisible by $k$, i.e. adds $+1$ to the answer; the number of such elements is $\Bigl\lfloor\dfrac{n}{k}\Bigr\rfloor$.

Next, every $k^2$-th element is divisible by $k^2$, i.e. adds another $+1$ to the answer (the first power of $k$ has already been counted in the previous paragraph). The number of such elements is $\Bigl\lfloor\dfrac{n}{k^2}\Bigr\rfloor$.

And so on, for every $i$ each $k^i$-th element adds another $+1$ to the answer, and there are $\Bigl\lfloor\dfrac{n}{k^i}\Bigr\rfloor$ such elements.

$$\Bigl\lfloor\dfrac{n}{k}\Bigr\rfloor + \Bigl\lfloor\dfrac{n}{k^2}\Bigr\rfloor + \ldots + \Bigl\lfloor\dfrac{n}{k^i}\Bigr\rfloor + \ldots$$

The sum is of course finite, since only approximately the first $\log_k n$ elements are not zeros. Thus, the runtime of this algorithm is $O(\log_k n)$.

### Implementation

int fact_pow (int n, int k) {
int res = 0;
while (n) {
n /= k;
res += n;
}
return res;
}

## Composite $k$

The same idea can't be applied directly. Instead we can factor $k$, representing it as $k = k_1^{p_1} \cdot \ldots \cdot k_m^{p_m}$. For each $k_i$, we find the number of times it is present in $n!$ using the algorithm described above - let's call this value $a_i$. The answer for composite $k$ will be

$$\min_ {i=1 \ldots m} \dfrac{a_i}{p_i}$$


text:
Back to Input Deck Cards

# UFD_BIOSPHERE¶

Specifies the Example Reference Biosphere (ERB) Model 1, which calculates annual dose to an individual drinking radioactive water from a well (IAEA 2003).

Under the SIMULATION block this process model is included by adding the UFD_BIOSPHERE block:

SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
UFD_BIOSPHERE <name_string>
/
SUBSURFACE_FLOW flow
MODE GENERAL
/
SUBSURFACE_TRANSPORT transport
GLOBAL_IMPLICIT
/
/
END

where <name_string> gives a user-defined name for the process model.

## Required Cards:¶

UFD_BIOSPHERE
Opens the UFD_BIOSPHERE block. Must have a matching END line.

### ERB_1A/B block¶

ERB_1A or ERB_1B <name_string>
Opens the ERB Model block with either the ERB_1A model type or the ERB_1B model type. Any number of ERB Models blocks may be specified. ERB_1A is used for water that is extracted by a sink (well) in the simulation. ERB_1B is used to approximate dose from a hypothetical well not explicitly modeled.

ERB_1A block required cards:

REGION <string>

Specifies the REGION associated with the ERB_1A model, where <string> indicates the name of that region. This region must be associated with a SOURCE_SINK which is located at the screened portion of a pumping water well.

REGION well-45b

INDIVIDUAL_CONSUMPTION_RATE <double> <unit_string>

Specifies the water consumption rate in units of volume/time that the typical person drinks out of the contaminated well.

INDIVIDUAL_CONSUMPTION_RATE 3.4 L/day

Optional Cards:

If this card is included, then the ERB_1A model will consider both supported and unsupported radionuclides in the calculation of dose. Unsupported radionuclides are those which are not explicitly modeled in the simulation (e.g. as a source term or in transport) due to short half-lives, but are considered in a dose calculation.

ERB_1B block required cards:

REGION <string>

Specifies the REGION associated with the ERB_1B model, where <string> indicates the name of that region. This region must NOT be associated with a SOURCE_SINK. Rather, the region is the location of a potential pumping water well, at the screened portion of the well.

REGION potential-well-22c

INDIVIDUAL_CONSUMPTION_RATE <double> <unit_string>

Specifies the water consumption rate in units of volume/time that the typical person drinks out of the contaminated well.

INDIVIDUAL_CONSUMPTION_RATE 2.1 L/day

DILUTION_FACTOR <double>

Specifies a unitless dilution factor that accounts for the dilution of radionuclides at the location of the potential well. The dilution factor is the user-estimated volumetric ratio of discharged well water to captured plume water. This number should be larger than or equal to 1. A value of 1 implies no dilution (e.g., discharged well water = captured plume water).

DILUTION_FACTOR 4.3d0

Optional Cards:

If this card is included, then the ERB_1B model will consider both supported and unsupported radionuclides in the calculation of dose. Unsupported radionuclides are those which are not explicitly modeled in the simulation (e.g. as a source or in transport) due to short half-lives, but are considered in a dose calculation.

Examples of the ERB_1A/B blocks:

ERB_1A A_model1
REGION well
INDIVIDUAL_CONSUMPTION_RATE 5.d0 L/day
/
ERB_1B B_model1
REGION potential-well-22c
DILUTION_FACTOR 2.d0
INDIVIDUAL_CONSUMPTION_RATE 2.d0 L/day
/
ERB_1B B_model2
REGION potential-well-67f
DILUTION_FACTOR 1.d0
INDIVIDUAL_CONSUMPTION_RATE 3.d0 L/day
/

Opens the SUPPORTED_RADIONUCLIDES block which indicates which supported radionuclides are considered in the dose calculations. The following cards or sub-blocks are required:

Specifies a sub-block for each supported radionuclide, where <string> indicates the name of the supported radionuclide. The radionuclide must be included as either a primary or secondary species in the CHEMISTRY block.

RADIONUCLIDE Tc-99
...
/

ELEMENT_KD <double>

Specifies the elemental Kd value within the material where the screen portion of the well resides. The units can be given in L-water/kg-solid or kg-water/m3-bulk, as long as the units are consistent for all supported and unsupported radionuclides.

ELEMENT_KD 3.5d5

DECAY_RATE <double> <unit_string>

Specifies the decay rate of the supported radionuclide in units of [1/time].

DECAY_RATE 1.29d-13 1/sec

INGESTION_DOSE_COEF <double>

Specifies the ingestion dose coefficient for the supported radionuclide in units of Sv/Bq.

INGESTION_DOSE_COEF 1.d-9

Examples of a full SUPPORTED_RADIONUCLIDES block:

SUPPORTED_RADIONUCLIDES
ELEMENT_KD 3.5d5  # at screened part of well
DECAY_RATE 1.4d-11 1/sec
INGESTION_DOSE_COEF 2.8d-7 # Sv/Bq
/
ELEMENT_KD 4.2d5  # at screened part of well
DECAY_RATE 1.04d-12 1/sec
INGESTION_DOSE_COEF 1.1d-7 # Sv/Bq
/
/

The concentration [Bq/L] of an unsupported radionuclide in secular equilibrium with its supported parent is the product of the sorption enhancement factor of the unsupported radionuclide, the emanation factor of the unsupported radionuclide, and the concentration [Bq/L] of the supported parent.

Opens the UNSUPPORTED_RADIONUCLIDES block which indicates which unsupported radionuclides are considered in the dose calculations. These are only included in the calculation for the ERB models that contain the card INCLUDE_UNSUPPORTED_RADS. The following cards or sub-blocks are required:

Specifies a sub-block for each unsupported radionuclide, where <string> indicates the name of the unsupported radionuclide. The radionuclide should NOT be included as either a primary or secondary species in the CHEMISTRY block because it is unsupported.

RADIONUCLIDE Bi-210
...
/

ELEMENT_KD <double>

Specifies the elemental Kd value within the material where the screen portion of the well resides. The units can be given in L-water/kg-solid or kg-water/m3-bulk, as long as the units are consistent for all supported and unsupported radionuclides. Kd is needed to calculate the sorption enhancement factor for each unsupported radionuclide. The sorption enhancement factor is the ratio of the retardation factor of the supported parent to that of the unsupported radionuclide.

ELEMENT_KD 6.5d5

DECAY_RATE <double> <unit_string>

Specifies the decay rate of the unsupported radionuclide in units of [1/time].

DECAY_RATE 8.9d-11 1/sec

INGESTION_DOSE_COEF <double>

Specifies the ingestion dose coefficient for the unsupported radionuclide in units of Sv/Bq.

INGESTION_DOSE_COEF 1.d-7

EMANATION_FACTOR <double>

Specifies the optional emanation factor in the dose calculation for an unsupported radionuclide. The default value is 1.0 if this keyword is omitted.

EMANATION_FACTOR 0.4

SUPPORTED_PARENT <string>

Indicates the name of the unsupported radionuclide’s supported parent radionuclide. The supported parent must be included as a supported radionuclide in the SUPPORTED_RADIONUCLIDES block.

SUPPORTED_PARENT Pb-210

Example of a full UNSUPPORTED_RADIONUCLIDES block:

UNSUPPORTED_RADIONUCLIDES
ELEMENT_KD 3.5d5  # at screened part of well
DECAY_RATE 1.d-14 1/sec
SUPPORTED_PARENT tracerA
INGESTION_DOSE_COEF 1.d-10 # Sv/Bq
EMANATION_FACTOR 0.4d0
/
ELEMENT_KD 4.2d5  # at screened part of well
DECAY_RATE 1.d-10 1/sec
SUPPORTED_PARENT tracerA
INGESTION_DOSE_COEF 1.d-9 # Sv/Bq
/
ELEMENT_KD 2.2d8  # at screened part of well
DECAY_RATE 1.d-9 1/sec
SUPPORTED_PARENT tracerB
INGESTION_DOSE_COEF 1.d-8 # Sv/Bq
EMANATION_FACTOR 0.4d0
/
/

### Output Control¶

At each time step, the process model will add dose calculations to an output file named *.bio. The output includes the total dose from all supported and unsupported* radionuclides, the annual dose from each radionuclide, as well as the annual dose from each supported radionuclide and its unsupported descendents. The average concentration of each radionuclide in the region associated with each ERB model is also reported.

OUTPUT_START_TIME <double> <unit_string>

If this card is included, the dose calculations will not be printed to the output.bio file until after this point in the simulation time. The default start time is 0.d0 seconds (e.g. starting at the first time step).

OUTPUT_START_TIME 5.d0 yr

## Examples:¶

UFD_BIOSPHERE

ERB_1A A_model1
REGION well
INDIVIDUAL_CONSUMPTION_RATE 5.d0 L/day
/
ERB_1B B_model1
REGION ERB1B_tester1
DILUTION_FACTOR 1.d0
INDIVIDUAL_CONSUMPTION_RATE 2.d0 L/day
/
ERB_1B B_model2
REGION ERB1B_tester2
DILUTION_FACTOR 5.d0
INDIVIDUAL_CONSUMPTION_RATE 3.d0 L/day
/

ELEMENT_KD 3.5d5  # screened part of well
DECAY_RATE 1.29d-15 1/sec
INGESTION_DOSE_COEF 1.d-10 # Sv/Bq
/
ELEMENT_KD 4.2d5  # screened part of well
DECAY_RATE 1.29d-15 1/sec
INGESTION_DOSE_COEF 1.d-9 # Sv/Bq
/
/

ELEMENT_KD 3.5d5  # at screened part of well
DECAY_RATE 1.d-14 1/sec
SUPPORTED_PARENT tracerA
INGESTION_DOSE_COEF 1.d-10 # Sv/Bq
EMANATION_FACTOR 0.4d0
/
ELEMENT_KD 4.2d5  # at screened part of well
DECAY_RATE 1.d-10 1/sec
SUPPORTED_PARENT tracerA
INGESTION_DOSE_COEF 1.d-9 # Sv/Bq
/


text:
# Toss two coins until two Heads and two Tails come up

You play a game where you toss two fair coins in the air. You always win \$1. However, if you have tossed 2 heads at least once, and 2 tails at least once, you surrender all winnings, and cannot play again. You may stop playing at anytime. What’s your strategy?

My thoughts were that this seems similar to the coupon collector. We have two bad events (2H and 2T). So after the occurence of the first bad event the second one will occur in an expected number of 4 turns. So my strategy is to stop after 3 tosses from the moment the first bad event occured. However, I can't prove it.

After the first bad event one has $$x\geq1$$ dollars in the pocket. The question is whether one shall go on playing. Denote by $$E(x)$$ the expected total win under the best strategy in this situation. If we decide to quit we have won $$x$$, and if we make a next turn we have (with probability $${3\over4}$$) won $$x+1$$ and the possibility of further play, or (with probability $${1\over4}$$) we have won $$0$$. This shows that $$E(x)=\max\left\{x, \ {3\over4}E(x+1)\right\}\ .\tag{1}$$ I'm assuming (without proof) that there will be an $$n$$ where we definitely won't play once more; hence $$E(n)=n$$. From this and $$(1)$$ we get $$E(n-1)=\max\left\{n-1,\ {3\over4}n\right\}=n-1\qquad(n\geq4)\ .$$ By downwards induction it follows that $$E(n)=n\qquad(n\geq 3)\ .$$ Using $$(1)$$ one then computes $$E(2)=\max\left\{2,\ {3\over4}E(3)\right\}={9\over4}\ ,\qquad E(1)=\max\left\{1,\ {3\over4}E(2)\right\}={27\over16}\ .$$ This means that one should play on when $$1\leq x\leq2$$, and quit otherwise.

Suppose you have won $$\x$$ by the first bad event, and plan to continue for $$y$$ more rounds. Then your expected winnings are $$\left(\frac34\right)^y(x+y)$$
For each value of $$x$$, write down the winnings for the various $$y$$, and pick the best $$y$$ for each $$x$$.

Each turn after the first bad event, you win another dollar with probability $$\frac{3}{4}$$ and loss everything with probability $$\frac{1}{4}$$.

Thus, assuming the first bad event happened at turn $$m$$, then at the $$n$$th turn after the bad event you have a $$(\frac{3}{4})^n$$ probability of earning $$n$$ dollars beyond the initial $$m$$, and $$1 - (\frac{3}{4})^n$$ probability of earning nothing. So in general you have an expected win of $$(m + n)(\frac{3}{4})^n$$ dollars. For $$n = 0$$ this is $$m$$, for $$n = 1$$ this is $$(m + 1)\frac{3}{4}$$, and for $$n = 2$$ this is $$(m + 2)\frac{9}{16}$$, for $$n = 3$$ this is $$(m + 3)\frac{27}{64}$$.

$$m \leq (m + 1)\frac{3}{4} \iff m \leq 3$$

$$(m + 1)\frac{3}{4} \leq (m + 2)\frac{9}{16} \iff 4m + 4 \leq 3m + 6 \iff m \leq 2$$

And in general,

$$(m + n)(\frac{3}{4})^n \leq (m + n + 1)(\frac{3}{4})^{n + 1} \iff 4m + 4n \leq 3m + 3n + 3 \iff m \leq 3 - n \iff n \leq 3- m$$

You can see that the optimal $$n$$ is then $$\max(0, 3 - m)$$

If while playing the game, the latest 3 tosses are HHT or TTH, then in the next move, you have a $$50\%$$ chance of losing. If you take the chance and manage to stay in the game, you again have a $$50\%$$ chance of losing. However the total probability of the game going in your favour will be $$25\%$$ by the second move itself. So, it somewhat depends on the amount of money you have already won. If the game goes something like HTHTHHTor HTHHHHHT and you have already won $$7-8\$$, then you should not take the $$50\%$$ chance. However, if you have won only $$3-4\$$, then you may take the chance as you have little to lose, but after that, it is best to stop playing as a $$25\%$$ chance clearly doesn't seem optimal.


text:
The semidirect product $(C_7\times C_{13})\rtimes C_3$

Following this page, in the classification of groups of order $273$, the product of the Sylow group $S:=C_7C_{13}\simeq C_{7}\times C_{13}$ is normal, hence can be acted on by the Sylow $C_3$ to produce semidirect products.

Since $\operatorname{Aut}(S)\simeq C_6\times C_{12}$, one only needs to find an order 3 element there to produce the homomorphism $C_3\rightarrow \operatorname{Aut}(S)$.

In particular, if $C_6=\langle x\rangle$ and $C_{12}=\langle y\rangle$, up to choice of generators of $C_3$, the group is classified by the 5 actions corresponding to: $(1,1), (x^2,1), (1,y^4), (x^2,y^4), (x^{-2},y^4)$.

It is easy to see that the center of the corresponding groups has order $273, 13, 7, 1, 1$ respectively. So my question is: how do we show that the last two groups are actually non-isomorphic?

P.S. In a similar example with $(C_7\times C_7)\rtimes C_3$ on the same page, the last two groups can be distinguished since every subgroup of order $7$ can be shown to be normal in one but not the other.

Edit: I didn't realize you were talking about $x$ and $y$ being in the automorphism group, and not the group itself. I've edited below to correct this.

We can take $\alpha$ as a generator of $C_7$, and then $x(\alpha)=\alpha^5$. So $x^2(\alpha) = \alpha^4$ and $x^{-2}(\alpha) = \alpha^2$.

We can take $\beta$ as a generator of $C_{13}$, and then $y(\beta)=\beta^2$, and so $y^4(\beta) = \beta^3$.

Consider your groups as $G=C_{91}\rtimes C_3$ and $H=C_{91}\rtimes C_3$. We can write $z$ as the generator of $C_{91}$, and $g$ and $h$ as the generator of the $C_3$ subgroups.

Then in $G$, we have $g^{-1}zg=z^{-10}$, since $-10$ is $4\pmod{7}$ and $3\pmod{13}$.

In $H$, we have $h^{-1}zh=z^{16}$, since $16$ is $2\pmod{7}$ and $3\pmod{13}$.

If $f:G\rightarrow H$ was an isomorphism, then we would have $f(z)=z^k$ for some integer $k$. We would also have $f(g)$ is some conjugate of $h$ or $h^2$, which thus acts the same as $h$ or $h^2$ on $z$.

So then

\begin{align*} z^{-10k} &= f(z^{-10})\\ &= f(g^{-1}zg)\\ &= f(g^{-1})z^kf(g)\\ &= z^{16k}\text{ or }z^{74k} \end{align*}

This implies $z^{26k}=1$ or $z^{84k}=1$, so that $7$ or $13$ divides $k$ respectively, which means $z^k$ is not a generator of $C_{91}$ in $H$. This is the contradiction.

This same idea can be used whenever you have a semidirect product $A\rtimes B$, and $A$ is abelian and $B$ is cyclic. It can be used to differentiate different semidirect products coming from non-conjugate images of $B$ in $\textrm{Aut}(A)$.

• This is satisfying~ By the way I think $f(g)$ is some conjugate of $h$ or $h^2$? – Kamineko Sep 12 '18 at 15:24
• Also $z^{26k}=1$ already says $z^k$ has order 26 hence not a generator. – Kamineko Sep 12 '18 at 15:24
• @Kamineko: yes, I had blindly thought that $h$ and $h^2$ would be conjugate in $H$, but they are not, I'll edit. – Steve D Sep 12 '18 at 15:33