print.webpower <- function(x, ...) {
cat(x$method, "\n\n", sep = "")
note <- x$note
url <- x$url
x[c("method", "note", "url", "alternative", "family", "parameter")] <- NULL
x <- do.call(cbind, x)
rownames(x) <- rep(" ", nrow(x))
print(x)
if (!is.null(note))
cat("\nNOTE: ", note, sep = "")
cat("\nURL: ", url, sep = "")
cat("\n")
invisible(x)
}
plot.webpower <- function(x, xvar = NULL, yvar = NULL, xlab = NULL, ylab = NULL,
...) {
wp <- x
wp[c("method", "note", "url", "alternative", "family", "parameter")] <- NULL
wp.matrix <- do.call(cbind, wp)
check <- is.null(xvar) + is.null(yvar)
if (check == 2) {
if (!is.null(xlab))
xlab <- xlab else xlab <- "Sample size"
if (!is.null(ylab))
ylab <- ylab else ylab <- "Power"
plot(wp$n, wp$power, xlab = xlab, ylab = ylab, ...)
lines(wp$n, wp$power, ...)
} else if (check == 1) {
stop("Both x and y need to be specified!")
} else {
if (!is.null(xlab))
xlab <- xlab else xlab <- xvar
if (!is.null(ylab))
ylab <- ylab else ylab <- yvar
plot(wp.matrix[, xvar], wp.matrix[, yvar], xlab = xlab, ylab = ylab,
...)
lines(wp.matrix[, xvar], wp.matrix[, yvar], ...)
}
}
wp.anova <- function(k = NULL, n = NULL, f = NULL, alpha = 0.05, power = NULL,
type = c("overall", "two.sided", "greater", "less")) {
type <- type[1]
if (sum(sapply(list(k, n, f, power, alpha), is.null)) != 1)
stop("exactly one of k, n, f, power, and alpha must be NULL")
if (!is.null(f) && min(f) < 0)
stop("f must be positive")
if (!is.null(k) && min(k) < 2)
stop("number of groups must be at least 2")
if (!is.null(n) && min(n) < 2)
stop("number of observations in each group must be at least 2")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
if (type == "overall")
p.body <- quote({
lambda <- n * f^2
pf(qf(alpha, k - 1, n - k, lower = FALSE), k - 1, n - k, lambda,
lower = FALSE)
})
if (type == "two.sided")
p.body <- quote({
lambda <- n * f^2
pf(qf(alpha, 1, n - k, lower = FALSE), 1, n - k, lambda, lower = FALSE)
})
if (type == "greater")
p.body <- quote({
lambda <- sqrt(n) * f
pt(qt(alpha, n - k, lower = FALSE), n - k, lambda, lower = FALSE)
})
if (type == "less")
p.body <- quote({
lambda <- sqrt(n) * f
pt(qt(alpha, n - k), n - k, lambda)
})
if (is.null(power))
power <- eval(p.body) else if (is.null(k))
k <- uniroot(function(k) eval(p.body) - power, c(2 + 1e-10, 100))$root else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + k + 1e-10,
1e+05))$root else if (is.null(f))
f <- uniroot(function(f) eval(p.body) - power, c(1e-07, 1e+07))$root else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10,
1 - 1e-10))$root else stop("internal error")
NOTE <- "n is the total sample size"
if (type == "overall") {
NOTE <- paste(NOTE, " (overall)", sep = "")
} else {
NOTE <- paste(NOTE, " (contrast, ", type, ")", sep = "")
}
METHOD <- "Power for One-way ANOVA"
URL <- "http://psychstat.org/anova"
structure(list(k = k, n = n, f = f, alpha = alpha, power = power, note = NOTE,
method = METHOD, url = URL), class = "webpower")
}
wp.prop <- function(h = NULL, n1 = NULL, n2 = NULL, alpha = 0.05, power = NULL,
type = c("1p", "2p", "2p2n"), alternative = c("two.sided", "less",
"greater")) {
type <- type[1]
if (type == "1p")
return(wp.one.prop(h = h, n = n1, alpha = alpha, power = power,
alternative = alternative))
if (type == "2p")
return(wp.two.prop(h = h, n = n1, alpha = alpha, power = power,
alternative = alternative))
if (type == "2p2n")
return(wp.two.prop.two.n(h = h, n1 = n1, n2 = n2, alpha = alpha,
power = power, alternative = alternative))
}
wp.one.prop <- function(h = NULL, n = NULL, alpha = 0.05, power = NULL,
alternative = c("two.sided", "less", "greater")) {
if (sum(sapply(list(h, n, power, alpha), is.null)) != 1)
stop("exactly one of h, n, power, and alpha must be NULL")
if (!is.null(n) && min(n) < 1)
stop("number of observations in each group must be at least 1")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
alternative <- match.arg(alternative)
tside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
if (tside == 2 && !is.null(h))
h <- abs(h)
if (tside == 2) {
p.body <- quote({
pnorm(qnorm(alpha/2, lower = FALSE) - h * sqrt(n), lower = FALSE) +
pnorm(qnorm(alpha/2, lower = TRUE) - h * sqrt(n), lower = TRUE)
})
}
if (tside == 3) {
p.body <- quote({
pnorm(qnorm(alpha, lower = FALSE) - h * sqrt(n), lower = FALSE)
})
}
if (tside == 1) {
p.body <- quote({
pnorm(qnorm(alpha, lower = TRUE) - h * sqrt(n), lower = TRUE)
})
}
if (is.null(power))
power <- eval(p.body) else if (is.null(h)) {
if (tside == 2) {
h <- uniroot(function(h) eval(p.body) - power, c(1e-10, 10))$root
}
if (tside == 1) {
h <- uniroot(function(h) eval(p.body) - power, c(-10, 5))$root
}
if (tside == 3) {
h <- uniroot(function(h) eval(p.body) - power, c(-5, 10))$root
}
} else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + 1e-10, 1e+05))$root else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10,
1 - 1e-10))$root else stop("internal error")
METHOD <- "Power for one-sample proportion test"
NOTE <- NULL
URL <- "http://psychstat.org/prop"
structure(list(h = h, n = n, alpha = alpha, power = power, url = URL,
method = METHOD, note = NOTE), class = "webpower")
}
wp.two.prop <- function(h = NULL, n = NULL, alpha = 0.05, power = NULL,
alternative = c("two.sided", "less", "greater")) {
if (sum(sapply(list(h, n, power, alpha), is.null)) != 1)
stop("exactly one of h, n, power, and alpha must be NULL")
if (!is.null(n) && min(n) < 1)
stop("number of observations in each group must be at least 1")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
alternative <- match.arg(alternative)
tside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
if (tside == 2 && !is.null(h))
h <- abs(h)
if (tside == 3) {
p.body <- quote({
pnorm(qnorm(alpha, lower = FALSE) - h * sqrt(n/2), lower = FALSE)
})
}
if (tside == 2) {
p.body <- quote({
pnorm(qnorm(alpha/2, lower = FALSE) - h * sqrt(n/2), lower = FALSE) +
pnorm(qnorm(alpha/2, lower = TRUE) - h * sqrt(n/2), lower = TRUE)
})
}
if (tside == 1) {
p.body <- quote({
pnorm(qnorm(alpha, lower = TRUE) - h * sqrt(n/2), lower = TRUE)
})
}
if (is.null(power))
power <- eval(p.body) else if (is.null(h)) {
if (tside == 2) {
h <- uniroot(function(h) eval(p.body) - power, c(1e-10, 10))$root
}
if (tside == 1) {
h <- uniroot(function(h) eval(p.body) - power, c(-10, 5))$root
}
if (tside == 3) {
h <- uniroot(function(h) eval(p.body) - power, c(-5, 10))$root
}
} else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + 1e-10, 1e+05))$root else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10,
1 - 1e-10))$root else stop("internal error")
NOTE <- "Sample sizes for EACH group"
URL <- "http://psychstat.org/prop2p"
METHOD <- "Power for two-sample proportion (equal n)"
structure(list(h = h, n = n, alpha = alpha, power = power, url = URL,
method = METHOD, note = NOTE), class = "webpower")
}
wp.two.prop.two.n <- function(h = NULL, n1 = NULL, n2 = NULL, alpha = 0.05,
power = NULL, alternative = c("two.sided", "less", "greater")) {
if (sum(sapply(list(h, n1, n2, power, alpha), is.null)) != 1)
stop("exactly one of h, n1, n2, power, and alpha must be NULL")
if (!is.null(n1) && min(n1) < 2)
stop("number of observations in the first group must be at least 2")
if (!is.null(n2) && min(n2) < 2)
stop("number of observations in the second group must be at least 2")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
alternative <- match.arg(alternative)
tside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
if (tside == 2 && !is.null(h))
h <- abs(h)
if (tside == 3) {
p.body <- quote({
pnorm(qnorm(alpha, lower = FALSE) - h * sqrt((n1 * n2)/(n1 +
n2)), lower = FALSE)
})
}
if (tside == 1) {
p.body <- quote({
pnorm(qnorm(alpha, lower = TRUE) - h * sqrt((n1 * n2)/(n1 +
n2)), lower = TRUE)
})
}
if (tside == 2) {
p.body <- quote({
pnorm(qnorm(alpha/2, lower = FALSE) - h * sqrt((n1 * n2)/(n1 +
n2)), lower = FALSE) + pnorm(qnorm(alpha/2, lower = TRUE) -
h * sqrt((n1 * n2)/(n1 + n2)), lower = TRUE)
})
}
if (is.null(power))
power <- eval(p.body) else if (is.null(h)) {
if (tside == 2) {
h <- uniroot(function(h) eval(p.body) - power, c(1e-10, 10))$root
}
if (tside == 1) {
h <- uniroot(function(h) eval(p.body) - power, c(-10, 5))$root
}
if (tside == 3) {
h <- uniroot(function(h) eval(p.body) - power, c(-5, 10))$root
}
} else if (is.null(n1))
n1 <- uniroot(function(n1) eval(p.body) - power, c(2 + 1e-10, 1e+05))$root else if (is.null(n2))
n2 <- uniroot(function(n2) eval(p.body) - power, c(2 + 1e-10, 1e+05))$root else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10,
1 - 1e-10))$root else stop("internal error")
NOTE <- "Sample size for each group"
URL <- "http://psychstat.org/prop2p2n"
METHOD <- "Power for two-sample proportion (unequal n)"
structure(list(h = h, n1 = n1, n2 = n2, alpha = alpha, power = power,
url = URL, method = METHOD, note = NOTE), class = "webpower")
}
wp.t <- function(n1 = NULL, n2 = NULL, d = NULL, alpha = 0.05, power = NULL,
type = c("two.sample", "one.sample", "paired", "two.sample.2n"), alternative = c("two.sided",
"less", "greater"), tol = .Machine$double.eps^0.25) {
type <- type[1]
if (type == "two.sample.2n") {
return(wp.t2(n1 = n1, n2 = n2, d = d, alpha = alpha, power = power,
alternative = alternative, tol = tol))
} else {
return(wp.t1(n = n1, d = d, alpha = alpha, power = power, type = type,
alternative = alternative, tol = tol))
}
}
wp.t1 <- function(n = NULL, d = NULL, alpha = 0.05, power = NULL, type = c("two.sample",
"one.sample", "paired"), alternative = c("two.sided", "less", "greater"),
tol = .Machine$double.eps^0.25) {
if (sum(sapply(list(n, d, power, alpha), is.null)) != 1)
stop("exactly one of n, d, power, and alpha must be NULL")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
type <- match.arg(type)
alternative <- match.arg(alternative)
tsample <- switch(type, one.sample = 1, two.sample = 2, paired = 1)
ttside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
tside <- switch(alternative, less = 1, two.sided = 2, greater = 1)
if (tside == 2 && !is.null(d))
d <- abs(d)
if (ttside == 1) {
p.body <- quote({
nu <- (n - 1) * tsample
pt(qt(alpha/tside, nu, lower = TRUE), nu, ncp = sqrt(n/tsample) *
d, lower = TRUE)
})
}
if (ttside == 2) {
p.body <- quote({
nu <- (n - 1) * tsample
qu <- qt(alpha/tside, nu, lower = FALSE)
pt(qu, nu, ncp = sqrt(n/tsample) * d, lower = FALSE) + pt(-qu,
nu, ncp = sqrt(n/tsample) * d, lower = TRUE)
})
}
if (ttside == 3) {
p.body <- quote({
nu <- (n - 1) * tsample
pt(qt(alpha/tside, nu, lower = FALSE), nu, ncp = sqrt(n/tsample) *
d, lower = FALSE)
})
}
if (is.null(power))
power <- eval(p.body) else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + 1e-10, 1e+07),
tol = tol, extendInt = "upX")$root else if (is.null(d)) {
if (ttside == 2) {
d <- uniroot(function(d) eval(p.body) - power, c(1e-07, 10),
tol = tol, extendInt = "upX")$root
}
if (ttside == 1) {
d <- uniroot(function(d) eval(p.body) - power, c(-10, 5), tol = tol,
extendInt = "upX")$root
}
if (ttside == 3) {
d <- uniroot(function(d) eval(p.body) - power, c(-5, 10), tol = tol,
extendInt = "upX")$root
}
} else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10,
1 - 1e-10), tol = tol, extendInt = "upX")$root else stop("internal error", domain = NA)
NOTE <- switch(type, paired = "n is number of *pairs*", two.sample = "n is number in *each* group",
NULL)
METHOD <- paste(switch(type, one.sample = "One-sample", two.sample = "Two-sample",
paired = "Paired"), "t-test")
URL <- "http://psychstat.org/ttest"
structure(list(n = n, d = d, alpha = alpha, power = power, alternative = alternative,
note = NOTE, method = METHOD, url = URL), class = "webpower")
}
wp.t2 <- function(n1 = NULL, n2 = NULL, d = NULL, alpha = 0.05, power = NULL,
alternative = c("two.sided", "less", "greater"), tol = .Machine$double.eps^0.25) {
if (sum(sapply(list(n1, n2, d, power, alpha), is.null)) != 1)
stop("exactly one of n1, n2, d, power, and alpha must be NULL")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
if (!is.null(n1) && min(n1) < 2)
stop("number of observations in the first group must be at least 2")
if (!is.null(n2) && min(n2) < 2)
stop("number of observations in the second group must be at least 2")
alternative <- match.arg(alternative)
tsample <- 2
ttside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
tside <- switch(alternative, less = 1, two.sided = 2, greater = 1)
if (tside == 2 && !is.null(d))
d <- abs(d)
if (ttside == 1) {
p.body <- quote({
nu <- n1 + n2 - 2
pt(qt(alpha/tside, nu, lower = FALSE), nu, ncp = d * (1/sqrt(1/n1 +
1/n2)), lower = FALSE)
})
}
if (ttside == 2) {
p.body <- quote({
nu <- n1 + n2 - 2
qu <- qt(alpha/tside, nu, lower = FALSE)
pt(qu, nu, ncp = d * (1/sqrt(1/n1 + 1/n2)), lower = FALSE) +
pt(-qu, nu, ncp = d * (1/sqrt(1/n1 + 1/n2)), lower = TRUE)
})
}
if (ttside == 3) {
p.body <- quote({
nu <- n1 + n2 - 2
pt(qt(alpha/tside, nu, lower = FALSE), nu, ncp = d * (1/sqrt(1/n1 +
1/n2)), lower = FALSE)
})
}
if (is.null(power))
power <- eval(p.body) else if (is.null(n1))
n1 <- uniroot(function(n1) eval(p.body) - power, c(2 + 1e-10, 1e+07),
tol = tol, extendInt = "upX")$root else if (is.null(n2))
n2 <- uniroot(function(n2) eval(p.body) - power, c(2 + 1e-10, 1e+07),
tol = tol, extendInt = "upX")$root else if (is.null(d)) {
if (ttside == 2) {
d <- uniroot(function(d) eval(p.body) - power, c(1e-07, 10),
tol = tol, extendInt = "upX")$root
}
if (ttside == 1) {
d <- uniroot(function(d) eval(p.body) - power, c(-10, 5), tol = tol,
extendInt = "upX")$root
}
if (ttside == 3) {
d <- uniroot(function(d) eval(p.body) - power, c(-5, 10), tol = tol,
extendInt = "upX")$root
}
} else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10,
1 - 1e-10), tol = tol, extendInt = "upX")$root else stop("internal error", domain = NA)
NOTE <- "n1 and n2 are number in *each* group"
METHOD <- "Unbalanced two-sample t-test"
URL <- "http://psychstat.org/ttest2n"
structure(list(n1 = n1, n2 = n2, d = d, alpha = alpha, power = power,
alternative = alternative, note = NOTE, method = METHOD, url = URL),
class = "webpower")
}
## function for repeated measures anova analysis
wp.rmanova <- function(n = NULL, ng = NULL, nm = NULL, f = NULL, nscor = 1,
alpha = 0.05, power = NULL, type = 0) {
if (sum(sapply(list(n, ng, nm, f, nscor, power, alpha), is.null)) !=
1)
stop("exactly one of n, ng, nm, power, and alpha must be NULL")
if (!is.null(f) && min(f) < 0)
stop("Effect size must be positive")
if (!is.null(n) && min(n) < 1)
stop("Sample size has to be larger than 1")
if (!is.null(ng) && min(ng) < 1)
stop("Number of groups have to be at least 1")
if (!is.null(nm) && min(nm) < 2)
stop("number of measurements must be at least 1")
if (!is.null(nscor) && !is.numeric(nscor) || any(0 > nscor | nscor >
1))
stop("Nonsphericity correction must be numeric in [0, 1]")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
## function to evaluate
if (type == 0) {
p.body <- quote({
df1 <- ng - 1
df2 <- n - ng
lambda <- f^2 * n * nscor
pf(qf(alpha, df1, df2, lower = FALSE), df1, df2, lambda, lower = FALSE)
})
}
if (type == 1) {
p.body <- quote({
df1 <- (nm - 1) * nscor
df2 <- (n - ng) * df1
lambda <- f^2 * n * nscor
pf(qf(alpha, df1, df2, lower = FALSE), df1, df2, lambda, lower = FALSE)
})
}
if (type == 2) {
p.body <- quote({
df1 <- (ng - 1) * (nm - 1) * nscor
df2 <- (n - ng) * (nm - 1) * nscor
lambda <- f^2 * n * nscor
pf(qf(alpha, df1, df2, lower = FALSE), df1, df2, lambda, lower = FALSE)
})
}
if (is.null(power))
power <- eval(p.body) else if (is.null(n)) {
n <- uniroot(function(n) eval(p.body) - power, c(5 + ng, 1e+07))$root
} else if (is.null(ng)) {
ndf <- uniroot(function(ng) eval(p.body) - power, c(1 + 1e-10,
1e+05))$root
} else if (is.null(nm)) {
ng <- uniroot(function(nm) eval(p.body) - power, c(1e-07, 1e+07))$root
} else if (is.null(f)) {
f <- uniroot(function(f) eval(p.body) - power, c(1e-07, 1e+07))$root
} else if (is.null(alpha)) {
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10,
1 - 1e-10))$root
} else stop("internal error")
if (type == 0)
NOTE <- "Power analysis for between-effect test"
if (type == 1)
NOTE <- "Power analysis for within-effect test"
if (type == 2)
NOTE <- "Power analysis for interaction-effect test"
URL <- "http://psychstat.org/rmanova"
METHOD <- "Repeated-measures ANOVA analysis"
structure(list(n = n, f = f, ng = ng, nm = nm, nscor = nscor, alpha = alpha,
power = power, note = NOTE, method = METHOD, url = URL), class = "webpower")
}
## Effect size and ICC calculation for CRT with 2 arms based on
## empirical data The data has to be in text format where the first
## column of the data is the ID variable, the second column represents
## cluster, the third column is the outcome variable, and the fourth
## column is the condition variable (0 for control, 1 for condition).
## The first line of the data should be the variable names.
wp.effect.CRT2arm <- function(file) {
if (is.character(file)){
dat <- read.table(file, header = TRUE)
}else{
dat <- file
}
J <- length(unique(dat[, 2]))
n <- nrow(dat)/J
dat.trt <- dat[dat[, 4] == 1, ]
dat.ctl <- dat[dat[, 4] == 0, ]
mutj <- aggregate(dat.trt[, 3], by = list(dat.trt[, 2]), FUN = mean)
mucj <- aggregate(dat.ctl[, 3], by = list(dat.ctl[, 2]), FUN = mean)
mu_j <- rbind(mutj, mucj)
names(mu_j) <- c(names(dat)[2], "mu")
dat.mu <- merge(dat, mu_j, by = names(dat)[2])
SSW <- sum((dat.mu[, 3] - dat.mu[5])^2)/(n * J - J)
mut <- mean(mutj[, 2])
muc <- mean(mucj[, 2])
SSB <- (sum((mutj[, 2] - mut)^2) + sum((mucj[, 2] - muc)^2))/(J - 2)
f <- (mut - muc)/sqrt(SSB + (1 - 1/n) * SSW)
rho <- (SSB - SSW/n)/(SSB + (1 - 1/n) * SSW)
return(list(f = f, ICC = rho))
}
## Effect size and ICC calculation for CRT with 3 arms based on
## empirical data The data has to be in text format, where the first
## column of the data is the ID variable, the second column represents
## cluster, the third column is the outcome variable, and the fourth
## column is the condition variable (0 for control, 1 for treatment1, 2
## for treatment2). The first line of the data should be the variable
## names.
wp.effect.CRT3arm <- function(file) {
if (is.character(file)){
dat <- read.table(file, header = TRUE)
}else{
dat <- file
}
J <- length(unique(dat[, 2]))
n <- nrow(dat)/J
dat.ctl <- dat[dat[, 4] == 0, ]
dat.trt1 <- dat[dat[, 4] == 1, ]
dat.trt2 <- dat[dat[, 4] == 2, ]
mut1j <- aggregate(dat.trt1[, 3], by = list(dat.trt1[, 2]), FUN = mean)
mut2j <- aggregate(dat.trt2[, 3], by = list(dat.trt2[, 2]), FUN = mean)
mucj <- aggregate(dat.ctl[, 3], by = list(dat.ctl[, 2]), FUN = mean)
mu_j <- rbind(mut1j, mut2j, mucj)
names(mu_j) <- c(names(dat)[2], "mu")
dat.mu <- merge(dat, mu_j, by = names(dat)[2])
SSW <- sum((dat.mu[, 3] - dat.mu[5])^2)/(n * J - J)
mut1 <- mean(mut1j[, 2])
mut2 <- mean(mut2j[, 2])
muc <- mean(mucj[, 2])
SSB <- (sum((mut1j[, 2] - mut1)^2) + sum((mut2j[, 2] - mut2)^2) + sum((mucj[,
2] - muc)^2))/(J - 3)
tau_hat <- SSB - SSW/n
f1 <- (0.5 * mut1 + 0.5 * mut2 - muc)/sqrt(tau_hat + SSW)
f2 <- (mut1 - mut2)/sqrt(tau_hat + SSW)
f3 <- sqrt((((0.5 * (mut1 + mut2) - muc)^2)/4.5 + ((mut1 - mut2)^2)/6)/(tau_hat +
SSW))
rho <- tau_hat/(tau_hat + SSW)
return(list(f1 = f1, f2 = f2, f3 = f3, ICC = rho))
# f1: effect size of treatment main effect f2: effect size of comparing
# two treatments f3: effect size of omnibus test
}
## Effect size for MRT with 2 arms based on empirical data The data has
## to be in text format where the first column of the data is the ID
## variable, the second column represents cluster, the third column is
## the outcome variable, and the fourth column is the condition variable
## (0 for control, 1 for condition). The first line of the data should
## be the variable names.
wp.effect.MRT2arm <- function(file) {
if (is.character(file)){
dat <- read.table(file, header = TRUE)
}else{
dat <- file
}
J <- length(unique(dat[, 2]))
n <- nrow(dat)/J
mutcj <- aggregate(dat[, 3], by = list(dat[, 2], dat[, 4]), FUN = mean)
mutj <- mutcj[mutcj[, 2] == 1, ]
mucj <- mutcj[mutcj[, 2] == 0, ]
beta0 <- aggregate(mutcj[, 3], by = list(mutcj[, 1]), FUN = mean)
names(beta0) <- c(names(dat)[2], "beta0")
dat1 <- merge(dat, beta0, by = names(dat)[2]) #merge beta0
mut_c <- merge(mutj, mucj, names(mutj)[1])
mut_c$beta1 <- mut_c[, 3] - mut_c[, 5]
mut_c1 <- mut_c[, c(1, 6)]
names(mut_c1) <- c(names(dat)[2], "beta1")
dat2 <- merge(dat1, mut_c1, by = names(dat)[2]) #merge beta1
mud <- sum(mut_c[, 3] - mut_c[, 5])/J
Xij <- dat[, 4] - 0.5
sg2 <- sum((dat2[, 3] - dat2[, 5] - dat2[, 6] * Xij)^2)/(J * (n - 2))
f <- mud/sqrt(sg2)
return(list(f = f))
}
## Effect size for MRT with 3 arms based on empirical data The data has
## to be in text format, where the first column of the data is the ID
## variable, the second column represents cluster, the third column is
## the outcome variable, and the fourth column is the condition variable
## (0 for control, 1 for treatment1, 2 for treatment2). The first line
## of the data should be the variable names.
wp.effect.MRT3arm <- function(file) {
if (is.character(file)){
dat <- read.table(file, header = TRUE)
}else{
dat <- file
}
J <- length(unique(dat[, 2]))
n <- nrow(dat)/J
mutcj <- aggregate(dat[, 3], by = list(dat[, 2], dat[, 4]), FUN = mean)
mut1j <- mutcj[mutcj[, 2] == 1, ]
names(mut1j)[3] <- "y1"
mut2j <- mutcj[mutcj[, 2] == 2, ]
names(mut2j)[3] <- "y2"
mucj <- mutcj[mutcj[, 2] == 0, ]
names(mucj)[3] <- "y0"
beta0 <- aggregate(mutcj[, 3], by = list(mutcj[, 1]), FUN = mean)
names(beta0) <- c(names(dat)[2], "beta0")
dat1 <- merge(dat, beta0, by = names(dat)[2]) #merge beta0
mut_c1 <- merge(mut1j, mucj, names(mut1j)[1])
mut_c <- merge(mut2j, mut_c1, names(mut2j)[1])
mut_c <- mut_c[, c(1, 3, 5, 7)]
mut_c$beta1 <- (mut_c$y1 + mut_c$y2)/2 - mut_c$y0
mut_c$beta2 <- (mut_c$y1 - mut_c$y2)
names(mut_c)[1] <- names(dat)[2]
dat2 <- merge(dat1, mut_c, by = names(dat)[2]) #merge beta1 and beta2
mud1 <- sum(mut_c$y1 - mut_c$y0)/J
mud2 <- sum(mut_c$y2 - mut_c$y0)/J
X1ij <- rep(0, nrow(dat)) # 0.5 for trt1 and trt2, -1 for ctl
X2ij <- rep(0, nrow(dat)) # 1 for trt1, -1 for trt2, 0 for ctl
X1ij[dat[, 4] == 1] <- 0.5
X1ij[dat[, 4] == 2] <- 0.5
X1ij[dat[, 4] == 0] <- -1
X2ij[dat[, 4] == 1] <- 1
X2ij[dat[, 4] == 2] <- -1
X2ij[dat[, 4] == 0] <- 0
sg2 <- sum((dat2[, 3] - dat2$beta0 - dat2$beta1 * X1ij - dat2$beta2 *
X2ij)^2)/(J * (n - 3))
f1 <- (mud1 + mud2)/2/sqrt(sg2) #effect size of treatment main effect
f2 <- (mud1 - mud2)/sqrt(sg2) #effect size of comparing the two treatments
return(list(f1 = f1, f2 = f2))
}
nuniroot <- function(f, interval, maxlength = 100) {
if (length(interval)!=2) stop("Please provide an interval with two values such as c(0,1).")
x <- seq(min(interval), max(interval), length = maxlength)
f.out <- f(x)
if (min(f.out) * max(f.out) > 0)
stop("The specified parameters do not yield valid results. Please try to supply a different interval, e.g., using interval=c(0,1), for your parameter.") else {
low <- max(f.out[f.out < 0])
high <- min(f.out[f.out > 0])
interval <- c(x[f.out == low][1], x[f.out == high][1])
uniroot(f, interval)
}
}
####################################################### function for cluster randomized trials with 2 arms##
wp.crt2arm <- function(n = NULL, f = NULL, J = NULL, icc = NULL, power = NULL,
alpha = 0.05, alternative = c("two.sided", "one.sided"), interval = NULL) {
alternative <- alternative[1]
if (sum(sapply(list(J, n, f, icc, alpha, power), is.null)) != 1)
stop("exactly one of J, n, f, icc, and alpha,power must be NULL")
if (!is.null(f) && min(f) < 0)
stop("Effect size must be positive")
if (!is.null(n) && min(n) < 1)
stop("Sample size has to be larger than 1")
if (!is.null(J) && min(J) < 2)
stop("Number of clusters has to be at least 2")
if (!is.null(icc) && !is.numeric(icc) || any(0 > icc | icc > 1))
stop("Intraclass correlation must be numeric in [0, 1]")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop("alpha must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop("Power must be numeric in [0, 1]")
## function to evaluate two-sided
if (alternative == "two.sided") {
p.body <- quote({
df <- J - 2
lambda <- sqrt(J * f^2/(4 * icc + 4 * (1 - icc)/n))
1 - pt(qt(1 - alpha/2, df), df, lambda) + pt(-qt(1 - alpha/2,
df), df, lambda)
})
} else {
## one-sided
p.body <- quote({
df <- J - 2
lambda <- sqrt(J * f^2/(4 * icc + 4 * (1 - icc)/n))
1 - pt(qt(1 - alpha, J - 2), df, lambda)
})
}
####
if (is.null(power))
power <- eval(p.body) else if (is.null(J))
J <- nuniroot(function(J) eval(p.body) - power, ifelse(is.null(interval), c(2 + 1e-10, 1000), interval))$root else if (is.null(n))
n <- nuniroot(function(n) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1, 1e+06), interval))$root else if (is.null(f))
f <- nuniroot(function(f) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-07, 1e+07), interval))$root else if (is.null(icc))
icc <- nuniroot(function(icc) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(0, 1), interval))$root else if (is.null(alpha))
alpha <- nuniroot(function(alpha) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-10, 1 - 1e-10), interval))$root else stop("internal error")
NOTE <- "n is the number of subjects per cluster."
METHOD <- "Cluster randomized trials with 2 arms"
URL <- "http://psychstat.org/crt2arm"
structure(list(J = J, n = n, f = f, icc = icc, power = power, alpha = alpha, note = NOTE, method = METHOD, url = URL), class = "webpower")
}
####################################################### function for cluster randomized trials with 3 arms##
wp.crt3arm <- function(n = NULL, f = NULL, J = NULL, icc = NULL, power = NULL,
alpha = 0.05, alternative = c("two.sided", "one.sided"), type = c("main",
"treatment", "omnibus"), interval = NULL) {
side <- alternative[1]
type <- type[1]
if (sum(sapply(list(J, n, f, icc, alpha, power), is.null)) != 1)
stop("exactly one of J, n, f, icc, and alpha,power must be NULL")
if (!is.null(f) && min(f) < 0)
stop("Effect size must be positive")
if (!is.null(n) && min(n) < 1)
stop("Sample size has to be larger than 1")
if (!is.null(J) && min(J) < 3)
stop("Number of clusters has to be at least 3")
if (!is.null(icc) && !is.numeric(icc) || any(0 > icc | icc > 1))
stop("Intraclass correlation must be numeric in [0, 1]")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop("alpha must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop("Power must be numeric in [0, 1]")
## compare the average treatment with control
if (type == "main") {
if (side == "two.sided") {
## two-sided
p.body <- quote({
df <- J - 3
lambda1 <- sqrt(J) * f/sqrt(4.5 * (icc + (1 - icc)/n))
t0 <- qt(1 - alpha/2, df)
1 - pt(t0, df, lambda1) + pt(-t0, df, lambda1)
})
} else {
## one-sided
p.body <- quote({
df <- J - 3
lambda1 <- sqrt(J) * f/sqrt(4.5 * (icc + (1 - icc)/n))
t0 <- qt(1 - alpha, df)
1 - pt(t0, df, lambda1)
})
}
}
## compare the two treatments
if (type == "treatment") {
if (side == "two.sided") {
## two-sided
p.body <- quote({
df <- J - 3
lambda2 <- sqrt(J) * f/sqrt(6 * (icc + (1 - icc)/n))
t0 <- qt(1 - alpha/2, df)
1 - pt(t0, df, lambda2) + pt(-t0, df, lambda2)
})
} else {
## one-sided
p.body <- quote({
df <- J - 3
lambda2 <- sqrt(J) * f/sqrt(6 * (icc + (1 - icc)/n))
t0 <- qt(1 - alpha, df)
1 - pt(t0, df, lambda2)
})
}
}
## omnibus test
if (type == "omnibus") {
p.body <- quote({
df1 <- 2
df2 <- J - 3
lambda3 <- J * f^2/(icc + (1 - icc)/n)
f0 <- qf(1 - alpha, df1, df2)
1 - pf(f0, df1, df2, lambda3)
})
}
####
if (is.null(power))
power <- eval(p.body) else if (is.null(J))
J <- nuniroot(function(J) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(3 + 1e-10, 1000), interval))$root else if (is.null(n))
n <- nuniroot(function(n) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(2 + 1e-10, 1e+06), interval))$root else if (is.null(f))
f <- nuniroot(function(f) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-07, 1e+07), interval))$root else if (is.null(icc))
icc <- nuniroot(function(icc) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(0, 1), interval))$root else if (is.null(alpha))
alpha <- nuniroot(function(alpha) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-10, 1 - 1e-10), interval))$root else stop("internal error")
NOTE <- "n is the number of subjects per cluster."
METHOD <- "Cluster randomized trials with 3 arms"
URL <- "http://psychstat.org/crt3arm"
structure(list(J = J, n = n, f = f, icc = icc, power = power, alpha = alpha,
note = NOTE, method = METHOD, url = URL), class = "webpower")
}
######################################################### function for multisite randomized trials with 2 arms##
wp.mrt2arm <- function(n = NULL, f = NULL, J = NULL, tau00 = NULL, tau11 = NULL,
sg2 = NULL, power = NULL, alpha = 0.05, alternative = c("two.sided",
"one.sided"), type = c("main", "site", "variance"), interval = NULL) {
type <- type[1]
side <- alternative[1]
if (sum(sapply(list(f, n, J, power), is.null)) != 1 & type == 1)
stop("exactly one of f, J, n and power must be NULL")
if (sum(sapply(list(n, J, power), is.null)) != 1 & type != 1)
stop("exactly one of J, n and power must be NULL")
if (!is.null(f) && min(f) < 0)
stop("Effect size must be positive")
if (!is.null(n) && min(n) < 2)
stop("Sample size has to be larger than 2")
if (!is.null(J) && min(J) < 2)
stop("Number of sites has to be at least 2")
if (!is.null(tau00) && min(tau00) < 0)
stop("Variance of site means must be positive")
if (!is.null(tau11) && min(tau11) < 0)
stop("Variance of treatment main effects across sites must be positive")
if ((!is.null(sg2) && min(sg2) < 0) || is.null(sg2))
stop("Between-person variation must be a positive number")
if (is.null(alpha) || (!is.null(alpha) && !is.numeric(alpha) || any(0 >
alpha | alpha > 1)))
stop("alpha must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop("Power must be numeric in [0, 1]")
if (is.null(tau11) && (type == 1 || type == 3))
stop("For this type of test, variance of treatment main effects across sites must be specified")
if (is.null(tau00) && type == 2)
stop("For this type of test, variance of site means must be specified")
## test treatment main effect #no tau00 needed
if (type == "main") {
if (side == "two.sided") {
## two-sided
p.body <- quote({
df <- J - 1
lambda1 <- sqrt(J) * f/sqrt(4/n + tau11/sg2)
t0 <- qt(1 - alpha/2, df)
1 - pt(t0, df, lambda1) + pt(-t0, df, lambda1)
})
} else {
## one-sided #no tau00 needed
p.body <- quote({
df <- J - 1
lambda1 <- sqrt(J) * f/sqrt(4/n + tau11/sg2)
t0 <- qt(1 - alpha, df)
1 - pt(t0, df, lambda1)
})
}
}
## test site variability #no tau11 and f needed
if (type == "site") {
p.body <- quote({
df1 <- J - 1
df2 <- J * (n - 2)
f0 <- qf(1 - alpha, df1, df2)
1 - pf(f0/(n * tau00/sg2 + 1), df1, df2)
})
}
## test variance of treatment effects #no tau00 and f needed
if (type == "variance") {
p.body <- quote({
df1 <- J - 1
df2 <- J * (n - 2)
f0 <- qf(1 - alpha, df1, df2)
1 - pf(f0/(n * tau11/sg2/4 + 1), df1, df2)
})
}
####
if (is.null(power))
power <- eval(p.body) else if (is.null(J))
J <- nuniroot(function(J) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1 + 1e-10, 1000), interval))$root else if (is.null(n))
n <- nuniroot(function(n) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(3 - 1e-10, 1e+06), interval))$root else if (is.null(f) & type == 1)
f <- nuniroot(function(f) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-07, 1e+07), interval))$root else stop("internal error")
NOTE <- "n is the number of subjects per cluster"
METHOD <- "Multisite randomized trials with 2 arms"
URL <- "http://psychstat.org/mrt2arm"
structure(list(J = J, n = n, f = f, tau00 = tau00, tau11 = tau11, sg2 = sg2,
power = power, alpha = alpha, note = NOTE, method = METHOD, url = URL),
class = "webpower")
}
######################################################### function for multisite randomized trials with 3 arms##
wp.mrt3arm <- function(n = NULL, f1 = NULL, f2 = NULL, J = NULL, tau = NULL,
sg2 = NULL, power = NULL, alpha = 0.05, alternative = c("two.sided",
"one.sided"), type = c("main", "treatment", "omnibus"), interval = NULL) {
type <- type[1]
side <- alternative[1]
if (type == "main")
if (sum(sapply(list(f1, n, J, power, alpha), is.null)) != 1)
stop("exactly one of f1, J, n, power and alpha must be NULL")
if (type == "treatment")
if (sum(sapply(list(f2, n, J, power, alpha), is.null)) != 1)
stop("exactly one of f2, J, n, power and alpha must be NULL")
if (type == "omnibus") {
if (sum(sapply(list(n, J, power, alpha), is.null)) != 1)
stop("exactly one of J, n, power and alpha must be NULL")
if (is.null(f1) || is.null(f2))
stop("f1 & f2 must be given")
}
if (!is.null(f1) && min(f1) < 0)
stop("Effect size must be positive")
if (!is.null(f2) && min(f2) < 0)
stop("Effect size must be positive")
if (!is.null(n) && min(n) < 3)
stop("Sample size has to be larger than 3")
if (!is.null(J) && min(J) < 2)
stop("Number of sites has to be at least 2")
if (!is.null(tau) && min(tau) < 0)
stop("Variance of treatment main effects across sites must be positive")
if ((!is.null(sg2) && min(sg2) < 0) || is.null(sg2))
stop("Between-person variation must be a positive number")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop("alpha must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop("Power must be numeric in [0, 1]")
## Test treatment main effect
if (type == "main") {
if (side == "two.sided") {
## two-sided
p.body <- quote({
df <- J - 1
lambda1 <- sqrt(J) * f1/sqrt(4.5/n + 1.5 * tau/sg2)
t0 <- qt(1 - alpha/2, df)
1 - pt(t0, df, lambda1) + pt(-t0, df, lambda1)
})
} else {
## one-sided
p.body <- quote({
df <- J - 1
lambda1 <- sqrt(J) * f1/sqrt(4.5/n + 1.5 * tau/sg2)
t0 <- qt(1 - alpha, df)
1 - pt(t0, df, lambda1)
})
}
}
## Compare two treatments
if (type == "treatment") {
if (side == "two.sided") {
## two-sided
p.body <- quote({
df <- J - 1
lambda2 <- sqrt(J) * f2/sqrt(6/n + 2 * tau/sg2)
t0 <- qt(1 - alpha/2, df)
1 - pt(t0, df, lambda2) + pt(-t0, df, lambda2)
})
} else {
## one-sided
p.body <- quote({
df <- J - 1
lambda2 <- sqrt(J) * f2/sqrt(6/n + 2 * tau/sg2)
t0 <- qt(1 - alpha, df)
1 - pt(t0, df, lambda2)
})
}
}
## Omnibus test
if (type == "omnibus") {
p.body <- quote({
df1 <- 2
df2 <- 2 * (J - 1)
lambda1 <- sqrt(J) * f1/sqrt(4.5/n + 1.5 * tau/sg2)
lambda2 <- sqrt(J) * f2/sqrt(6/n + 2 * tau/sg2)
lambda3 <- lambda1^2 + lambda2^2
f0 <- qf(1 - alpha, df1, df2)
1 - pf(f0, df1, df2, lambda3)
})
}
####
if (is.null(power))
power <- eval(p.body) else if (is.null(J))
J <- nuniroot(function(J) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(2 - 1e-10, 1000), interval))$root else if (is.null(n))
n <- nuniroot(function(n) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(3 - 1e-10, 1e+07), interval))$root else if (is.null(f1) && type != 2)
f1 <- nuniroot(function(f1) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-07, 1e+07), interval))$root else if (is.null(f2) && type != 1)
f2 <- nuniroot(function(f2) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-07, 1e+07), interval))$root else stop("internal error")
NOTE <- "n is the number of subjects per cluster"
METHOD <- "Multisite randomized trials with 3 arms"
URL <- "http://psychstat.org/mrt3arm"
structure(list(J = J, n = n, f1 = f1, f2 = f2, tau = tau, sg2 = sg2,
power = power, alpha = alpha, note = NOTE, method = METHOD, url = URL),
class = "webpower")
}
wp.regression <- function (n = NULL, p1 = NULL, p2 = 0,
f2 = NULL, alpha = 0.05, power = NULL,
type=c("regular", "Cohen")){
if (sum(sapply(list(n, f2, power, alpha), is.null)) != 1)
stop("exactly one of n, f2, power, and alpha must be NULL")
if (!is.null(f2) && min(f2) < 0)
stop("f2 must be positive")
if (!is.null(n) && min(n) < 5)
stop("sample size has to be at least 5")
if (!is.null(p1) && min(p1) < 1)
stop("number of predictor in the full model has to be larger than 1")
if (!is.null(p2) && p1 < p2)
stop("number of predictor in the full model has to be larger than that in the reduced model")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha |
alpha > 1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power |
power > 1))
stop(sQuote("power"), " must be numeric in [0, 1]")
p.body <- quote({
u <- p1 - p2
v <- n - p1 - 1
if (type[1]=="Cohen"){
lambda <- f2 * (u + v + 1)
}else{
lambda <- f2 * n
}
pf(qf(alpha, u, v, lower = FALSE), u, v, lambda, lower = FALSE)
})
if (is.null(power))
power <- eval(p.body)
else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(5 +
p1 + 1e-10, 1e+05))$root
else if (is.null(f2))
f2 <- uniroot(function(f2) eval(p.body) - power, c(1e-07,
1e+07))$root
else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power,
c(1e-10, 1 - 1e-10))$root
else stop("internal error")
METHOD <- "Power for multiple regression"
URL <- "http://psychstat.org/regression"
structure(list(n = n, p1 = p1, p2 = p2, f2 = f2, alpha = alpha,
power = power, method = METHOD, url = URL), class = "webpower")
}
wp.logistic <- function(n = NULL, p0 = NULL, p1 = NULL, alpha = 0.05, power = NULL,
alternative = c("two.sided", "less", "greater"), family = c("Bernoulli",
"exponential", "lognormal", "normal", "Poisson", "uniform"), parameter = NULL) {
sig.level <- alpha
if (sum(sapply(list(n, power), is.null)) != 1)
stop("exactly one of n, power, and alpha must be NULL")
if (is.null(p0) || !is.null(p0) && !is.numeric(p0) || any(0 > p0 |
p0 > 1))
stop(sQuote("p0"), " must be numeric in (0,1)")
if (is.null(p1) || !is.null(p1) && !is.numeric(p1) || any(0 > p1 |
p1 > 1))
stop(sQuote("p1"), " must be numeric in (0,1)")
if (!is.null(n) && min(n) < 5)
stop("number of observations must be at least 5")
if (is.null(alpha) || !is.null(alpha) & !is.numeric(alpha) || any(alpha <
0 | alpha > 1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) & !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
alternative <- match.arg(alternative)
tside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
p.body <- quote({
s * pnorm(-qnorm(1 - alpha) - sqrt(n)/sqrt(g * v0 + (1 - g) * v1) *
beta1) + t * pnorm(-qnorm(1 - alpha) + sqrt(n)/sqrt(g * v0 +
(1 - g) * v1) * beta1)
})
if (family == "Bernoulli") {
## binomial predictor
if (is.null(parameter)) {
B <- 0.5
} else {
B <- parameter
}
g <- 0
odds <- p1/(1 - p1)/(p0/(1 - p0))
beta1 <- log(odds) ##beta1 under H_A in the population
beta0 <- log(p0/(1 - p0))
d <- B * p1 * (1 - p1) + (1 - B) * p0 * (1 - p0)
e <- B * p1 * (1 - p1)
f <- B * p1 * (1 - p1)
v1 <- d/(d * f - e^2)
### v0
mu1 <- B * p1 + (1 - B) * p0
i00 <- log(mu1/(1 - mu1))
pn <- 1/(1 + exp(-i00))
a <- pn * (1 - pn) #00
b <- B * pn * (1 - pn) #11
c <- B * pn * (1 - pn) #01
v0 <- b/(a * b - c^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power))
power <- eval(p.body)
if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + 1e-10,
1e+07))$root
}
if (family == "exponential") {
## Exponential predictor
if (is.null(parameter)) {
lambda <- 1
} else {
lambda <- parameter
}
g <- 0
beta0 <- log(p0/(1 - p0))
odds <- p1/(1 - p1)/(p0/(1 - p0))
beta1 <- log(odds) ##beta1 under H_A in the population##beta1 under H_A in the population
d <- integrate(function(x) (1 - (1 + exp(-beta0 - beta1 * x))^(-1)) *
(1 + exp(-beta0 - beta1 * x))^(-1) * dexp(x, lambda), 0, Inf,
subdivisions = 100L)$value
e <- integrate(function(x) x * (1 - (1 + exp(-beta0 - beta1 * x))^(-1)) *
(1 + exp(-beta0 - beta1 * x))^(-1) * dexp(x, lambda), 0, Inf,
subdivisions = 100L)$value
f <- integrate(function(x) x^2 * (1 - (1 + exp(-beta0 - beta1 *
x))^(-1)) * (1 + exp(-beta0 - beta1 * x))^(-1) * dexp(x, lambda),
0, Inf, subdivisions = 100L)$value
v1 <- d/(d * f - e^2)
mu1 <- integrate(function(x) 1/(1 + exp(-beta0 - beta1 * x)) *
dexp(x, lambda), 0, Inf, subdivisions = 100L)$value
i00 <- log(mu1/(1 - mu1))
pn <- 1/(1 + exp(-i00))
a <- pn * (1 - pn) #00
b <- 2 * lambda^(-2) * pn * (1 - pn) #11
c <- lambda^(-1) * pn * (1 - pn) #01
v0 <- a/(a * b - c^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power))
power <- eval(p.body)
if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + 1e-10,
1e+07))$root
}
if (family == "lognormal") {
g <- 0
if (is.null(parameter)) {
mu <- 0
sigma <- 1
} else {
if (length(parameter) != 2)
stop("Both mean and standard deviation of the log-normal distribution have to be provided as a vector.")
mu <- parameter[1]
sigma <- parameter[2]
}
beta0 <- log(p0/(1 - p0))
odds <- p1/(1 - p1)/(p0/(1 - p0))
beta1 <- log(odds) ##beta1 under H_A in the population
d <- integrate(function(x) (1 - (1 + exp(-beta0 - beta1 * x))^(-1)) *
(1 + exp(-beta0 - beta1 * x))^(-1) * dlnorm(x, mu, sigma),
0, Inf, subdivisions = 100L)$value
e <- integrate(function(x) x * (1 - (1 + exp(-beta0 - beta1 * x))^(-1)) *
(1 + exp(-beta0 - beta1 * x))^(-1) * dlnorm(x, mu, sigma),
0, Inf, subdivisions = 100L)$value
f <- integrate(function(x) x^2 * (1 - (1 + exp(-beta0 - beta1 *
x))^(-1)) * (1 + exp(-beta0 - beta1 * x))^(-1) * dlnorm(x,
mu, sigma), 0, Inf, subdivisions = 100L)$value
v1 <- d/(d * f - e^2)
mu1 <- integrate(function(x) 1/(1 + exp(-beta0 - beta1 * x)) *
dlnorm(x, mu, sigma), 0, Inf, subdivisions = 100L)$value
i00 <- log(mu1/(1 - mu1))
pn <- 1/(1 + exp(-i00))
a <- pn * (1 - pn) #00
b <- (exp(sigma * sigma) - 1) * exp(2 * mu + sigma^2) * pn * (1 -
pn) #11
c <- exp(mu + 0.5 * sigma^2) * pn * (1 - pn) #01
v0 <- a/(a * b - c^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power))
power <- eval(p.body)
if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + 1e-10,
1e+07))$root
}
if (family == "normal") {
g <- 0
if (is.null(parameter)) {
mu <- 0
sigma <- 1
} else {
if (length(parameter) != 2)
stop("Both mean and standard deviation of the normal distribution have to be provided as a vector.")
mu <- parameter[1]
sigma <- parameter[2]
}
beta0 <- log(p0/(1 - p0))
odds <- p1/(1 - p1)/(p0/(1 - p0))
beta1 <- log(odds) ##beta1 under H_A in the population
d <- integrate(function(x) (1 - (1 + exp(-beta0 - beta1 * x))^(-1)) *
(1 + exp(-beta0 - beta1 * x))^(-1) * dnorm(x, mu, sigma), -Inf,
Inf, subdivisions = 100L)$value
e <- integrate(function(x) x * (1 - (1 + exp(-beta0 - beta1 * x))^(-1)) *
(1 + exp(-beta0 - beta1 * x))^(-1) * dnorm(x, mu, sigma), -Inf,
Inf, subdivisions = 100L)$value
f <- integrate(function(x) x^2 * (1 - (1 + exp(-beta0 - beta1 *
x))^(-1)) * (1 + exp(-beta0 - beta1 * x))^(-1) * dnorm(x, mu,
sigma), -Inf, Inf, subdivisions = 100L)$value
v1 <- d/(d * f - e^2)
mu1 <- integrate(function(x) 1/(1 + exp(-beta0 - beta1 * x)) *
dnorm(x, mu, sigma), -Inf, Inf, subdivisions = 100L)$value
i00 <- log(mu1/(1 - mu1))
pn <- 1/(1 + exp(-i00))
a <- pn * (1 - pn) #00
b <- sigma^2 * pn * (1 - pn) #11
c <- mu * pn * (1 - pn) #01
v0 <- a/(a * b - c^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power))
power <- eval(p.body)
if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + 1e-10,
1e+07))$root
}
if (family == "Poisson") {
g <- 0
if (is.null(parameter)) {
lambda <- 1
} else {
lambda <- parameter
}
beta0 <- log(p0/(1 - p0))
odds <- p1/(1 - p1)/(p0/(1 - p0))
beta1 <- log(odds) ##beta1 under H_A in the population
d <- sum(sapply(0:1e+05, function(x) (1 - (1 + exp(-beta0 - beta1 *
x))^(-1)) * (1 + exp(-beta0 - beta1 * x))^(-1) * dpois(x, lambda)))
e <- sum(sapply(0:1e+05, function(x) x * (1 - (1 + exp(-beta0 -
beta1 * x))^(-1)) * (1 + exp(-beta0 - beta1 * x))^(-1) * dpois(x,
lambda)))
f <- sum(sapply(0:1e+05, function(x) x^2 * (1 - (1 + exp(-beta0 -
beta1 * x))^(-1)) * (1 + exp(-beta0 - beta1 * x))^(-1) * dpois(x,
lambda)))
v1 <- d/(d * f - e^2)
mu1 <- sum(sapply(0:1e+05, function(x) 1/(1 + exp(-beta0 - beta1 *
x)) * dpois(x, lambda)))
i00 <- log(mu1/(1 - mu1))
pn <- 1/(1 + exp(-i00))
a <- pn * (1 - pn) #00
b <- lambda * pn * (1 - pn) #11
c <- lambda * pn * (1 - pn) #01
v0 <- a/(a * b - c^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power))
power <- eval(p.body)
if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + 1e-10,
1e+07))$root
}
if (family == "uniform") {
g <- 0
if (is.null(parameter)) {
L <- 0
R <- 1
} else {
if (length(parameter) != 2)
stop("The lower and upper bounds have to be provided as a vector")
L <- parameter[1]
R <- parameter[2]
}
beta0 <- log(p0/(1 - p0))
odds <- p1/(1 - p1)/(p0/(1 - p0))
beta1 <- log(odds) ##beta1 under H_A in the population
d <- integrate(function(x) (1 - (1 + exp(-beta0 - beta1 * x))^(-1)) *
(1 + exp(-beta0 - beta1 * x))^(-1)/(R - L), L, R, subdivisions = 100L)$value
e <- integrate(function(x) x * (1 - (1 + exp(-beta0 - beta1 * x))^(-1)) *
(1 + exp(-beta0 - beta1 * x))^(-1)/(R - L), L, R, subdivisions = 100L)$value
f <- integrate(function(x) x^2 * (1 - (1 + exp(-beta0 - beta1 *
x))^(-1)) * (1 + exp(-beta0 - beta1 * x))^(-1)/(R - L), L,
R, subdivisions = 100L)$value
v1 <- d/(d * f - e^2)
mu1 <- integrate(function(x) 1/(1 + exp(-beta0 - beta1 * x))/(R -
L), L, R, subdivisions = 100L)$value
i00 <- log(mu1/(1 - mu1))
pn <- 1/(1 + exp(-i00))
a <- pn * (1 - pn) #00
b <- (R - L)^2/12 * pn * (1 - pn) #11
c <- (R + L)/2 * pn * (1 - pn) #01
v0 <- a/(a * b - c^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power))
power <- eval(p.body)
if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + 1e-10,
1e+07))$root
}
METHOD <- "Power for logistic regression"
URL <- "http://psychstat.org/logistic"
structure(list(p0 = p0, p1 = p1, beta0 = beta0, beta1 = beta1, n = n,
alpha = sig.level, power = power, alternative = alternative, family = family,
parameter = parameter, method = METHOD, url = URL), class = "webpower")
}
wp.mediation <- function(n = NULL, power = NULL, a = 0.5, b = 0.5, varx = 1,
vary = 1, varm = 1, alpha = 0.05, interval = NULL) {
if (sum(sapply(list(n, a, b, varx, varm, vary, power, alpha), is.null)) !=
1)
stop("exactly one of n, a, b, varx, varm, vary, power, and alpha must be NULL")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
## power formulae
p.body <- quote({
numer <- sqrt(n) * a * b
denom <- sqrt(a^2 * vary/(varm - a^2 * varx) + b^2 * (varm - a^2 *
varx)/varx)
delta <- numer/denom
alpha2 <- alpha/2
za2 <- qnorm(1 - alpha2)
1 - pnorm(za2 - delta) + pnorm(-za2 - delta)
})
if (is.null(power))
power <- eval(p.body) else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-07, 1e+07), interval))$root else if (is.null(a)) {
astart <- varm/varx
alow <- -sqrt(astart) + 1e-06
aup <- sqrt(astart) - 1e-06
a <- nuniroot(function(a) eval(p.body) - power, c(alow, aup))$root
} else if (is.null(b))
b <- nuniroot(function(b) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(-10, 10), interval))$root else if (is.null(varx))
varx <- nuniroot(function(varx) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-10, 1e+07), interval))$root else if (is.null(vary))
vary <- nuniroot(function(vary) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-10, 1e+07), interval))$root else if (is.null(varm))
varm <- nuniroot(function(varm) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-10, 1e+07), interval))$root else if (is.null(alpha))
alpha <- nuniroot(function(alpha) eval(p.body) - power, ifelse(rep(is.null(interval), 2), c(1e-10, 1 - 1e-10), interval))$root else stop("internal error")
METHOD <- "Power for simple mediation"
URL <- "http://psychstat.org/mediation"
structure(list(n = n, power = power, a = a, b = b, varx = varx, varm = varm,
vary = vary, alpha = alpha, method = METHOD, url = URL), class = "webpower")
}
estCRT2arm <- function(file) {
dat <- read.table(file, header = TRUE)
X <- dat[, 4] - 0.5
newdata <- data.frame(dat, X)
model <- summary(lmer(score ~ X + (1 | cluster), data = newdata))
parest <- model$coefficients
J <- model$ngrps
p.value <- (1 - pt(abs(parest[, 3]), J - 2)) * 2
test <- cbind(parest, p.value)
cat("Fixed effects (X: Treatment main effect):\n")
print(test)
}
estCRT3arm <- function(file) {
dat <- read.table(file, header = T)
X1 <- X2 <- rep(0, dim(dat)[1])
X1[dat[, 4] == 1] <- X1[dat[, 4] == 2] <- 1/3
X1[dat[, 4] == 0] <- -2/3
X2[dat[, 4] == 1] <- 1/2
X2[dat[, 4] == 2] <- -1/2
newdata <- data.frame(dat, X1, X2)
model <- summary(lmer(score ~ X1 + X2 + (1 | cluster), data = newdata))
parest <- model$coefficients
rownames(parest)[2] <- "X"
rownames(parest)[3] <- "X1-X2"
J <- model$ngrps
p.value <- (1 - pt(abs(parest[, 3]), J - 3)) * 2
test <- cbind(parest, p.value)
# Omnibus test
model1 <- lmer(score ~ X1 + X2 + (1 | cluster), data = newdata, REML = FALSE)
model2 <- lmer(score ~ (1 | cluster), data = newdata, REML = FALSE)
omnibus <- anova(model1, model2)$"Pr(>Chisq)"[2]
cat("Fixed effects (X: Treatment main effect; X1-X2: Comparing two treatments):\n")
print(test)
cat("Omnibust test: p-value = ")
cat(omnibus)
}
estMRT2arm <- function(file) {
dat <- read.table(file, header = TRUE)
X <- dat[, 4] - 0.5
newdata <- data.frame(dat, X)
model <- summary(lmer(score ~ X + (X | cluster), data = newdata))
parest <- model$coefficients
J <- model$ngrps
p.value <- (1 - pt(abs(parest[, 3]), J - 1)) * 2
test <- cbind(parest, p.value)
cat("Fixed effects (X: Treatment main effect):\n")
print(test)
}
estMRT3arm <- function(file) {
dat <- read.table(file, header = T)
X1 <- X2 <- rep(0, dim(dat)[1])
X1[dat[, 4] == 1] <- X1[dat[, 4] == 2] <- 1/3
X1[dat[, 4] == 0] <- -2/3
X2[dat[, 4] == 1] <- 1/2
X2[dat[, 4] == 2] <- -1/2
newdata <- data.frame(dat, X1, X2)
model <- summary(lmer(score ~ X1 + X2 + (X1 + X2 | cluster), data = newdata))
parest <- model$coefficients
rownames(parest)[2] <- "X"
rownames(parest)[3] <- "X1-X2"
J <- model$ngrps
p.value <- (1 - pt(abs(parest[, 3]), J - 1)) * 2
test <- cbind(parest, p.value)
# Omnibus test
model1 <- lmer(score ~ X1 + X2 + (X1 + X2 | cluster), data = newdata,
REML = FALSE)
model2 <- lmer(score ~ (X1 + X2 | cluster), data = newdata, REML = FALSE)
omnibus <- anova(model1, model2)$"Pr(>Chisq)"[2]
cat("Fixed effects (X: Treatment main effect; X1-X2: Comparing two treatments):\n")
print(test)
cat("Omnibust test: p-value = ")
cat(omnibus)
}
wp.correlation <- function(n = NULL, r = NULL, power = NULL, p = 0, rho0 = 0,
alpha = 0.05, alternative = c("two.sided", "less", "greater")) {
if (sum(sapply(list(n, r, power, alpha), is.null)) != 1)
stop("exactly one of n, r, power, and alpha must be NULL")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha >
1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power >
1))
stop(sQuote("power"), " must be numeric in [0, 1]")
if (!is.null(n) && min(n) < 4)
stop("number of observations must be at least 4")
alternative <- match.arg(alternative)
tside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
if (tside == 2 && !is.null(r))
r <- abs(r)
if (tside == 3) {
p.body <- quote({
delta <- sqrt(n - 3 - p) * (log((1 + r)/(1 - r))/2 + r/(n -
1 - p)/2 * (1 + (5 + r^2)/(n - 1 - p)/4 + (11 + 2 * r^2 +
3 * r^4)/(n - 1 - p)^2/8) - log((1 + rho0)/(1 - rho0))/2 -
rho0/(n - 1 - p)/2)
v <- (n - 3 - p)/(n - 1 - p) * (1 + (4 - r^2)/(n - 1 - p)/2 +
(22 - 6 * r^2 - 3 * r^4)/(n - 1 - p)^2/6)
zalpha <- qnorm(1 - alpha)
pnorm((delta - zalpha)/sqrt(v))
})
}
if (tside == 1) {
p.body <- quote({
delta <- sqrt(n - 3 - p) * (log((1 + r)/(1 - r))/2 + r/(n -
1 - p)/2 * (1 + (5 + r^2)/(n - 1 - p)/4 + (11 + 2 * r^2 +
3 * r^4)/(n - 1 - p)^2/8) - log((1 + rho0)/(1 - rho0))/2 -
rho0/(n - 1 - p)/2)
v <- (n - 3 - p)/(n - 1 - p) * (1 + (4 - r^2)/(n - 1 - p)/2 +
(22 - 6 * r^2 - 3 * r^4)/(n - 1 - p)^2/6)
zalpha <- qnorm(1 - alpha)
pnorm((-delta - zalpha)/sqrt(v))
})
}
if (tside == 2) {
p.body <- quote({
delta <- sqrt(n - 3 - p) * (log((1 + r)/(1 - r))/2 + r/(n -
1 - p)/2 * (1 + (5 + r^2)/(n - 1 - p)/4 + (11 + 2 * r^2 +
3 * r^4)/(n - 1 - p)^2/8) - log((1 + rho0)/(1 - rho0))/2 -
rho0/(n - 1 - p)/2)
v <- (n - 3 - p)/(n - 1 - p) * (1 + (4 - r^2)/(n - 1 - p)/2 +
(22 - 6 * r^2 - 3 * r^4)/(n - 1 - p)^2/6)
zalpha <- qnorm(1 - alpha/2)
pnorm((delta - zalpha)/sqrt(v)) + pnorm((-delta - zalpha)/sqrt(v))
})
}
if (is.null(power))
power <- eval(p.body) else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(4 + p + 1e-10,
1e+07))$root else if (is.null(r)) {
if (tside == 2) {
r <- uniroot(function(r) eval(p.body) - power, c(1e-10, 1 -
1e-10))$root
} else {
r <- uniroot(function(r) eval(p.body) - power, c(-1 + 1e-10,
1 - 1e-10))$root
}
} else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10,
1 - 1e-10))$root else stop("internal error")
METHOD <- "Power for correlation"
URL <- "http://psychstat.org/correlation"
structure(list(n = n, r = r, alpha = alpha, power = power, alternative = alternative,
method = METHOD, url = URL), class = "webpower")
}
wp.poisson <- function(n = NULL, exp0 = NULL, exp1 = NULL, alpha = 0.05, power = NULL,
alternative = c("two.sided", "less", "greater"),
family = c("Bernoulli", "exponential", "lognormal",
"normal", "Poisson", "uniform"),
parameter = NULL, subdivisions=200L,
i.method=c("numerical", "MC"), mc.iter=20000)
{
sig.level <- alpha
if(sum(sapply(list(n, power), is.null)) != 1)
stop("exactly one of n, power, andalpha must be NULL")
if(is.null(exp0) || !is.null(exp0) && !is.numeric(exp0) ||
exp0 <= 0)
stop(sQuote("exp0"), " must be numeric in (0,Infinity)")
if(is.null(exp1) || !is.null(exp1) && !is.numeric(exp1) ||
exp1 <= 0)
stop(sQuote("exp1"), " must be numeric in (0,Infinity)")
if(!is.null(n) && min(n) < 5)
stop("number of observations must be at least 5")
if(is.null(alpha) || !is.null(alpha) & !is.numeric(alpha) ||
any(alpha < 0 | alpha > 1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if(!is.null(power) & !is.numeric(power) || any(0 > power |
power > 1))
stop(sQuote("power"), " must be numeric in [0, 1]")
alternative <- match.arg(alternative)
tside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
p.body <- quote({
s * pnorm(-qnorm(1 - alpha) - sqrt(n)/sqrt(v1) * beta1) +
t * pnorm(-qnorm(1 - alpha) + sqrt(n)/sqrt(v1) *
beta1)
})
if (family == "Bernoulli") {
if (is.null(parameter)) {
B <- 0.5
}
else {
B <- parameter
}
beta1 <- log(exp1)
beta0 <- log(exp0)
a <- (1 - B) * exp(beta0) + B * exp(beta0 + beta1)
b <- B * exp(beta0 + beta1)
c <- B * exp(beta0 + beta1)
v1 <- a/(a * c - b^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power)) {
power <- eval(p.body)
}
if (is.null(n)) {
n <- uniroot(function(n) eval(p.body) - power, c(2 +
1e-10, 1e+07))$root
}
}
if (family == "exponential") {
if (is.null(parameter)) {
lambda <- 1
}
else {
lambda <- parameter
}
beta1 <- log(exp1); beta0 <- log(exp0)
if (i.method[1] == "MC"){
x.s=rexp(mc.iter,lambda)
d=mean(exp0*exp(log(exp1)*x.s))
e=mean(x.s*exp0*exp(log(exp1)*x.s))
f=mean(x.s^2*exp0*exp(log(exp1)*x.s))
}else{
d.integ<-function(x,beta0,beta1,lamba){#beta1!=lambda
lambda*exp(beta0)/(beta1-lambda)*exp((beta1-lambda)*x)
}
e.integ<-function(x,beta0,beta1,lambda){
lambda*exp(beta0)/(beta1-lambda)*x*exp((beta1-lambda)*x) -1/(beta1-lambda)*d.integ(x,beta0,beta1,lambda)
}
f.integ<-function(x,beta0,beta1,lambda){
x^2*d.integ(x,beta0,beta1,lambda)-2/(beta1-lambda)*e.integ(x,beta0,beta1,lambda)
}
if(beta1<lambda){
d=d.integ(100/(lambda-beta1),beta0,beta1,lambda)-d.integ(0,beta0,beta1,lambda)
e=e.integ(100/(lambda-beta1),beta0,beta1,lambda)-e.integ(0,beta0,beta1,lambda)
f=f.integ(100/(lambda-beta1), beta0,beta1,lambda)-f.integ(0,beta0,beta1,lambda)
v1=d/(d*f-e^2)}
if(beta1>lambda){
d=d.integ(100/abs(lambda-beta1),beta0,beta1,lambda)-d.integ(0,beta0,beta1,lambda)
e=e.integ(100/abs(lambda-beta1),beta0,beta1,lambda)-e.integ(0,beta0,beta1,lambda)
f=f.integ(100/abs(lambda-beta1), beta0,beta1, lambda)-f.integ(0,beta0,beta1,lambda)
v1=d/(d*f-e^2)
}else{
v1=0.00001#the actual value is 0
}
}
v1 <- d/(d * f - e^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power)) {
power <- eval(p.body)
}
if (is.null(n)) {
n <- uniroot(function(n) eval(p.body) - power, c(2 +
1e-10, 1e+07))$root
}
}
if (family == "lognormal") {
if (is.null(parameter)) {
mu <- 0
sigma <- 1
}
else {
if (length(parameter) != 2)
stop("Both mean and standard deviation of the log-normal distribution have to be provided as a vector.")
mu <- parameter[1]
sigma <- parameter[2]
}
beta1 <- log(exp1)
beta0 <- log(exp0)
if (i.method[1] == "MC"){
x.s=rlnorm(mc.iter,mu,sigma)
d=mean(exp0*exp(log(exp1)*x.s))
e=mean(x.s*exp0*exp(log(exp1)*x.s))
f=mean(x.s^2*exp0*exp(log(exp1)*x.s))
}else{
d <- integrate(function(x) exp(beta0 + beta1 * x) * dlnorm(x, mu, sigma), 0, Inf, subdivisions = subdivisions)$value
e <- integrate(function(x) x * exp(beta0 + beta1 * x) * dlnorm(x, mu, sigma), 0, Inf, subdivisions = subdivisions)$value
f <- integrate(function(x) x^2 * exp(beta0 + beta1 * x) * dlnorm(x, mu, sigma), 0, Inf, subdivisions = subdivisions)$value
}
v1 <- d/(d * f - e^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power)) {
power <- eval(p.body)
}
if (is.null(n)) {
n <- uniroot(function(n) eval(p.body) - power, c(2 +
1e-10, 1e+07))$root
}
}
if (family == "normal") {
if (is.null(parameter)) {
mu <- 0
sigma <- 1
}
else {
if (length(parameter) != 2)
stop("Both mean and standard deviation of the normal distribution have to be provided as a vector.")
mu <- parameter[1]
sigma <- parameter[2]
}
beta1 <- log(exp1)
beta0 <- log(exp0)
d=exp0*exp(beta1^2*sigma^2/2+beta1*mu)
e=exp0*exp(beta1^2*sigma^2/2+beta1*mu)*(mu+beta1*sigma^2)
f=exp0*exp(beta1^2*sigma^2/2+beta1*mu)*((mu+beta1*sigma^2)^2+sigma^2)
v1 <- d/(d * f - e^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power)) {
power <- eval(p.body)
}
if (is.null(n)) {
n <- uniroot(function(n) eval(p.body) - power, c(2 +
1e-10, 1e+07))$root
}
}
if (family == "Poisson") {
if (is.null(parameter)) {
lambda <- 1
}
else {
lambda <- parameter
}
beta1 <- log(exp1)
beta0 <- log(exp0)
d <- sum(sapply(0:1e+05, function(x) exp(beta0 + beta1 *
x) * dpois(x, lambda)))
e <- sum(sapply(0:1e+05, function(x) x * exp(beta0 +
beta1 * x) * dpois(x, lambda)))
f <- sum(sapply(0:1e+05, function(x) x^2 * exp(beta0 +
beta1 * x) * dpois(x, lambda)))
v1 <- d/(d * f - e^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power)) {
power <- eval(p.body)
}
if (is.null(n)) {
n <- uniroot(function(n) eval(p.body) - power, c(2 +
1e-10, 1e+07))$root
}
}
if (family == "uniform") {
if (is.null(parameter)) {
L <- 0
R <- 1
}
else {
if (length(parameter) != 2)
stop("The lower and upper bounds have to be provided as a vector")
L <- parameter[1]
R <- parameter[2]
}
beta1 <- log(exp1)
beta0 <- log(exp0)
d <- integrate(function(x) exp(beta0 + beta1 * x)/(R -
L), L, R, subdivisions = subdivisions)$value
e <- integrate(function(x) x * exp(beta0 + beta1 * x)/(R -
L), L, R, subdivisions = subdivisions)$value
f <- integrate(function(x) x^2 * exp(beta0 + beta1 *
x)/(R - L), L, R, subdivisions = subdivisions)$value
v1 <- d/(d * f - e^2)
if (tside == 1) {
s <- 1
t <- 0
alpha <- alpha
}
if (tside == 2) {
s <- 1
t <- 1
alpha <- alpha/2
}
if (tside == 3) {
s <- 0
t <- 1
alpha <- alpha
}
if (is.null(power)) {
power <- eval(p.body)
}
if (is.null(n)) {
n <- uniroot(function(n) eval(p.body) - power, c(2 +
1e-10, 1e+07))$root
}
}
METHOD <- "Power for Poisson regression"
URL <- "http://psychstat.org/poisson"
structure(list(n = n, power = power, alpha = sig.level, exp0 = exp0,
exp1 = exp1, beta0 = beta0, beta1 = beta1, alternative = alternative,
family = family, method = METHOD,
url = URL), class = "webpower")
}
wp.sem.chisq <-
function (n = NULL, df = NULL, effect = NULL, power = NULL, alpha = 0.05)
{
if (sum(sapply(list(n, df, power, alpha, effect), is.null)) != 1)
stop("exactly one of sample size, degrees of freedom, effect size, alpha and power must be NULL")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 >
alpha | alpha > 1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power |
power > 1))
stop(sQuote("power"), " must be numeric in [0, 1]")
p.body <- quote({
ncp<-(n-1)*effect
#ncp<-n*effect
calpha<-qchisq(1-alpha, df)
1-pchisq(calpha, df, ncp)
})
if (is.null(power))
power <- eval(p.body)
else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(10 +
1e-10, 1e+08))$root
else if (is.null(df))
df <- uniroot(function(df) eval(p.body) - power, c(1,
10000))$root
else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) -
power, c(1e-10, 1 - 1e-10))$root
else if (is.null(effect))
effect <- uniroot(function(effect) eval(p.body) - power, c(0, 1))$root
else stop("internal error")
METHOD <- "Power for SEM (Satorra & Saris, 1985)"
URL <- "http://psychstat.org/semchisq"
structure(list(n = n, df = df, effect = effect, power = power, alpha = alpha, url=URL, method=METHOD),
class = "webpower")
}
wp.sem.rmsea <-
function (n = NULL, df = NULL, rmsea0 = NULL, rmsea1 = NULL, power = NULL, alpha = 0.05, type=c('close','notclose')) {
if (sum(sapply(list(n, df, power, alpha, rmsea0, rmsea1), is.null)) != 1)
stop("exactly one of sample size, degrees of freedom, rmsea0, rmsea1, alpha and power must be NULL")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 >
alpha | alpha > 1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power |
power > 1))
stop(sQuote("power"), " must be numeric in [0, 1]")
type <- match.arg(type)
htype <- switch(type, close = 1, notclose = 2)
if (htype==1){
p.body <- quote({
ncp0<-(n-1)*df*rmsea0^2
ncp1<-(n-1)*df*rmsea1^2
calpha<-qchisq(1-alpha, df, ncp0)
1-pchisq(calpha, df, ncp1)
})
}else{
p.body <- quote({
ncp0<-(n-1)*df*rmsea0^2
ncp1<-(n-1)*df*rmsea1^2
calpha<-qchisq(alpha, df, ncp0)
pchisq(calpha, df, ncp1)
})
}
if (is.null(power))
power <- eval(p.body)
else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 +
1e-10, 1e+05))$root
else if (is.null(df))
df <- uniroot(function(df) eval(p.body) - power, c(1,
10000))$root
else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) -
power, c(1e-10, 1 - 1e-10))$root
else if (is.null(rmsea0))
rmsea0 <- uniroot(function(rmsea0) eval(p.body) -
power, c(0, 1))$root
else if (is.null(rmsea1))
rmsea1 <- uniroot(function(rmsea1) eval(p.body) -
power, c(0, 1))$root
else stop("internal error")
METHOD <- "Power for SEM based on RMSEA"
URL <- "http://psychstat.org/rmsea"
structure(list(n = n, df = df, rmsea0 = rmsea0, rmsea1 = rmsea1, power = power, alpha = alpha, method=METHOD, url=URL),
class = "webpower")
}
wp.kanova<-function (n = NULL, ndf = NULL, f = NULL, ng=NULL, alpha = 0.05, power = NULL)
{
if (sum(sapply(list(n, ndf, f, ng, power, alpha), is.null)) !=
1)
stop("exactly one of n, ndf, f, ng, power, alpha must be NULL")
if (!is.null(f) && min(f) < 0)
stop("f must be positive")
if (!is.null(n) && min(n) < 1)
stop("degree of freedom u for numerator must be at least 1")
if (!is.null(ndf) && min(ndf) < 1)
stop("degree of freedom v for denominator must be at least 1")
if (!is.null(ng) && min(ng) < 1)
stop("number of groups must be at least 1")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 >
alpha | alpha > 1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power |
power > 1))
stop(sQuote("power"), " must be numeric in [0, 1]")
p.body <- quote({
lambda <- f^2 * n
ddf <- n - ng
pf(qf(alpha, ndf, ddf, lower = FALSE), ndf, ddf, lambda,
lower = FALSE)
})
if (is.null(power))
power <- eval(p.body)
else if (is.null(n)){
n <- uniroot(function(n) eval(p.body) - power, c(1 +
ng, 1e+07))$root
ddf <- n - ng
}
else if (is.null(ndf)){
ndf <- uniroot(function(ndf) eval(p.body) - power, c(1 +
1e-10, 1e+05))$root
ddf <- n -ng
}
else if (is.null(ng)) {
ng <- uniroot(function(ng) eval(p.body) - power, c(1e-07,
1e+07))$root
ddf <- n -ng
}
else if (is.null(f)){
f <- uniroot(function(f) eval(p.body) - power, c(1e-07,
1e+07))$root
ddf <- n -ng
}
else if (is.null(alpha)){
alpha <- uniroot(function(alpha) eval(p.body) -
power, c(1e-10, 1 - 1e-10))$root
ddf <- n -ng
}
else stop("internal error")
METHOD <- "Multiple way ANOVA analysis"
URL <- "http://psychstat.org/kanova"
NOTE <- "Sample size is the total sample size"
structure(list(n = n, ndf = ndf, ddf=ddf, f = f, ng=ng, alpha = alpha,
power = power, method = METHOD, url=URL, note=NOTE), class = "webpower")
}
## ANOVA with binary data
wp.anova.binary<-function(k = NULL, n = NULL, V = NULL, alpha = 0.05, power = NULL){
if (sum(sapply(list(k, n, V, power, alpha), is.null)) != 1)
stop("exactly one of k, n, V, power, and alpha must be NULL")
if (!is.null(V) && min(V) < 0)
stop("V must be positive")
if (!is.null(k) && min(k) < 2)
stop("number of groups must be at least 2")
if (!is.null(n) && min(n) < 2)
stop("number of observations in each group must be at least 2")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha > 1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power > 1))
stop(sQuote("power"), " must be numeric in [0, 1]")
p.body <- quote({
chi <- (V)^2*n*(k-1)
df <- k - 1
crit.val <- qchisq(p=1-alpha,df=df,ncp=0)
1-pchisq(crit.val,df=df,ncp=chi)
})
if (is.null(power))
power <- eval(p.body)
else if (is.null(k))
k <- uniroot(function(k) eval(p.body) - power, c(2 + 1e-10, 100))$root
else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + k + 1e-10, 1e+05))$root
else if (is.null(V))
V <- uniroot(function(V) eval(p.body) - power, c(1e-07, 1e+07))$root
else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10, 1 - 1e-10))$root
else stop("internal error")
NOTE <- "n is the total sample size"
METHOD <- "One-way Analogous ANOVA with Binary Data"
URL <- "http://psychstat.org/anovabinary"
structure(list(k = k, n = n, V = V, alpha = alpha, power = power,
note = NOTE, method = METHOD, url=URL), class = "webpower")
}
wp.anova.count<-function(k = NULL, n = NULL, V = NULL, alpha = 0.05, power = NULL){
if (sum(sapply(list(k, n, V, power, alpha), is.null)) != 1)
stop("exactly one of k, n, V, power, and alpha must be NULL")
if (!is.null(V) && min(V) < 0)
stop("V must be positive")
if (!is.null(k) && min(k) < 2)
stop("number of groups must be at least 2")
if (!is.null(n) && min(n) < 2)
stop("number of observations in each group must be at least 2")
if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha > 1))
stop(sQuote("alpha"), " must be numeric in [0, 1]")
if (!is.null(power) && !is.numeric(power) || any(0 > power | power > 1))
stop(sQuote("power"), " must be numeric in [0, 1]")
p.body <- quote({
chi <- (V)^2*n*(k-1)
df <- k - 1
crit.val <- qchisq(p=1-alpha,df=df,ncp=0)
1-pchisq(crit.val,df=df,ncp=chi)
})
if (is.null(power))
power <- eval(p.body)
else if (is.null(k))
k <- uniroot(function(k) eval(p.body) - power, c(2 + 1e-10, 100))$root
else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2 + k + 1e-10, 1e+05))$root
else if (is.null(V))
V <- uniroot(function(V) eval(p.body) - power, c(1e-07, 1e+07))$root
else if (is.null(alpha))
alpha <- uniroot(function(alpha) eval(p.body) - power, c(1e-10, 1 - 1e-10))$root
else stop("internal error")
NOTE <- "n is the total sample size"
METHOD <- "One-way Analogous ANOVA with Count Data"
URL <- "http://psychstat.org/anovacount"
structure(list(k = k, n = n, V = V, alpha = alpha, power = power, note = NOTE, method = METHOD, url=URL), class = "webpower")
}
### Univariate Latent Change Score Models
wp.lcsm<-function(N=100, T=5, R=1000,
betay=0, my0=0, mys=0, varey=1, vary0=1, varys=1, vary0ys=0, alpha=0.05, ...){
#if (sum(N < 2*T)>0) stop("The sample size has to be at least 2 times of the number of occasions")
pop.model <- function(T){
## latent y
## Intercept
model<-"y0 =~ 1*y1\n"
## path from y(t-1) to y(t) with path 1
for (i in 2:T){
model<-paste(model, "y",i,"~1*y",(i-1),"\n", sep="")
}
## loading from dy(t) to y(t) with path 1
for (i in 2:T){
model<-paste(model, "dy",i,"=~1*y",i,"\n", sep="")
}
## path from y(t) to dy(t+1) with path betay
for (i in 2:T){
model<-paste(model, "dy",i,"~", betay, "*y", (i-1), "\n", sep="")
}
## latent slope ys factor model
for (i in 2:T){
model<-paste(model, "ys=~1*dy", i, "\n", sep="")
}
## variance for dy constraints to 0
for (i in 2:T){
model<-paste(model, "dy",i,"~~0*dy",i,"\n", sep="")
}
## variance for y constraints to 0
for (i in 1:T){
model<-paste(model, "y",i,"~~0*y",i,"\n", sep="")
}
## variance and covariance for intercept and slope
model<-paste(model, "ys~~", vary0ys, "*y0\n", sep="")
model<-paste(model, "y0~~", vary0, "*y0\n", sep="")
model<-paste(model, "ys~~", varys, "*ys\n", sep="")
model<-paste(model, "ys~", mys, "*1\n", sep="")
model<-paste(model, "y0~", my0, "*1\n", sep="")
## constrain means of y and dy to be zero
for (i in 1:T){
model<-paste(model, "y",i,"~0*1\n", sep="")
}
for (i in 2:T){
model<-paste(model, "dy",i,"~0*1\n", sep="")
}
## for observed data part
## y(t) to Y(t)
for (i in 1:T){
model<-paste(model, "y",i,"=~1*", "Y",i, "\n", sep="")
}
## set means of Y to be zero
for (i in 1:T){
model<-paste(model, "Y",i, "~0*1\n", sep="")
}
## set the variance for Y
for (i in 1:T){
model<-paste(model, "Y",i, "~~", varey, "*", "Y",i, "\n", sep="")
}
model
}
fit.model <- function(T){
## latent y
## Intercept
model<-"y0 =~ 1*y1\n"
## path from y(t-1) to y(t) with path 1
for (i in 2:T){
model<-paste(model, "y",i,"~1*y",(i-1),"\n", sep="")
}
## loading from dy(t) to y(t) with path 1
for (i in 2:T){
model<-paste(model, "dy",i,"=~1*y",i,"\n", sep="")
}
## path from y(t) to dy(t+1) with path betay
for (i in 2:T){
model<-paste(model, "dy",i,"~start(", betay, ")*y", (i-1)," + betay*y", (i-1), "\n", sep="")
}
## latent slope ys factor model
for (i in 2:T){
model<-paste(model, "ys=~1*dy", i, "\n", sep="")
}
## variance for dy constraints to 0
for (i in 2:T){
model<-paste(model, "dy",i,"~~0*dy",i,"\n", sep="")
}
## variance for y constraints to 0
for (i in 1:T){
model<-paste(model, "y",i,"~~0*y",i,"\n", sep="")
}
## variance and covariance for intercept and slope
model<-paste(model, "ys~~start(", vary0ys, ")*y0 + vary0ys*y0\n", sep="")
model<-paste(model, "y0~~start(", vary0, ")*y0 + vary0*y0\n", sep="")
model<-paste(model, "ys~~start(", varys, ")*ys + v