## Inference procedure from a Simon's two-stage design. ## Common inputs are # DES contains the design parameter (n1, nt, r1, rt) # p0 ... P[success] under H0 # p1 ... P[success] under H1 # n1 ... Stage I sample size # r1 ... Stage I critical value # nt ... Total sample size # rt ... Critical value at the end of stage II ## Some functions additionally require # xt ... Total number of responses # alp ... type I error rate of a design (for computing a confidence interval) # p ... any p ## Inference procedure from Simon's two-stage design. ## Some main functions # all(DES, p0) # For a given design (DES) under p0, this computes # pvalues (correct and naive) and confidence interval for each outcome (nt). # find.p(DES) # For a given design (DES), this computes four estimates of p for each outcome (nt). # find.Epall(DES, p) # For a given design (DES) and the assumed true p, this computes four estimates of p and their bias, variance and MSE. ## -- ## -- ## -- ## px <- function(DES, p, xt){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 temp1 <- dbinom(xt, n1, p) pos.x1 <- seq(r1+1, xt) cor.x2 <- xt - pos.x1 temp2 <- sum(dbinom(pos.x1, n1, p) * dbinom(cor.x2, n2, p)) c(temp1, temp2)[1+(xt > r1)]} ## -- ## -- ## -- ## ## -- ## -- ## -- ## px.all <- function(DES, p){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 xtv <- pxt <- seq(0, nt) ; L <- length(xtv) for(i in 1:L){pxt[i] <- px(DES, p, xtv[i])} out <- data.frame(cbind(xtv, pxt)) names(out) <- c('Total X', 'P[X = x]') out} ## -- ## -- ## -- ## ## -- ## -- ## -- ## pval <- function(DES, p0, xt){ temp <- px.all(DES, p0) 1-sum(temp[0:xt, 2])} ## -- ## -- ## -- ## ## -- ## -- ## -- ## pval.all <- function(DES, p0){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 xtv <- pv <- seq(0, nt) ; L <- length(xtv) for(i in 1:L){pv[i] <- pval(DES, p0, xtv[i])} out <- data.frame(cbind(xtv, pv)) names(out) <- c('Total X', 'pvalue') out} ## -- ## -- ## -- ## ## -- ## -- ## -- ## pval.naive <- function(DES, p0){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 xtv <- seq(0, nt) pvn <- 1 - pbinom(xtv-1, n1 + n2* (xtv > r1), p0) cbind(xtv, pvn)} ## -- ## -- ## -- ## ## -- ## -- ## -- ## pval.both <- function(DES, p0, truncate=F){ temp1 <- pval.all(DES, p0) temp2 <- pval.naive(DES, p0) x <- temp1[,1] ; pv <- temp1[,2] ; pvn <- temp2[,2] ok <- rep(T, length(x)) if(truncate){ ok1 <- pmin(pv, pvn) > 1e-5 ok2 <- (1 - pmax(pv, pvn)) > 1e-5 ok <- ok1 & ok2} out <- cbind(x, formatC(cbind(pv, pvn), format='f'))[ok, ] data.frame(out)} ## -- ## -- ## -- ## ## -- ## -- ## -- ## conf <- function(DES, xt, alp=.05){ LB <- UB <- 0 if(xt!=0){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 f.L <- function(p0t, DESt, xtt, alpt){pval(DESt, p0t, xtt) - alpt} f.U <- function(p0t, DESt, xtt, alpt){pval(DESt, p0t, xtt) - (1-alpt)} LB <- uniroot(f.L, c(0, 1), DESt=DES, xtt=xt, alpt=alp)$root UB <- uniroot(f.U, c(0, 1), DESt=DES, xtt=xt, alpt=alp)$root } cbind('lower'=LB, 'upper'=UB)} ## -- ## -- ## -- ## ## -- ## -- ## -- ## conf.all <- function(DES, alp=.05){ p0 <- p01[1] ; p1 <- p01[2] n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 xt <- seq(0, nt) ; L <- nt + 1 conf.m <- matrix(0, ncol=2, nrow=L) for(i in 1:L){conf.m[i,] <- conf(DES, xt[i], alp)} out <- data.frame(conf.m) names(out) <- c('Lower', 'Upper') out} ## -- ## -- ## -- ## ## -- ## -- ## -- ## all <- function(DES, p0, alp=.05){ px <- formatC(px.all(DES, p0)[,2], format='f') temp1 <- pval.both(DES, p0) temp2 <- round(conf.all(DES, alp), 5) temp3 <- data.frame(cbind(temp1[,1], px, temp1[,c(2,3)], temp2)) names(temp3) <- c('xt', 'px', 'correct', 'naive', 'lower', 'upper') temp3} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.phat <- function(DES){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 xtv <- seq(0, nt) phat <- xtv / (n1 + n2* (xtv > r1)) cbind(xtv, phat)} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.phat1 <- function(DES, xt){ temp <- find.phat(DES) temp[,2][temp[,1] == xt]} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.Ep <- function(DES, p){ temp <- find.phat(DES) pxx <- px.all(DES, p) sum(temp[,2] * pxx[,2])} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.pcup1 <- function(DES, xt){ insidepcup <- function(pt, DESt, xtt){ one <- find.Ep(DES, pt) two <- find.phat(DES)[xtt+1, 2] one - two} n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 out <- 0 if( (xt!=0) & (xt<nt)) { out <- uniroot(insidepcup, c(0, 1), DESt = DES, xtt = xt)$root} if( xt>=nt) { out <- 1} out} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.pcup <- function(DES){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 xtv <- pcup <- seq(0, nt) ; L <- length(pcup) for(i in 1:L){pcup[i] <- find.pcup1(DES, xtv[i])} cbind(xtv, pcup)} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.pube1 <- function(DES, xt){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 x1 <- seq( max(r1+1, xt-n2), min(n1, xt)) num <- sum( choose(n1-1, x1-1) * choose(n2, xt - x1) ) den <- sum( choose(n1, x1) * choose(n2, xt - x1) ) num/den} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.pube <- function(DES){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 xtv <- pube <- seq(0, nt) ; L <- length(pube) for(i in 1:L){pube[i] <- find.pube1(DES, xtv[i])} cbind(xtv, pube)} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.pchu1 <- function(DES, xt){ phatv <- find.phat(DES) phat <- phatv[,2][phatv[,1] == xt] fphat <- find.Ep(DES, phat) 2 * phat - fphat} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.pchu <- function(DES){ n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1 xtv <- pchun <- seq(0, nt) ; L <- length(pchun) for(i in 1:L){pchun[i] <- find.pchu1(DES, xtv[i])} cbind(xtv, pchun)} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.p <- function(DES){ phat <- find.phat(DES) xt <- phat[,1] ; phat <- phat[,2] pcup <- find.pcup(DES)[,2] pube <- find.pube(DES)[,2] pchu <- find.pchu(DES)[,2] data.frame(xt, phat, pcup, pube, pchu)} ## -- ## -- ## -- ## ## -- ## -- ## -- ## find.Epall <- function(DES, p){ phat <- find.phat(DES)[,2] pcup <- find.pcup(DES)[,2] pube <- find.pube(DES)[,2] pchu <- find.pchu(DES)[,2] pxx <- px.all(DES, p)[,2] e.phat <- sum( pxx * phat ) e.pcup <- sum( pxx * pcup ) e.pube <- sum( pxx * pube ) e.pchu <- sum( pxx * pchu ) esv <- c(e.phat, e.pcup, e.pube, e.pchu) biasv <- p - esv pbv <- biasv / p * 100 mse.phat <- sum((phat - p)^2 * pxx) mse.pcup <- sum((pcup - p)^2 * pxx) mse.pube <- sum((pube - p)^2 * pxx) mse.pchu <- sum((pchu - p)^2 * pxx) msev <- c(mse.phat, mse.pcup, mse.pube, mse.pchu) var.phat <- mse.phat - biasv[1]^2 var.pcup <- mse.pcup - biasv[2]^2 var.pube <- mse.pube - biasv[3]^2 var.pchu <- mse.pchu - biasv[4]^2 varv <- c(var.phat, var.pcup, var.pube, var.pchu) out <- data.frame(rbind(esv, biasv, pbv, msev, varv), row.names=c('Estimate', 'Bias', 'Bias%', 'MSE', 'VAR')) names(out)=c('mle', 'wh', 'ub', 'chun') out} ## -- ## -- ## -- ## ## -- ## -- ## -- ## p0 <- .2 ; p1 <- .4 n1 <- 24 ; r1 <- 5 ; nt <- 45 ; rt <- 13 x1 <- 8 ; x2 <- 5 ; xt <- x1 + x2 p01 <- c(p0, p1) DES <- c(n1, nt, r1, rt)