# R/powDT.R In DunnettTests: Software implementation of step-down and step-up Dunnett test procedures

#### Documented in powDT

```powDT <- function(r,k,mu,mu0,n,n0,contrast,sigma=NA,df=Inf,alpha=0.05,mcs=1e+05,testcall){
#step 1: calculate the standardized effect size
theta <- mu-mu0
if(contrast=="props" & is.na(sigma)){
pooled.mu <- (mu*n*k+mu0*n0)/(n*k+n0)
sigma <- sqrt(pooled.mu*(1-pooled.mu))
}
corr <- n/(n0+n)
delta <- theta/(sigma*sqrt(1/n+1/n0))

#step 2: calculate the critical values
if(testcall=="SD"){
cvSet <- cvSDDT(k=k, alpha=alpha, alternative="U", df=df, corr=corr)
}
if(testcall=="SU"){
cvSet <- cvSUDT(k=k, alpha=alpha, alternative="U", df=df, corr=corr)
}

#step 3: numerically approximate the power by monte carlo.
z0 <- rnorm(mcs)
if(df == Inf){u <- 1}else{u <- sqrt(rchisq(mcs,df=df)/df)}
dvSet <- sapply(cvSet,function(x){(x*u+sqrt(corr)*z0)/sqrt(1-corr)})
# with columns corresponds to c1,c2,...,ck
list.J.Fun <- list(J2.fun,J3.fun,J4.fun,J5.fun,J6.fun,J7.fun,J8.fun,J9.fun,J10.fun,J11.fun,J12.fun,J13.fun,J14.fun,J15.fun,J16.fun)
if(testcall=="SU"){
#initialize the pow as the probability of rejecting all the k nulls
pow <- mean((1-pnorm(dvSet[,1],mean=delta/sqrt(1-corr)))^k)
if(r<k){
#add the probabilities of rejecting k-s out of k nulls, s=1,...,k-r
s <- 1
p.rej <- (1-pnorm(dvSet[,s+1],mean=delta/sqrt(1-corr)))^(k-s)
Psi.d1 <- pnorm(dvSet[,s],mean=delta/sqrt(1-corr))
J1 <- Psi.d1
list.J <- c(list(1),list(J1))
list.Psi <- list(Psi.d1)
pow <- pow+choose(k,s)*mean(p.rej*J1)
if(r<(k-1)){
for(s in 2:(k-r)){
p.rej <- (1-pnorm(dvSet[,s+1],mean=delta/sqrt(1-corr)))^(k-s)
Psi.ds <- pnorm(dvSet[,s],mean=delta/sqrt(1-corr))
Js <- list.J.Fun[[s-1]](Psi.ds,list.J,list.Psi)
list.J <- c(list.J,list(Js))
list.Psi <- c(list.Psi,list(Psi.ds))
pow <- pow+choose(k,s)*mean(p.rej*Js)
}
}
}
}

if(testcall=="SD"){
powSet <- NULL #with elements correspondng to reject 1, 2,...,k in sequence
Psi.d1 <- 1-pnorm(dvSet[,k],mean=delta/sqrt(1-corr))
J1 <- Psi.d1
list.J <- c(list(1),list(J1))
list.Psi <- list(Psi.d1)
p.acc <- pnorm(dvSet[,k-1],mean=delta/sqrt(1-corr))^(k-1)
powSet <- c(powSet,choose(k,k-1)*mean(p.acc*J1)) #reject only 1
if(k>2){
for(s in (k-2):0){
if(s==0){p.acc <- 1}else{p.acc <- pnorm(dvSet[,s],mean=delta/sqrt(1-corr))^s}
Psi.ds <- 1-pnorm(dvSet[,s+1],mean=delta/sqrt(1-corr))
Js <- list.J.Fun[[k-s-1]](Psi.ds,list.J,list.Psi)
powSet <- c(powSet,choose(k,s)*mean(p.acc*Js)) #reject k-s
list.J <- c(list.J,list(Js))
list.Psi <- c(list.Psi,list(Psi.ds))
}
}
pow <- sum(powSet[r:k])
}
pow
}
```

## Try the DunnettTests package in your browser

Any scripts or data that you put into this service are public.

DunnettTests documentation built on May 2, 2019, 9:13 a.m.