Nothing
#' Sample size calculations for linear mixed models of rate of change based on
#' lmer, lme, or gee "placebo" pilot estimates.
#'
#' These functions compute sample size for linear mixed models based on the
#' formula due to Diggle (2002) or Liu and Liang (1997). These formulae are
#' expressed in terms of marginal model or Generalized Estimating Equations
#' (GEE) parameters. These functions translate pilot mixed effect model
#' parameters (e.g. random intercept and/or slope, fixed effects, etc.) into
#' marginal model parameters so that either formula can be applied to
#' equivalent affect. Pilot estimates are assumed to be from an appropriate
#' "placebo" group and the parameter of interest is assumed to be the rate of
#' change over time of the outcome.
#'
#' Any parameters not explicitly stated are extracted from the fitted
#' \code{object}.
#'
#' @name lmmpower
#' @aliases lmmpower-methods lmmpower,ANY-method lmmpower,merMod-method
#' lmmpower.default lmmpower.lme lmmpower.gee lmmpower.numeric lmmpower.double
#' @docType methods
#' @param object an object returned by lme4
#' @param n sample size per group
#' of a mixed-effects model object to placebo data assumed to have either a
#' random intercept, or a random intercept and random effect for time (slope);
#' and fixed effect representing the rate of change in a placebo group.
#' @param parameter the name or position
#' of the rate of change parameter of interest, e.g. (\code{"time"},
#' \code{"t"}, or \code{2} if it is the second specified fixed effect).
#' @param pct.change the percent change
#' in the pilot estimate of the parameter of interest (\code{beta}, the
#' placebo/null effect)
#' @param delta the change in the pilot estimate
#' of the parameter of interest, computed from \code{pct.change} if left
#' missing.
#' @param t vector of time points
#' @param sig.level Type I error
#' @param power power
#' @param alternative \code{"two.sided"} or \code{"one.sided"}
#' @param beta pilot estimate of the placebo
#' effect (slope or rate of change in the outcome)
#' @param beta.CI 95\% confidence limits of
#' the pilot estimate of beta
#' @param delta.CI 95\% confidence limits of
#' the effect size
#' @param sig2.i pilot estimate of variance
#' of random intercept
#' @param sig2.s pilot estimate of variance
#' of random slope
#' @param sig2.e pilot estimate of residual
#' variance
#' @param cov.s.i pilot estimate of
#' covariance of random slope and intercept
#' @param cor.s.i pilot estimate of
#' correlation of random slope and intercept
#' @param R pilot estimate of a marginal
#' model working correlation matrix
#' @param p proportion vector for both groups; if i indexes visits, p[i] = the
#' proportion whose last visit was at visit i (p sums to 1)
#' @param method the formula to use. Defaults
#' to \code{"diggle"} for Diggle et al (2002). Alternatively \code{"liuliang"}
#' can be selected for Liu & Liang (1997), \code{"edland"} for Ard & Edland (2011),
#' or \code{"hu"} for Hu, Mackey & Thomas (2021).
#' @param tol numerical tolerance used in root finding.
#' @param ... other arguments
#' @return An object of class \code{power.htest} giving the calculated sample
#' size, N, per group and other parameters.
#' @author Michael C. Donohue
#' @seealso \code{\link{liu.liang.linear.power}},
#' \code{\link{diggle.linear.power}}, \code{\link{edland.linear.power}},
#' \code{\link{hu.mackey.thomas.linear.power}}
#' @references Diggle P.J., Heagerty P.J., Liang K., Zeger S.L. (2002)
#' \emph{Analysis of longitudinal data}. Second Edition. Oxford Statistical
#' Science Series.
#'
#' Liu, G., and Liang, K. Y. (1997) Sample size calculations for studies with
#' correlated observations. \emph{Biometrics}, 53(3), 937-47.
#'
#' Ard, C. and Edland, S.D. (2011) Power calculations for clinical trials in Alzheimer's disease.
#' \emph{Journal of Alzheimer's Disease.} 21:369-377.
#'
#' Hu, N., Mackey, H., & Thomas, R. (2021). Power and sample size
#' for random coefficient regression models in randomized experiments with
#' monotone missing data. \emph{Biometrical Journal}, 63(4), 806-824.
#'
#' @keywords power sample size mixed effects random effects marginal model
#' methods
#' @examples
#'
#' \dontrun{
#' browseVignettes(package = "longpower")
#' }
#'
#' lmmpower(delta=1.5, t = seq(0,1.5,0.25),
#' sig2.i = 55, sig2.s = 24, sig2.e = 10, cov.s.i=0.8*sqrt(55)*sqrt(24), power = 0.80)
#' lmmpower(n=208, t = seq(0,1.5,0.25),
#' sig2.i = 55, sig2.s = 24, sig2.e = 10, cov.s.i=0.8*sqrt(55)*sqrt(24), power = 0.80)
#' lmmpower(beta = 5, pct.change = 0.30, t = seq(0,1.5,0.25),
#' sig2.i = 55, sig2.s = 24, sig2.e = 10, cov.s.i=0.8*sqrt(55)*sqrt(24), power = 0.80)
#'
#' \dontrun{
#' library(lme4)
#' fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
#' lmmpower(fm1, pct.change = 0.30, t = seq(0,9,1), power = 0.80)
#'
#' library(nlme)
#' fm2 <- lme(Reaction ~ Days, random=~Days|Subject, sleepstudy)
#' lmmpower(fm2, pct.change = 0.30, t = seq(0,9,1), power = 0.80)
#'
#' # random intercept only
#' fm3 <- lme(Reaction ~ Days, random=~1|Subject, sleepstudy)
#' lmmpower(fm3, pct.change = 0.30, t = seq(0,9,1), power = 0.80)
#'
#' library(gee)
#' fm4 <- gee(Reaction ~ Days, id = Subject,
#' data = sleepstudy,
#' corstr = "exchangeable")
#' lmmpower(fm4, pct.change = 0.30, t = seq(0,9,1), power = 0.80)
#' }
#'
#' @method lmmpower default
#' @export
lmmpower.default <- function(object=NULL,
n=NULL,
parameter = 2,
pct.change = NULL,
delta = NULL,
t = NULL,
sig.level = 0.05,
power = NULL,
alternative = c("two.sided", "one.sided"),
beta=NULL,
beta.CI=NULL,
delta.CI=NULL,
sig2.i=NULL,
sig2.s=NULL,
sig2.e=NULL,
cov.s.i=NULL,
cor.s.i=NULL,
R=NULL,
p=NULL,
method = c("diggle", "liuliang", "edland", "hu"),
tol = .Machine$double.eps^2,
...)
{
if(sum(!sapply(list(delta, pct.change), is.null))==2)
stop("Only one of delta and pct.change must be specified.")
if(is.null(delta)&!is.null(beta)&!is.null(pct.change))
delta<-pct.change*beta
if (sum(sapply(list(n, delta, power, sig.level), is.null)) != 1)
stop("exactly one of 'n', 'delta', 'power', and 'sig.level' must be NULL")
if (!is.null(sig.level) && !is.numeric(sig.level) || any(0 >
sig.level | sig.level > 1))
stop("'sig.level' must be numeric in [0, 1]")
alternative <- match.arg(alternative)
method <- match.arg(method)
m <- length(t)
if(is.null(R) & method %in% c("diggle", "liuliang")){
D <- matrix(c(sig2.i, cov.s.i, cov.s.i, sig2.s), nrow=2)
R <- cbind(1,t)%*%D%*%rbind(1,t)
R <- R + diag(sig2.e, m, m)
}
if(method == "liuliang"){
u <- list(u1 = t, u2 = rep(0,m))
v <- list(v1 = cbind(1,1,t),
v2 = cbind(1,0,t))
if(!is.null(n)) N <- n*2 else N <- NULL
}
results <- switch(method,
hu = hu.mackey.thomas.linear.power(n=n, delta=delta, t=t,
sig2.i=sig2.i, cor.s.i=cor.s.i, sig2.s=sig2.s, sig2.e=sig2.e,
p=p, sig.level=sig.level, power=power,
alternative=alternative,tol=tol,...),
edland = edland.linear.power(n=n, delta=delta, t=t,
sig2.int=sig2.i, sig2.s=sig2.s, sig2.e=sig2.e,
sig.level=sig.level,
power=power,
alternative=alternative,tol=tol,...),
diggle = diggle.linear.power(n=n, delta=delta, t=t, R=R,
sig.level=sig.level,
power=power,
alternative=alternative,tol=tol,...),
liuliang = liu.liang.linear.power(N=N, delta=delta, u=u, v=v, R=R,
sig.level=sig.level,
power=power,
alternative=alternative,tol=tol,...))
if(is.null(delta.CI)&!is.null(beta.CI)) results$delta.CI <- (results$delta/beta)*beta.CI
if(!is.null(beta)) results$beta <- beta
if(!is.null(beta.CI)) results$beta.CI <- beta.CI
if(!is.null(results$delta.CI)){
n.upper <- switch(method,
hu = hu.mackey.thomas.linear.power(n=NULL, results$delta.CI[1], t=t,
sig2.i=sig2.i, cor.s.i=cor.s.i, sig2.s=sig2.s, sig2.e=sig2.e,
p=p, sig.level=sig.level, power=power,
alternative=alternative, tol=tol,...)$n,
edland = edland.linear.power(n=NULL, results$delta.CI[1], t=t,
sig2.int=sig2.i, sig2.s=sig2.s, sig2.e=sig2.e,
sig.level=sig.level,
power=power,
alternative=alternative,tol=tol,...)$n,
diggle = diggle.linear.power(n=NULL, results$delta.CI[1], t=t, R=R,
sig.level=sig.level,
power=power,
alternative=alternative,tol=tol,...)$n,
liuliang = liu.liang.linear.power(N=NULL, results$delta.CI[1], u=u, v=v, R=R,
sig.level=sig.level,
power=power,tol=tol,...)$n)
n.lower <- switch(method,
hu = hu.mackey.thomas.linear.power(n=NULL, results$delta.CI[2], t=t,
sig2.i=sig2.i, cor.s.i=cor.s.i, sig2.s=sig2.s, sig2.e=sig2.e,
p=p, sig.level=sig.level, power=power,
alternative=alternative,tol=tol,...)$n,
edland = edland.linear.power(n=NULL, results$delta.CI[2], t=t,
sig2.int=sig2.i, sig2.s=sig2.s, sig2.e=sig2.e,
sig.level=sig.level, power=power,
alternative=alternative,tol=tol,...)$n,
diggle = diggle.linear.power(n=NULL, results$delta.CI[2], t=t, R=R,
sig.level=sig.level,
power=power,tol=tol,...)$n,
liuliang = liu.liang.linear.power(N=NULL, results$delta.CI[2], u=u, v=v, R=R,
sig.level=sig.level,
power=power,tol=tol,...)$n)
n.CI <- c(n.lower, n.upper)
if(n.CI[1]>n.CI[2]) n.CI <- n.CI[2:1]
results$n.CI <- n.CI
}
if(is.character(parameter)){
names(results)[names(results) == "beta"] <- parameter
names(results)[names(results) == "beta.CI"] <- paste(parameter, "CI")
}
results <- results[unlist(lapply(results, function(x) !is.null(x)))]
structure(results, class = "power.longtest")
}
#' @export
#' @method lmmpower double
lmmpower.double <- lmmpower.default
#' @export
#' @method lmmpower numeric
lmmpower.numeric <- lmmpower.default
#' @export
lmmpower <- function(object, ...) UseMethod("lmmpower")
setGeneric("lmmpower")
#' @importFrom nlme getVarCov
#' @method lmmpower lme
#' @export
lmmpower.lme <- function(object,
n = NULL,
parameter = 2,
pct.change = NULL,
delta = NULL,
t = NULL,
sig.level = 0.05,
power = NULL,
alternative = c("two.sided", "one.sided"),
beta=NULL,
beta.CI=NULL,
delta.CI=NULL,
sig2.i=NULL,
sig2.s=NULL,
sig2.e=NULL,
cov.s.i=NULL,
cor.s.i=NULL,
p=NULL,
method = c("diggle", "liuliang", "edland", "hu"),
tol = .Machine$double.eps^2,
...)
{
alternative <- match.arg(alternative)
method <- match.arg(method)
if(is.numeric(parameter)) parameter <- rownames(summary(object)$tTable)[parameter]
tab <- getVarCov(object)
if(nrow(tab)>2) stop("Too many random effects. Function is
equipped to handle at most a random intercept and slope.")
if(is.null(beta))
beta = summary(object)$tTable[parameter,'Value']
if(is.null(beta.CI))
beta.CI = rep(summary(object)$tTable[parameter,'Value'],2) +
c(-1,1)*qnorm(0.025)*summary(object)$tTable[parameter,'Std.Error']
if(beta.CI[1] > beta.CI[2]) beta.CI <- beta.CI[2:1]
# var of random intercept
if(is.null(sig2.i))
sig2.i = tab["(Intercept)", "(Intercept)"]
# var of random slope
if(is.null(sig2.s))
sig2.s = ifelse(nrow(tab)==1, 0,
ifelse(nrow(tab)==2, tab[2, 2], NA))
# residual var
if(is.null(sig2.e))
sig2.e = object$sigma^2
# covariance of slope and intercep
if(is.null(cov.s.i))
cov.s.i = ifelse(nrow(tab)==1, 0,
ifelse(nrow(tab)==2, tab[1, 2], NA))
# cor of slope and intercep
if(is.null(cor.s.i))
cor.s.i = ifelse(nrow(tab)==1, 0,
ifelse(nrow(tab)==2, tab[1, 2]/(sqrt(tab[1,1]) * sqrt(tab[2,2])), NA))
lmmpower(object=NULL,
n = n,
parameter = parameter,
pct.change = pct.change,
delta = delta,
t = t,
sig.level = sig.level,
power = power,
alternative = alternative,
beta = beta,
beta.CI = beta.CI,
delta.CI = delta.CI,
sig2.i = sig2.i,
sig2.s = sig2.s,
sig2.e = sig2.e,
cov.s.i = cov.s.i,
cor.s.i = cor.s.i,
p = p,
method = method,
tol=tol, ...)
}
#' @method lmmpower gee
#' @export
lmmpower.gee <- function(object,
n = NULL,
parameter = 2,
pct.change = NULL,
delta = NULL,
t = NULL,
sig.level = 0.05,
power = NULL,
alternative = c("two.sided", "one.sided"),
beta=NULL,
beta.CI=NULL,
delta.CI=NULL,
p = p,
method = c("diggle", "liuliang", "edland", "hu"),
tol = .Machine$double.eps^2,
...)
{
alternative <- match.arg(alternative)
method <- match.arg(method)
if(is.numeric(parameter)) parameter <- rownames(summary(object)$coefficients)[parameter]
if(is.null(beta))
beta = summary(object)$coefficients[parameter,'Estimate']
if(is.null(beta.CI))
beta.CI = rep(summary(object)$coefficients[parameter,'Estimate'],2) +
c(-1,1)*qnorm(0.025)*summary(object)$coefficients[parameter,'Robust S.E.']
if(beta.CI[1] > beta.CI[2]) beta.CI <- beta.CI[2:1]
R <- summary(object)$working.correlation * summary(object)$scale
lmmpower(object=NULL,
n = n,
parameter = parameter,
pct.change = pct.change,
sig.level = sig.level,
power = power,
alternative = alternative,
sig2.e=1,
sig2.s=0,
beta=beta,
beta.CI=beta.CI,
delta.CI=delta.CI,
R=R,
t=t,
p=p,
method=method,
tol=tol, ...)
}
#' @importFrom lme4 VarCorr fixef getME
#' @export
setMethod("lmmpower", signature(object = "merMod"),
function(object,
n = NULL,
parameter = 2,
pct.change = NULL,
delta = NULL,
t = NULL,
sig.level = 0.05,
power = NULL,
alternative = c("two.sided", "one.sided"),
beta=NULL,
beta.CI=NULL,
delta.CI=NULL,
sig2.i=NULL,
sig2.s=NULL,
sig2.e=NULL,
cov.s.i=NULL,
cor.s.i=NULL,
p=NULL,
method = c("diggle", "liuliang", "edland", "hu"),
tol = .Machine$double.eps^2,
...)
{
if (!(sum(sapply(list(n, delta, power, sig.level), is.null)) == 1 |
sum(sapply(list(n, pct.change, power, sig.level), is.null)) == 1))
stop("exactly one of 'n', 'delta' (or 'pct.change'), 'power', and 'sig.level' must be NULL")
if (!is.null(sig.level) && !is.numeric(sig.level) || any(0 >
sig.level | sig.level > 1))
stop("'sig.level' must be numeric in [0, 1]")
alternative <- match.arg(alternative)
method <- match.arg(method)
if(is.numeric(parameter)) parameter <- rownames(coef(summary(object)))[parameter]
tab <- VarCorr(object)
if(length(tab)>1) stop("Too many grouping levels. Function is
equipped to handle at most one grouping level.")
tab <- tab[[1]]
if(nrow(tab)>3) stop("Too many random effects. Function is
equipped to handle at most a random intercept and slope.")
m <- length(t)
if(is.null(beta))
beta = fixef(object)[parameter][[1]]
if(is.null(beta.CI))
beta.CI = rep(coef(summary(object))[parameter,"Estimate"],2) +
c(-1,1)*qnorm(0.025)*coef(summary(object))[parameter,"Std. Error"]
if(beta.CI[1] > beta.CI[2]) beta.CI <- beta.CI[2:1]
# var of random intercept
if(is.null(sig2.i))
sig2.i = tab[1, 1]
# var of random slope
if(is.null(sig2.s))
sig2.s = ifelse(nrow(tab)==1, 0,
ifelse(nrow(tab)==2, tab[2 ,2], NA))
# residual var
if(is.null(sig2.e))
sig2.e = getME(object, "sigma")^2
# covariance of slope and intercept
if(is.null(cov.s.i))
cov.s.i = ifelse(nrow(tab)==1, 0,
ifelse(nrow(tab)==2, tab[2, 1], NA))
# cor of slope and intercept
if(is.null(cor.s.i))
cor.s.i = ifelse(nrow(tab)==1, 0,
ifelse(nrow(tab)==2, tab[2, 1]/(sqrt(tab[1, 1]) * sqrt(tab[2, 2])), NA))
lmmpower(n=n, object=NULL,
parameter = parameter,
pct.change = pct.change,
delta = delta,
t = t,
sig.level = sig.level,
power = power,
alternative = alternative,
beta=beta,
beta.CI=beta.CI,
delta.CI=delta.CI,
sig2.i=sig2.i,
sig2.s=sig2.s,
sig2.e=sig2.e,
cov.s.i=cov.s.i,
cor.s.i=cor.s.i,
p=p,
method = method,
tol=tol, ...)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.