R/webpower.R

Defines functions print.webpower plot.webpower wp.anova wp.prop wp.one.prop wp.two.prop wp.two.prop.two.n wp.t wp.t1 wp.t2 wp.rmanova wp.effect.CRT2arm wp.effect.CRT3arm wp.effect.MRT2arm wp.effect.MRT3arm nuniroot wp.crt2arm wp.crt3arm wp.mrt2arm wp.mrt3arm wp.regression wp.logistic wp.mediation estCRT2arm estCRT3arm estMRT2arm estMRT3arm wp.correlation wp.poisson wp.anova.binary wp.anova.count wp.lcsm plot.lcs.power wp.blcsm summary.power wp.popPar wp.mc.sem.basic wp.mc.sem.boot wp.mc.sem.power.curve wp.mc.t

Documented in estCRT2arm estCRT3arm estMRT2arm estMRT3arm nuniroot plot.lcs.power plot.webpower print.webpower summary.power wp.anova wp.anova.binary wp.anova.count wp.blcsm wp.correlation wp.crt2arm wp.crt3arm wp.effect.CRT2arm wp.effect.CRT3arm wp.effect.MRT2arm wp.effect.MRT3arm wp.lcsm wp.logistic wp.mc.sem.basic wp.mc.sem.boot wp.mc.sem.power.curve wp.mc.t wp.mediation wp.mrt2arm wp.mrt3arm wp.poisson wp.popPar wp.prop wp.regression wp.rmanova wp.t

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) {
    dat <- read.table(file, header = TRUE)
    
    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) {
    dat <- read.table(file, header = TRUE)
    
    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) {
    dat <- read.table(file, header = TRUE)
    
    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) {
    dat <- read.table(file, header = TRUE)
    
    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) {
    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.") 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")) {
    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, c(2 + 1e-10, 1000))$root else if (is.null(n)) 
        n <- nuniroot(function(n) eval(p.body) - power, c(1, 1e+06))$root else if (is.null(f)) 
        f <- nuniroot(function(f) eval(p.body) - power, c(1e-07, 1e+07))$root else if (is.null(icc)) 
        icc <- nuniroot(function(icc) eval(p.body) - power, c(0, 1))$root else if (is.null(alpha)) 
        alpha <- nuniroot(function(alpha) eval(p.body) - power, c(1e-10, 
            1 - 1e-10))$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")) {
    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, c(3 + 1e-10, 1000))$root else if (is.null(n)) 
        n <- nuniroot(function(n) eval(p.body) - power, c(2 + 1e-10, 1e+06))$root else if (is.null(f)) 
        f <- nuniroot(function(f) eval(p.body) - power, c(1e-07, 1e+07))$root else if (is.null(icc)) 
        icc <- nuniroot(function(icc) eval(p.body) - power, c(0, 1))$root else if (is.null(alpha)) 
        alpha <- nuniroot(function(alpha) eval(p.body) - power, c(1e-10, 
            1 - 1e-10))$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")) {
    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, c(1 + 1e-10, 1000))$root else if (is.null(n)) 
        n <- nuniroot(function(n) eval(p.body) - power, c(3 - 1e-10, 1e+06))$root else if (is.null(f) & type == 1) 
        f <- nuniroot(function(f) eval(p.body) - power, c(1e-07, 1e+07))$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")) {
    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, c(2 - 1e-10, 1000))$root else if (is.null(n)) 
        n <- nuniroot(function(n) eval(p.body) - power, c(3 - 1e-10, 1e+07))$root else if (is.null(f1) && type != 2) 
        f1 <- nuniroot(function(f1) eval(p.body) - power, c(1e-07, 1e+07))$root else if (is.null(f2) && type != 1) 
        f2 <- nuniroot(function(f2) eval(p.body) - power, c(1e-07, 1e+07))$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) {
    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
        lambda <- f2 * (u + v + 1)
        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) {
    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, c(1e-10, 1e+07))$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, c(-10, 10))$root else if (is.null(varx)) 
        varx <- nuniroot(function(varx) eval(p.body) - power, c(1e-10, 
            1e+07))$root else if (is.null(vary)) 
        vary <- nuniroot(function(vary) eval(p.body) - power, c(1e-10, 
            1e+07))$root else if (is.null(varm)) 
        varm <- nuniroot(function(varm) eval(p.body) - power, c(1e-10, 
            1e+07))$root else if (is.null(alpha)) 
        alpha <- nuniroot(function(alpha) eval(p.body) - power, c(1e-10, 
            1 - 1e-10))$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) {
		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") {
        # bernoulli
        
        
        if (is.null(parameter)) {
            B <- 0.5
        } else {
            B <- parameter
        }
        
        
        beta1 <- log(exp1)  ##beta1 under H_A in the population 
        beta0 <- log(exp0)
        a <- (1 - B) * exp(beta0) + B * exp(beta0 + beta1)  ##00
        b <- B * exp(beta0 + beta1)  #10
        c <- B * exp(beta0 + beta1)  #11
        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") {
        ## exponential predictor
        if (is.null(parameter)) {
            lambda <- 1
        } else {
            lambda <- parameter
        }
        
        beta1 <- log(exp1)  ##beta1 under H_A in the population 
        beta0 <- log(exp0)
        d <- integrate(function(x) exp(beta0 + (beta1 - lambda) * x) * 
            lambda, 0, Inf, subdivisions = 100L)$value  ##00
        e <- integrate(function(x) x * exp(beta0 + (beta1 - lambda) * x) * 
            lambda, 0, Inf, subdivisions = 100L)$value  #01
        f <- integrate(function(x) x^2 * exp(beta0 + (beta1 - lambda) * 
            x) * lambda, 0, Inf, subdivisions = 100L)$value  #11
        
        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)  ##beta1 under H_A in the population 
        beta0 <- log(exp0)
        d <- integrate(function(x) exp(beta0 + beta1 * x) * dlnorm(x, mu, 
            sigma), 0, Inf, subdivisions = 100L)$value
        e <- integrate(function(x) x * exp(beta0 + beta1 * x) * dlnorm(x, 
            mu, sigma), 0, Inf, subdivisions = 100L)$value
        f <- integrate(function(x) x^2 * exp(beta0 + beta1 * x) * dlnorm(x, 
            mu, sigma), 0, Inf, subdivisions = 100L)$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)  ##beta1 under H_A in the population 
        beta0 <- log(exp0)
        
        d <- integrate(function(x) exp(beta0 + beta1 * x) * dnorm(x, mu, 
            sigma), -Inf, Inf, subdivisions = 100L)$value
        e <- integrate(function(x) x * exp(beta0 + beta1 * x) * dnorm(x, 
            mu, sigma), -Inf, Inf, subdivisions = 100L)$value
        f <- integrate(function(x) x^2 * exp(beta0 + beta1 * x) * dnorm(x, 
            mu, sigma), -Inf, Inf, subdivisions = 100L)$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 == "Poisson") {
        
        if (is.null(parameter)) {
            lambda <- 1
        } else {
            lambda <- parameter
        }
        
        beta1 <- log(exp1)  ##beta1 under H_A in the population 
        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)  ##beta1 under H_A in the population 
        beta0 <- log(exp0)
        
        
        d <- integrate(function(x) exp(beta0 + beta1 * x)/(R - L), L, R, 
            subdivisions = 100L)$value
        e <- integrate(function(x) x * exp(beta0 + beta1 * x)/(R - L), 
            L, R, subdivisions = 100L)$value
        f <- integrate(function(x) x^2 * exp(beta0 + beta1 * x)/(R - L), 
            L, R, subdivisions = 100L)$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, paremeter = parameter, 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 + varys*ys\n", sep="")
		
		model<-paste(model, "ys~start(", mys, ")*1 + label('mys')*1\n", sep="")	
		model<-paste(model, "y0~start(", my0, ")*1 + label('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, "~~start(", varey, ")*", "Y", i, " + varey*Y",i, "\n", sep="")		
		}
		model
	}
	
	sem.est <- function(model, data){
		temp.res <- sem(model=model, data=data)
		label <- temp.res@ParTable$label
		c(temp.res@ParTable$est[label!=""], temp.res@ParTable$se[label!=""])
	}
	## do it once for a given N and T
	fit.once <- function(N, T){
		## generate data
		pop.model.T <- pop.model(T)
		pop.model.T.res <- sem(pop.model.T, do.fit=FALSE)
		pop.model.T.cov <- inspect(pop.model.T.res, "cov.ov")
		pop.model.T.mean <- inspect(pop.model.T.res, "mean.ov")
		ynames <- row.names(pop.model.T.cov)
		gen.data <- lapply(1:R, mvrnorm, n=N, mu=pop.model.T.mean, Sigma=pop.model.T.cov)
		
		## conduct the analysis
		fit.model.T <- fit.model(T)
		fit.res <- lapply(gen.data, sem.est, model=fit.model.T)
		
		## run once to get the model information
		model.info.res <- sem(fit.model.T, gen.data[[1]])
		label <- model.info.res@ParTable$label
		label <- label[label!=""]
		label.unique <- !duplicated(label)
		label <- label[label.unique]
		npar <- length(label)
		## get the parameter estimates, sd, se, power, CI of power
		all.res <- do.call(rbind, fit.res)
		all.res <- all.res[, c(label.unique, label.unique)]
		all.res <- na.omit(all.res)
		mc.est <- colMeans(all.res[, 1:npar], na.rm=TRUE)
		mc.se <- apply(all.res[, (npar+1):(2*npar)], 2, mean, na.rm=TRUE)
		mc.sd <- apply(all.res[, 1:npar], 2, sd, na.rm=TRUE)
		
		mc.z.score <- all.res[, 1:npar]/all.res[, (npar+1):(2*npar)]
		mc.z.score.check <- abs(mc.z.score) >= qnorm(1-alpha/2)
		
		mc.power <- colMeans(mc.z.score.check, na.rm=TRUE)
		
		pop.par <- unlist(lapply(label, function(x){eval(parse(text=x))}))
		mc.output <- cbind(pop.par, mc.est, mc.sd, mc.se, mc.power, N, T)
		row.names(mc.output) <- label
		label.sort <- sort(label)
		mc.output[label.sort, ]
	}
	
	if (length(N)>1 | length(T)>1){
		all.output <- list()
		for (i in N){
			for (j in T){
				all.output [[paste('N',i,'-T',j, sep="")]]<- fit.once(i,j)
			}
		}
	}else{
		all.output <- fit.once(N,T)
	}
	class(all.output) <- "lcs.power"
	all.output 	
}

plot.lcs.power <- function(x, parameter, ...){
	## x is the output from power analysis
	power.mat <- do.call('rbind', x)
	power.par <- power.mat[rownames(power.mat)==parameter, ]
	unique.N <- unique(power.par[ ,6])
	unique.T <- unique(power.par[ ,7])
	
	if (length(unique.N)==1 & length(unique.T)==1) stop("Multiple N or T is needed for power plot.")
	
	if (length(unique.N)==1){
		## plot the power along T
		plot(power.par[, 7], power.par[, 5], type='l', xlab='Number of Occasions', ylab=paste('Power of ',parameter, sep=""), ylim=c(0,1))
		points(power.par[, 7], power.par[, 5])
	}
	
	if (