You are here:
Vanderbilt Biostatistics Wiki
>
Main Web
>
TWikiUsers
>
TatsukiKoyama
>
TwoStageInference
>
TwoStageInferenceR
(05 Aug 2021,
TatsukiKoyama
)
(raw view)
E
dit
A
ttach
The contents of this page have been moved to [[TwoStageInference][TwoStageInference]]. <!-- <verbatim> ## 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 <- pKoyama <- seq(0, nt) ; L <- length(pKoyama) for(i in 1:L){pKoyama[i] <- find.pchu1(DES, xtv[i])} cbind(xtv, pKoyama)} ## -- ## -- ## -- ## ## -- ## -- ## -- ## 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', 'Koyama') 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) </verbatim> -->
E
dit
|
A
ttach
|
P
rint version
|
H
istory
: r5
<
r4
<
r3
<
r2
|
B
acklinks
|
V
iew topic
|
Edit
w
iki text
|
M
ore topic actions
Topic revision: r5 - 05 Aug 2021,
TatsukiKoyama
Main
Department Home Page
Biostatistics Graduate Program
Vanderbilt University Medical Center
Main Web
Main Web Home
Search
Recent Changes
Changes
Topic list
Biostatistics Webs
Archive
Main
Sandbox
System
Register
|
Log In
Copyright © 2013-2022 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding Vanderbilt Biostatistics Wiki?
Send feedback