fe.power <- function(d, n1, n2, p1, alpha = 0.05) {
nmx <- n1 + n2 + 1
lgnmx <- lgamma(1:nmx)
power <- 0
zz1 <- .Fortran("ferej",
as.integer(nmx),
as.integer(n1),
as.integer(n2),
as.double(alpha),
fcl = as.integer(matrix(0,nmx,2)),
as.double(lgnmx),
PACKAGE="clinfun")
zz3 <- sapply(d,function(d,nmx,n1,n2,p1,fcl,lgnmx) {
p2 <- p1+d
zz2 <- .Fortran("fepow",
as.integer(nmx),
as.integer(n1),
as.integer(n2),
as.double(p1),
as.double(p2),
as.integer(fcl),
as.double(lgnmx),
power = as.double(double(1)),
PACKAGE="clinfun")
c(p2,zz2$power)
},nmx,n1,n2,p1,zz1$fcl,lgnmx)
zz3 <- t(zz3)
zz3 <- as.data.frame(zz3)
names(zz3) <- c("p2","power")
zz3
}
fe.mdor <- function(ncase,ncontrol,pcontrol,alpha=0.05,power=0.8) {
za <- qnorm(1-alpha/2)
zb <- qnorm(power)
cpc <- ncontrol/ncase
A1 <- (za + zb)^2
B1 <- 1 + 2 * cpc * pcontrol
C1 <- 2 * cpc * pcontrol * (ncase * (1 - pcontrol) - A1 * pcontrol)
if (C1 <= 0) {
warning("Small Sample Sizes and Large pcontrol")
smdor <- 1e8
}
else {
temp <- sqrt((A1 * B1)^2 + 2 * (1 + cpc) * A1 * C1)
r1 <- 1 + (A1 * B1 - temp)/C1
r2 <- 1 + (A1 * B1 + temp)/C1
smdor <- max(r1, r2)
}
omdor <- matrix(0,3,2)
omdor[1,1] <- smdor
omdor[2,1] <- orCPS(ncase,ncontrol,pcontrol,alpha,power)
n1 <- ncontrol
n2 <- ncase
nmx <- n1 + n2 + 1
lgnmx <- lgamma(1:nmx)
zz2 <- .Fortran("femdor",
as.integer(nmx),
as.integer(ncontrol),
as.integer(ncase),
as.double(pcontrol),
as.double(alpha),
as.double(power),
as.integer(matrix(0,nmx,2)),
as.double(lgnmx),
omdor = as.double(omdor),
PACKAGE="clinfun")
xx <- matrix(zz2$omdor,3,2)
dimnames(xx) <- list(c("Schlesselman","CPS","Fisher Exact"),c("Odds Ratio", "Exact Power"))
xx
}
fe.ssize <- function(p1,p2,alpha=0.05,power=0.8,r=1,npm=5,mmax=1000) {
m <- mCPS(p1,p2,alpha,power,r)
if (m <= mmax) {
ossiz <- matrix(0,2,3)
ossiz[1, 1] <- m
ossiz[1, 2] <- m*r
nmx <- ceiling((m+npm)*(1+r)) + 2
lgnmx <- lgamma(1:nmx)
zz2 <- .Fortran("fessiz",
as.integer(nmx),
as.double(p1),
as.double(p2),
as.double(r),
as.double(alpha),
as.double(power),
as.integer(npm),
as.integer(matrix(0,nmx,2)),
as.double(lgnmx),
ossiz = as.double(ossiz),
PACKAGE="clinfun")
xx <- matrix(zz2$ossiz,2,3)
xx[1,1] <- ceiling(xx[1,1])
xx[1,2] <- ceiling(xx[1,2])
xx[2,1] <- ceiling(xx[2,1])
xx[2,2] <- ceiling(xx[2,2])
dimnames(xx) <- list(c("CPS","Fisher Exact"),c("Group 1", "Group 2", "Exact Power"))
} else {
message("\n Exact power is not computed for sample size greater than ", mmax, "\n")
xx <- matrix(0,1,3)
xx[1,1] <- ceiling(m)
xx[1,2] <- ceiling(m*r)
xx[1,3] <- power
dimnames(xx) <- list(c("CPS"),c("Group 1", "Group 2", "Approx Power"))
}
xx
}
mCPS <- function(p1,p2,alpha,power,r) {
za <- qnorm(1-alpha/2)
zb <- qnorm(power)
p <- (p1+r*p2)/(1+r)
m1 <- (za*sqrt((1+r)*p*(1-p)) + zb*sqrt(r*p1*(1-p1)+p2*(1-p2)))^2/(r*(p1-p2)^2)
(m1*(1+sqrt(1+2*(r+1)/(m1*r*abs(p1-p2))))^2)/4
}
CPS.ssize <- function(p1,p2,alpha=0.05,power=0.8,r=1) {
m <- mCPS(p1,p2,alpha=alpha,power=power,r=r)
xx <- matrix(0, 1, 3)
xx[1,1] <- m
xx[1,2] <- m*r
xx[1,1] <- ceiling(xx[1,1])
xx[1,2] <- ceiling(xx[1,2])
xx[1,3] <- power
dimnames(xx) <- list(c("CPS"),c("Group 1", "Group 2", "Approx Power"))
xx
}
orCPS <- function(ncase,ncontrol,pcontrol,alpha=0.05,power=0.8) {
r <- ncase/ncontrol
ff <- function(d, r, ncontrol, pcontrol, alpha, power) {
za <- qnorm(1-alpha/2)
zb <- qnorm(power)
p1 <- pcontrol
p2 <- pcontrol + d
p <- (p1+r*p2)/(1+r)
A <- (za*sqrt((1+r)*p*(1-p)) + zb*sqrt(r*p1*(1-p1)+p2*(1-p2)))^2
m1 <- A/(r*(p1-p2)^2)
(m1*(1+sqrt(1+2*(r+1)/(m1*r*abs(p1-p2))))^2)/4 - ncontrol
}
d0 <- uniroot(ff, lower=0.0001, upper=1-pcontrol-0.0001, r=r, ncontrol=ncontrol, pcontrol=pcontrol, alpha=alpha, power=power)$root
(pcontrol+d0)*(1-pcontrol)/(pcontrol*(1-pcontrol-d0))
}
or2pcase <- function(pcontrol, OR) {OR*pcontrol/(1 + (OR -1)*pcontrol)}