In floatofmath/adaperm: Permutation tests for adaptive designs
knitr::opts_chunk$set(eval=FALSE,echo=FALSE,warning=FALSE,results='hide',fig.width=7,fig.height=7,message=FALSE)
library(magrittr)
library(plyr)
library(dplyr)
library(microbenchmark)
library(pander)
library(parallel)
library(bt88.03.704)
library(adaperm)
#devtools::install_github('floatofmath/adaperm')
library(ldbounds)
options(mc.cores=detectCores()-1)
Adjusted p-values
Let $T_{(1)},...,T_{(|G|)}$ the ordered resampling test statistics of
the pre-planned test. Then the $p$-value of the test following an
adaptive interim analysis is given by $\tilde{p} = 1-\frac{r^*}{|G|}$
where
$$
r^*=\min { r : p_E \leq A(T_{(r)})}
$$
and $A(t)$ is defined as the conditional survival function of the
preplanned test given the first stage treatment assignments:
This is easily computed:
sim_pvals <- function(B,delta,nulldist=rnorm,...){
replicate(B,{
x1 <- nulldist(10,...)+rep(0:1,5)*delta
x2 <- nulldist(10,...)+rep(0:1,5)*delta
xE <- nulldist(6,...)+rep(0:1,3)*delta
pt <- t.test(c(x1,x2,xE)~rep(0:1,13),alternative='less',var.equal=T)$p.value
pdist <- perm_dist(x1,x2,rep(0:1,5),rep(0:1,5),ranksum,B=10000)
pw <- wilcox.test(c(x1,x2,xE)~rep(0:1,13),alternative='less')$p.value
cdist <- adaperm:::cond_dist(x1,x2,rep(0:1,5),rep(0:1,5),ranksum,B=10000)
edist <- perm_dist(x2,xE,rep(0:1,5),rep(0:1,3),ranksum,B=10000)
t2 <- ranksum(c(x2,xE),rep(0:1,8))
c(sum(pdist>=quantile(cdist,sum(edist<t2)/length(edist)))/length(pdist),pt,pw)
})
}
## pval_dist <- sim_pvals(1000,1,rnorm_cont,csd=10)
## head(t(pval_dist))
## cor(t(pval_dist)[,c(1,3)])
## plot(t(pval_dist)[,c(1,3)])
## par(mfrow=c(3,1))
## plot(sort(pval_dist[1,]),seq(0,1,length.out=1000))
## plot(sort(pval_dist[2,]),seq(0,1,length.out=1000))
## plot(sort(pval_dist[3,]),seq(0,1,length.out=1000))
## lookin good!
## rowSums(pval_dist<.05)/1000
## hist(pval_dist[1,])
adaptive_permtest_ts <- function(x1,x2,xE,g1,g2,gE,stat,permutations=10000){
pdist <- perm_dist(x1,x2,g1,g2,stat,B=permutations)
cdist <- adaperm:::cond_dist(x1,x2,g1,g2,stat,B=permutations)
edist <- perm_dist(x2,xE,g2,gE,stat,B=permutations)
t2 <- stat(c(x2,xE),c(g2,gE))
pval <- sum(pdist>=quantile(cdist,sum(edist<t2)/length(edist)))/length(pdist)
names(t2) <- deparse(substitute(stat))
out <- list(statistic = t2,
parameter = c('permutations'=permutations),
p.value = pval,
alternative = "one.sided",
null.value = c("distribution of treated samples"='different from control samples'),
data.name = paste(deparse(substitute(c(x1,x2,xE,g1,g2,gE))),collapse=', '),
method = "Two-sample permutation test for adaptive designs")
class(out) <- 'htest'
out
}
str(t.test(rnorm(10)))
t.test(rnorm(10))
adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,maritzm)
adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,ranksum)
str(adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,maritzm))
Robus m-Tests based on maritz' m-estimator
Maritz et al. and Rosenbaum suggest the use of an $m$-estimator
(Huber) for use with a rerandomization test. In adaperm we have
implemented such a test statistic in the function maritzm.
compare_adaptive_tests <- function(n1,n,rule,rdist,test_statistic,mpfr=F,resam=100,perms=100,...){
x <- rdist(n,...)
ne <- rule(x[1:n1])
if(ne>n){
x <- c(x,rdist(ne-n,...))
} else {
ne <- n
}
list(ne = ne,
permtestT = adaptive_permtest_os(x,n1,n,ne,diffmean,perms=perms),
invnormT = adaptive_invnorm_ttest_os(x,n1,n,ne),
permtestW = adaptive_permtest_os(x,n1,n,ne,signedranks,perms=perms),
invnormW = adaptive_invnorm_wilcoxtest_os(x,n1,n,ne),
permtestM = adaptive_permtest_os(x,n1,n,ne,maritzm,perms=perms))
}
##' Compare operating characteristics for a number of one-sample tests
##'
##' \code{test_funs} needs to be a list of functions that take arguments \code{x,n1,n,ne,permutations} and return \code{TRUE} of \code{FALSE} depending on whether the corresponding decision rule rejects the null hypothesis or not.
##'
##' \code{rule} needs to be a function that takes as argument the vector of first stage observations and return an interger number, which (if larger than \code{n}) will be used as the new total sample size.
##'
##' @title Compare one-sample test procedures
##' @param test_funs list of test functions (see Details)
##' @param n1 first stage sample size
##' @param n pre-planned total sample size
##' @param rule sample size reassessment rule (see Details)
##' @param rdist distribution from which to draw observations
##' @param ... additional arguments to rdist
##' @return let's see
##' @author Florian Klinglmueller
compare_tests_onesample <- function(B,test_funs,n1,n,rule,rdist,...){
run <- function(test_funs,n1,n,rule,rdist,...){
x <- rdist(n,...)
ne <- rule(x[1:n1])
if(ne>n){
x <- c(x,rdist(ne-n,...))
} else {
ne <- n
}
c(n1=n1,n=n,ne=ne,lapply(test_funs,do.call,list(x,n1,n,ne)))
}
results <- mclapply2(1:B,function(i) run(test_funs,n1,n,rule,rdist,...))
results %>% bind_rows -> results
eval(expr=parse(text=paste0("summarize(results,ASN=mean(ne),maxSN=max(ne),sdSN=sd(ne),",paste0(names(tests),"=","mean(",names(tests),")",collapse=','),")")))
}
tests <- list(
permtestT = function(x,n1,n,ne) adaptive_permtest_os(x,n1,n,ne,diffmean,perms=10000),
invnormT = function(x,n1,n,ne) adaptive_invnorm_ttest_os(x,n1,n,ne),
permtestW = function(x,n1,n,ne) adaptive_permtest_os(x,n1,n,ne,signedranks,perms=10000),
invnormW = function(x,n1,n,ne) adaptive_invnorm_wilcoxtest_os(x,n1,n,ne),
permtestM = function(x,n1,n,ne) adaptive_permtest_os(x,n1,n,ne,maritzm,perms=10000))
## Type I error
## compare_tests_onesample(1000,tests,10,20,cond_power_rule_norm,rnorm_cont,mean=0,sd=1.2,cshift=0,csd=3)
## ## Power
## compare_tests_onesample(1000,tests,10,20,cond_power_rule_norm,rnorm_cont,mean=.5,sd=1,cshift=.5,csd=5)
delta <- 0
x1 <- rnorm(10)+rep(0:1,5)*delta
x2 <- rnorm(10)+rep(0:1,5)*delta
xE <- rnorm(6)+rep(0:1,3)*delta
g1 <- sign(x1);g2 <- sign(x2);gE <- sign(xE)
pt <- t.test(c(x1,x2,xE),alternative='less')$p.value
pdist <- perm_dist(x1,x2,g1,g2,maritzm,B=10000,restricted=F)
pw <- wilcox.test(c(x1,x2,xE),alternative='less')$p.value
cdist <- adaperm:::cond_dist(x1,x2,g1,g2,maritzm,B=10000,restricted=F)
edist <- perm_dist(x2,xE,g2,gE,maritzm,B=10000,restricted=F)
t2 <- maritzm(c(x2,xE),rep(0:1,8))
c(sum(pdist>=quantile(cdist,sum(edist<t2)/length(edist)))/length(pdist),pt,pw)
adaptive_permtest_os(c(x1,x2,xE),10,20,26,maritzm,perms=10000)
permutation_CER(abs(x1),g1>0,abs(x2),maritzm,permutations=10000,restricted=F)
perm_test(abs(x2),abs(xE),g2,gE,maritzm,B=10000,restricted=F)
par(mfrow=c(3,1))
xl <- min(c(pdist,cdist,edist))
xu <- max(c(pdist,cdist,edist))
hist(pdist,xlim=c(xl,xu))
hist(cdist,xlim=c(xl,xu))
hist(edist,xlim=c(xl,xu))
abline(v=t2)
Point estimates and confidence intervals
##' Permutation test p-value
##'
##' Type "non-randomized" will give the standard p-value, which is define by the proportion of sample permutations with larger or equal test statistic than the observed one. "midp" will give the mid-p value which is the proportion of sample permutations with larger test statistic plus one half the proportion of permutations with equal test statistic. "randomized" will return the proportion of permutations with larger test statistic plus a random number that is uniformly distributed between zero and the proportion of permutations with test statistic equal to the observed.
##'
##' @title Permutation p-value
##' @param t
##' @param dist
##' @param type
##' @return p-value
##' @author Florian Klinglmueller
t2p <- function(t,dist,type=c('non-randomized','midp','randomized')){
switch(type[1],
'non-randomized'=mean(dist>=t),
'midp'=mean(dist>t) + .5*mean(dist==t),
'randomized'=mean(dist>t) + runif(1)*mean(dist==t))
}
##' Level alpha critical boundary of permutation test
##'
##' For type 'non-randomized' this returns the smallest unique value of the permutation distribution such that the proportion of statistics larger is smaller than alpha. Type 'randomized' is currently not working but should give the critila value of a randomized test with exact level p.
##'
##' @title p to t
##' @param p p-value/alpha level
##' @param dist permutation distribution
##' @param type see details
##' @return critical value
##' @author Florian Klinglmueller
p2t <- function(p,dist,type=c('non-randomized')){#,'randomized')){
n <- length(dist)
ta <- max(dist[rank(-dist,ties='min') >= p*length(dist)])
switch(type[1],
'non-randomized'= ta)#,
# 'randomized' = ta - (p-floor(p*n)/n)/2
}
perm_DR <- function (x1, x2, g1, g2, stat, B, x3 = NULL, g3 = NULL,
restricted,
type=c('non-randomized','randomized','midp'),alpha)
{
perm_test(x1,x2,g1,g2,stat,B,x3,g3,restricted,type)<=alpha
}
perm_test <- function (x1, x2, g1, g2, stat, B, x3 = NULL, g3 = NULL,
restricted,
type=c('non-randomized','randomized','midp'))
{
dist <- perm_dist(x1, x2, g1, g2, stat, B, x3, g3, restricted = restricted)
t <- stat(c(x1, x2, x3), c(g1, g2, g3))
switch(type[1],
'non-randomized'=t2p(t,dist,type="non-randomized"),
'randomized'= t2p(t,dist,type="randomized"),
'midp'=t2p(t,dist,type="midp"))
}
adaptive_permtest_quick <- function(x1,x2,xE,g1,g2,gE,stat,permutations=10000,restricted){
pdist <- perm_dist(x1,x2,g1,g2,stat,B=permutations,restricted=restricted)
cdist <- adaperm:::cond_dist(x1,x2,g1,g2,stat,B=permutations,restricted=restricted)
edist <- perm_dist(x2,xE,g2,gE,stat,B=permutations,restricted=restricted)
t2 <- stat(c(x2,xE),c(g2,gE))
pval <- mean(pdist >= quantile(cdist, mean(edist<t2),type=1))
pval
}
bisect <- function(err,start,end,tolerance,maxit){
l <- err(start)
u <- err(end)
if(sign(l*u)>0) stop('err() values at end points not of opposite sign')
if(maxit<1) {
warning('Bisection reached maximum number of iterations')
return(err((end-start)/2))
}
m <- err((end-start)/2)
if(m==0|abs(m) < tolerance) return(start-end)
if(sign(u*m)>0){
bisect(err,start,(end-start)/2,tolerance,maxit-1)
} else {
bisect(err,(end-start)/2,end,tolerance,maxit-1)
}
}
adaptive_permtest_ts <- function(x1,x2,xE,g1,g2,gE,stat,permutations=10000,conf_level=.05){
pval <- adaptive_permtest_quick(x1,x2,xE,g1,g2,gE,stat,permutations=10000,restricted=T)
err <- function(d,conf_level=.5){
x1 <- x1-d*g1
x2 <- x2-d*g2
xE <- xE-d*gE
adaptive_permtest_quick(x1,x2,xE,g1,g2,gE,stat,permutations=10000,restricted=T)-conf_level
}
t2 <- stat(c(x2,xE),c(g2,gE))
names(t2) <- deparse(substitute(stat))
## rms
guesstimate <- 2*sqrt(mean(c(x1,x2,xE)^2))*sign(.5-pval)
estimate <- c("constant additive effect"=uniroot(err,c(0,guesstimate))$root)
conf.int <- c(uniroot(err,c(0,guesstimate*sign(conf_level-pval)),conf_level=conf_level)$root,Inf)
attr(conf.int,"conf.level") <- 1-conf_level
out <- list(statistic = t2,
parameter = c('permutations'=permutations),
p.value = pval,
estimate = estimate,
conf.int=conf.int,
alternative = "one.sided",
null.value = c("distribution of treated samples"='different from control samples'),
data.name = paste(deparse(substitute(c(x1,x2,xE,g1,g2,gE))),collapse=', '),
method = "Two-sample permutation test for adaptive designs")
class(out) <- 'htest'
out
}
g1 <- rep(0:1,5)
g2 <- rep(0:1,5)
gE <- rep(0:1,3)
x1 <- rnorm(10)+g1
x2 <- rnorm(10)+g2
xE <- rnorm(6)+gE
adaptive_permtest_quick(x1,x2,xE,g1,g2,gE,maritzm,perm=10000,restricted=T)
t.test(-c(x1,x2,xE)~c(g1,g2,gE),alternative='greater')
adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,maritzm,perm=10000)
adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,ranksum,perm=10000)
adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,meandiff,perm=10000)
Case study
Synthetic data
print(n0 <- power.z.test(delta=6,sd=6,power=.8))
set.seed(5449219)
x <- rnorm_cont(16,25,8,cshift=0,csd=24)
y <- rnorm_cont(16,19,8,cshift=0,csd=24)
boxplot(cbind(x,y),ylim=c(0,56))
## using the normal distribution conditional power rule suggested overall sample size is 48
cond_power_rule_norm_ts(x[1:8],y[1:8],delta=6,target=.8,rob_var=F)
## using the robust variance estiamte the suggested sample size rule is 32
print(nE <- cond_power_rule_norm_ts(x[1:8],y[1:8],delta=6,target=.8,rob_var=T))
xE <- rnorm_cont(nE-16,25,8,cshift=0,csd=24)
yE <- rnorm_cont(nE-16,19,8,cshift=0,csd=24)
adaperm(c(y,yE,x,xE),g=rep(0:1,each=nE),8,16,test_statistic=ranksum,cer_type='randomized',atest_type='non-randomized')
adaperm(c(y,yE,x,xE),g=rep(0:1,each=nE),8,16,test_statistic=tstat,cer_type='randomized',atest_type='non-randomized')
adaptive_invnormtest_2s(c(y,yE),c(x,xE),8,16,nE)
## invnorm estimates and CI
dec <- function(delta,alpha){
2*adaptive_invnormtest_2s(c(y,yE)+delta,c(x,xE),8,16,nE,alpha=alpha)-1
}
bisect(-10,10,function(x) dec(x,alpha=.025))
bisect(-10,10,function(x) dec(x,alpha=.5))
bisect(-15,15,function(x) dec(x,alpha=.975))
adaptive_invnorm_wilcoxtest_2s(c(y,yE),c(x,xE),8,16,nE)
dec <- function(delta,alpha){
2*adaptive_invnorm_wilcoxtest_2s(c(y,yE)+delta,c(x,xE),8,16,nE,alpha=alpha)-1
}
bisect(-10,10,function(x) dec(x,alpha=.025))
bisect(-10,10,function(x) dec(x,alpha=.5))
bisect(-15,15,function(x) dec(x,alpha=.975))
## almost identical - need to compute the wilcox.test twice to get matching p-values CIs
wilcox.test(x,y,conf.int=T,alternative='greater',correct=FALSE,exact=TRUE,conf.level=0.975)
wilcox.test(x,y,conf.int=T,alternative='two.sided',correct=FALSE,exact=TRUE,conf.level=0.95)
## may take a bit due to many permutations
adaperm(c(y,x),g=rep(0:1,each=16),8,16,test_statistic=ranksum,permutations=10^5,cer_type='randomized')
adaperm(c(y,x),g=rep(0:1,each=16),8,16,test_statistic=ranksum,permutations=10^5,cer_type='non-randomized')
delta <- 0
x1 <- c(y[1:8],x[1:8])
x2 <- c(y[9:16],x[9:16])
x3 <- c(yE,xE)
g1 <- rep(0:1,each=8)
g2 <- rep(0:1,each=8)
g3 <- rep(0:1,each=nE-16)
pt <- t.test(c(x1,x2,x3)~c(g1,g2,g3),alternative='less')$p.value
pdist <- perm_dist(x1,x2,g1,g2,ranksum,B=10000,restricted=T)
pw <- wilcox.test(c(x1,x2,x3)~c(g1,g2,g3),alternative='less')$p.value
cdist <- adaperm:::cond_dist(x1,x2,g1,g2,ranksum,B=10000,restricted=T)
edist <- perm_dist(x2,x3,g2,g3,ranksum,B=10000,restricted=T)
t2 <- ranksum(c(x2,x3),c(g2,g3))
c(sum(pdist>=quantile(cdist,sum(edist<t2)/length(edist)))/length(pdist),pt,pw)
par(mfrow=c(3,1))
xl <- min(c(pdist,cdist,edist))
xu <- max(c(pdist,cdist,edist))
hist(pdist,xlim=c(xl,xu),main='Permutation distribution of rank-sum statistic')
hist(cdist,xlim=c(xl,xu),main='Conditional permutation distribution given first stage treatment assignments')
hist(edist,xlim=c(xl,xu),main='Permutation distribution of extended second stage data')
abline(v=t2)
## Group sequential stuff
alpha <- bounds(c(.5,1),alpha=.025)$exit.pr
adaperm(c(y,x),g=rep(0:1,each=16),8,16,alpha0=alpha[1],test_statistic=ranksum,permutations=10^5)
Real data
library(coin)
library(gMCP)
library(ggplot2)
library(nparcomp)
data(rotarod)
## not useful
ggplot(rotarod,aes(group,time))+geom_boxplot()
data(neuropathy)
## heteroskedastic
ggplot(neuropathy,aes(group,pain))+geom_boxplot()
library(multcomp)
data(recovery)
head(recovery)
ggplot(recovery,aes(blanket,minutes))+geom_boxplot()
recovery %>% group_by(blanket) %>% group_size
#maybe(b0 vs b3)
data(mtept)
head(mtept)
ggplot(mtept,aes(treatment,E4))+geom_boxplot()
mtept %>% group_by(blanket) %>% group_size
library(MASS)
data(anorexia)
head(anorexia)
ggplot(anorexia,aes(Treat,Prewt-Postwt))+geom_boxplot()
anorexia %>% group_by(Treat) %>% group_size
power.t.test(delta=.25,sd=.25,power=.8,sig.level=.025,alternative="one.sided")
s <- sample(10000:99999,1)
set.seed(s)
stage <- sample(rep(1:3,c(9,9,12)))
neuropathy$stage <- c(stage[-which(stage==3)[(1:2)]],stage)
x1 <- neuropathy[neuropathy$stage==1,'pain']
g1 <- rep(0:1,each=9)
x2 <- neuropathy[neuropathy$stage==2,'pain']
g2 <- rep(0:1,each=9)
xE <- neuropathy[neuropathy$stage==3,][-(19:20),'pain']
gE <- rep(0:1,each=10)
cond_power_rule_t_ts(x1[g1==0],x1[g1==1],delta=1/4,target=.8)
73541::22
92203::27
38478::24
neuropathy %>% group_by(group) %>% group_size()
neuropathy %>% group_by(group) %>% mutate(stage = sample
Group sequential procedures
## Obrien Fleming type boundaries
## inputs
alpha <- bounds(c(.5,1),alpha=.025)$exit.pr
g1 <- rep(0:1,10)
g2 <- rep(0:1,10)
x1 <- rnorm(20)
x2 <- rnorm(20)
## function
gsd2_permtest_ts <- function(x1,x2,g1,g2,alpha,stat,permutations){
G <- omega(g1,g2,B= permutations)
dist1 <- stat(x1,G[1:length(g1),])
dist2 <- stat(c(x1,x2),G)
to1 <- stat(x1,g1)
to2 <- stat(c(x1,x2),c(g1,g2))
B <- length(dist1)
R1 <- rank(dist1,ties.method="min")>= B*(1-alpha[1])
R2 <- rep(FALSE,B)
R2[!R1] <- rank(dist2[!R1],ties.method="min")>= (B*(1-alpha[2]) + alpha[2]*sum(R1))
Ta1 = min(dist1[R1])
Ta2 = min(dist2[R2])
## Joint rejection probability
list(jrp = sum(R1 | R2)/B,
Ta1 = Ta1,
Ta2 = Ta2,
to1 = to1,
to2 = to2,
decision = (to1>=Ta1|to2>=Ta2))
}
gsd2_permtest_ts(x1,x2,g1,x2,alpha,meandiff,10^6)
Adaptive group sequential procedures
## Obrien Fleming type boundaries
## inputs
alpha <- bounds(c(.5,1),alpha=.025)$exit.pr
g1 <- rep(0:1,10)
g2 <- rep(0:1,10)
gE <- rep(0:1,5)
x1 <- rnorm(20)
x2 <- rnorm(20)
xE <- rnorm(10)
permutations <- 10^4
stat <- maritzm
gsd2_permtest_ts <- function(x1,x2,xE,g1,g2,gE,alpha,stat,permutations){
G <- omega(g1,g2,B= permutations)
dist1 <- stat(x1,G[1:length(g1),])
dist2 <- stat(c(x1,x2),G)
to1 <- stat(x1,g1)
to2 <- stat(c(x1,x2),c(g1,g2))
B <- length(dist1)
R1 <- rank(dist1,ties.method="min")>= B*(1-alpha[1])
R2 <- rep(FALSE,B)
R2[!R1] <- rank(dist2[!R1],ties.method="min")>= (B*(1-alpha[2]) + alpha[2]*sum(R1))
Ta1 = min(dist1[R1])
Ta2 = min(dist2[R2])
jrp = sum(R1 | R2)/B
cdist <- adaperm:::cond_dist(x1,x2,g1,g2,stat,permutations)
## reuse the alpha trimmings - this may be controversial if A=0 even though the test would control the error rate
A <- sum(cdist>=Ta2)/length(cdist)+max(0,alpha[2]-jrp)
GE <- omega(c(g2,gE),B=permutations)
## dirty trick to get non-finite cricital values if condiitonal error is exactly 0
distE <- c(stat(c(x2,xE),GE),Inf)
toE <- stat(c(x2,xE),c(g1,gE))
## The ranks are not changed by adding Inf, but now the vector won't be empty if A = 0
TaE <- min(distE[rank(distE,ties.method='min') > (1-A)*(length(distE)-1)])
## Joint rejection probability
list(a1 = sum(R1)/B,
jrp = jrp,
Ta1 = Ta1,
TaE = TaE,
to1 = to1,
toE = toE,
decision = (to1>=Ta1|toE>=TaE))
}
gsd2_permtest_ts(x1+g1,x2+g2,xE+gE,g1,g2,gE,alpha,meandiff,10^5)
Efficiency considerations
different ways to subset some subsetting
select1 <- function(G,sg1){
tG <- t(G)
colnames(tG) <- paste0('V',1:ncol(tG))
cond <- paste0('V',1:length(sg1),'==sg1[',1:length(sg1),']',collapse='&')
tG[with(data.frame(tG),eval(parse(text=cond))),]
}
select2 <- function(G,sg1){
tG <- t(G)
merge(matrix(sg1,nr=1),tG,by=1:length(sg1))
}
select3 <- function(G,sg1){
tG <- as.data.table(t(G))
sg1 <- as.data.table(as.data.frame(matrix(sg1,nr=1)))
setkeyv(tG,paste0('V',1:length(sg1)))
setkeyv(sg1,paste0('V',1:length(sg1)))
tG[sg1,]
}
select4 <- function(G,sg1){
inner_join(t(G) %>% as.data.frame %>% tbl_df,
matrix(sg1,nr=1) %>% as.data.frame %>% tbl_df,by=paste0('V',1:length(sg1)))
}
select5 <- function(G,sg1){
## this may be dangerous and does not save time!!
suppressMessages(inner_join(t(G) %>% as.data.frame %>% tbl_df,
matrix(sg1,nr=1) %>% as.data.frame %>% tbl_df))
}
n1=10
g1 <- rep(0:1,n1)
g2 <- rep(0:1,n1)
G <- omega(g1,g2,B=10^6)
dim(G)
sg1 <- sample(g1)
dim(select1(G,sg1))
microbenchmark(select1(G,sg1),
# select2(G,sg1),
# select3(G,sg1),
select4(G,sg1),
select5(G,sg1))
Conditional error of the randomized test
Typically the randomized decision rule of the permutation test is
defined as
$$
\phi = \left{\begin{array}{ll} 1 & \mbox{if}\; T_o > T_\alpha \
\gamma & \mbox{if}\; T_o = T_\alpha \
0 & \mbox{otherwise.}
\end{array}\right.
$$
where $T_\alpha$ is the $\lceil (1-\alpha) * |G|\rceil$ order
statistic (using the minimum rank for ties) of the orbit ${T(x): x
\in S(x)}$ and $\gamma$ is defined as
$$
\frac{\alpha - P(T > T_\alpha)}{P(T=T_\alpha)}.
$$
The corresponding conditional error function is then
$$
A_R'=A' + \gamma P(T = T_\alpha| x_{|n_1} ).
$$
If either the $T_\alpha$ or the conditional permutation distribution
given the first stage observations are estimated using Monte-Carlo
permutation algorithm, it may be convenient to compute the conditional
error function of a slightly different randomized decision rule, which
flips a coin for all values of the test statistic falling within
$\epsilon$ of the critical value:
$$
\phi_\epsilon = \gamma \mbox{if}\; T_o \in [T_\alpha - \epsilon,T_\alpha).
$$
In this case $\gamma$ is given by
$$
\gamma = \frac{\alpha - P(T > T_\alpha)}{P(T \in [T_\alpha-\epsilon,T_\alpha))}
$$
which can be easily computed even for a Monte Carlo sample from the permutation space.
And the conditional error function is:
$$
A'\epsilon = A' + P(T \in [T\alpha-\epsilon,T_\alpha)|x_{|n_1})\frac{\alpha - P(T > T_\alpha)}{P(T \in [T_\alpha-\epsilon,T_\alpha))}
$$
For $\epsilon \rightarrow \infty$ one gets:
$$
A'\infty = A' + (1-A')\frac{\alpha - P(T >T\alpha)}{1-P(T > T_\alpha)}
$$
Alternatively one may use the conditional error function:
$$
A_t'= A' + (\alpha - P(T> T_\alpha))
$$
which is also valid.
It can be easily seen that
$$
E((1-A')\frac{\alpha - P(T >T_\alpha)}{1-P(T > T_\alpha)}) = (\alpha - P(T> T_\alpha)).
$$
More interesting one sees that $\phi_t$ uses a constant randomization probability for all $T_o \leq T_\alpha$ whereas $\phi_t$ adds relatively more to the unrandomized conditional error for outcomes with larger conditional error rates, and less otherwise. Consequently using $A_t'$ makes more sense; using $A_R'$ makes probably most sense.
x <- rep(1,4)
n1 <- 2
n <- 4
ne <- length(x)
perms=1000
alpha=2/16
test_statistic=sumdiff
options(mc.cores=detectCores()-2)
exact <- {perms>=2^n}
obs <- adaperm:::split_sample_os(x,n1,n,ne)
xs <- obs[[1]]
gs <- obs[[2]]
pdist <- adaperm:::perm_dist(xs[[1]],xs[[2]],gs[[1]],gs[[2]],test_statistic,perms,restricted=FALSE)
table(pdist)
cdist <- adaperm:::cond_dist(xs[[1]],xs[[2]],gs[[1]],gs[[2]],test_statistic,perms,restricted=FALSE)
table(cdist)
t_alpha <- max(pdist[rank(-pdist,ties='min') >= alpha*length(pdist)])
t_eps <- max(pdist[rank(-pdist,ties='min') >= alpha*(length(pdist)-1)])
trimmings <- alpha - mean(pdist > t_alpha)
aprime <- mean(cdist>t_alpha)
aR <- aprime + mean(cdist==t_alpha)*trimmings/mean(pdist == t_alpha)
aprime <- permutation_CER(xs[[1]],gs[[1]],xs[[2]],test_statistic,perms,alpha,gs[[2]],restricted=FALSE)
permutation_cer_os <- function(x,n1,n,test_statistic,perms=50000,alpha=0.025,verbose=FALSE,epsilon=0){
exact <- {perms>=2^n}
ne <- length(x)
obs <- adaperm:::split_sample_os(x,n1,n,ne)
xs <- obs[[1]]
gs <- obs[[2]]
pdist <- adaperm:::perm_dist(xs[[1]],xs[[2]],gs[[1]],gs[[2]],test_statistic,perms,restricted=FALSE)
cdist <- adaperm:::cond_dist(xs[[1]],xs[[2]],gs[[1]],gs[[2]],test_statistic,perms,restricted=FALSE)
t_alpha <- max(pdist[rank(-pdist,ties='max') >= alpha*length(pdist)])
if(epsilon==Inf){
t_eps <- -Inf
} else {
t_eps <- sort(unique(pdist[pdist<t_alpha]),decreasing=TRUE)[epsilon]
}
trimmings <- alpha - mean(pdist > t_alpha)
aprime <- mean(cdist>t_alpha)
if(exact&epsilon==0){
## a_r
A <- aprime + mean(cdist==t_alpha)*trimmings/mean(pdist == t_alpha)
} else {
## a_epsilon
A <- aprime + mean(cdist <= t_alpha & cdist >= t_eps)*trimmings/mean(pdist <= t_alpha & pdist >= t_eps)
}
if(verbose){
return(c(A,aprime+trimmings,aprime))
} else {
return(A)
}
}
adist <- function(x,n1,n,test_statistic,perms=50000,alpha=0.025,epsilon=0){
x <- abs(x)
x2 <- x[(n1+1):n]
G <- (omega(rep(1,n1),restricted=F,B=perms)*2-1)
s <- test_statistic(G,rep(1,n1))
rbind(s,apply(G,2,function(i) permutation_cer_os(c(i[1:n1],x2),n1,n,test_statistic,perms,alpha,verbose=TRUE,epsilon=0)))
}
foo <- adist(rep(1,4),2,4,sumdiff,,alpha=1/8,epsilon=0)
apply(foo[-1,],1,function(b) colSums(t(matrix(rep(pnorm(rnorm(10,mean=2),lower=F),4),ncol=4)) <= b))
n <- 2*5
G <- (e1071::bincombinations(n)*2)-1
rowMeans(Adist <- apply(G,1,function(x) permutation_cer_os(x,n/2,n,sumdiff,alpha=.025,verbose=T,epsilon=0)))
n <- 2*8
G <- (e1071::bincombinations(n)*2)-1
system.time(rowMeans(Adist <- simplify2array(mclapply2(1:nrow(G),function(i) permutation_cer_os(G[i,],n/2,n,sumdiff,alpha=.025,verbose=T,epsilon=0,perms=10^6)))))
rowMeans(Adist)
Gnorm <- t(t(G)*abs(rnorm(n)))
rowMeans(Adist_norm <- simplify2array(mclapply2(1:nrow(G),function(i) permutation_cer_os(Gnorm[i,],n/2,n,sumdiff,alpha=.025,verbose=T,epsilon=0,perms=10^6))))
The non-randomized conditional error function looses a lot in size
when applied to discrete outcomes. For continous data, the loss with 8
observations per group is tiny.
Expected/Mid-p values
midp_test <- function(x1,x2,g1,g2,stat,B,x3=NULL,g3=NULL,restricted){
cdist <- perm_dist(x1,x2,g1,g2,stat=test_statistic,B=B,x3=x3,g3=g3,restricted=restricted)
tobs <- test_statistic(c(x1,x2,x3),c(g1,g2,g3))
mean(cdist>tobs) + .5 * mean(cdist == tobs)
}
x <- rep(1,4)
midp_test_os(x,sumdiff)
x1 <- c(1,1)
x2 <- c(1,1)
xE <- c(1,1)
g1 <- c(1,1)
g2 <- c(1,1)
gE <- c(1,-1)
stat <- sumdiff
restricted=F
permutations = 10000
adaptive_midptest_quick <- function(x1,x2,xE,g1,g2,gE,stat,permutations=10000,restricted=T){
## permutation distribution of preplanned test
pdist <- perm_dist(x1,x2,g1,g2,stat,B=permutations,restricted=restricted)
f_pdist <- plyr::count(pdist) %>%
mutate(dens=freq/length(pdist),cdf=1-cumsum(dens))
## conditional distribution of preplanned test given the first stage data
cdist <- adaperm:::cond_dist(x1,x2,g1,g2,stat,B=permutations,restricted=restricted)
f_cdist <- plyr::count(cdist) %>%
mutate(cdens=freq/length(cdist),ccdf=1-cumsum(cdens))
## permutation distribution of adapted test
edist <- perm_dist(x2,xE,g2,gE,stat,B=permutations,restricted=restricted)
f_edist <- plyr::count(edist) %>%
mutate(edens=freq/length(edist),ecdf=1-cumsum(edens))
## observed test statistic of adapted test
t2 <- stat(c(x2,xE),c(g2,gE))
## adaptive p-value
mid_p2 <- mean(edist>t2)+1/2*mean(edist==t2)
df <- plyr::join(f_pdist,f_cdist,by=x,type='left')
## smallest level where one would reject for sure
is <- with(df,min(cdf[!is.na(ccdf) & ccdf>=mid_p2]))
## smallest level where one would not reject anymore
iss <- with(df,max(cdf[!is.na(ccdf) & ccdf<=mid_p2]))
## mid-p value
pval <- (is+iss)/2
pval
}
adaptive_midptest_os <- function (x, n1, n, ne, test_statistic, alpha = 0.025)
{
if(n>25){
stop("Exact distribution not available for n>25 observations")
}
perms = 2^n+1
obs <- adaperm:::split_sample_os(x, n1, n, ne)
if (ne > n) {
xs <- obs[[1]]
gs <- obs[[2]]
return(adaptive_midptest_quick(xs[[1]], xs[[2]], xs[[3]],
gs[[1]], gs[[2]], gs[[3]], test_statistic, restricted = FALSE,
perms))
}
else {
xs <- obs[[1]]
gs <- obs[[2]]
return(midp_test(xs[[1]], xs[[2]], gs[[1]], gs[[2]],
test_statistic, restricted = FALSE, B = perms))
}
}
## test
B <- 10000
X <- matrix(sample(c(-1,1),20*B,rep=T),ncol=B)
Xnorm <- matrix(rnorm(20*B),ncol=B)
pvals_midp <- unlist(mclapply2(1:B,function(i) adaptive_midptest_os(X[,i],n1=8,n=16,ne=20,maritzm,alpha=.05)))
pvals_perm <- unlist(mclapply2(1:B,function(i) adaptive_permtest_os(X[,i],n1=8,n=16,ne=20,maritzm,alpha=.05)))
pvals_midp_norm <- unlist(mclapply2(1:B,function(i) adaptive_midptest_os(Xnorm[,i],n1=8,n=16,ne=20,maritzm,alpha=.05)))
pvals_perm_norm <- unlist(mclapply2(1:B,function(i) adaptive_permtest_os(Xnorm[,i],n1=8,n=16,ne=20,maritzm,alpha=.05,perms=10^6)))
mean(pvals_midp<.05)
mean(pvals_perm<.05)
mean(pvals_midp_norm > pvals_perm_norm)
mean(pvals_midp_norm < pvals_perm_norm)
mean(pvals_midp_norm<.05)
mean(pvals_perm_norm<.05)
out <- adaperm:::split_sample_os(Xnorm[,3],8,16,20)
pdist <- perm_dist(out[[1]][[1]],out[[1]][[2]],out[[2]][[1]],out[[2]][[1]],maritzm,restricted=F,B=10^6)
cdist <- adaperm:::cond_dist(out[[1]][[1]],out[[1]][[2]],out[[2]][[1]],out[[2]][[1]],maritzm,restricted=F,B=10^6)
edist <- perm_dist(out[[1]][[2]],out[[1]][[3]],out[[2]][[2]],out[[2]][[3]],maritzm,restricted=F,B=10^6)
stat1 <- maritzm(c(out[[1]][[1]],out[[1]][[2]]),c(out[[2]][[1]],out[[2]][[2]]))
stat2 <- maritzm(c(out[[1]][[2]],out[[1]][[3]]),c(out[[2]][[2]],out[[2]][[3]]))
c(t2p(stat1,pdist)>t2p(stat1,pdist,'midp'))
c(t2p(stat2,pdist)>t2p(stat2,pdist,'midp'))
permutation_cer_os(Xnorm[,3],8,16,maritzm,perms=10^6)
permutation_CER(out[[1]][[2]],out[[2]][[1]],out[[1]][[2]],stat=maritzm,g2=out[[2]][[2]],restricted=F,perm=10^6)
## reject if
mean(pdist >= quantile(cdist, mean(edist < stat2), type = 1))
mean(edist < stat2)
1-mean(edist>=stat2)
mean(edist>=stat2)
t2p(stat2,edist,'midp')
t2p(stat2,edist)
mean(pdist >= quantile(cdist, mean(edist < stat2),type=1))
mean(pdist > quantile(cdist, mean(edist <= stat2) - .5*mean(edist==stat2),type=1)) +
.5* mean(pdist == quantile(cdist, mean(edist <= stat2) - .5*mean(edist==stat2),type=1))
adaptive_permtest_os(Xnorm[,3],n1=8,n=16,ne=20,maritzm,alpha=.05,perms=10^6)
mean(pdist >= quantile(cdist, mean(edist > stat2) + .5 * mean(edist==stat2), type = 1))
mean(pdist > quantile(cdist, mean(edist > stat2) + .5 * mean(edist==stat2), type = 1))
adaptive_midptest_os(Xnorm[,3],n1=8,n=16,ne=20,maritzm,alpha=.05)
The midp midp value of the adaptive permutation test is still
off. Here is a conjecture
mean(pdist > quantile(cdist, mean(edist <= stat2) - .5*mean(edist==stat2),type=1)) +
.5* mean(pdist == quantile(cdist, mean(edist <= stat2) - .5*mean(edist==stat2),type=1))
Variance of Ta
Most of the time we estimate Ta from a subset of all
permutations. Then we use that as a critical value when estimating the
conditional error. However, we estimate the conditional distribution
using another subset of permutations.
ta <- function(x,B,alpha,stat){
pdist <- stat(x,omega(rep(0:1,length.out=length(x)),B=B))
c(Ta=min(pdist[rank(pdist,ties='min') >= (1-alpha)*B]),Tap=min(pdist[rank(pdist,ties='min') >= (1-alpha+.01)*B]),Tam=min(pdist[rank(pdist,ties='min') >= (1-alpha-.01)*B]),min=min(pdist),max=max(pdist),q1=quantile(pdist,.25),q2=quantile(pdist,.75))
}
x <- rnorm(20)
ta3 <- replicate(1000,ta(x,10^3,.025,ranksum))
ta4 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^4,.025,ranksum)))
ta5 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^5,.025,ranksum)))
sd(ta3[1,])
mean(ta3[2,]-ta3[3,])
sd(ta4[1,])
mean(ta4[2,]-ta4[3,])
sd(ta5[1,])
mean(ta5[2,]-ta5[3,])
ta3 <- replicate(1000,ta(x,10^3,.025,maritzm))
ta4 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^4,.025,maritzm)))
ta5 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^5,.025,maritzm)))
sd(ta3[1,])
mean(ta3[2,]-ta3[3,])
sd(ta4[1,])
mean(ta4[2,]-ta4[3,])
sd(ta5[1,])
mean(ta5[2,]-ta5[3,])
ta3 <- replicate(1000,ta(x,10^3,.025,meandiff))
ta4 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^4,.025,meandiff)))
ta5 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^5,.025,meandiff)))
sd(ta3[1,])
mean(ta3[2,]-ta3[3,])
sd(ta4[1,])
mean(ta4[2,]-ta4[3,])
sd(ta5[1,])
mean(ta5[2,]-ta5[3,])
Conclusion, using $10^4$ ($10^5$) permutations to estimate $T_\alpha$ puts
about 10 (50) standard deviations between $T_{\alpha+1\%}$ and $T_{\alpha-1\%}$.
Is it faster to perform the test or to compute the p-value
library(microbenchmark)
microbenchmark(adaptive_permtest_os(rnorm(20),10,15,20,maritzm),adaptive_permtest2_os(rnorm(20),10,15,20,maritzm))
## hm interestingly the p-value based test is slower by about 30%, I
## don't understand why though
Fix cpp combinations
library(Rcpp)
sourceCpp("../src/combinations.cpp")
library(microbenchmark)
## about 90 fold time improvement
microbenchmark(all_combinations_cpp(16,8),gtools::combinations(16,8))
microbenchmark(all_reassignments_cpp(16,8),gtools::combinations(16,8))
null <- all_combinations(24,12)
## stress test
replicate(1000,{null <- all_combinations(20,10);NULL})
for(n in 1:24){
for(k in 1:n){
if(!all(t(all_combinations(n,k)) == gtools::combinations(n,k))){
stop(paste('error at n =',n,'k =',k))
}
}
}
for(n in 1:24){
for(k in 1:n){
if(!all(all_reassignments_cpp(n,k) == adaperm:::all_reassignments(n,k))){
stop(paste('error at n =',n,'k =',k))
}
}
}
faster cpp permutations
library(Rcpp)
sourceCpp("../src/random_reassignments_cpp.cpp")
library(microbenchmark)
## about 90 fold time improvement
microbenchmark(all_combinations_cpp(16,8),gtools::combinations(16,8))
microbenchmark(all_reassignments_cpp(16,8),gtools::combinations(16,8))
null <- all_combinations(24,12)
## stress test
replicate(1000,{null <- all_combinations(20,10);NULL})
Robust scale estimates
The normal distribution based unbiased estimate of the variance is
unreliable as a basis for sample size reassessment when observations
are not normally distributed. We therefore investigated other
estimates of scale for the use in sample size reassessment.
We consider 3 robust scale estimates
1. The median absolute deviation
2. $S_n$ as suggested in @rousseeuw1993alternatives
3. $Q_n$ as suggested in @rousseeuw1993alternatives
In a first step we look at the bias and standard deviation of the
three estimates for a single sample of size between 2 and 20 as well
as 50 and 100.
library(bt88.03.704)
options(mc.cores=20)
library(parallel)
library(matrixStats)
set.seed('5381')
est <- list()
sd <- list()
sn <- function(x,factor=1.1926){
n <- length(x)
factor*sort(matrixStats::rowOrderStats(abs(matrix(x,n,n) - matrix(x,n,n,byrow=T)),which=floor(n/2)+1))[floor((n+1)/2)]
}
qn <- function(x,factor=2.2219){
n <- length(x)
xm <- abs(matrix(x,n,n) - matrix(x,n,n,byrow=T))
factor*sort(xm[upper.tri(xm)])[choose(floor(n/2)+1,2)]
}
qnp <- function(x,y,factor=2.2219){
n1 <- length(x)
n2 <- length(x)
if(n1 != n2) stop('Unequal sample sizes not yet implemented')
xm <- abs(matrix(x,n1,n1) - matrix(x,n1,n1,byrow=T))
ym <- abs(matrix(y,n2,n2) - matrix(y,n2,n2,byrow=T))
factor*sort(c(xm[upper.tri(xm)],ym[upper.tri(ym)]))[choose(floor(n1/2)+1,2)+choose(floor(n2/2)+1,2)]
}
snp <- function(x,y,factor=1.1926){
n1 <- length(x)
n2 <- length(y)
if(n1 != n2) stop('Unequal sample sizes not yet implemented')
xm <- matrixStats::rowOrderStats(abs(matrix(x,n1,n1) - matrix(x,n1,n1,byrow=T)),which=floor(n1/2)+1)
ym <- matrixStats::rowOrderStats(abs(matrix(y,n2,n2) - matrix(y,n2,n2,byrow=T)),which=floor(n2/2)+1)
factor*sort(c(xm,ym))[n1]
}
gc()
B <- 10^6
square <- list()
for(i in c(2:21,50,51,100,101)){
square[[i]] <- rowMeans(replicate2(B,{
x <- rnorm(i)
y <- rnorm(i)+1
c(qn=qn(x)^2,
qnp=qnp(x,y)^2,
sn=sn(x)^2,
snp=snp(x,y)^2,
sd = var(x),
mad = mad(x)^2,
iqr = (IQR(x)/1.349)^2)}))
}
root <- list()
for(i in c(2:21,50,51,100,101)){
root[[i]] <- rowMeans(replicate2(B,{
x <- rnorm(i)
y <- rnorm(i)+1
c(qn=qn(x),
qnp=qnp(x,y),
sn=sn(x),
snp=snp(x,y),
sd = sd(x),
mad = mad(x),
iqr = (IQR(x)/1.349))}))
}
}
save(sim,file=vfile('scale_estim'))
facs <- do.call(cbind,est)
effs <- do.call(cbind,sd)
save(facs,effs,file=vfile('fin_sample_const'))
## Estimated finite sample efficiencies for qn and sn
load('fin_sample_const_node6_160609.Rd')
efficiencies <- effs[5,]/t(effs[c(1,3,6,7),])
colnames(efficiencies) <- c(expression(Q(n)),expression(S[n]),'MAD','IQR')
pander(efficiencies)
biases <- t(facs[c(1,3,5,6,7),]) - 1
colnames(biases) <- c(expression(Q(n)),expression(S[n]),'VAR','MAD','IQR')
pander(round(biases,3))
## It seems that sn is mor efficient for small sample sizes!
##
sfac <- facs[4,]
qfac <- facs[2,]
plot(2:15,(1+qfac)/qfac)
plot(2:15,(1+sfac)/sfac)
n1 <- 2:15
cn <- ifelse(n1%%2<1,
-1 - 0.276 * n1,
-1.894 - 0.639 * n1)
load('facs_mm_node6_160608.Rd') #facs <- do.call(cbind,est)
sf <- sqrt(facs[2,])
sfo <- sf[1+2*(1:7)]
o <- 1+2*(1:7)
sfe <- sf[2*(1:7)]
e <- 2*(1:7)
lm(sfaco ~ 0 +
summary(mo <- lm(I(sfo/(1-sfo))~o))
summary(me <- lm(I(sfe/(1-sfe))~e))
plot(o,sfo/(1-sfo))
abline(mo)
plot(e,sfe/(1-sfe))
abline(me)
save(facs,file=vfile('facs_mm'))
floatofmath/adaperm documentation built on May 16, 2019, 1:18 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(eval=FALSE,echo=FALSE,warning=FALSE,results='hide',fig.width=7,fig.height=7,message=FALSE)
library(magrittr) library(plyr) library(dplyr) library(microbenchmark) library(pander) library(parallel) library(bt88.03.704) library(adaperm) #devtools::install_github('floatofmath/adaperm') library(ldbounds) options(mc.cores=detectCores()-1)
Adjusted p-values
Let $T_{(1)},...,T_{(|G|)}$ the ordered resampling test statistics of the pre-planned test. Then the $p$-value of the test following an adaptive interim analysis is given by $\tilde{p} = 1-\frac{r^*}{|G|}$ where
$$ r^*=\min { r : p_E \leq A(T_{(r)})} $$
and $A(t)$ is defined as the conditional survival function of the preplanned test given the first stage treatment assignments:
This is easily computed:
sim_pvals <- function(B,delta,nulldist=rnorm,...){ replicate(B,{ x1 <- nulldist(10,...)+rep(0:1,5)*delta x2 <- nulldist(10,...)+rep(0:1,5)*delta xE <- nulldist(6,...)+rep(0:1,3)*delta pt <- t.test(c(x1,x2,xE)~rep(0:1,13),alternative='less',var.equal=T)$p.value pdist <- perm_dist(x1,x2,rep(0:1,5),rep(0:1,5),ranksum,B=10000) pw <- wilcox.test(c(x1,x2,xE)~rep(0:1,13),alternative='less')$p.value cdist <- adaperm:::cond_dist(x1,x2,rep(0:1,5),rep(0:1,5),ranksum,B=10000) edist <- perm_dist(x2,xE,rep(0:1,5),rep(0:1,3),ranksum,B=10000) t2 <- ranksum(c(x2,xE),rep(0:1,8)) c(sum(pdist>=quantile(cdist,sum(edist<t2)/length(edist)))/length(pdist),pt,pw) }) } ## pval_dist <- sim_pvals(1000,1,rnorm_cont,csd=10) ## head(t(pval_dist)) ## cor(t(pval_dist)[,c(1,3)]) ## plot(t(pval_dist)[,c(1,3)]) ## par(mfrow=c(3,1)) ## plot(sort(pval_dist[1,]),seq(0,1,length.out=1000)) ## plot(sort(pval_dist[2,]),seq(0,1,length.out=1000)) ## plot(sort(pval_dist[3,]),seq(0,1,length.out=1000)) ## lookin good! ## rowSums(pval_dist<.05)/1000 ## hist(pval_dist[1,]) adaptive_permtest_ts <- function(x1,x2,xE,g1,g2,gE,stat,permutations=10000){ pdist <- perm_dist(x1,x2,g1,g2,stat,B=permutations) cdist <- adaperm:::cond_dist(x1,x2,g1,g2,stat,B=permutations) edist <- perm_dist(x2,xE,g2,gE,stat,B=permutations) t2 <- stat(c(x2,xE),c(g2,gE)) pval <- sum(pdist>=quantile(cdist,sum(edist<t2)/length(edist)))/length(pdist) names(t2) <- deparse(substitute(stat)) out <- list(statistic = t2, parameter = c('permutations'=permutations), p.value = pval, alternative = "one.sided", null.value = c("distribution of treated samples"='different from control samples'), data.name = paste(deparse(substitute(c(x1,x2,xE,g1,g2,gE))),collapse=', '), method = "Two-sample permutation test for adaptive designs") class(out) <- 'htest' out } str(t.test(rnorm(10))) t.test(rnorm(10)) adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,maritzm) adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,ranksum) str(adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,maritzm))
Robus m-Tests based on maritz' m-estimator
Maritz et al. and Rosenbaum suggest the use of an $m$-estimator
(Huber) for use with a rerandomization test. In adaperm we have
implemented such a test statistic in the function maritzm.
compare_adaptive_tests <- function(n1,n,rule,rdist,test_statistic,mpfr=F,resam=100,perms=100,...){ x <- rdist(n,...) ne <- rule(x[1:n1]) if(ne>n){ x <- c(x,rdist(ne-n,...)) } else { ne <- n } list(ne = ne, permtestT = adaptive_permtest_os(x,n1,n,ne,diffmean,perms=perms), invnormT = adaptive_invnorm_ttest_os(x,n1,n,ne), permtestW = adaptive_permtest_os(x,n1,n,ne,signedranks,perms=perms), invnormW = adaptive_invnorm_wilcoxtest_os(x,n1,n,ne), permtestM = adaptive_permtest_os(x,n1,n,ne,maritzm,perms=perms)) } ##' Compare operating characteristics for a number of one-sample tests ##' ##' \code{test_funs} needs to be a list of functions that take arguments \code{x,n1,n,ne,permutations} and return \code{TRUE} of \code{FALSE} depending on whether the corresponding decision rule rejects the null hypothesis or not. ##' ##' \code{rule} needs to be a function that takes as argument the vector of first stage observations and return an interger number, which (if larger than \code{n}) will be used as the new total sample size. ##' ##' @title Compare one-sample test procedures ##' @param test_funs list of test functions (see Details) ##' @param n1 first stage sample size ##' @param n pre-planned total sample size ##' @param rule sample size reassessment rule (see Details) ##' @param rdist distribution from which to draw observations ##' @param ... additional arguments to rdist ##' @return let's see ##' @author Florian Klinglmueller compare_tests_onesample <- function(B,test_funs,n1,n,rule,rdist,...){ run <- function(test_funs,n1,n,rule,rdist,...){ x <- rdist(n,...) ne <- rule(x[1:n1]) if(ne>n){ x <- c(x,rdist(ne-n,...)) } else { ne <- n } c(n1=n1,n=n,ne=ne,lapply(test_funs,do.call,list(x,n1,n,ne))) } results <- mclapply2(1:B,function(i) run(test_funs,n1,n,rule,rdist,...)) results %>% bind_rows -> results eval(expr=parse(text=paste0("summarize(results,ASN=mean(ne),maxSN=max(ne),sdSN=sd(ne),",paste0(names(tests),"=","mean(",names(tests),")",collapse=','),")"))) } tests <- list( permtestT = function(x,n1,n,ne) adaptive_permtest_os(x,n1,n,ne,diffmean,perms=10000), invnormT = function(x,n1,n,ne) adaptive_invnorm_ttest_os(x,n1,n,ne), permtestW = function(x,n1,n,ne) adaptive_permtest_os(x,n1,n,ne,signedranks,perms=10000), invnormW = function(x,n1,n,ne) adaptive_invnorm_wilcoxtest_os(x,n1,n,ne), permtestM = function(x,n1,n,ne) adaptive_permtest_os(x,n1,n,ne,maritzm,perms=10000)) ## Type I error ## compare_tests_onesample(1000,tests,10,20,cond_power_rule_norm,rnorm_cont,mean=0,sd=1.2,cshift=0,csd=3) ## ## Power ## compare_tests_onesample(1000,tests,10,20,cond_power_rule_norm,rnorm_cont,mean=.5,sd=1,cshift=.5,csd=5) delta <- 0 x1 <- rnorm(10)+rep(0:1,5)*delta x2 <- rnorm(10)+rep(0:1,5)*delta xE <- rnorm(6)+rep(0:1,3)*delta g1 <- sign(x1);g2 <- sign(x2);gE <- sign(xE) pt <- t.test(c(x1,x2,xE),alternative='less')$p.value pdist <- perm_dist(x1,x2,g1,g2,maritzm,B=10000,restricted=F) pw <- wilcox.test(c(x1,x2,xE),alternative='less')$p.value cdist <- adaperm:::cond_dist(x1,x2,g1,g2,maritzm,B=10000,restricted=F) edist <- perm_dist(x2,xE,g2,gE,maritzm,B=10000,restricted=F) t2 <- maritzm(c(x2,xE),rep(0:1,8)) c(sum(pdist>=quantile(cdist,sum(edist<t2)/length(edist)))/length(pdist),pt,pw) adaptive_permtest_os(c(x1,x2,xE),10,20,26,maritzm,perms=10000) permutation_CER(abs(x1),g1>0,abs(x2),maritzm,permutations=10000,restricted=F) perm_test(abs(x2),abs(xE),g2,gE,maritzm,B=10000,restricted=F) par(mfrow=c(3,1)) xl <- min(c(pdist,cdist,edist)) xu <- max(c(pdist,cdist,edist)) hist(pdist,xlim=c(xl,xu)) hist(cdist,xlim=c(xl,xu)) hist(edist,xlim=c(xl,xu)) abline(v=t2)
Point estimates and confidence intervals
##' Permutation test p-value ##' ##' Type "non-randomized" will give the standard p-value, which is define by the proportion of sample permutations with larger or equal test statistic than the observed one. "midp" will give the mid-p value which is the proportion of sample permutations with larger test statistic plus one half the proportion of permutations with equal test statistic. "randomized" will return the proportion of permutations with larger test statistic plus a random number that is uniformly distributed between zero and the proportion of permutations with test statistic equal to the observed. ##' ##' @title Permutation p-value ##' @param t ##' @param dist ##' @param type ##' @return p-value ##' @author Florian Klinglmueller t2p <- function(t,dist,type=c('non-randomized','midp','randomized')){ switch(type[1], 'non-randomized'=mean(dist>=t), 'midp'=mean(dist>t) + .5*mean(dist==t), 'randomized'=mean(dist>t) + runif(1)*mean(dist==t)) } ##' Level alpha critical boundary of permutation test ##' ##' For type 'non-randomized' this returns the smallest unique value of the permutation distribution such that the proportion of statistics larger is smaller than alpha. Type 'randomized' is currently not working but should give the critila value of a randomized test with exact level p. ##' ##' @title p to t ##' @param p p-value/alpha level ##' @param dist permutation distribution ##' @param type see details ##' @return critical value ##' @author Florian Klinglmueller p2t <- function(p,dist,type=c('non-randomized')){#,'randomized')){ n <- length(dist) ta <- max(dist[rank(-dist,ties='min') >= p*length(dist)]) switch(type[1], 'non-randomized'= ta)#, # 'randomized' = ta - (p-floor(p*n)/n)/2 } perm_DR <- function (x1, x2, g1, g2, stat, B, x3 = NULL, g3 = NULL, restricted, type=c('non-randomized','randomized','midp'),alpha) { perm_test(x1,x2,g1,g2,stat,B,x3,g3,restricted,type)<=alpha } perm_test <- function (x1, x2, g1, g2, stat, B, x3 = NULL, g3 = NULL, restricted, type=c('non-randomized','randomized','midp')) { dist <- perm_dist(x1, x2, g1, g2, stat, B, x3, g3, restricted = restricted) t <- stat(c(x1, x2, x3), c(g1, g2, g3)) switch(type[1], 'non-randomized'=t2p(t,dist,type="non-randomized"), 'randomized'= t2p(t,dist,type="randomized"), 'midp'=t2p(t,dist,type="midp")) } adaptive_permtest_quick <- function(x1,x2,xE,g1,g2,gE,stat,permutations=10000,restricted){ pdist <- perm_dist(x1,x2,g1,g2,stat,B=permutations,restricted=restricted) cdist <- adaperm:::cond_dist(x1,x2,g1,g2,stat,B=permutations,restricted=restricted) edist <- perm_dist(x2,xE,g2,gE,stat,B=permutations,restricted=restricted) t2 <- stat(c(x2,xE),c(g2,gE)) pval <- mean(pdist >= quantile(cdist, mean(edist<t2),type=1)) pval } bisect <- function(err,start,end,tolerance,maxit){ l <- err(start) u <- err(end) if(sign(l*u)>0) stop('err() values at end points not of opposite sign') if(maxit<1) { warning('Bisection reached maximum number of iterations') return(err((end-start)/2)) } m <- err((end-start)/2) if(m==0|abs(m) < tolerance) return(start-end) if(sign(u*m)>0){ bisect(err,start,(end-start)/2,tolerance,maxit-1) } else { bisect(err,(end-start)/2,end,tolerance,maxit-1) } } adaptive_permtest_ts <- function(x1,x2,xE,g1,g2,gE,stat,permutations=10000,conf_level=.05){ pval <- adaptive_permtest_quick(x1,x2,xE,g1,g2,gE,stat,permutations=10000,restricted=T) err <- function(d,conf_level=.5){ x1 <- x1-d*g1 x2 <- x2-d*g2 xE <- xE-d*gE adaptive_permtest_quick(x1,x2,xE,g1,g2,gE,stat,permutations=10000,restricted=T)-conf_level } t2 <- stat(c(x2,xE),c(g2,gE)) names(t2) <- deparse(substitute(stat)) ## rms guesstimate <- 2*sqrt(mean(c(x1,x2,xE)^2))*sign(.5-pval) estimate <- c("constant additive effect"=uniroot(err,c(0,guesstimate))$root) conf.int <- c(uniroot(err,c(0,guesstimate*sign(conf_level-pval)),conf_level=conf_level)$root,Inf) attr(conf.int,"conf.level") <- 1-conf_level out <- list(statistic = t2, parameter = c('permutations'=permutations), p.value = pval, estimate = estimate, conf.int=conf.int, alternative = "one.sided", null.value = c("distribution of treated samples"='different from control samples'), data.name = paste(deparse(substitute(c(x1,x2,xE,g1,g2,gE))),collapse=', '), method = "Two-sample permutation test for adaptive designs") class(out) <- 'htest' out } g1 <- rep(0:1,5) g2 <- rep(0:1,5) gE <- rep(0:1,3) x1 <- rnorm(10)+g1 x2 <- rnorm(10)+g2 xE <- rnorm(6)+gE adaptive_permtest_quick(x1,x2,xE,g1,g2,gE,maritzm,perm=10000,restricted=T) t.test(-c(x1,x2,xE)~c(g1,g2,gE),alternative='greater') adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,maritzm,perm=10000) adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,ranksum,perm=10000) adaptive_permtest_ts(x1,x2,xE,g1,g2,gE,meandiff,perm=10000)
Case study
Synthetic data
print(n0 <- power.z.test(delta=6,sd=6,power=.8)) set.seed(5449219) x <- rnorm_cont(16,25,8,cshift=0,csd=24) y <- rnorm_cont(16,19,8,cshift=0,csd=24) boxplot(cbind(x,y),ylim=c(0,56)) ## using the normal distribution conditional power rule suggested overall sample size is 48 cond_power_rule_norm_ts(x[1:8],y[1:8],delta=6,target=.8,rob_var=F) ## using the robust variance estiamte the suggested sample size rule is 32 print(nE <- cond_power_rule_norm_ts(x[1:8],y[1:8],delta=6,target=.8,rob_var=T)) xE <- rnorm_cont(nE-16,25,8,cshift=0,csd=24) yE <- rnorm_cont(nE-16,19,8,cshift=0,csd=24) adaperm(c(y,yE,x,xE),g=rep(0:1,each=nE),8,16,test_statistic=ranksum,cer_type='randomized',atest_type='non-randomized') adaperm(c(y,yE,x,xE),g=rep(0:1,each=nE),8,16,test_statistic=tstat,cer_type='randomized',atest_type='non-randomized') adaptive_invnormtest_2s(c(y,yE),c(x,xE),8,16,nE) ## invnorm estimates and CI dec <- function(delta,alpha){ 2*adaptive_invnormtest_2s(c(y,yE)+delta,c(x,xE),8,16,nE,alpha=alpha)-1 } bisect(-10,10,function(x) dec(x,alpha=.025)) bisect(-10,10,function(x) dec(x,alpha=.5)) bisect(-15,15,function(x) dec(x,alpha=.975)) adaptive_invnorm_wilcoxtest_2s(c(y,yE),c(x,xE),8,16,nE) dec <- function(delta,alpha){ 2*adaptive_invnorm_wilcoxtest_2s(c(y,yE)+delta,c(x,xE),8,16,nE,alpha=alpha)-1 } bisect(-10,10,function(x) dec(x,alpha=.025)) bisect(-10,10,function(x) dec(x,alpha=.5)) bisect(-15,15,function(x) dec(x,alpha=.975)) ## almost identical - need to compute the wilcox.test twice to get matching p-values CIs wilcox.test(x,y,conf.int=T,alternative='greater',correct=FALSE,exact=TRUE,conf.level=0.975) wilcox.test(x,y,conf.int=T,alternative='two.sided',correct=FALSE,exact=TRUE,conf.level=0.95) ## may take a bit due to many permutations adaperm(c(y,x),g=rep(0:1,each=16),8,16,test_statistic=ranksum,permutations=10^5,cer_type='randomized') adaperm(c(y,x),g=rep(0:1,each=16),8,16,test_statistic=ranksum,permutations=10^5,cer_type='non-randomized') delta <- 0 x1 <- c(y[1:8],x[1:8]) x2 <- c(y[9:16],x[9:16]) x3 <- c(yE,xE) g1 <- rep(0:1,each=8) g2 <- rep(0:1,each=8) g3 <- rep(0:1,each=nE-16) pt <- t.test(c(x1,x2,x3)~c(g1,g2,g3),alternative='less')$p.value pdist <- perm_dist(x1,x2,g1,g2,ranksum,B=10000,restricted=T) pw <- wilcox.test(c(x1,x2,x3)~c(g1,g2,g3),alternative='less')$p.value cdist <- adaperm:::cond_dist(x1,x2,g1,g2,ranksum,B=10000,restricted=T) edist <- perm_dist(x2,x3,g2,g3,ranksum,B=10000,restricted=T) t2 <- ranksum(c(x2,x3),c(g2,g3)) c(sum(pdist>=quantile(cdist,sum(edist<t2)/length(edist)))/length(pdist),pt,pw) par(mfrow=c(3,1)) xl <- min(c(pdist,cdist,edist)) xu <- max(c(pdist,cdist,edist)) hist(pdist,xlim=c(xl,xu),main='Permutation distribution of rank-sum statistic') hist(cdist,xlim=c(xl,xu),main='Conditional permutation distribution given first stage treatment assignments') hist(edist,xlim=c(xl,xu),main='Permutation distribution of extended second stage data') abline(v=t2) ## Group sequential stuff alpha <- bounds(c(.5,1),alpha=.025)$exit.pr adaperm(c(y,x),g=rep(0:1,each=16),8,16,alpha0=alpha[1],test_statistic=ranksum,permutations=10^5)
Real data
library(coin) library(gMCP) library(ggplot2) library(nparcomp) data(rotarod) ## not useful ggplot(rotarod,aes(group,time))+geom_boxplot() data(neuropathy) ## heteroskedastic ggplot(neuropathy,aes(group,pain))+geom_boxplot() library(multcomp) data(recovery) head(recovery) ggplot(recovery,aes(blanket,minutes))+geom_boxplot() recovery %>% group_by(blanket) %>% group_size #maybe(b0 vs b3) data(mtept) head(mtept) ggplot(mtept,aes(treatment,E4))+geom_boxplot() mtept %>% group_by(blanket) %>% group_size library(MASS) data(anorexia) head(anorexia) ggplot(anorexia,aes(Treat,Prewt-Postwt))+geom_boxplot() anorexia %>% group_by(Treat) %>% group_size power.t.test(delta=.25,sd=.25,power=.8,sig.level=.025,alternative="one.sided") s <- sample(10000:99999,1) set.seed(s) stage <- sample(rep(1:3,c(9,9,12))) neuropathy$stage <- c(stage[-which(stage==3)[(1:2)]],stage) x1 <- neuropathy[neuropathy$stage==1,'pain'] g1 <- rep(0:1,each=9) x2 <- neuropathy[neuropathy$stage==2,'pain'] g2 <- rep(0:1,each=9) xE <- neuropathy[neuropathy$stage==3,][-(19:20),'pain'] gE <- rep(0:1,each=10) cond_power_rule_t_ts(x1[g1==0],x1[g1==1],delta=1/4,target=.8) 73541::22 92203::27 38478::24 neuropathy %>% group_by(group) %>% group_size() neuropathy %>% group_by(group) %>% mutate(stage = sample
Group sequential procedures
## Obrien Fleming type boundaries ## inputs alpha <- bounds(c(.5,1),alpha=.025)$exit.pr g1 <- rep(0:1,10) g2 <- rep(0:1,10) x1 <- rnorm(20) x2 <- rnorm(20) ## function gsd2_permtest_ts <- function(x1,x2,g1,g2,alpha,stat,permutations){ G <- omega(g1,g2,B= permutations) dist1 <- stat(x1,G[1:length(g1),]) dist2 <- stat(c(x1,x2),G) to1 <- stat(x1,g1) to2 <- stat(c(x1,x2),c(g1,g2)) B <- length(dist1) R1 <- rank(dist1,ties.method="min")>= B*(1-alpha[1]) R2 <- rep(FALSE,B) R2[!R1] <- rank(dist2[!R1],ties.method="min")>= (B*(1-alpha[2]) + alpha[2]*sum(R1)) Ta1 = min(dist1[R1]) Ta2 = min(dist2[R2]) ## Joint rejection probability list(jrp = sum(R1 | R2)/B, Ta1 = Ta1, Ta2 = Ta2, to1 = to1, to2 = to2, decision = (to1>=Ta1|to2>=Ta2)) } gsd2_permtest_ts(x1,x2,g1,x2,alpha,meandiff,10^6)
Adaptive group sequential procedures
## Obrien Fleming type boundaries ## inputs alpha <- bounds(c(.5,1),alpha=.025)$exit.pr g1 <- rep(0:1,10) g2 <- rep(0:1,10) gE <- rep(0:1,5) x1 <- rnorm(20) x2 <- rnorm(20) xE <- rnorm(10) permutations <- 10^4 stat <- maritzm gsd2_permtest_ts <- function(x1,x2,xE,g1,g2,gE,alpha,stat,permutations){ G <- omega(g1,g2,B= permutations) dist1 <- stat(x1,G[1:length(g1),]) dist2 <- stat(c(x1,x2),G) to1 <- stat(x1,g1) to2 <- stat(c(x1,x2),c(g1,g2)) B <- length(dist1) R1 <- rank(dist1,ties.method="min")>= B*(1-alpha[1]) R2 <- rep(FALSE,B) R2[!R1] <- rank(dist2[!R1],ties.method="min")>= (B*(1-alpha[2]) + alpha[2]*sum(R1)) Ta1 = min(dist1[R1]) Ta2 = min(dist2[R2]) jrp = sum(R1 | R2)/B cdist <- adaperm:::cond_dist(x1,x2,g1,g2,stat,permutations) ## reuse the alpha trimmings - this may be controversial if A=0 even though the test would control the error rate A <- sum(cdist>=Ta2)/length(cdist)+max(0,alpha[2]-jrp) GE <- omega(c(g2,gE),B=permutations) ## dirty trick to get non-finite cricital values if condiitonal error is exactly 0 distE <- c(stat(c(x2,xE),GE),Inf) toE <- stat(c(x2,xE),c(g1,gE)) ## The ranks are not changed by adding Inf, but now the vector won't be empty if A = 0 TaE <- min(distE[rank(distE,ties.method='min') > (1-A)*(length(distE)-1)]) ## Joint rejection probability list(a1 = sum(R1)/B, jrp = jrp, Ta1 = Ta1, TaE = TaE, to1 = to1, toE = toE, decision = (to1>=Ta1|toE>=TaE)) } gsd2_permtest_ts(x1+g1,x2+g2,xE+gE,g1,g2,gE,alpha,meandiff,10^5)
Efficiency considerations
different ways to subset some subsetting
select1 <- function(G,sg1){ tG <- t(G) colnames(tG) <- paste0('V',1:ncol(tG)) cond <- paste0('V',1:length(sg1),'==sg1[',1:length(sg1),']',collapse='&') tG[with(data.frame(tG),eval(parse(text=cond))),] } select2 <- function(G,sg1){ tG <- t(G) merge(matrix(sg1,nr=1),tG,by=1:length(sg1)) } select3 <- function(G,sg1){ tG <- as.data.table(t(G)) sg1 <- as.data.table(as.data.frame(matrix(sg1,nr=1))) setkeyv(tG,paste0('V',1:length(sg1))) setkeyv(sg1,paste0('V',1:length(sg1))) tG[sg1,] } select4 <- function(G,sg1){ inner_join(t(G) %>% as.data.frame %>% tbl_df, matrix(sg1,nr=1) %>% as.data.frame %>% tbl_df,by=paste0('V',1:length(sg1))) } select5 <- function(G,sg1){ ## this may be dangerous and does not save time!! suppressMessages(inner_join(t(G) %>% as.data.frame %>% tbl_df, matrix(sg1,nr=1) %>% as.data.frame %>% tbl_df)) } n1=10 g1 <- rep(0:1,n1) g2 <- rep(0:1,n1) G <- omega(g1,g2,B=10^6) dim(G) sg1 <- sample(g1) dim(select1(G,sg1)) microbenchmark(select1(G,sg1), # select2(G,sg1), # select3(G,sg1), select4(G,sg1), select5(G,sg1))
Conditional error of the randomized test
Typically the randomized decision rule of the permutation test is defined as
$$ \phi = \left{\begin{array}{ll} 1 & \mbox{if}\; T_o > T_\alpha \ \gamma & \mbox{if}\; T_o = T_\alpha \ 0 & \mbox{otherwise.} \end{array}\right. $$
where $T_\alpha$ is the $\lceil (1-\alpha) * |G|\rceil$ order statistic (using the minimum rank for ties) of the orbit ${T(x): x \in S(x)}$ and $\gamma$ is defined as
$$ \frac{\alpha - P(T > T_\alpha)}{P(T=T_\alpha)}. $$
The corresponding conditional error function is then
$$ A_R'=A' + \gamma P(T = T_\alpha| x_{|n_1} ). $$
If either the $T_\alpha$ or the conditional permutation distribution given the first stage observations are estimated using Monte-Carlo permutation algorithm, it may be convenient to compute the conditional error function of a slightly different randomized decision rule, which flips a coin for all values of the test statistic falling within $\epsilon$ of the critical value:
$$ \phi_\epsilon = \gamma \mbox{if}\; T_o \in [T_\alpha - \epsilon,T_\alpha). $$
In this case $\gamma$ is given by
$$ \gamma = \frac{\alpha - P(T > T_\alpha)}{P(T \in [T_\alpha-\epsilon,T_\alpha))} $$
which can be easily computed even for a Monte Carlo sample from the permutation space. And the conditional error function is:
$$ A'\epsilon = A' + P(T \in [T\alpha-\epsilon,T_\alpha)|x_{|n_1})\frac{\alpha - P(T > T_\alpha)}{P(T \in [T_\alpha-\epsilon,T_\alpha))} $$
For $\epsilon \rightarrow \infty$ one gets:
$$ A'\infty = A' + (1-A')\frac{\alpha - P(T >T\alpha)}{1-P(T > T_\alpha)} $$
Alternatively one may use the conditional error function:
$$ A_t'= A' + (\alpha - P(T> T_\alpha)) $$
which is also valid.
It can be easily seen that
$$ E((1-A')\frac{\alpha - P(T >T_\alpha)}{1-P(T > T_\alpha)}) = (\alpha - P(T> T_\alpha)). $$
More interesting one sees that $\phi_t$ uses a constant randomization probability for all $T_o \leq T_\alpha$ whereas $\phi_t$ adds relatively more to the unrandomized conditional error for outcomes with larger conditional error rates, and less otherwise. Consequently using $A_t'$ makes more sense; using $A_R'$ makes probably most sense.
x <- rep(1,4) n1 <- 2 n <- 4 ne <- length(x) perms=1000 alpha=2/16 test_statistic=sumdiff options(mc.cores=detectCores()-2) exact <- {perms>=2^n} obs <- adaperm:::split_sample_os(x,n1,n,ne) xs <- obs[[1]] gs <- obs[[2]] pdist <- adaperm:::perm_dist(xs[[1]],xs[[2]],gs[[1]],gs[[2]],test_statistic,perms,restricted=FALSE) table(pdist) cdist <- adaperm:::cond_dist(xs[[1]],xs[[2]],gs[[1]],gs[[2]],test_statistic,perms,restricted=FALSE) table(cdist) t_alpha <- max(pdist[rank(-pdist,ties='min') >= alpha*length(pdist)]) t_eps <- max(pdist[rank(-pdist,ties='min') >= alpha*(length(pdist)-1)]) trimmings <- alpha - mean(pdist > t_alpha) aprime <- mean(cdist>t_alpha) aR <- aprime + mean(cdist==t_alpha)*trimmings/mean(pdist == t_alpha) aprime <- permutation_CER(xs[[1]],gs[[1]],xs[[2]],test_statistic,perms,alpha,gs[[2]],restricted=FALSE) permutation_cer_os <- function(x,n1,n,test_statistic,perms=50000,alpha=0.025,verbose=FALSE,epsilon=0){ exact <- {perms>=2^n} ne <- length(x) obs <- adaperm:::split_sample_os(x,n1,n,ne) xs <- obs[[1]] gs <- obs[[2]] pdist <- adaperm:::perm_dist(xs[[1]],xs[[2]],gs[[1]],gs[[2]],test_statistic,perms,restricted=FALSE) cdist <- adaperm:::cond_dist(xs[[1]],xs[[2]],gs[[1]],gs[[2]],test_statistic,perms,restricted=FALSE) t_alpha <- max(pdist[rank(-pdist,ties='max') >= alpha*length(pdist)]) if(epsilon==Inf){ t_eps <- -Inf } else { t_eps <- sort(unique(pdist[pdist<t_alpha]),decreasing=TRUE)[epsilon] } trimmings <- alpha - mean(pdist > t_alpha) aprime <- mean(cdist>t_alpha) if(exact&epsilon==0){ ## a_r A <- aprime + mean(cdist==t_alpha)*trimmings/mean(pdist == t_alpha) } else { ## a_epsilon A <- aprime + mean(cdist <= t_alpha & cdist >= t_eps)*trimmings/mean(pdist <= t_alpha & pdist >= t_eps) } if(verbose){ return(c(A,aprime+trimmings,aprime)) } else { return(A) } } adist <- function(x,n1,n,test_statistic,perms=50000,alpha=0.025,epsilon=0){ x <- abs(x) x2 <- x[(n1+1):n] G <- (omega(rep(1,n1),restricted=F,B=perms)*2-1) s <- test_statistic(G,rep(1,n1)) rbind(s,apply(G,2,function(i) permutation_cer_os(c(i[1:n1],x2),n1,n,test_statistic,perms,alpha,verbose=TRUE,epsilon=0))) } foo <- adist(rep(1,4),2,4,sumdiff,,alpha=1/8,epsilon=0) apply(foo[-1,],1,function(b) colSums(t(matrix(rep(pnorm(rnorm(10,mean=2),lower=F),4),ncol=4)) <= b)) n <- 2*5 G <- (e1071::bincombinations(n)*2)-1 rowMeans(Adist <- apply(G,1,function(x) permutation_cer_os(x,n/2,n,sumdiff,alpha=.025,verbose=T,epsilon=0))) n <- 2*8 G <- (e1071::bincombinations(n)*2)-1 system.time(rowMeans(Adist <- simplify2array(mclapply2(1:nrow(G),function(i) permutation_cer_os(G[i,],n/2,n,sumdiff,alpha=.025,verbose=T,epsilon=0,perms=10^6))))) rowMeans(Adist) Gnorm <- t(t(G)*abs(rnorm(n))) rowMeans(Adist_norm <- simplify2array(mclapply2(1:nrow(G),function(i) permutation_cer_os(Gnorm[i,],n/2,n,sumdiff,alpha=.025,verbose=T,epsilon=0,perms=10^6))))
The non-randomized conditional error function looses a lot in size when applied to discrete outcomes. For continous data, the loss with 8 observations per group is tiny.
Expected/Mid-p values
midp_test <- function(x1,x2,g1,g2,stat,B,x3=NULL,g3=NULL,restricted){ cdist <- perm_dist(x1,x2,g1,g2,stat=test_statistic,B=B,x3=x3,g3=g3,restricted=restricted) tobs <- test_statistic(c(x1,x2,x3),c(g1,g2,g3)) mean(cdist>tobs) + .5 * mean(cdist == tobs) } x <- rep(1,4) midp_test_os(x,sumdiff) x1 <- c(1,1) x2 <- c(1,1) xE <- c(1,1) g1 <- c(1,1) g2 <- c(1,1) gE <- c(1,-1) stat <- sumdiff restricted=F permutations = 10000 adaptive_midptest_quick <- function(x1,x2,xE,g1,g2,gE,stat,permutations=10000,restricted=T){ ## permutation distribution of preplanned test pdist <- perm_dist(x1,x2,g1,g2,stat,B=permutations,restricted=restricted) f_pdist <- plyr::count(pdist) %>% mutate(dens=freq/length(pdist),cdf=1-cumsum(dens)) ## conditional distribution of preplanned test given the first stage data cdist <- adaperm:::cond_dist(x1,x2,g1,g2,stat,B=permutations,restricted=restricted) f_cdist <- plyr::count(cdist) %>% mutate(cdens=freq/length(cdist),ccdf=1-cumsum(cdens)) ## permutation distribution of adapted test edist <- perm_dist(x2,xE,g2,gE,stat,B=permutations,restricted=restricted) f_edist <- plyr::count(edist) %>% mutate(edens=freq/length(edist),ecdf=1-cumsum(edens)) ## observed test statistic of adapted test t2 <- stat(c(x2,xE),c(g2,gE)) ## adaptive p-value mid_p2 <- mean(edist>t2)+1/2*mean(edist==t2) df <- plyr::join(f_pdist,f_cdist,by=x,type='left') ## smallest level where one would reject for sure is <- with(df,min(cdf[!is.na(ccdf) & ccdf>=mid_p2])) ## smallest level where one would not reject anymore iss <- with(df,max(cdf[!is.na(ccdf) & ccdf<=mid_p2])) ## mid-p value pval <- (is+iss)/2 pval } adaptive_midptest_os <- function (x, n1, n, ne, test_statistic, alpha = 0.025) { if(n>25){ stop("Exact distribution not available for n>25 observations") } perms = 2^n+1 obs <- adaperm:::split_sample_os(x, n1, n, ne) if (ne > n) { xs <- obs[[1]] gs <- obs[[2]] return(adaptive_midptest_quick(xs[[1]], xs[[2]], xs[[3]], gs[[1]], gs[[2]], gs[[3]], test_statistic, restricted = FALSE, perms)) } else { xs <- obs[[1]] gs <- obs[[2]] return(midp_test(xs[[1]], xs[[2]], gs[[1]], gs[[2]], test_statistic, restricted = FALSE, B = perms)) } } ## test B <- 10000 X <- matrix(sample(c(-1,1),20*B,rep=T),ncol=B) Xnorm <- matrix(rnorm(20*B),ncol=B) pvals_midp <- unlist(mclapply2(1:B,function(i) adaptive_midptest_os(X[,i],n1=8,n=16,ne=20,maritzm,alpha=.05))) pvals_perm <- unlist(mclapply2(1:B,function(i) adaptive_permtest_os(X[,i],n1=8,n=16,ne=20,maritzm,alpha=.05))) pvals_midp_norm <- unlist(mclapply2(1:B,function(i) adaptive_midptest_os(Xnorm[,i],n1=8,n=16,ne=20,maritzm,alpha=.05))) pvals_perm_norm <- unlist(mclapply2(1:B,function(i) adaptive_permtest_os(Xnorm[,i],n1=8,n=16,ne=20,maritzm,alpha=.05,perms=10^6))) mean(pvals_midp<.05) mean(pvals_perm<.05) mean(pvals_midp_norm > pvals_perm_norm) mean(pvals_midp_norm < pvals_perm_norm) mean(pvals_midp_norm<.05) mean(pvals_perm_norm<.05) out <- adaperm:::split_sample_os(Xnorm[,3],8,16,20) pdist <- perm_dist(out[[1]][[1]],out[[1]][[2]],out[[2]][[1]],out[[2]][[1]],maritzm,restricted=F,B=10^6) cdist <- adaperm:::cond_dist(out[[1]][[1]],out[[1]][[2]],out[[2]][[1]],out[[2]][[1]],maritzm,restricted=F,B=10^6) edist <- perm_dist(out[[1]][[2]],out[[1]][[3]],out[[2]][[2]],out[[2]][[3]],maritzm,restricted=F,B=10^6) stat1 <- maritzm(c(out[[1]][[1]],out[[1]][[2]]),c(out[[2]][[1]],out[[2]][[2]])) stat2 <- maritzm(c(out[[1]][[2]],out[[1]][[3]]),c(out[[2]][[2]],out[[2]][[3]])) c(t2p(stat1,pdist)>t2p(stat1,pdist,'midp')) c(t2p(stat2,pdist)>t2p(stat2,pdist,'midp')) permutation_cer_os(Xnorm[,3],8,16,maritzm,perms=10^6) permutation_CER(out[[1]][[2]],out[[2]][[1]],out[[1]][[2]],stat=maritzm,g2=out[[2]][[2]],restricted=F,perm=10^6) ## reject if mean(pdist >= quantile(cdist, mean(edist < stat2), type = 1)) mean(edist < stat2) 1-mean(edist>=stat2) mean(edist>=stat2) t2p(stat2,edist,'midp') t2p(stat2,edist) mean(pdist >= quantile(cdist, mean(edist < stat2),type=1)) mean(pdist > quantile(cdist, mean(edist <= stat2) - .5*mean(edist==stat2),type=1)) + .5* mean(pdist == quantile(cdist, mean(edist <= stat2) - .5*mean(edist==stat2),type=1)) adaptive_permtest_os(Xnorm[,3],n1=8,n=16,ne=20,maritzm,alpha=.05,perms=10^6) mean(pdist >= quantile(cdist, mean(edist > stat2) + .5 * mean(edist==stat2), type = 1)) mean(pdist > quantile(cdist, mean(edist > stat2) + .5 * mean(edist==stat2), type = 1)) adaptive_midptest_os(Xnorm[,3],n1=8,n=16,ne=20,maritzm,alpha=.05)
The midp midp value of the adaptive permutation test is still off. Here is a conjecture
mean(pdist > quantile(cdist, mean(edist <= stat2) - .5*mean(edist==stat2),type=1)) +
.5* mean(pdist == quantile(cdist, mean(edist <= stat2) - .5*mean(edist==stat2),type=1))
Variance of Ta
Most of the time we estimate Ta from a subset of all permutations. Then we use that as a critical value when estimating the conditional error. However, we estimate the conditional distribution using another subset of permutations.
ta <- function(x,B,alpha,stat){ pdist <- stat(x,omega(rep(0:1,length.out=length(x)),B=B)) c(Ta=min(pdist[rank(pdist,ties='min') >= (1-alpha)*B]),Tap=min(pdist[rank(pdist,ties='min') >= (1-alpha+.01)*B]),Tam=min(pdist[rank(pdist,ties='min') >= (1-alpha-.01)*B]),min=min(pdist),max=max(pdist),q1=quantile(pdist,.25),q2=quantile(pdist,.75)) } x <- rnorm(20) ta3 <- replicate(1000,ta(x,10^3,.025,ranksum)) ta4 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^4,.025,ranksum))) ta5 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^5,.025,ranksum))) sd(ta3[1,]) mean(ta3[2,]-ta3[3,]) sd(ta4[1,]) mean(ta4[2,]-ta4[3,]) sd(ta5[1,]) mean(ta5[2,]-ta5[3,]) ta3 <- replicate(1000,ta(x,10^3,.025,maritzm)) ta4 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^4,.025,maritzm))) ta5 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^5,.025,maritzm))) sd(ta3[1,]) mean(ta3[2,]-ta3[3,]) sd(ta4[1,]) mean(ta4[2,]-ta4[3,]) sd(ta5[1,]) mean(ta5[2,]-ta5[3,]) ta3 <- replicate(1000,ta(x,10^3,.025,meandiff)) ta4 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^4,.025,meandiff))) ta5 <- simplify2array(mclapply(1:1000,function(i) ta(x,10^5,.025,meandiff))) sd(ta3[1,]) mean(ta3[2,]-ta3[3,]) sd(ta4[1,]) mean(ta4[2,]-ta4[3,]) sd(ta5[1,]) mean(ta5[2,]-ta5[3,])
Conclusion, using $10^4$ ($10^5$) permutations to estimate $T_\alpha$ puts about 10 (50) standard deviations between $T_{\alpha+1\%}$ and $T_{\alpha-1\%}$.
Is it faster to perform the test or to compute the p-value
library(microbenchmark) microbenchmark(adaptive_permtest_os(rnorm(20),10,15,20,maritzm),adaptive_permtest2_os(rnorm(20),10,15,20,maritzm)) ## hm interestingly the p-value based test is slower by about 30%, I ## don't understand why though
Fix cpp combinations
library(Rcpp) sourceCpp("../src/combinations.cpp") library(microbenchmark) ## about 90 fold time improvement microbenchmark(all_combinations_cpp(16,8),gtools::combinations(16,8)) microbenchmark(all_reassignments_cpp(16,8),gtools::combinations(16,8)) null <- all_combinations(24,12) ## stress test replicate(1000,{null <- all_combinations(20,10);NULL}) for(n in 1:24){ for(k in 1:n){ if(!all(t(all_combinations(n,k)) == gtools::combinations(n,k))){ stop(paste('error at n =',n,'k =',k)) } } } for(n in 1:24){ for(k in 1:n){ if(!all(all_reassignments_cpp(n,k) == adaperm:::all_reassignments(n,k))){ stop(paste('error at n =',n,'k =',k)) } } }
faster cpp permutations
library(Rcpp) sourceCpp("../src/random_reassignments_cpp.cpp") library(microbenchmark) ## about 90 fold time improvement microbenchmark(all_combinations_cpp(16,8),gtools::combinations(16,8)) microbenchmark(all_reassignments_cpp(16,8),gtools::combinations(16,8)) null <- all_combinations(24,12) ## stress test replicate(1000,{null <- all_combinations(20,10);NULL})
Robust scale estimates
The normal distribution based unbiased estimate of the variance is unreliable as a basis for sample size reassessment when observations are not normally distributed. We therefore investigated other estimates of scale for the use in sample size reassessment.
We consider 3 robust scale estimates 1. The median absolute deviation 2. $S_n$ as suggested in @rousseeuw1993alternatives 3. $Q_n$ as suggested in @rousseeuw1993alternatives
In a first step we look at the bias and standard deviation of the three estimates for a single sample of size between 2 and 20 as well as 50 and 100.
library(bt88.03.704) options(mc.cores=20) library(parallel) library(matrixStats) set.seed('5381') est <- list() sd <- list() sn <- function(x,factor=1.1926){ n <- length(x) factor*sort(matrixStats::rowOrderStats(abs(matrix(x,n,n) - matrix(x,n,n,byrow=T)),which=floor(n/2)+1))[floor((n+1)/2)] } qn <- function(x,factor=2.2219){ n <- length(x) xm <- abs(matrix(x,n,n) - matrix(x,n,n,byrow=T)) factor*sort(xm[upper.tri(xm)])[choose(floor(n/2)+1,2)] } qnp <- function(x,y,factor=2.2219){ n1 <- length(x) n2 <- length(x) if(n1 != n2) stop('Unequal sample sizes not yet implemented') xm <- abs(matrix(x,n1,n1) - matrix(x,n1,n1,byrow=T)) ym <- abs(matrix(y,n2,n2) - matrix(y,n2,n2,byrow=T)) factor*sort(c(xm[upper.tri(xm)],ym[upper.tri(ym)]))[choose(floor(n1/2)+1,2)+choose(floor(n2/2)+1,2)] } snp <- function(x,y,factor=1.1926){ n1 <- length(x) n2 <- length(y) if(n1 != n2) stop('Unequal sample sizes not yet implemented') xm <- matrixStats::rowOrderStats(abs(matrix(x,n1,n1) - matrix(x,n1,n1,byrow=T)),which=floor(n1/2)+1) ym <- matrixStats::rowOrderStats(abs(matrix(y,n2,n2) - matrix(y,n2,n2,byrow=T)),which=floor(n2/2)+1) factor*sort(c(xm,ym))[n1] } gc() B <- 10^6 square <- list() for(i in c(2:21,50,51,100,101)){ square[[i]] <- rowMeans(replicate2(B,{ x <- rnorm(i) y <- rnorm(i)+1 c(qn=qn(x)^2, qnp=qnp(x,y)^2, sn=sn(x)^2, snp=snp(x,y)^2, sd = var(x), mad = mad(x)^2, iqr = (IQR(x)/1.349)^2)})) } root <- list() for(i in c(2:21,50,51,100,101)){ root[[i]] <- rowMeans(replicate2(B,{ x <- rnorm(i) y <- rnorm(i)+1 c(qn=qn(x), qnp=qnp(x,y), sn=sn(x), snp=snp(x,y), sd = sd(x), mad = mad(x), iqr = (IQR(x)/1.349))})) } } save(sim,file=vfile('scale_estim')) facs <- do.call(cbind,est) effs <- do.call(cbind,sd) save(facs,effs,file=vfile('fin_sample_const')) ## Estimated finite sample efficiencies for qn and sn load('fin_sample_const_node6_160609.Rd') efficiencies <- effs[5,]/t(effs[c(1,3,6,7),]) colnames(efficiencies) <- c(expression(Q(n)),expression(S[n]),'MAD','IQR') pander(efficiencies) biases <- t(facs[c(1,3,5,6,7),]) - 1 colnames(biases) <- c(expression(Q(n)),expression(S[n]),'VAR','MAD','IQR') pander(round(biases,3)) ## It seems that sn is mor efficient for small sample sizes! ## sfac <- facs[4,] qfac <- facs[2,] plot(2:15,(1+qfac)/qfac) plot(2:15,(1+sfac)/sfac) n1 <- 2:15 cn <- ifelse(n1%%2<1, -1 - 0.276 * n1, -1.894 - 0.639 * n1) load('facs_mm_node6_160608.Rd') #facs <- do.call(cbind,est) sf <- sqrt(facs[2,]) sfo <- sf[1+2*(1:7)] o <- 1+2*(1:7) sfe <- sf[2*(1:7)] e <- 2*(1:7) lm(sfaco ~ 0 + summary(mo <- lm(I(sfo/(1-sfo))~o)) summary(me <- lm(I(sfe/(1-sfe))~e)) plot(o,sfo/(1-sfo)) abline(mo) plot(e,sfe/(1-sfe)) abline(me) save(facs,file=vfile('facs_mm'))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.