power.mmrm | R Documentation |
This function performs the sample size calculation for a mixed model of repeated measures with general correlation structure. See Lu, Luo, & Chen (2008) for parameter definitions and other details. This function executes Formula (3) on page 4.
power.mmrm(
N = NULL,
Ra = NULL,
ra = NULL,
sigmaa = NULL,
Rb = NULL,
rb = NULL,
sigmab = NULL,
lambda = 1,
delta = NULL,
sig.level = 0.05,
power = NULL,
alternative = c("two.sided", "one.sided"),
tol = .Machine$double.eps^2
)
N |
total sample size |
Ra |
correlation matrix for group a |
ra |
retention in group a |
sigmaa |
standard deviation of observation of interest in group a |
Rb |
correlation matrix for group a |
rb |
retention in group b |
sigmab |
standard deviation of observation of interest in group b. If
NULL, |
lambda |
allocation ratio |
delta |
effect size |
sig.level |
type one error |
power |
power |
alternative |
one- or two-sided test |
tol |
numerical tolerance used in root finding. |
See Lu, Luo, & Chen (2008).
The number of subject required per arm to attain the specified
power
given sig.level
and the other parameter estimates.
Michael C. Donohue
Lu, K., Luo, X., Chen, P.-Y. (2008) Sample size estimation for repeated measures analysis in randomized clinical trials with missing data. International Journal of Biostatistics, 4, (1)
power.mmrm.ar1
, lmmpower
,
diggle.linear.power
# reproduce Table 1 from Lu, Luo, & Chen (2008)
phi1 <- c(rep(1, 6), 2, 2)
phi2 <- c(1, 1, rep(2, 6))
lambda <- c(1, 2, sqrt(1/2), 1/2, 1, 2, 1, 2)
ztest <- ttest1 <- c()
for(i in 1:8){
Na <- (phi1[i] + lambda[i] * phi2[i])*(qnorm(0.05/2) + qnorm(1-0.90))^2*(0.5^-2)
Nb <- Na/lambda[i]
ztest <- c(ztest, Na + Nb)
v <- Na + Nb - 2
Na <- (phi1[i] + lambda[i] * phi2[i])*(qt(0.05/2, df = v) + qt(1-0.90, df = v))^2*(0.5^-2)
Nb <- Na/lambda[i]
ttest1 <- c(ttest1, Na + Nb)
}
data.frame(phi1, phi2, lambda, ztest, ttest1)
Ra <- matrix(0.25, nrow = 4, ncol = 4)
diag(Ra) <- 1
ra <- c(1, 0.90, 0.80, 0.70)
sigmaa <- 1
power.mmrm(Ra = Ra, ra = ra, sigmaa = sigmaa, delta = 0.5, power = 0.80)
power.mmrm(N = 174, Ra = Ra, ra = ra, sigmaa = sigmaa, delta = 0.5)
power.mmrm(N = 174, Ra = Ra, ra = ra, sigmaa = sigmaa, power = 0.80)
power.mmrm(Ra = Ra, ra = ra, sigmaa = sigmaa, delta = 0.5, power = 0.80, lambda = 2)
power.mmrm(N = 174, Ra = Ra, ra = ra, sigmaa = sigmaa, delta = 0.5, lambda = 2)
power.mmrm(N = 174, Ra = Ra, ra = ra, sigmaa = sigmaa, power = 0.80, lambda = 2)
# Extracting paramaters from gls objects with general correlation
# Create time index:
Orthodont$t.index <- as.numeric(factor(Orthodont$age, levels = c(8, 10, 12, 14)))
with(Orthodont, table(t.index, age))
fmOrth.corSym <- gls( distance ~ Sex * I(age - 11),
Orthodont,
correlation = corSymm(form = ~ t.index | Subject),
weights = varIdent(form = ~ 1 | age) )
summary(fmOrth.corSym)$tTable
C <- corMatrix(fmOrth.corSym$modelStruct$corStruct)[[1]]
sigmaa <- fmOrth.corSym$sigma *
coef(fmOrth.corSym$modelStruct$varStruct, unconstrained = FALSE)['14']
ra <- seq(1,0.80,length=nrow(C))
power.mmrm(N=100, Ra = C, ra = ra, sigmaa = sigmaa, power = 0.80)
# Extracting paramaters from gls objects with compound symmetric correlation
fmOrth.corCompSymm <- gls( distance ~ Sex * I(age - 11),
Orthodont,
correlation = corCompSymm(form = ~ t.index | Subject),
weights = varIdent(form = ~ 1 | age) )
summary(fmOrth.corCompSymm)$tTable
C <- corMatrix(fmOrth.corCompSymm$modelStruct$corStruct)[[1]]
sigmaa <- fmOrth.corCompSymm$sigma *
coef(fmOrth.corCompSymm$modelStruct$varStruct, unconstrained = FALSE)['14']
ra <- seq(1,0.80,length=nrow(C))
power.mmrm(N=100, Ra = C, ra = ra, sigmaa = sigmaa, power = 0.80)
# Extracting paramaters from gls objects with AR1 correlation
fmOrth.corAR1 <- gls( distance ~ Sex * I(age - 11),
Orthodont,
correlation = corAR1(form = ~ t.index | Subject),
weights = varIdent(form = ~ 1 | age) )
summary(fmOrth.corAR1)$tTable
C <- corMatrix(fmOrth.corAR1$modelStruct$corStruct)[[1]]
sigmaa <- fmOrth.corAR1$sigma *
coef(fmOrth.corAR1$modelStruct$varStruct, unconstrained = FALSE)['14']
ra <- seq(1,0.80,length=nrow(C))
power.mmrm(N=100, Ra = C, ra = ra, sigmaa = sigmaa, power = 0.80)
power.mmrm.ar1(N=100, rho = C[1,2], ra = ra, sigmaa = sigmaa, power = 0.80)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.