# This is the empirical example for the paper "Subsampling Inference on Quantile Regression Processes" by V. Chernozhukov and I. Fernandez-Val #NB: these programs are specialized to my empirical example # Updated on 02/14/2014 to remove some commands and packages that have been superseded in R version 3.0.1 ############ table.rq.fn<- function (formula, taus = c(0.25, 0.5, 0.75), method = "fn", data=data, ...) { m <- length(taus) tab <- NULL for (i in 1:m) { fit <- rq(formula, taus[i], method = "fn", data=data) tab <- rbind(tab, coefficients(fit)) } invisible(tab) } table.rq.res<- function (formula, taus = c(0.25, 0.5, 0.75), method = "fn", data=data, ...) { m <- length(taus) tab <- NULL setab<- NULL for (i in 1:m) { fit <- rq(formula, taus[i], method = "fn", data=data) tab <- rbind(tab, coefficients(fit)) setab <- rbind(setab, coef(summary(fit))[,2] ) } return(list(tab = tab, setab = setab)) } ############## computes the inference processes for i-th subset of data ############################## estimates<- function(i,ss, data, taus, formula){ datas<- data[ss[,i], ] if (det(t(as.matrix(datas)) %*% as.matrix(datas))>0.1) {; #sol<- t(table.rq(formula, taus=taus, method="br", data=datas)$a[,,"coefs"]) ; sol<- table.rq.fn(formula=formula, taus=taus, method="fn", data=datas); solqr<- sol[,2] ; solols<- lm(formula=formula,data=datas)$coef[-1] ; inference.process.shift<- solqr-solols[1] ; inference.process.dom<- (-(solqr > 0) + (solqr < 0)) * solqr ; inference.process.noeffect<- solqr; infer<- as.numeric(c( inference.process.shift, inference.process.dom, inference.process.noeffect)); write.table(t(as.matrix(infer)), file="h:\\IQR\\application\\results\\BQR.em.res.txt", sep=" ", append=T, col.names = F, na = "NA"); }; } ############## subsampling function that computes the processes i \leq B bootstrap subsamples ############################ subsample<- function( data, nn=nn, b=b, B, formula) { ss<- matrix(sample(nn, b * B, replace = T), b, B); outpt<- lapply (1:B, FUN=estimates, ss, data, taus, formula=formula) return(outpt) } ########################################################################################################################## # Part II: Implementation library(quantreg); # library(limma); Penn46<- as.data.frame(read.table("h:\\IQR\\application\\data\\Penn46_ascii.txt", header=T )); nn<- dim(Penn46)[1]; formula <- log(duration)~treatment+female+black+ndependents+factor(quarter)+ recall+young+old+durable+lusd #ols results z <- summary(lm(formula, data = Penn46)); ols <- z$coef; vars <- names(z$coef[, 1]); p <- length(z$coefficients[, 1]); Jn <- solve(z$cov.unscaled); taus<- c(3:17)/20; Jt <- length(taus); res <- table.rq.fn(formula, taus = taus, method = "fn", data = Penn46, hs=F ) ; solqr<- res[,2] ; solols<- ols["treatment",1] ; res.to.plot<- table.rq.res ( formula, taus=taus, method="fn", data = Penn46, hs=F ) postscript("h:\\IQR\\application\\results\\QTE.ps",horizontal=F, pointsize=15,width=7.0,height=7.0) par(mfrow=c(1,1)) b<- res.to.plot$tab[,"treatment"] ; b.p<- b + 1.69*res.to.plot$setab[,"treatment"] ; b.m<- b - 1.69*res.to.plot$setab[,"treatment"] ; plot( c(0,1 ),range(c( 0.12, b.m,b.p, -.12 )), xlim =c(0, 1), type="n",xlab="Quantile Index", ylab="Coefficient", main="Quantile Treatment Effect for Unemployment Duration"); polygon(c(taus,rev(taus)),c(b.p,rev(b.m)),density=-3, border=T, col=3); points(taus,b); lines(taus,b); abline(h=0); dev.off() inference.processN.shift <- solqr -solols; inference.processN.noeffect <- solqr; inference.processN.dom <- (-(solqr > 0) + (solqr < 0)) * solqr; # subsampling stage to determine critical value set.seed(runif(1)*100); set.seed(88); B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) B<- 20; b<- 3000; output<- subsample(data=Penn46, nn=nn, b, B, formula=formula) # CENTERED SUBSAMPLING out<- read.table("h:\\IQR\\application\\results\\BQR.em.res.txt"); B <- min(10000, dim(out)[1]) out<- out[1:B,-1] # out<- matrix(unlist(output), B, 3*Jt, byrow=T) out.shift<- out[,1:Jt] out.dom<- out[,(Jt+1): (2*Jt) ] ; out.noeffect<- out[,(2*Jt+1): (3*Jt) ] ; varr.shift<- b*apply(out.shift, 2, var); varr.dom<- b*apply(out.dom, 2, var); varr.noeffect<- b*apply(out.noeffect, 2, var); vb.shift<- sqrt( b*(out.shift -matrix(inference.processN.shift, B, Jt, byrow=T))^2/matrix(varr.shift, B, Jt, byrow=T) ); vb.noeffect<- sqrt( b*(out.noeffect -matrix(inference.processN.noeffect, B, Jt, byrow=T))^2/matrix(varr.noeffect, B, Jt, byrow=T)); vb.dom<- sqrt(b)* (out.dom -matrix(inference.processN.dom, B, Jt, byrow=T))/sqrt(matrix(varr.noeffect, B, Jt, byrow=T)) ; vb.shiftN<- max(sqrt( nn*(inference.processN.shift)^2/varr.shift )); vb.noeffectN<- max(sqrt( nn*(inference.processN.noeffect)^2 / varr.noeffect )); vb.domN<- max(sqrt(nn)* (inference.processN.dom) /sqrt(varr.dom)) ; max0<- function(y) { max(y, 0) } vb.dom2<- apply(vb.dom, c(1,2), max0)^2; # KSM.shift<- apply(vb.shift, 1, mean); KS.shift<- sqrt(nn*mean(inference.processN.shift^2/varr) ); # KSM.noeffect<- apply(vb.noeffect, 1, mean); KS.noeffect<- sqrt(nn*mean(inference.processN.noeffect^2/varr) ); # KSM.dom<- apply(vb.dom2, 1, mean); KS.dom<- sqrt(nn*mean(vN.dom2/varr) ); KS.shift<- apply(vb.shift, 1, max); KS.noeffect<- apply(vb.noeffect, 1, max); KS.dom<- apply(vb.dom, 1, max); KSM.shift<- apply(vb.shift, 1, mean); KSM.noeffect<- apply(vb.noeffect, 1, mean); KSM.dom<- apply(vb.dom, 1, mean); # Unadjusted critical values critical.shift<- quantile(KS.shift, .95); critical.noeffect<- quantile(KS.noeffect, .95); critical.dom<- quantile(KS.dom, .95); # record the KS, critical value, and the decision reject.shift<- ifelse(vb.shiftN>critical.shift, "Reject", "Accept" ) reject.noeffect<- ifelse(vb.noeffectN>critical.noeffect, "Reject", "Accept" ) reject.dom<- ifelse(vb.domN>critical.dom, "Reject", "Accept" ) result.table <- rbind( c(vb.noeffectN, critical.noeffect, reject.noeffect ), c(vb.shiftN, critical.shift, reject.shift ), c(vb.domN, critical.dom, reject.dom ) ) result.table # Adjusted critical values alp <- 0.05; critical.shift<- quantile(KS.shift, .95) + qnorm(alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KS.shift),quantile(KS.shift,.95))[1]); critical.noeffect<- quantile(KS.noeffect, .95) + qnorm(alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KS.noeffect),quantile(KS.noeffect,.95))[1]); critical.dom<- quantile(KS.dom, .95) + qnorm(alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KS.dom),quantile(KS.dom,.95))[1]); # record the KS, critical value, and the decision reject.shift<- ifelse(vb.shiftN>critical.shift, "Reject", "Accept" ) reject.noeffect<- ifelse(vb.noeffectN>critical.noeffect, "Reject", "Accept" ) reject.dom<- ifelse(vb.domN>critical.dom, "Reject", "Accept" ) result.table <- rbind( c(vb.noeffectN, critical.noeffect, reject.noeffect ), c(vb.shiftN, critical.shift, reject.shift ), c(vb.domN, critical.dom, reject.dom ) ) result.table # UNCENTERED (CANONICAL) SUBSAMPLING out<- read.table("h:\\IQR\\application\\results\\BQR.em.res.txt"); B <- min(10000, dim(out)[1]) out<- out[1:B,-1] # out<- matrix(unlist(output), B, 3*Jt, byrow=T) out.shift<- out[,1:Jt] out.dom<- out[,(Jt+1): (2*Jt) ] ; out.noeffect<- out[,(2*Jt+1): (3*Jt) ] ; varr.shift<- b*apply(out.shift, 2, var); varr.dom<- b*apply(out.dom, 2, var); varr.noeffect<- b*apply(out.noeffect, 2, var); vb.shift<- sqrt( b*(out.shift)^2/matrix(varr.shift, B, Jt, byrow=T) ); vb.noeffect<- sqrt( b*(out.noeffect)^2/matrix(varr.noeffect, B, Jt, byrow=T)); vb.dom<- sqrt(b)* (out.dom)/sqrt(matrix(varr.noeffect, B, Jt, byrow=T)) ; vb.shiftN<- max(sqrt( nn*(inference.processN.shift)^2/varr.shift )); vb.noeffectN<- max(sqrt( nn*(inference.processN.noeffect)^2 / varr.noeffect )); vb.domN<- max(sqrt(nn)* (inference.processN.dom) /sqrt(varr.dom)) ; max0<- function(y) { max(y, 0) } vb.dom2<- apply(vb.dom, c(1,2), max0)^2; # KSM.shift<- apply(vb.shift, 1, mean); KS.shift<- sqrt(nn*mean(inference.processN.shift^2/varr) ); # KSM.noeffect<- apply(vb.noeffect, 1, mean); KS.noeffect<- sqrt(nn*mean(inference.processN.noeffect^2/varr) ); # KSM.dom<- apply(vb.dom2, 1, mean); KS.dom<- sqrt(nn*mean(vN.dom2/varr) ); KS.shift<- apply(vb.shift, 1, max); KS.noeffect<- apply(vb.noeffect, 1, max); KS.dom<- apply(vb.dom, 1, max); KSM.shift<- apply(vb.shift, 1, mean); KSM.noeffect<- apply(vb.noeffect, 1, mean); KSM.dom<- apply(vb.dom, 1, mean); # Unadjusted critical values critical.shift<- quantile(KS.shift, .95); critical.noeffect<- quantile(KS.noeffect, .95); critical.dom<- quantile(KS.dom, .95); # record the KS, critical value, and the decision reject.shift<- ifelse(vb.shiftN>critical.shift, "Reject", "Accept" ) reject.noeffect<- ifelse(vb.noeffectN>critical.noeffect, "Reject", "Accept" ) reject.dom<- ifelse(vb.domN>critical.dom, "Reject", "Accept" ) result.table <- rbind( c(vb.noeffectN, critical.noeffect, reject.noeffect ), c(vb.shiftN, critical.shift, reject.shift ), c(vb.domN, critical.dom, reject.dom ) ) result.table # Adjusted critical values alp <- 0.05; critical.shift<- quantile(KS.shift, .95) + qnorm(alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KS.shift),quantile(KS.shift,.95))[1]); critical.noeffect<- quantile(KS.noeffect, .95) + qnorm(alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KS.noeffect),quantile(KS.noeffect,.95))[1]); critical.dom<- quantile(KS.dom, .95) + qnorm(alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KS.dom),quantile(KS.dom,.95))[1]); # record the KS, critical value, and the decision reject.shift<- ifelse(vb.shiftN>critical.shift, "Reject", "Accept" ) reject.noeffect<- ifelse(vb.noeffectN>critical.noeffect, "Reject", "Accept" ) reject.dom<- ifelse(vb.domN>critical.dom, "Reject", "Accept" ) result.table <- rbind( c(vb.noeffectN, critical.noeffect, reject.noeffect ), c(vb.shiftN, critical.shift, reject.shift ), c(vb.domN, critical.dom, reject.dom ) ) result.table