Nothing
binding <-
function(l, u, theta, sigma, n1, n2, t.vec, type1, type2, gamma = rep(-4,2),
crange=c(-10,10), drange=c(-10,10), force=TRUE,
plot=TRUE, ll=3, ul=6, n.sim=1e4, seed=NULL)
{
HSD <- function(t,error,gamma) error*(1-exp(-gamma*t))/(1-exp(-gamma))
I.vec <- HSD(t.vec,type1, gamma[1])
II.vec <- HSD(t.vec,type2, gamma[2])
K <- length(t.vec)
nn1 <- ceiling(n1 * t.vec); nn1[K] <- n1
nn2 <- ceiling(n2 * t.vec); nn2[K] <- n2
if(!is.null(seed)) set.seed(seed)
simdataL <- data.frame(matrix(0,n.sim*K,3))
colnames(simdataL)<- c("stage", "t.L", "t.U")
x1 <- matrix(rnorm(n1*n.sim, 0, sigma), n.sim, n1)
x2 <- matrix(rnorm(n2*n.sim, l, sigma), n.sim, n2)
for(k in 1:K){
x1.stage <- x1[,1:nn1[k]]
mean1.temp <- apply(x1.stage, 1, mean)
v1.temp <- apply(x1.stage, 1, var)
x2.stage <- x2[,1:nn2[k]]
mean2.temp <- apply(x2.stage, 1, mean)
v2.temp <- apply(x2.stage, 1, var)
v.pool <- (v1.temp*(nn1[k]-1)+ v2.temp*(nn2[k]-1))/(nn1[k]+ nn2[k]-2)
se.pool <- sqrt(v.pool*(1/nn1[k]+1/nn2[k]))
diff <- mean2.temp - mean1.temp
t.L <- (diff - l) / se.pool
t.U <- (diff - u) / se.pool
simdataL[(k-1)*n.sim + 1:n.sim,] <- cbind(k, t.L, t.U)
}
simdata0 <- data.frame(matrix(0,n.sim*K,3))
colnames(simdata0)<- c("stage", "t.L", "t.U")
x1 <- matrix(rnorm(n1*n.sim, 0, sigma), n.sim, n1)
x2 <- matrix(rnorm(n2*n.sim, theta, sigma), n.sim, n2)
for(k in 1:K){
x1.stage <- x1[,1:nn1[k]]
mean1.temp <- apply(x1.stage, 1, mean)
v1.temp <- apply(x1.stage, 1, var)
x2.stage <- x2[,1:nn2[k]]
mean2.temp <- apply(x2.stage, 1, mean)
v2.temp <- apply(x2.stage, 1, var)
v.pool <- (v1.temp*(nn1[k]-1)+ v2.temp*(nn2[k]-1))/(nn1[k]+ nn2[k]-2)
se.pool <- sqrt(v.pool*(1/nn1[k]+1/nn2[k]))
diff <- mean2.temp - mean1.temp
t.L <- (diff - l) / se.pool
t.U <- (diff - u) / se.pool
simdata0[(k-1)*n.sim + 1:n.sim,] <- cbind(k, t.L, t.U)
}
ct2.L <- rep(99, K); ct2.U <- rep(-99, K)
dt2.L <- rep(999, K); dt2.U <- rep(-999, K)
fL <- function(x) {
temp1 <- simdataL[(simdataL$stage == 1),]$t.L > x
temp2 <- simdataL[(simdataL$stage == 1),]$t.U < -x
result <- sum( temp1 & temp2) / n.sim - I.vec[1]
return(result)
}
ct2.L[1] <- uniroot(fL, crange, tol=0.001)$root
ct2.U[1] <- -ct2.L[1]
f0 <- function(x) {
temp1 <- simdata0[(simdata0$stage == 1),]$t.L <= x
temp2 <- simdata0[(simdata0$stage == 1),]$t.U >= -x
result <- sum( temp1 | temp2) / n.sim - II.vec[1]
return(result)
}
dt2.L[1] <- uniroot(f0, drange, tol=0.001)$root
dt2.U[1] <- -dt2.L[1]
stage1.data <- simdataL[(simdataL$stage == 1),]
temp1 <- (stage1.data$t.L > ct2.L[1])
temp2 <- (stage1.data$t.U < ct2.U[1])
temp3 <- (stage1.data$t.L <= dt2.L[1])
temp4 <- (stage1.data$t.U >= dt2.U[1])
outL <- (temp1 & temp2) | (temp3 | temp4)
stage1.data <- simdata0[(simdata0$stage == 1),]
temp1 <- (stage1.data$t.L > ct2.L[1])
temp2 <- (stage1.data$t.U < ct2.U[1])
temp3 <- (stage1.data$t.L <= dt2.L[1])
temp4 <- (stage1.data$t.U >= dt2.U[1])
out0 <- (temp1 & temp2) | (temp3 | temp4)
for (k in 2:K) {
stagealpha <- I.vec[k] - I.vec[k-1]
stagek.dataL <- simdataL[(simdataL$stage == k),]
if(all(outL)){ct2.L<- ct2.L[1:(k-1)]; ct2.U<- ct2.U[1:(k-1)]; break}
else{
temp <- stagek.dataL[!outL,]
f <- function(x) {
temp1 <- (temp$t.L > x)
temp2 <- (temp$t.U < -x)
result<- sum( temp1 & temp2) / n.sim - stagealpha
return(result)
}
if(nrow(temp) < n.sim*stagealpha)
stop(paste("No solutions for the equivalence boundaries at stage ", as.character(k)))
else{
ct2.L[k] <- uniroot(f, crange, tol=0.001)$root
ct2.U[k] <- -ct2.L[k]
}
}
stagebeta <- II.vec[k] - II.vec[k-1]
stagek.data0 <- simdata0[(simdata0$stage == k),]
if(all(out0)){ dt2.L<- dt2.L[1:(k-1)]; dt2.U<- dt2.U[1:(k-1)]; break}
else{
temp <- stagek.data0[!out0,]
f <- function(x) {
temp1 <- temp$t.L <= x
temp2 <- temp$t.U >= -x
result <- sum( temp1 | temp2) / n.sim- stagebeta
return(result)
}
if(nrow(temp) >= n.sim*stagebeta ){
dt2.L[k] <- uniroot(f, drange, tol=0.001)$root
dt2.U[k] <- -dt2.L[k]
}
}
temp1 <- (stagek.dataL$t.L > ct2.L[k])
temp2 <- (stagek.dataL$t.U < ct2.U[k])
temp3 <- (stagek.dataL$t.L <= dt2.L[k])
temp4 <- (stagek.dataL$t.U >= dt2.U[k])
out <- (temp1 & temp2) | (temp3 | temp4)
outL <- (outL|out)
temp1 <- (stagek.data0$t.L > ct2.L[k])
temp2 <- (stagek.data0$t.U < ct2.U[k])
temp3 <- (stagek.data0$t.L <= dt2.L[k])
temp4 <- (stagek.data0$t.U >= dt2.U[k])
out <- (temp1 & temp2) | (temp3 | temp4)
out0 <- (out0|out)
}
if(force) { dt2.L[K] <- ct2.L[K]; dt2.U[K] <- ct2.U[K]}
result<- list(typeI= I.vec, typeII= II.vec,
equivL=ct2.L,equivU=ct2.U, futilL=dt2.L,futilU=dt2.U)
if(plot){
if(K<=3) par(mfrow=c(1,K), mar=c(4,2,1,1), ask=T)
else if(K==4) par(mfrow=c(2,2), mar=c(4,2,1,1), ask=T)
else if(K>=5 & K<=6) par(mfrow=c(2,3), mar=c(4,2,1,1), ask=T)
figureEF(result, K, ll, ul)
}
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.