In floatofmath/blindConfidence: Simulation of estimation after blinded interim analysis

knitr::opts_chunk$set(echo=FALSE,warning=FALSE,results='hide',fig.width=5,fig.height=5,message=FALSE,fig.align="center")

## For running the simulations and plotting the results the following
## packages have to be loaded:
library(parallel)
## devtools::install_github('floatofmath/blindConfidence')
library(blindConfidence)
library(ggplot2)
library(reshape2)
library(plyr)
library(pander)
options(mc.cores=min(10,detectCores()-1))

\newcommand{\nmax}{n_{2{\max}}}

In this document we present the numerical results shown in the article ``Estimation after blinded sample size reassessment'' as well as additional results that - for the sake of brevity - have not been included in the main article. Written in Rmarkdown this document serves both as supplementary material to the article but also as a program that permits a complete reproduction of all numerical results and graphical presentations of the results. This document is distributed as a vignette to the r-package blindConfidence. Therefore, all results shown in the paper as well as in this document can be replicated by compiling the corresponding source file in the package. Simply compiling the document - as happens on installation of the package - will recompute all figures using stored simulation results. To rerun the simulations the corresponding R chunks can be reconfigured by setting \texttt{eval=TRUE} so that they will be evaluated on compilation. Note that this may take a very long time and require a considerable amount of main memory. Computation of the publication results took over two weeks on a simulation server with $40$ processor cores and 64 GB of shared memory. Note that the code has been optimized to run on Linux, so while the package should work fine on Windows in theory, it has not been tested.

In Section 1 we outline the setup and reproduce the results of the simulation study presented in Section 4 of the article. Additionally in Section 2 of this document show simulation results concerning the inflation of the non-coverage for the adjusted sample size rule, and the (squared) standard error of the mean estimate. In Section 3 we investigate the maximum absolute bias of the variance and mean estimates for given first stage sample sizes, as well as the maximum of the Type I error of the $t$-test following blinded sample size reassesment. In Section 4 we reproduce the simulations shown in the case study of our paper. In Section 5 we present a simulation study investigating the bias and coverage probabilities of estimaten and confidence intervals following blinded sample size reassessment when the second stage sample size is restricted to not exceed twice the preplanned second stage sample size (which in our example is equal to the first-stage sample size).

In this document we present 22 Figures (15 of which are not shown in the Article) which take up considerable space. To keep the text readable, we have pushed all figures to the end of document. The following table gives an index of which Figures in the original paper correspond to which in this document. Note that Figure 4 in the article combines three figures from this document into one panel.

figure_table <- data.frame("Article Figure"=paste('Figure',c(1,2,3,4,5)),
                           "Supplementary Figure"=paste('Figure',c(3,4,2,"7, 13, 14",16)))

pander(figure_table)

1. Coverage probabilities and bias of the mean and variance estimates

In a first simulation study we investigated the coverage probabilities of conventional t-test confidence intervals following blinded interim analysis, as well as the bias of the estimates of the mean and variance. The second stage sample size was computed based on the blinded first stage data using either the adjusted

$$ n^a_2(S^2_{1,OS})=\min\left{\nmax,\max\left{n_{2\min},2(z_{1-\alpha} + z_{1-\beta})^2\left( \frac{S^2_{1,OS}}{\delta_0^2}-\frac{n_1}{4n_1-2}\right)-n_1+1 \right}\right}.\label{eq:sssa} $$

or the unadjusted

$$ n^u_2(S^2_{1,OS})=\min\left{\nmax,\max\left{n_{2\min},2(z_{1-\alpha} + z_{1-\beta})^2 \frac{S^2_{1,OS}}{\delta_0^2}-n_1+1 \right}\right},\label{eq:ssss} $$

reassessment rules.

To optimize the computational efficiency of the simulation algorithm, the joint distribution of the stage-wise t-statistics was simulated by sampling the stage-wise means and sample variances from independent normal and chi-squared distributions (instead of sampling individual observations). After submission of our manuscript we became aware of a paper by \cite{lu2016distribution} on sample size reassessment for non-inferiority trials, where the same approach is proposed.

We assume that the preplanned sample size of the trial was planned to reach a target power of 80% (using the normal approximation) assuming an effect size $\delta_0=1$, standard deviations $\sigma_0$ between $0.5$ and 2 in steps of $.5$ and a one-sided significance level $\alpha=0.025$. The first stage sample size $n_1$ was set to halve of the preplanned sample size. For this simulation the maximum and minimum second stage sample sizes where not restricted (i.e. $\nmax=\infty$, $n_{2\min}=0$)The simulations were performed in R with $5*10^7$ simulation runs per scenario. The data used in the final manuscript version, is distributed with this package and can be loaded with the command data(gridsim).

Difference between actual and nominal coverage probabilities

In Figure \ref{fig:coverage} we plot how far the coverage probabilities of one and two-sided confidence intervals deviate from their nominal coverage if sample sizes are reassessed using the adjusted variance estimator. In Figure \ref{fig:uc.coverage} we show the difference between actual and nominal coverage probabilities, for designs that use the unadjusted lumped variance at interim for sample size reassessment (i.e. $n_2^u$).

Bias of the effect size estimate

Next we look at estimates of the mean difference between treatment groups. In Figure \ref{fig:bias} solid lines show the bias in the mean estimate for designs where sample sizes are adjusted based on the adjusted interim variance estimate ($n^a_2$), dashed lines for designs where sample sizes are adjusted based on the unadjusted interim variance estimate ($n^u_2$). The dotted lines show the bounds for the bias of the mean estimate that can be attained under any (blinded) samplesize reassessment rule.

Bias of the variance estimate

In Figure \ref{fig:sd} we show the bias when estimating the variance at the end of the adapted trial. Solid lines show the bias of the variance estimate for designs where sample sizes are adjusted based on the adjusted interim variance estimate ($n^a_2$), dashed lines for designs where sample sizes are adjusted based on the unadjusted interim variance estimate ($n^u_2$). The dotted lines show the bias of the mean estimate using the samplesize reassessment that maximizes the negative or positive bias, respectively.

Estimates of the variance of the mean estimate (squared standard error)

To compare the actual variance of the mean to the estimated variance, we computed for each simulated trial the variance estimate $\hat S^2_e= \frac{2S^2}{n_1+n_2}$ as well as the actual variance $\sigma^2_{f}=\frac{2\sigma^2}{n_1+n_2}$ of a design with fixed sample sizes $n_1, n_2$ (thus, ignoring that $n_2$ is dependent on the first stage sample). Both quantities were then averaged over the Monte-Carlo samples. Finally, we compute the variance of the mean estimate across all Monte Carlo samples, which is an unbiased estimate of the true variance of the mean ($\sigma^2_e$). The results show that the true variance of the mean estimate of an adaptive design with blinded sample size reassessment is smaller than the average variance of a fixed sample design with the same sample sizes. The bias of the estimate of the variance of the mean estimate is close to zero around the null hypothesis. This holds both for designs using the adjusted and unadjusted interim variance estimates to reassess the sample size and gives some intuition why there is no relevant inflation of the coverage probability of confidence intervals under the null hypotheses, even though the variance estimate has a considerable negative bias. We show $\hat S^2_e$, $\sigma^2_{e}$ and $\sigma^2_{f}$ for trials using the adjusted (Figure \ref{fig:varmeanBias}) and the adjusted sample size reassessment rules (Figure \ref{fig:uc.varmeanBias}).

2. Maximum bias of the mean estimate and worst case inflation of the non-coverage

We also maximize the different absolute bias of mean and variance estimates as well as the (absolute) difference between actual and nominal coverage probabilities over $\delta$ and $\sigma$, for fixed first stage sample sizes. As the different quantities can only be evaluated with some simulation error, we implemented a three-stage optimization procedure.

1. We evaluate bias, and coverage probabilities using an initial simulation, for all per group first stage sample sizes between 2 and 50 over grid of parameter values $\delta \in [0,4]$ and $\sigma \in [.5,4]$ in steps of $.1$.
2. For each first stage sample size and quantity (mean bias, variance bias, actual - nominal coverage probabily) we identify the maximizing parameter settings based on the results of the simulation study in step one. We add additional candidate scenarios in the vicinity of each of the identified parameter settings (in a resolution that is higher then the grid in step 1) and perform another set of simulations (using an increased number of simulation runs).
3. Again, for each first stage sample size and quantity (mean bias, variance bias, actual - nominal coverage probabily) we identify the maximizing parameter settings based on the results of the simulation study in step two. Finally we run a last set of simulations using the identified scenarios to avoid any selection bias.

In the Figures \ref{fig:maxmean}-\ref{fig:maxnull} we present the results from the final simulation run in Step 3. For the bias of the variance we performed some additional simulations as our results suggested that the optimizing scenarios are given for $\delta=0$ and for the undadjusted sample size rule also $\sigma=\infty$.

Step 1: Simulation over a grid

We simulate adaptive trials for true differences of means between 0 and 4, true standard deviations between .5 and 4 (both in steps of $.1$) and per group first stage sample sizes between 2 and 50. For each scenario we performed $10^6$ simulation runs. Steps 2 and 3 of the maximizing procedure were then performed for each quantity seperately.

Step 2 and 3: Maximum bias of the mean for given $n_1$

To improve our estimate of the maximum bias of the mean estimate we refine the grid around the optima found in the simulation performed in Step 1 runs by adding points every .01 in a cubic vicinity - with a width of .1 - of the maximizing $\delta$ and $\sigma$ and simulate again using $4*10^6$ runs. We select the parameter values leading to the maximal bias and resimulate using $10^7$ runs to remove any selection bias. The results for the maximum bias of the mean estimate are shown in Figure \ref{fig:maxmean}

Step 2 and 3: Maximum bias of the variance estimate

To improve our estimate of the maximum bias of the variance estimate we refine the grid around the optima found in the simulation performed in Step 1 runs by adding points every .01 in a cubic vicinity - with a width of .1 - of the maximizing $\delta$ and $\sigma$ and simulate again using $4*10^6$ runs. We select the parameter values leading to the maximal bias and resimulate using $10^7$ runs to remove any selection bias. The results of this simulation can be seen in Figure \ref{fig:maxvar}.

It appears that the absolute bias of the variance estimate (using both sample size rules) is maximised for $\delta=0$. Above figure shows that in the simulation study the maximum is consistenly attained for $\delta$ close to zero. Also for the unadjusted reassessment rule, values of $\sigma$ that maximise the absolute variance bias are consistently found at the upper border of the simulation grid. It appears that the maximum is obtained for large values of $\sigma$. Consequently we additionally ran the following simulations:

1. Maximize the absolute variance bias using a fixed value of $\delta=0$, and values of $\sigma$ ranging from $0.5$ to $4$
2. Estimate the absolute variance using simulations setting $\delta=0$ and $\sigma=4$ (i.e. at the border of the chosen parameter grid)
3. Estimate the absolute variance using simulations setting $\delta=0$ and $\sigma=50$ (i.e. approaching infinity)

Comparison of simulated optima for the unadjusted reassessment rule

In Figure \ref{fig:unadjusted0} we compare the maximum of the absolute variance bias among scenarios where $\delta = 0$ is fixed to the maximum absolute bias among scenarios where $\delta$ was allowed to range from $0$ to $4$. It seems that the difference in the bias curves is due only to Monte Carlo error. We therefore conclude that the maximum is attained at $\delta=0$.

In Figure \ref{fig:sigma50} we compare the maximum of the absolute variance bias over scenarios where $\sigma$ may range betwenn $.5$ and $4$, where $\sigma=4$ is fixed, and where $\sigma$ is set to $50$. We observe that - especially for smaller sample sizes - the absolute bias for scenarios with $\sigma=50$ exceeds the maximum found for scenarios with $\sigma\leq 4$. The difference is consistent across different first stage sample sizes. Also, fixing $\sigma=4$ does not result in a decrease of the absolute bias compared to settings where $\sigma$ may range between $.5$ and $4$ - but rather in a small but consistent increase\footnote{This may seem as surprising as the $\sigma = 4$ is included in the scenario where $\sigma$ is not fixed. However, if the local optimum is attained at $\sigma=4$ optimizing over a range of values is likely to select a different (suboptimal) value due to Monte Carlo error, thus resulting in a smaller estimate of the maximum (absolute) bias.}. In consequence we conclude that the maximum absolute bias of the variance (using the unadjusted sample size reassessment rule) is increasing in $\sigma$. In the remainder we will use $\sigma=4$, which even though not optimal (in terms of maximising the bias) compared to $\sigma=50$ (or even larger values) reflects a value that may be of practical relevance.

Comparison of simulated optima for the adjusted reassessment rule

Similar to the unadjusted sample size reassessment rule, the results shown in Figure \ref{fig:adjusted0} suggest that the difference between the maximum absolute variance bias found for scenarios where $\delta = 0$ fixed, and scenarios where $\delta$ was allowed to range from $0$ to $4$ is due only to Monte Carlo error. We therefore conclude that the maximum is attained at $\delta=0$.

In Figure \ref{fig:adjusted50}, we show that, for the adjusted sample size rule, setting $\delta = 0$ and $\sigma = 50$ decreases the absolute variance bias compared to the maximum values found when $\delta=0$ and $\sigma$ ranged from $.5$ to $4$; and also, that setting $\sigma=4$ does decrease the bias compared to the maximum values found when $\delta=0$ and $\sigma$ ranged from $.5$ to $4$, but not as much as for $\sigma=50$. The differences are consistent across different first stage sample sizes and appear larger than the Monte Carlo Error. We therefore conclude that the the (absolute) bias is maximized for values of $\sigma$ lying within the parameter grid.

Final results included in the paper

Based on above results we concluded to present the results that show the maximum absolute variance bias for $\delta=0$ and $\sigma$ ranging from $.5$ to $4$ - when using the adjusted sample size rule; and $\delta = 0$ and $\sigma=4$ - when using the unadjusted sample size rule. Even though not a global maximum we opted for the latter scenario as it reflects the maximum absolute variance bias within a range of parameters that seem of practical concern. The corresponding results are shown in Figure \ref{fig:maxvar}.

Step 2 and 3: Maximum difference between nominal and actual coverage probability

To improve our estimate of the absolute difference between nominal and actual coverage probability we refine the grid around the optima found in the simulation performed in Step 1 runs by adding points every .005 in a cubic vicinity - with a width of .4 - of the maximizing $\delta$ and $\sigma$ and simulate again using $4*10^6$ runs. We select the parameter values leading to the maximal bias and resimulate using $10^7$ runs to remove any selection bias. The results of this simulation can be seen in Figure \ref{fig:maxcov}.

Maximimizing the inflation of the non-coverage probability under the Null

We also performed a simulation based optimization of the inflation of the non-coverage probability under the null hypotheses (i.e. $\delta=0$). This corresponds to the inflation of the Type I error rate of the associated test procedure. The results show that there is minor inflation of the Type I error rate for first stage per group sample sizes below 10 (.5 to .01 percentage points). For larger sample sizes the inflation is less than the simulation error. See Figure \ref{fig:maxnull} for the corresponding results.

3. Case Study

alpha = 0.025
beta = 0.2
d = 5.5
s = 8

n1 <- 15
s1o <- 6
s1a <- 5.3

set.seed(659376)
G <- expand.grid(delta = round(seq(-2*d, 2*d,d/20),2),
                 sigma = round(seq(1,20,1)),
                 d=5.5,
                 s=8,n1=15)

runs=10^7

casesim <- simulate_batch(G,runs)
runs <- 10^5
bsim <- rbindlist(lapply(1:nrow(G),
                             function(i) {c(G[i,],
                                            simMBIA(delta=G[i,]$delta,
                                                    sigma=G[i,]$sigma,
                                                    n1=G[i,]$n1,
                                                    runs=runs))}))

#load('fallback.rda')
#bsim <- as.data.frame(apply(bsim,2,unlist),row.names=NA)

pN <- 2*ceiling(zss(s,d,alpha,beta))

gsim <- casesim

gsim[,mean.upper.bias := bsim$m.bias]
gsim[,mean.lower.bias := bsim$m.bias.n]
casesim <- gsim

print(fname <- paste('casesim_',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep=''))
save(casesim,file=fname)

As an example @Kieser03 present the trial reported in @Malsch01, which is a randomized placebo controlled trial of the efficacy of the kava-kava special extract "WS$^{\textregistered}$ 1490" for the treatment of anxiety. The primary endpoint of the study was the change in the Hamilton Anxiety Scale (HAMA) between baseline and end of treatment. Assuming a mean difference between treatment and control of $\delta_0 = r d$ and standard deviation of $\sigma_0 = 8$ (i.e. a variance of $\sigma_0^2=64$) results in a sample size of $34$ patients per group to provide a power of $80\%$.

variance.bias <- lowerBound(n1,d,alpha,beta)

data('casesim')


max.mean.bias <- which.max(casesim$mean.bias)
uc.max.mean.bias <- which.max(casesim$uc.mean.bias)
min.mean.bias <- which.min(casesim$mean.bias)
uc.min.mean.bias <- which.min(casesim$uc.mean.bias)

max.var.bias <- which.min(casesim$variance.bias)
uc.max.var.bias <- which.min(casesim$uc.variance.bias)

max.nc.upper <- which.max(casesim$upper.prob)
max.nc.lower <- which.max(casesim$lower.prob)
max.nc.total <- which.max(casesim$total.prob)
uc.max.nc.upper <- which.max(casesim$uc.upper.prob)
uc.max.nc.lower <- which.max(casesim$uc.lower.prob)
uc.max.nc.total <- which.max(casesim$uc.total.prob)

@Kieser03 assume that sample size reassessment based on blinded data is performed after $15$ patients per group and consider that the interim estimate of the standard deviation is $S_{1,OS} = 6$ and adjusting for $\delta_0 = 5.5$ gives $S_{1,OS,\delta_0}=5.3$. The corresponding updated sample sizes are then given by $n^u_2 = 4.7$ and $n^a_2 = 0.6$ patients (per group), respectively. Consequently, using the unadjusted sample size rule 5 patients per group would have been recruited in the second stage; using the adjusted rule only 1 patient per group would have been recruited in the second stage.

If the sample size reassessment is based on the unadjusted interim variance estimate $S_{1,OS}$, the variance estimate from the completed trial will be negatively biased. Theorem 4 gives a lower bound for the magnitude of that bias which in this case is $r round(variance.bias,2)$. We also performed a simulation study, similar to that of Section 4 in the paper, where we fix $\delta_0=5.5$, $\sigma_0=8$ and $n_1=15$ - but not $S_{1,OS}$, $S_{1,OS,\delta_0}$ - and let $\delta$ vary from $-2\delta_0=-11$ to $2\delta_0=11$ in steps of $0.05$ and $\sigma$ from $1$ to $20$ in steps of $1$. Detailed results are shown in Figure \ref{fig:case} where we plot the coverage of the confidence intervals and biases of the mean, and the variance estimate. The variance bias gets as low as $r round(casesim[uc.max.var.bias,uc.variance.bias],2)$ if the true mean effect is $\delta = r round(casesim[uc.max.var.bias,delta],2)$ and the true standard deviation $\sigma = r round(casesim[uc.max.var.bias,sigma],2)$ (i.e. it reaches the theoretical boundary within simulation error). For smaller values of the true standard deviation the bias is smaller in magnitude. The mean bias goes up to $r round(casesim[uc.max.mean.bias,uc.mean.bias],2)$ for negative values of the true effect size and as low as $r round(casesim[uc.min.mean.bias,uc.mean.bias],2)$ for positive values of the true effect size. The absolute bias reaches its maximum for a true mean differences of $\delta = \pm r round(casesim[uc.min.mean.bias,delta],2)$) and standard deviation of $\sigma = r round(casesim[uc.min.mean.bias,sigma])$). The maximum inflation of the coverage probabilities of one-sided $97.5\%$ and two-sided $95\%$ confidence intervals is $r 100*round(casesim[uc.max.nc.upper,uc.upper.prob]-.025,3)$ and $r 100*round(casesim[uc.max.nc.total,uc.total.prob]-.05,3)$ percentage points, respectively. The actual coverage probabilities are smallest for large absolute values of the true mean difference and standard deviations of around $5$.

If the sample size reassessment is based on the adjusted variance estimate, the absolute bias of the variance and mean estimate will be even larger taking values up to $r abs(round(casesim[max.var.bias,variance.bias],2))$ for the variance and up to $r round(casesim[max.mean.bias,mean.bias],2)$ for the mean, respectively. The inflation of the coverage probabilities goes up to $r 100*round(casesim[max.nc.upper,upper.prob]-.025,3)$ percentage points for the one-sided confidence intervals and $r 100*round(casesim[max.nc.total,total.prob]-.05,3)$ percentage points for the two-sided intervals.

4. Restricted sample size rules

Additionally we investigate the effect of blinded sample size reassessment on the coverage probabilities of conventional t-test confidence intervals, as well as the bias of the estimates of the mean and variance, when the second stage sample size is limited by a predetermined upper bound. Specifically we compute the second stage sample size using the rules $n^a_2$ and $n^u_2$ as defined above, setting the maximum second stage sample size to $\nmax = 2n_1$, but leaving the minimum sample size at zero (i.e. $n_{2\min}= 0$). Because the (blinded) interim analysis is performed after one half of the preplanned total (per-group) sample size is observed, setting $\nmax = 2n_1$ is equivalent to limiting the second stage sample size to be at most twice its preplanned value.

As in Section 1 we assume that the preplanned sample size of the trial was planned to reach a target power of 80% (using the normal approximation) assuming an effect size $\delta_0=1$, standard deviations $\sigma_0$ between $0.5$ and 2 in steps of $.5$ and a one-sided significance level $\alpha=0.025$. The upper bound for the second stage per group sample size $\nmax$ was set to twice the preplanned second stage sample size (i.e. $\nmax= 2n_1$). The simulations were performed in R with $10^7$ simulation runs per scenario. The data used in the final manuscript version, is distributed with this package and can be loaded with the command data(restsim).

Overall we observe that for $\sigma_0 < \sigma$ the bias of effect and variance estimates (Figures \ref{fig:rest.bias} and \ref{fig:rest.sd}, respectively) as well as the difference (in absolute terms) between nominal and actual coverage probabilities (Figures \ref{fig:rest.coverage} and \ref{fig:rest.uc.coverage}) is substantially reduced. This is not surprising as for $\sigma_0 << \sigma$ the adjusted second stage sample size is essentially fixed at $\nmax$ resulting in unbiased estimates and correct coverage probabilities. Note, however, that the corresponding trials will not control the Type II error at the desired level anymore. This can also be seen as the variance of the mean estimate (Figures \ref{fig:rest.varmeanBias} and \ref{fig:rest.uc.varmeanBias}) is substantially larger in these cases and does not decrease for increasing absolute effect sizes, compared to the unrestricted samplesize reassessment rules.

For $\sigma_0 \geq \sigma$ we observe that the shape of the bias of estimates and the difference between nominal and actual coverage probabilities as functions of the true effect size is similar to when unrestricted sample size rules are applied. While the (absolute) bias of the effect estimate is reduced by about a factor $\nmax/(n_1+\nmax) = 2/3$ (following from Theorem 2 in the article) the bias of the variance estimate is comparable in magnitude to when unrestricted sample size reassessment rules are applied. The effect on the difference between nominal and actual coverage probabilities is not as clear cut. For $\sigma$ close to $\sigma_0$ there appears to be a reduction (to about 80%) compared to designs without sample size restrictions, for $\sigma < \sigma_0$ the magnitude of the differences is similar.

\clearpage

Figures

#This little function takes care of annotating the plots.
lverbose  <- function(labels,from=c("anc","t1","mmb","mvb"),to=c("actual - nominal coverage [%]","Type I error inflation","maximum mean bias","absolute bias"),multi_line=TRUE){
    variable <- names(labels)
    my_labels <- list()
    if(variable == "s"){
        my_labels[[1]] <- sapply(labels[[1]],function(val) substitute(paste(sigma[0]," = ",foo,", ",n[1]," = ",n1,sep=''),
                                              list(foo=val,n1=ceiling(1/2*zss(val,1,.025,.2)))))
    } else {
        my_labels[[1]] <- sapply(labels[[1]],function(val) substitute(paste(sigma," = ",foo,sep=""),list(foo=val)))
    }
    return(my_labels)
}

## The next chunk generates the matrix of parameters for which
## simulations are run.
G <- expand.grid(delta=round(seq(-4,4,.1),1),
                 sigma=round(seq(.5,2,.5),2),
                 d = 1,
                 s = round(seq(.5,2,.5),2))
G$n1 <- ceiling(1/2*zss(G$s,G$d,.025,.2))
G <- subset(G,n1 < 65 & n1 > 1)

## The next chunk runs the simulations by applying the simulation
## function, in parallel, across rows of the parameter grid $G$, which
## contains the parameter settings for each scenario. In order to rerun
## them it has to be reconfigured for evaluation. Note that running the
## simulations may take quite long.
set.seed(50014056)
runs <- 5*10^7
gridsim <- t(simplify2array(mclapply(1:nrow(G),
                                     function(i) {c(G[i,],
                                                    simVBIA(n1=G[i,]$n1,
                                                            delta=G[i,]$delta,
                                                            sigma=G[i,]$sigma,
                                                            d=G[i,]$d,
                                                            s=G[i,]$s,
                                                            cf=1,
                                                            runs=runs,
                                                            alpha=.025,
                                                            beta=.2))})))


gsim <- as.data.frame(t(apply(gridsim,1,unlist)))

## compute relative bias
gsim$var.rbias <- gsim$variance.bias/gsim$sigma^2

## compute the theoretical lower bound for the variance bias
gsim$bound <- lowerBound(gsim$n1,gsim$d)

runs <- 10^5
bsim <- mclapply(1:nrow(gsim),function(i) {c(gsim[i,],
                                             simMBIA(delta=gsim[i,]$delta,
                                                     sigma=gsim[i,]$sigma,
                                                     n1=gsim[i,]$n1,
                                                     runs=runs))})

bsim <- do.call('cbind',bsim)
bsim <- as.data.frame(apply(bsim,1,unlist),row.names=NA)

gridsim <- bsim
gridsim$brannath <- .4*sqrt(2)*gridsim$sigma/sqrt(gridsim$n1)

fname <- paste('gridsim_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(gridsim,file=fname)

data(gridsim)

full.sim <- gridsim

gridsim <- subset(full.sim,d==1 & sigma>=.5 & n1 > 2)
gridsim$ylow <- 0

```r",fig.pos="hb!"}

In this and the following chunks we set up the figures using ggplot2

objects.

tupd <- ggplot(gridsim,mapping=aes(x=delta,y=100(.025-upper.prob))) + geom_line(y=0,colour='gray') + geom_path(lty=2) + geom_path(aes(y=100(.025-lower.prob)),lty=3) + geom_path(aes(y=100*(.05-total.prob))) + facet_grid(s~sigma,labeller=lverbose,scales="free_y") + theme_bw() + #ylim(-.025,.01) + xlab(expression(delta)) + ylab('Actual - nominal coverage [%]') print(tupd)

```r$: upper confidence bound (dashed line),  lower confidence bound (dotted line) and  two-sided interval (solid line). Rows refer to the a priori assumed standard deviations $\\sigma_0$ determining the first stage sample size $n_1$. The columns correspond to actual standard deviations. The x-axis in each graph denotes the true treatment effect $\\delta$, the y-axis shows the difference of actual to nominal coverage probability such that negative values indicate settings where the confidence bound is not valid.\\label{fig:uc.coverage}"}

tupd.uc <- ggplot(gridsim,mapping=aes(x=delta,y=100*(.025-uc.upper.prob))) +
    geom_line(y=0,colour='gray') +
    geom_path(lty=2) + 
    geom_path(aes(y=100*(.025-uc.lower.prob)),lty=3) + 
    geom_path(aes(y=100*(.05-uc.total.prob))) + 
    facet_grid(s~sigma,labeller=lverbose) + theme_bw()+ #ylim(-.025,.01) +
    xlab(expression(delta)) + ylab('Actual - nominal coverage [%]')
print(tupd.uc)

```rBias of the mean under blinded sample size reassessment using the unadjusted (solid line) and adjusted (dashed line) interim variance estimate. The dotted lines show maximum negative and positive bias that can be attained under any blinded sample size reassessment rule according to Theorem 2. The dashed gray line shows the maximum bias that can be attained under any (unblinded) sample size reassessment rule. The treatment effect used for planning is set to $\delta_0=1$. Rows refer to the a priori assumed standard deviations $\sigma_0$ determining the first stage sample size $n_1$. The columns correspond to actual standard deviations. The x-axis in each graph denotes the true treatment effect $\delta$, the y-axis shows the bias."}

tupd.mean.bias <- ggplot(gridsim,mapping=aes(x=delta,y=mean.bias)) + geom_path(lty=2) + geom_line(y=0,colour='gray') + geom_path(aes(y=m.bias),lty=3,col='black',data=subset(gridsim,m.bias<=0.2 & n1 == 2)) + geom_path(aes(y=m.bias.n),lty=3,col='black',data=subset(gridsim,m.bias.n>=-0.2 & n1 == 2)) + geom_path(aes(y=m.bias),lty=3,col='black',data=subset(gridsim,m.bias<=0.1 & n1 == 8)) + geom_path(aes(y=m.bias.n),lty=3,col='black',data=subset(gridsim,m.bias.n>=-0.1 & n1 == 8)) + geom_path(aes(y=m.bias),lty=3,col='black',data=subset(gridsim,m.bias<=0.1 & n1 == 18)) + geom_path(aes(y=m.bias.n),lty=3,col='black',data=subset(gridsim,m.bias.n>=-0.1 & n1 == 18)) + geom_path(aes(y=m.bias),lty=3,col='black',data=subset(gridsim,m.bias<=0.1 & n1 == 32)) + geom_path(aes(y=m.bias.n),lty=3,col='black',data=subset(gridsim,m.bias.n>=-0.1 & n1 == 32)) + geom_path(aes(y=uc.mean.bias),lty=1) + geom_path(aes(y=brannath),lty=2,col='darkgray',data=subset(gridsim,brannath <= 0.1)) + geom_path(aes(y=-brannath),lty=2,col='darkgray',data=subset(gridsim,brannath <= 0.1)) +

facet_grid(s~sigma,labeller=lverbose) + theme_bw() +
xlab(expression(delta)) + ylab('Bias of the mean') #+ ylim(-.4,.4)

print(tupd.mean.bias)

```rBias of the variance under blinded sample size reassessment using the unadjusted (solid line) and adjusted (dashed line) interim variance estimate. The gray line gives the lower bound from Theorem 4 for the bias under sample size reassessment based on the unadjusted variance estimate. The treatment effect used for planning is set to $\\delta_0=1$. Rows refer to the a priori assumed standard deviations $\\sigma_0$ determining the first stage sample size $n_1$. The columns correspond to actual standard deviations. The x-axis in each graph denotes the true treatment effect $\\delta$, the y-axis shows the bias."}

tupd.var.bias <- ggplot(gridsim,mapping=aes(x=delta,y=variance.bias)) +
#    geom_line(y=0,colour='gray',lty=2) +
        geom_path(aes(y=bound),color="darkgray") +
    geom_path(lty=2) + geom_path(aes(y=uc.variance.bias),lty=1) +
    facet_grid(s~sigma,labeller=lverbose) + theme_bw() + #ylim(-0.3,0.1) +
    xlab(expression(delta)) + ylab('Bias of the variance')
print(tupd.var.bias)

\clearpage

```r Adjusted sample size reassessment rule: actual variance of the mean estimate $\sigma^2_e$ (solid line), average of the estimated variances of the mean $S^2_e$ (dashed line), average of the variances $\sigma^2_f$ of the mean of fixed sample trials with the same sample sizes."}

tupd.vm.bias <- ggplot(gridsim,mapping=aes(x=delta,y=vm)) +geom_line(lty=1) + geom_path(aes(y=ev),lty=2) + geom_path(aes(y=exv),lty=3) + facet_grid(s~sigma,labeller=lverbose) + theme_bw() + geom_blank(aes(x=delta,y=ylow)) + xlab(expression(delta)) + ylab('Variance of the mean estimate')

print(tupd.vm.bias)

```r Variance of the effect estimates (i.e., square of the standard errors) for the unadjusted sample size reassessment rule:  actual variance $\\sigma^2_e$ of the mean  estimate (solid line), average of the estimated variances $S^2_e$ of the mean  (dashed line), average of the variance $\\sigma^2_f$ of the mean of corresponding fixed sample trials (dotted line)."}

tupd.uc.vm.bias <- ggplot(gridsim,mapping=aes(x=delta,y=uc.vm)) +
    geom_line(lty=1) + geom_path(aes(y=uc.ev),lty=2) +
    geom_path(aes(y=uc.exv),lty=3) + facet_grid(s~sigma,labeller=lverbose) +
    theme_bw() + xlab(expression(delta)) + geom_blank(aes(x=delta,y=ylow)) +  ylab('Variance of the mean estimate')

print(tupd.uc.vm.bias)

## this chunk produces the files that go into the paper
setwd('~/repos/overleaf/blind_confidence/graphs/')

cairo_ps('coverage.eps',height=5)
print(tupd)
dev.off()
cairo_ps('uc_coverage.eps',height=5)
print(tupd.uc)
dev.off()
cairo_ps('varmeanBias.eps',height=5)
print(tupd.vm.bias)
dev.off()
cairo_ps('uc_varmeanBias.eps',height=5)
print(tupd.uc.vm.bias)
dev.off()
cairo_ps('mean_bias.eps',height=5)
print(tupd.mean.bias)
dev.off()
cairo_ps('var_bias.eps',height=5)
print(tupd.var.bias)
dev.off()

G <- expand.grid(delta=round(seq(0,4,.1),2),
                 sigma=round(seq(.5,4,.1),2),
                 d = c(1),
                 n1 = 2:50)
G$s <- zsd(G$n1,G$d,.5,.025,.2)

set.seed(50014056)
options(mc.cores=min(38,detectCores()-1))
maxsim <- simulate_batch(G,10^6)

fname <- paste('maxsim_',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(maxsim,file=fname)

data(maxsim)

full.sim <- maxsim
maxsim$ylow <- 0

                                        #This little function takes care of annotating the plots.
lverbose  <- function(labels,from=c("anc","t1","mmb","mvb"),to=c("AC-NC [%]","Inflation [%]","max. mean bias","abs. variance bias"),multi_line=TRUE){
    str_labels <- label_value(labels,multi_line)
    par_labels <- label_parsed(labels,multi_line)
    strings <- grep("str_",str_labels[[1]])
    str_labels[[1]] <- sub("str_","",str_labels[[1]])
    if(!is.null(from)){
        str_labels[[1]] <- mapvalues(str_labels[[1]],from,to)
    }
    par_labels[[1]][strings] <- str_labels[[1]][strings]
    return(par_labels)
}

maxmeanopt1 <- select_results(maxsim,c('mean.bias','uc.mean.bias'),base_columns=c('n1','tn1','delta','sigma','d','s'))

G2mean <- add_epsilon(maxmeanopt1,.05,.01,c('delta','sigma'))
plot_grid(maxsim,maxmeanopt1,G2mean,'delta')
options(mc.cores=min(20,detectCores()-1))
maxmeansim2 <- simulate_batch(G2mean,4*10^6)
## reshape and resim
maxmeanopt2 <- select_results(maxmeansim2,c('mean.bias','uc.mean.bias'),base_columns=c('n1','tn1','delta','sigma','d','s'))
maxmeansim3 <- simulate_batch(maxmeanopt2,10^7)
fname <- paste('resim_mean',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(maxmeansim3,file=fname)

data(maxmeansim3)

max.bias2 <- ddply(maxmeansim3,"n1",transform,
                  max.mean.bias= mean.bias,
                  uc.max.mean.bias= uc.mean.bias,
                  method=c("mean.bias"="adjusted","uc.mean.bias"="unadjusted")[what.max]
                  )

max.bias2 <- max.bias2[,c('method','n1','delta','sigma','max.mean.bias','uc.max.mean.bias')]

```r Maximum absolute bias of the mean estimate for given per group first stage sample sizes $n_1$ between $2$ and $50$. Left column shows the results for the unadjusted sample size reassessment rule, right column for the adjusted sample size reassessment rule. The first row shows the value of the bias, the second the effect size $\delta$ and the third the standard deviation $\sigma$ at which the specific value is attained. The gray line denotes a Loess smoothed estimate.",fig.width=4,fig.height=6,fig.align='center'}

df <- melt(max.bias2,id=c('n1','method')) df <- subset(df,!(method=='adjusted' & variable == 'uc.max.mean.bias')) df <- subset(df,!(method=='unadjusted' & variable == 'max.mean.bias'))

df <- within(df,scale <- c("absolute","relative")[grepl("rel",variable)+1]) df <- within(df,variable <- sub("rel\.","",variable)) df <- within(df,variable <- sub("uc\.","",variable)) df <- df[df$variable %in% c('max.mean.bias','delta','sigma'),] df <- within(df,variable <- factor(variable)) df <- within(df,variable <- revalue(variable,c('max.mean.bias'='str_mmb'))) df <- within(df,variable <- relevel(variable,'delta')) df <- within(df,variable <- relevel(variable,'str_mmb')) df <- within(df,method <- relevel(method,'unadjusted'))

dfa <- df[df$scale=="absolute",]#"relative",]# scalehelper <- data.frame(n1=rep(-10,3),variable=factor(1:3,label=c("str_mmb",expression(delta),expression(sigma))),value=rep(0,3),scale=rep("absolute",3),method=rep("adjusted",3)) dfa <- rbind(dfa,scalehelper)

mplot <- ggplot(dfa)+ geom_line(aes(n1,-value),data=subset(dfa,variable=="str_mmb")) + geom_smooth(aes(n1,value),data=subset(dfa,variable!="str_mmb"),col='darkgray',fill='gray') + geom_point(aes(n1,value),data=subset(dfa,variable!="str_mmb"),size=.2) + facet_grid(variable~method,scale='free_y',labeller=lverbose) +xlab(expression(n[1]))+ylab('')+theme_bw()+xlim(c(2,50))

print(mplot)

relativer bias (constallation die den bias,bzw,type I error, maximiert)

```r
maxvarianceopt1 <- select_results(maxsim,c('variance.bias','uc.variance.bias'),base_columns=c('n1','tn1','delta','sigma','d','s'))


G2variance <- add_epsilon(maxvarianceopt1,.05,.01,c('delta','sigma'))
plot_grid(maxsim,maxvarianceopt1,G2variance,'sigma')
options(mc.cores=min(20,detectCores()-1))
maxvariancesim2 <- simulate_batch(G2variance,4*10^6)
## reshape and resim
maxvarianceopt2 <- select_results(maxvariancesim2,c('variance.bias','uc.variance.bias'),base_columns=c('n1','tn1','delta','sigma','d','s'))
plot_grid(G2variance,maxvarianceopt1,maxvarianceopt2,'delta')
maxvariancesim3 <- simulate_batch(maxvarianceopt2,10^7)
fname <- paste('resim_variance',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(maxvariancesim3,file=fname)

data(maxvariancesim3)

max.varbias2 <- ddply(maxvariancesim3,"n1",transform,
                  max.variance.bias= variance.bias,
                  uc.max.variance.bias= uc.variance.bias,
                      method=c("variance.bias"="adjusted","uc.variance.bias"="unadjusted")[what.max],
                  theoretical = bound
                  )

max.varbias2 <- max.varbias2[,c('method','n1','delta','sigma','max.variance.bias','uc.max.variance.bias','theoretical')]

```r Maximum absolute bias of the variance estimate for given per group first stage sample sizes $n_1$ between $2$ and $50$. Left column shows the results for the unadjusted sample size reassessment rule, right column for the adjusted sample size reassessment rule. The first row shows the value of the bias, the second the effect size $\delta$ and the third the standard deviation $\sigma$ at which the specific value is attained. The gray line denotes a Loess smoothed estimate.",fig.align='center'}

df <- melt(max.varbias2,id=c('n1','method','theoretical')) df <- subset(df,!(method=='adjusted' & variable == 'uc.max.variance.bias')) df <- subset(df,!(method=='unadjusted' & variable == 'max.variance.bias'))

df <- within(df,scale <- c("absolute","relative")[grepl("rel",variable)+1]) df <- within(df,variable <- sub("rel\.","",variable)) df <- within(df,variable <- sub("uc\.","",variable)) df <- df[df$variable %in% c('max.variance.bias','delta','sigma'),] df <- within(df,variable <- factor(variable)) df <- within(df,variable <- revalue(variable,c('max.variance.bias'='str_mvb'))) df <- within(df,variable <- relevel(variable,'delta')) df <- within(df,variable <- relevel(variable,'str_mvb')) df <- within(df,method <- relevel(method,'unadjusted'))

dfa <- df[df$scale=="absolute",]#"relative",]# scalehelper <- data.frame(n1=rep(-10,3),variable=factor(1:3,label=c("str_mvb",expression(delta),expression(sigma))),value=rep(0,3),scale=rep("absolute",3),method=rep("adjusted",3),theoretical=3) dfa <- rbind(dfa,scalehelper)

vplot <- ggplot(dfa,aes(n1,value))+ geom_line(aes(n1,-value),data=subset(dfa,variable=="str_mvb")) + geom_line(aes(n1,-theoretical),data=subset(dfa,variable=="str_mvb" & method=='unadjusted'),linetype='dotted') + geom_smooth(data=subset(dfa,variable==levels(variable)[3]),col='darkgray',fill='gray') + geom_point(data=subset(dfa,variable!="str_mvb"),size=.2) + geom_smooth(data=subset(dfa,variable==levels(variable)[2]),col='darkgray',fill='gray') + facet_grid(variable~method,scale='free_y',labeller=lverbose) +xlab(expression(n[1]))+ylab('')+theme_bw()+xlim(c(2,50)) print(vplot)

```r

maxvarianceoptR1 <- select_results(subset(maxsim,delta==0),c('variance.bias','uc.variance.bias'),base_columns=c('n1','tn1','delta','sigma','d','s'))


G2varianceR <- add_epsilon(maxvarianceoptR1,.1,.01,c('sigma'))
plot_grid(maxsim,maxvarianceoptR1,G2varianceR,'sigma')
options(mc.cores=min(20,detectCores()-1))
maxvariancesimR2 <- simulate_batch(data.frame(G2varianceR),10^7,use=T)
## reshape and resim
maxvarianceoptR2 <- select_results(maxvariancesimR2,c('variance.bias','uc.variance.bias'),base_columns=c('n1','tn1','delta','sigma','d','s'))

## simulate with delta set to zero!
## this will not be plotte if you want to plot it
## modify the code below
plot_grid(G2varianceR,maxvarianceoptR1,maxvarianceoptR2,'sigma')
maxvariancesimR3 <- simulate_batch(maxvarianceoptR2,10^7)
fname <- paste('resim_variance_d0',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(maxvariancesimR3,file=fname)

## simulate with delta = 0, and large sigma=50, we need to up the sim runs this takes up quite a bit of ram! 
## this will not be plotte if you want to plot it
## modify the code below
options(mc.cores=min(2,detectCores()-1))
maxvarianceoptRL2 <- maxvarianceoptR2
maxvarianceoptRL2[what.max=='uc.variance.bias',sigma:=50]
maxvariancesimRL3 <- simulate_batch(maxvarianceoptRL2,10^8,use=T)
fnameL <- paste('resim_variance_d0sl',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(maxvariancesimRL3,file=fnameL)

## simulate with delta = 0, and large sigma=50, we need to up the sim runs this takes up quite a bit of ram! 
## this will not be plotte if you want to plot it
## modify the code below
options(mc.cores=min(2,detectCores()-1))
maxvarianceoptR502 <- maxvarianceoptR2
maxvarianceoptR502[,sigma:=50]
maxvariancesimR503 <- simulate_batch(data.frame(maxvarianceoptR502),10^8,use=T)
fnameL <- paste('resim_variance_d0s50',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(maxvariancesimR503,file=fnameL)

## simulate with delta = 0, and large sigma=50, we need to up the sim runs this takes up quite a bit of ram! 
## this will not be plotte if you want to plot it
## modify the code below
options(mc.cores=min(2,detectCores()-1))
maxvarianceoptR42 <- maxvarianceoptR2
maxvarianceoptR42[,sigma:=4]
maxvariancesimR43 <- simulate_batch(data.frame(maxvarianceoptR42),10^8,use=T)
fnameL <- paste('resim_variance_d0s4',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(maxvariancesimR43,file=fnameL)

```rComparison of maximal absolute variance bias for different parameter settings: Solid line - $\delta=0$, $\sigma \in [.5,4]$; Dashed line - $\delta \in [0,4]$, $\sigma=[.5,4]$; transparent red line - theoretical bound."} data(maxvariancesimR3) data(maxvariancesimR43) data(maxvariancesimRL3) data(maxvariancesimR503)

ggplot(subset(maxvariancesimR3,what.max=='uc.variance.bias'),aes(n1,-uc.variance.bias))+geom_line()+ geom_line(data=subset(maxvariancesim3,what.max=='uc.variance.bias'),linetype=2) + geom_line(data=subset(maxvariancesim3,what.max=='uc.variance.bias'),aes(n1,-bound),col='red',size=3,alpha=.5)+ylab("Maximum of the absolute variance bias") + theme_bw()

```rComparison of maximal absolute variance bias for different parameter settings: Solid line - $\\delta=0$, $\\sigma=50$; Dashed line - $\\delta=0$, $\\sigma=4$; Dotted line $\\delta=0$, $\\sigma \\in [.5,4]$; transparent red line - theoretical bound."}
ggplot(subset(maxvariancesimRL3,what.max=='uc.variance.bias'),aes(n1,-uc.variance.bias))+geom_line()+
    geom_line(data=subset(maxvariancesimR43,what.max=='uc.variance.bias'),linetype=2) +
            geom_line(data=subset(maxvariancesimR3,what.max=='uc.variance.bias'),linetype=3) +
            geom_line(data=subset(maxvariancesim3,what.max=='uc.variance.bias'),aes(n1,-bound),col='red',size=3,alpha=.5)+ylab("Maximum of the absolute variance bias") + theme_bw()

```rComparison of maximal absolute variance bias for different parameter settings: Solid line - $\delta=0$, $\sigma \in [.5,4]$; Dashed line - $\delta \in [0,4]$, $\sigma=[.5,4]$."} ggplot(subset(maxvariancesimRL3,what.max=='variance.bias'),aes(n1,-variance.bias))+geom_line()+ geom_line(data=subset(maxvariancesim3,what.max=='variance.bias'),linetype=2)+ylab("Maximum of the absolute variance bias")+ theme_bw()

```rComparison of maximal absolute variance bias for different parameter settings: Solid line - $\\delta=0$, $\\sigma=50$; Dashed line - $\\delta=0$, $\\sigma=4$; Dotted line $\\delta=0$, $\\sigma \\in [.5,4]$.",fig.align="center"}
ggplot(subset(maxvariancesimRL3,what.max=='variance.bias'),aes(n1,-variance.bias))+geom_line(linetype=3)+
    geom_line(data=subset(maxvariancesimR43,what.max=='variance.bias'),linetype=2)+
        geom_line(data=subset(maxvariancesimR503,what.max=='variance.bias'),linetype=1)+ylab("Maximum of the absolute variance bias")+ theme_bw()

set(maxvariancesimR3,j=which(duplicated(names(maxvariancesimR3))),value=NULL)
set(maxvariancesimR43,j=which(duplicated(names(maxvariancesimR43))),value=NULL)
set(maxvariancesimRL3,j=which(duplicated(names(maxvariancesimRL3))),value=NULL)

maxvariancesimR43[what.max=='variance.bias',] <- maxvariancesimRL3[what.max=='variance.bias',]

max.varbias2 <- ddply(maxvariancesimR43,"n1",transform,
                  max.variance.bias= variance.bias,
                  uc.max.variance.bias= uc.variance.bias,
                      method=c("variance.bias"="adjusted","uc.variance.bias"="unadjusted")[what.max],
                  theoretical = bound
                  )

max.varbias2 <- max.varbias2[,c('method','n1','delta','sigma','max.variance.bias','uc.max.variance.bias','theoretical')]

```r Maximum absolute bias of the variance estimate for given per group first stage sample sizes $n_1$ between $2$ and $50$ using fixed values $\delta=0$ and for the unadjusted reassessment rule also $\sigma=4$. Left column shows the results for the unadjusted sample size reassessment rule, right column for the adjusted sample size reassessment rule. The first row shows the value of the bias, the second the effect size $\delta$ and the third the standard deviation $\sigma$ at which the specific value is attained. The gray line denotes a Loess smoothed estimate.",fig.align='center'}

df <- melt(max.varbias2,id=c('n1','method','theoretical')) df <- subset(df,!(method=='adjusted' & variable == 'uc.max.variance.bias')) df <- subset(df,!(method=='unadjusted' & variable == 'max.variance.bias'))

df <- within(df,scale <- c("absolute","relative")[grepl("rel",variable)+1]) df <- within(df,variable <- sub("rel\.","",variable)) df <- within(df,variable <- sub("uc\.","",variable)) df <- df[df$variable %in% c('max.variance.bias','delta','sigma'),] df <- within(df,variable <- factor(variable)) df <- within(df,variable <- revalue(variable,c('max.variance.bias'='str_mvb'))) df <- within(df,variable <- relevel(variable,'delta')) df <- within(df,variable <- relevel(variable,'str_mvb')) df <- within(df,method <- relevel(method,'unadjusted'))

dfa <- df[df$scale=="absolute",]#"relative",]# scalehelper <- data.frame(n1=rep(-10,3),variable=factor(1:3,label=c("str_mvb",expression(delta),expression(sigma))),value=rep(0,3),scale=rep("absolute",3),method=rep("adjusted",3),theoretical=3) dfa <- rbind(dfa,scalehelper)

vplot <- ggplot(dfa,aes(n1,value))+ geom_line(aes(n1,-value),data=subset(dfa,variable=="str_mvb")) + geom_line(aes(n1,-theoretical),data=subset(dfa,variable=="str_mvb" & method=='unadjusted'),size=1.5,alpha=.2) + geom_point(data=subset(dfa,variable==levels(variable)[2]),size=.2) + geom_point(data=subset(dfa,variable==levels(variable)[3] & value < 5),size=.2) + geom_smooth(data=subset(dfa,variable==levels(variable)[3] & value <4),col='darkgray',fill='gray') +

geom_text(label="paste(sigma,'=',50)",parse=T,data=data.frame(n1=25,method="unadjusted",variable="sigma",value=1.25),size=8) +

geom_line(data=subset(dfa,variable==levels(variable)[2])) +

facet_grid(variable~method,scale='free',labeller=lverbose) +xlab(expression(n[1]))+ylab('')+theme_bw()+xlim(c(2,50))

print(vplot)

```r
## max.coverage <- ddply(maxsim,"n1",summarise,
##                       max.coverage = -(max(total.prob)-.05),
##                       max.coverage.bias.delta = delta[which.max(total.prob)],
##                       max.coverage.bias.sigma = sigma[which.max(total.prob)],
##                       uc.max.coverage= -(max(uc.total.prob)-.05),
##                       uc.max.coverage.bias.delta = delta[which.max(uc.total.prob)],
##                       uc.max.coverage.bias.sigma = sigma[which.max(uc.total.prob)],
##                       max.onesided.coverage = -(max(upper.prob)-.025),
##                       max.onesided.coverage.bias.delta = delta[which.max(upper.prob)],
##                       max.onesided.coverage.bias.sigma = sigma[which.max(upper.prob)],
##                       uc.max.onesided.coverage= -(max(uc.upper.prob)-.025),
##                       uc.max.onesided.coverage.bias.delta = delta[which.max(uc.upper.prob)],
##                       uc.max.onesided.coverage.bias.sigma = sigma[which.max(uc.upper.prob)]
##                       )

maxcoverageopt1 <- select_results(maxsim,c('total.prob','upper.prob','uc.total.prob','uc.upper.prob'),base_columns=c('n1','tn1','delta','sigma','d','s'),functional=which.max)

G2coverage <- add_epsilon(maxcoverageopt1,.2,.005,c('delta','sigma'))
plot_grid(maxsim,maxcoverageopt1,G2coverage,'sigma')
options(mc.cores=min(20,detectCores()-1))
maxcoveragesim2 <- simulate_batch(G2coverage,4*10^6)
## reshape and resim
maxcoverageopt2 <- select_results(maxcoveragesim2,c('total.prob','upper.prob','uc.total.prob','uc.upper.prob'),base_columns=c('n1','tn1','delta','sigma','d','s'),functional=which.max)
options(mc.cores=min(16,detectCores()-1))
maxcoveragesim3 <- simulate_batch(maxcoverageopt2,10^7)
fname <- paste('resim_coverage',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(maxcoveragesim3,file=fname)

data(maxcoveragesim3)

max.coverage2 <- ddply(maxcoveragesim3,"n1",transform,
                       max.twosided= -(total.prob-.05),
                       uc.max.twosided = -(uc.total.prob-.05),
                       max.onesided= -(upper.prob-.025),
                       uc.max.onesided = -(uc.upper.prob-.025),
                       method=c("upper.prob"="adjusted","uc.upper.prob"="unadjusted","total.prob"="adjusted","uc.total.prob"="unadjusted")[what.max],
                       side=c("total.prob"="two-sided","uc.total.prob"="two-sided","upper.prob"="one-sided","uc.upper.prob"="one-sided")[what.max]
                  )

max.coverage2 <- max.coverage2[,c('method','side','n1','delta','sigma','max.twosided','uc.max.twosided','max.onesided','uc.max.onesided')]

```r Maximum (absolute) difference between actual and nominal coverage probabilities for a given per group first stage sample sizes $n_1$ between $2$ and $50$. Left column shows the results for the unadjusted sample size reassessment rule, right column for the adjusted sample size reassessment rule. The first row shows the difference between actual and nominal coverage probabilities, the second the effect size $\delta$ and the third the standard deviation $\sigma$ at which the specific value is attained. The gray line denotes a Loess smoothed estimate.",fig.align='center'} df <- melt(max.coverage2,id=c('n1','method','side')) df <- subset(df,!(method=='adjusted' & variable %in% c('uc.max.onesided','uc.max.twosided'))) df <- subset(df,!(method=='unadjusted' & variable %in% c('max.onesided','max.twosided'))) df <- subset(df,!(side=='one-sided' & variable %in% c('uc.max.twosided','max.twosided'))) df <- subset(df,!(side=='two-sided' & variable %in% c('max.onesided','uc.max.onesided'))) df <- within(df,scale <- c("absolute","relative")[grepl("rel",variable)+1]) df <- within(df,variable <- sub("rel\.","",variable)) df <- within(df,variable <- sub("uc\.","",variable)) df <- df[df$variable %in% c('max.onesided','max.twosided','delta','sigma'),] df <- within(df,variable <- sub("(onesided)|(twosided)","coverage",variable)) df <- within(df,variable <- factor(variable)) df <- within(df,variable <- revalue(variable,c('max.coverage'='str_anc','delta'=expression(delta),'sigma'=expression(sigma)))) df <- within(df,variable <- relevel(variable,'delta')) df <- within(df,variable <- relevel(variable,'str_anc')) df <- within(df,method <- relevel(method,'unadjusted'))

dfa <- df[df$scale=="absolute",]#"relative",]# scalehelper <- data.frame(n1=rep(-10,6),variable=factor(1:3,label=c("str_anc","delta","sigma")),value=rep(0,6),scale=rep("absolute",3),method=rep("adjusted",3),side=rep(c('one-sided','two-sided'),each=3)) dfa <- rbind(dfa,scalehelper)

cplot <- ggplot(dfa) + geom_line(aes(n1,-100*value,linetype=side,shape=side),data=subset(dfa,variable=="str_anc")) + geom_smooth(aes(n1,value),data=subset(dfa,variable!="str_anc"),col='darkgray',fill='gray') + geom_point(aes(n1,value,shape=side),data=subset(dfa,variable!="str_anc"),size=.5) + facet_grid(variable~method,scale='free_y',labeller=lverbose)+xlab(expression(n[1]))+ylab('')+theme_bw()+xlim(c(2,50))

print(cplot)

```r
maxcoverageopt0 <- select_results(subset(maxsim,delta==0),c('total.prob','upper.prob','uc.total.prob','uc.upper.prob'),base_columns=c('n1','tn1','delta','sigma','d','s'),functional=which.max)

G2coverage <- add_epsilon(maxcoverageopt0,.2,.005,c('sigma'))
plot_grid(maxsim,maxcoverageopt0,G2coverage,'sigma')
options(mc.cores=min(20,detectCores()-1))
maxcoveragesim02 <- simulate_batch(data.frame(G2coverage),4*10^6,use=T)
##save(maxcoveragesim02,file="nullsim.rda")
## reshape and resim
maxcoverageopt02 <- select_results(maxcoveragesim02,c('total.prob','upper.prob','uc.total.prob','uc.upper.prob'),base_columns=c('n1','tn1','delta','sigma','d','s'),functional=which.max)
options(mc.cores=min(16,detectCores()-1))
maxcoveragesim03 <- simulate_batch(data.frame(maxcoverageopt02),10^7,use=T)
fname <- paste('~/newBlindedResults/resim_coverage0',Sys.info()['nodename'],'_',format(Sys.time(),"%y%m%d"),'.Rd',sep='')
save(maxcoveragesim03,file=fname)

data(nullsim)


max.coverage2 <- ddply(nullsim,"n1",transform,
                       max.twosided= -(total.prob-.05),
                       uc.max.twosided = -(uc.total.prob-.05),
                       max.onesided= -(upper.prob-.025),
                       uc.max.onesided = -(uc.upper.prob-.025),
                       method=c("upper.prob"="adjusted","uc.upper.prob"="unadjusted","total.prob"="adjusted","uc.total.prob"="unadjusted")[what.max],
                       side=c("total.prob"="two-sided","uc.total.prob"="two-sided","upper.prob"="one-sided","uc.upper.prob"="one-sided")[what.max]
                  )

max.coverage2 <- max.coverage2[,c('method','side','n1','delta','sigma','max.twosided','uc.max.twosided','max.onesided','uc.max.onesided')]

```r Maximum Type I error rate inflation of the $t$-test following blinded sample size reassessment based on blinded interim data for given per group first stage sample sizes $n_1$ between $2$ and $50$. Left column shows the results for the unadjusted sample size reassessment rule, right column for the adjusted sample size reassessment rule. The first row shows the difference between actual and nominal coverage probabilities, the second the effect size $\delta$ and the third the standard deviation $\sigma$ at which the specific value is attained. The gray line denotes a Loess smoothed estimate.",fig.align='center'} df <- melt(max.coverage2,id=c('n1','method','side')) df <- subset(df,!(method=='adjusted' & variable %in% c('uc.max.onesided','uc.max.twosided'))) df <- subset(df,!(method=='unadjusted' & variable %in% c('max.onesided','max.twosided'))) df <- subset(df,!(side=='one-sided' & variable %in% c('uc.max.twosided','max.twosided'))) df <- subset(df,!(side=='two-sided' & variable %in% c('max.onesided','uc.max.onesided'))) df <- within(df,scale <- c("absolute","relative")[grepl("rel",variable)+1]) df <- within(df,variable <- sub("rel\.","",variable)) df <- within(df,variable <- sub("uc\.","",variable)) df <- df[df$variable %in% c('max.onesided','max.twosided','delta','sigma'),] df <- within(df,variable <- sub("(onesided)|(twosided)","coverage",variable)) df <- within(df,variable <- factor(variable)) df <- within(df,variable <- revalue(variable,c('max.coverage'='str_t1','delta'=expression(delta),'sigma'=expression(sigma)))) df <- within(df,variable <- relevel(variable,'delta')) df <- within(df,variable <- relevel(variable,'str_t1')) df <- within(df,method <- relevel(method,'unadjusted'))

dfa <- df[df$scale=="absolute",]#"relative",]# scalehelper <- data.frame(n1=rep(-10,6),variable=factor(1:3,label=c("str_t1","delta","sigma")),value=rep(0,6),scale=rep("absolute",3),method=rep("adjusted",3),side=rep(c('one-sided','two-sided'),each=3)) dfa <- rbind(dfa,scalehelper)

t1plot <- ggplot(dfa) + geom_line(aes(n1,-100*value,linetype=side,shape=side),data=subset(dfa,variable=="str_t1")) + geom_smooth(aes(n1,value),data=subset(dfa,variable!="str_t1"),col='darkgray',fill='gray') + geom_point(aes(n1,value,shape=side),data=subset(dfa,variable!="str_t1"),size=.5) + facet_grid(variable~method,scale='free_y',labeller=lverbose)+xlab(expression(n[1]))+ylab('')+theme_bw()+xlim(c(2,50))

print(t1plot)

```r
library(bt88.03.704)

mplot <- mplot+theme(strip.text.y = element_text(size=rel(.8)),axis.text=element_text(size=rel(.6)),plot.margin=unit(c(.7,-.2,-.2,-1),'lines'))
vplot <- vplot+theme(strip.text.y = element_text(size=rel(.8)),axis.text=element_text(size=rel(.6)),plot.margin=unit(c(.7,-.2,-.2,-1),'lines'))
cplot <- cplot+theme(strip.text.y = element_text(size=rel(.8)),axis.text=element_text(size=rel(.6)),plot.margin=unit(c(.7,-1,-.2,-1),'lines'),legend.text=element_text(size=rel(.6)))
setwd('~/repos/overleaf/blind_confidence/graphs/')
cairo_ps('blinded-max.eps',width=8.5,height=4.4,pointsize=8)
multiplot(mplot,vplot,cplot,cols=3,widths=unit(c(.3,.3,.4),'npc'),titles = c("Mean bias","Variance bias","Actual - nominal coverage"))
dev.off()

lverbose = function(labels,multi_line=TRUE){
    variable <- names(labels)
    my_labels <- list()
    if(variable == "sigma"){
        my_labels[[1]] <- sapply(labels[[1]],function(val) substitute(paste(sigma," = ",foo,sep=""),list(foo=val)))#label_parsed(variable,label_both(variable,value))
        return(my_labels)
    } else {
        label_value(labels,multi_line)
    }

}

```r Coverage probabilities and bias of mean and variance estimates for the case study. The first panels shows actual - nominal coverage probabilities (AC-NC) for the $97.5\%$ upper (dashed line), lower (dotted line) and the $95\%$ two-sided confidence intervals, for the unadjusted reassessment rule $n^u_2$. The second panel shows the bias of the mean estimate if $n^u_2$ (solid line) or $n^a_2$ (dashed line) is used. The dotted line shows upper and lower bounds for the bias for general blinded sample size reassessment rules based on $S^2_{1,OS}$. The dashed gray line shows the bounds for the bias that can be attaind with a general unblinded sample size reassessment rule. The third panel shows the bias of the variance estimate if either $n_2^u$ (solid line) or $n^a_2$ (dashed line) is used. The red line shows the theoretical boundary for the bias give in Theorem 4."}

tcasesim <- casesim tcasesim <- subset(tcasesim, sigma %in% c(1,5,10,15,20))

df <- rbind(tcasesim,tcasesim,tcasesim) df$variable <- rep(c("AC-NC [%]","Bias of the mean","Bias of the variance"),each=nrow(tcasesim))

caseplot <- ggplot(df) + facet_grid(variable~sigma,scales='free_y',labeller=lverbose) + geom_path(aes(delta,uc.mean.bias),data=subset(df,variable == "Bias of the mean")) + geom_path(aes(delta,mean.bias),lty=2,data=subset(df,variable == "Bias of the mean")) + geom_path(aes(delta,mean.upper.bias),lty=3,data=subset(df,variable == "Bias of the mean" & mean.upper.bias < .8))+ geom_path(aes(delta,-mean.upper.bias),lty=3,data=subset(df,variable == "Bias of the mean" & mean.upper.bias < .8))+ geom_path(aes(delta,brannath),lty=2,col='darkgray',data=subset(df,variable == "Bias of the mean" & brannath < 1.2))+ geom_path(aes(delta,-brannath),lty=2,col='darkgray',data=subset(df,variable == "Bias of the mean" & brannath < 1.2))+ geom_path(aes(delta,uc.variance.bias),data=subset(df,variable == "Bias of the variance")) + geom_path(aes(delta,variance.bias),lty=2,data=subset(df,variable == "Bias of the variance")) + geom_path(aes(delta,bound),col='darkgray',data=subset(df,variable == "Bias of the variance")) + geom_path(aes(x=delta,y=100(.025-uc.upper.prob)),lty=2,data=subset(df,variable=="AC-NC [%]" & abs(uc.upper.prob < .035))) + geom_line(aes(x=delta,y=0),colour='gray',data=subset(df,variable=="AC-NC [%]" & abs(uc.upper.prob < .035))) + geom_path(aes(x=delta,y=100(.025-uc.lower.prob)),lty=3,data=subset(df,variable=="AC-NC [%]" & abs(uc.lower.prob < .035))) + geom_path(aes(x=delta,y=100*(.05-uc.total.prob)),data=subset(df,variable=="AC-NC [%]" & abs(uc.total.prob < .06))) + geom_blank(aes(x=-5,y=-0.7),data=subset(df,variable=="AC-NC [%]" & abs(uc.upper.prob < .035))) + geom_blank(aes(x=-5,y=0.05),data=subset(df,variable=="AC-NC [%]" & abs(uc.upper.prob < .035))) + ylab('')+xlab(expression(delta))+theme_bw()

print(caseplot)

```r
## this chunk produces the file that goes into the paper

cairo_ps('case.eps',height=5)
print(caseplot)
dev.off()

## The next chunk runs the simulations by applying the simulation
## function, in parallel, across rows of the parameter grid $G$, which
## contains the parameter settings for each scenario. In order to rerun
## them it has to be reconfigured for evaluation. Note that running the
## simulations may take quite long.
G <- expand.grid(delta=round(seq(-4,4,.1),1),
                 sigma=round(seq(.5,2,.5),2),
p                 d = 1,
                 s = round(seq(.5,2,.5),2))
G$n1 <- ceiling(1/2*zss(G$s,G$d,.025,.2))
G$n2max <- 2*G$n1

G <- subset(G,n1 < 65 & n1 > 1)
G <- data.table(G)
set.seed(693325)
runs <- 10^7
restsim <- simulate_batch(G,runs,T,T)

full.sim <- restsim

restsim <- subset(full.sim,d==1 & sigma>=.5 & n1 > 2)
restsim$ylow <- 0

gsim <- as.data.frame(t(apply(restsim,1,unlist)))

## compute relative bias
gsim$var.rbias <- gsim$variance.bias/gsim$sigma^2

## compute the theoretical lower bound for the variance bias
gsim$bound <- lowerBound(gsim$n1,gsim$d)

runs <- 10^5
bsim <- mclapply2(1:nrow(gsim),function(i) {c(gsim[i,],
                                             simMBIA(delta=gsim[i,]$delta,
                                                     sigma=gsim[i,]$sigma,
                                                     n1=gsim[i,]$n1,
                                                     n2max=gsim[i,]$n2max,
                                                     runs=runs))})

bsim <- do.call('cbind',bsim)
bsim <- as.data.frame(apply(bsim,1,unlist),row.names=NA)

restsim <- bsim

restsim$brannath <- with(restsim,n2max/(n1 +n2max)*.4*sqrt(2)*sigma/sqrt(n1))




save(restsim,file=bt88.03.704::vfile('restsim'))

## The next chunk runs the simulations by applying the simulation
## function, in parallel, across rows of the parameter grid $G$, which
## contains the parameter settings for each scenario. In order to rerun
## them it has to be reconfigured for evaluation. Note that running the
## simulations may take quite long.
G <- expand.grid(delta=round(seq(-4,4,.1),1),
                 sigma=round(seq(.5,2,.5),2),
                 d = 1,
                 s = round(seq(.5,2,.5),2))
G$n1 <- ceiling(1/2*zss(G$s,G$d,.025,.2))
G$n2max <- 2*G$n1
G$n2min <- G$n1/2

G <- subset(G,n1 < 65 & n1 > 1)
G <- data.table(G)
set.seed(693325)
runs <- 10^7
restsim <- simulate_batch(G,runs,T,T)

full.sim <- restsim
save(restsim,file=vfile('furest_sim_nomaxbias'))

restsim <- subset(full.sim,d==1 & sigma>=.5 & n1 > 2)
restsim$ylow <- 0

gsim <- as.data.frame(t(apply(restsim,1,unlist)))

## compute relative bias
gsim$var.rbias <- gsim$variance.bias/gsim$sigma^2

## compute the theoretical lower bound for the variance bias
gsim$bound <- lowerBound(gsim$n1,gsim$d)

runs <- 10^5
bsim <- mclapply2(1:nrow(gsim),function(i) {c(gsim[i,],
                                             simMBIA(delta=gsim[i,]$delta,
                                                     sigma=gsim[i,]$sigma,
                                                     n1=gsim[i,]$n1,
                                                     n2max=gsim[i,]$n2max,
                                                     runs=runs))})

bsim <- do.call('cbind',bsim)
bsim <- as.data.frame(apply(bsim,1,unlist),row.names=NA)

restsim <- bsim

restsim$brannath <- with(restsim,n2max/(n1 +n2max)*.4*sqrt(2)*sigma/sqrt(n1))




save(restsim,file=bt88.03.704::vfile('restsim'))

data(restsim)

```r"}

foo <- plot_coverage(restsim,adjusted=T) foo + geom_point(aes(x,y),data=data.frame(x=0,y=-0.78,sigma=1,s=1.5),alpha=0) + geom_point(aes(x,y),data=data.frame(x=0,y=-0.48,sigma=1,s=2),alpha=0)

```r$  with a maximal second stage sample size limited to $2n_1$: upper confidence bound (dashed line),  lower confidence bound (dotted line) and  two-sided interval (solid line). Rows refer to the a priori assumed standard deviations $\\sigma_0$ determining the first stage sample size $n_1$. The columns correspond to actual standard deviations. The x-axis in each graph denotes the true treatment effect $\\delta$, the y-axis shows the difference of actual to nominal coverage probability such that negative values indicate settings where the confidence bound is not valid.\\label{fig:rest.uc.coverage}"}
foo  <- plot_coverage(restsim,adjusted=T)
foo + geom_point(aes(x,y),data=data.frame(x=0,y=-0.99,sigma=1,s=1.5),alpha=0) + geom_point(aes(x,y),data=data.frame(x=0,y=-0.99,sigma=1,s=2),alpha=0)

```rBias of the mean under blinded sample size reassessment with a maximal second stage sample size limited to $2n_1$ using the unadjusted (solid line) and adjusted (dashed line) interim variance estimate. The dotted lines show maximum negative and positive bias that according to Theorem 2 can be attained under any blinded sample size reassessment rule with a maximal second stage sample size limited to $2n_1$. The dashed gray line shows the maximum bias that can be attained under any (unblinded) sample size reassessment rule with a maximal second stage sample size limited to $2n_1$. The treatment effect used for planning is set to $\delta_0=1$. Rows refer to the a priori assumed standard deviations $\sigma_0$ determining the first stage sample size $n_1$. The columns correspond to actual standard deviations. The x-axis in each graph denotes the true treatment effect $\delta$, the y-axis shows the bias."}

foo <- mean_bias_plot(restsim) foo

```rBias of the variance under blinded sample size reassessment  with a maximal second stage sample size limited to $2n_1$ using the unadjusted (solid line) and adjusted (dashed line) interim variance estimate. The gray line gives the lower bound from Theorem 4 for the bias under sample size reassessment based on the unadjusted variance estimate. The treatment effect used for planning is set to $\\delta_0=1$. Rows refer to the a priori assumed standard deviations $\\sigma_0$ determining the first stage sample size $n_1$. The columns correspond to actual standard deviations. The x-axis in each graph denotes the true treatment effect $\\delta$, the y-axis shows the bias."}
foo <- variance_bias_plot(restsim)
foo + geom_point(aes(x,y),data=data.frame(x=0,y=-0.095,sigma=1,s=1),alpha=0) + geom_point(aes(x,y),data=data.frame(x=0,y=-0.095,sigma=1,s=1.5),alpha=0) + geom_point(aes(x,y),data=data.frame(x=0,y=-0.095,sigma=1,s=2),alpha=0)

```r Adjusted sample size reassessment rule with a maximal second stage sample size limited to $2n_1$: actual variance of the mean estimate $\sigma^2_e$ (solid line), average of the estimated variances of the mean $S^2_e$ (dashed line), average of the variances $\sigma^2_f$ of the mean of fixed sample trials with the same sample sizes."} foo <- sem_plot(restsim,adjusted=T) foo

```r Variance of the effect estimates (i.e., square of the standard errors) for the unadjusted sample size reassessment rule  with a maximal second stage sample size limited to $2n_1$:  actual variance $\\sigma^2_e$ of the mean  estimate (solid line), average of the estimated variances $S^2_e$ of the mean  (dashed line), average of the variance $\\sigma^2_f$ of the mean of corresponding fixed sample trials (dotted line)."}
foo <- sem_plot(restsim)
foo

\clearpage

References

floatofmath/blindConfidence documentation built on May 16, 2019, 1:20 p.m.