power.mmrm.ar1 | R Documentation |
This function performs the sample size calculation for a mixed model of repeated measures with AR(1) correlation structure. See Lu, Luo, & Chen (2008) for parameter definitions and other details.
power.mmrm.ar1(
N = NULL,
rho = NULL,
ra = NULL,
sigmaa = NULL,
rb = NULL,
sigmab = NULL,
lambda = 1,
times = 1:length(ra),
delta = NULL,
sig.level = 0.05,
power = NULL,
alternative = c("two.sided", "one.sided"),
tol = .Machine$double.eps^2
)
N |
total sample size |
rho |
AR(1) correlation parameter |
ra |
retention in group a |
sigmaa |
standard deviation of observation of interest in group a |
rb |
retention in group a (assumed same as |
sigmab |
standard deviation of observation of interest in group b. If
NULL, |
lambda |
allocation ratio |
times |
observation times |
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
, lmmpower
,
diggle.linear.power
# reproduce Table 2 from Lu, Luo, & Chen (2008)
tab <- c()
for(J in c(2,4))
for(aJ in (1:4)/10)
for(p1J in c(0, c(1, 3, 5, 7, 9)/10)){
rJ <- 1-aJ
r <- seq(1, rJ, length = J)
# p1J = p^(J-1)
tab <- c(tab, power.mmrm.ar1(rho = p1J^(1/(J-1)), ra = r, sigmaa = 1,
lambda = 1, times = 1:J,
delta = 1, sig.level = 0.05, power = 0.80)$phi1)
}
matrix(tab, ncol = 6, byrow = TRUE)
# approximate simulation results from Table 5 from Lu, Luo, & Chen (2008)
ra <- c(100, 76, 63, 52)/100
rb <- c(100, 87, 81, 78)/100
power.mmrm.ar1(rho=0.6, ra=ra, sigmaa=1, rb = rb,
lambda = sqrt(1.25/1.75), power = 0.904, delta = 0.9)
power.mmrm.ar1(rho=0.6, ra=ra, sigmaa=1, rb = rb,
lambda = 1.25/1.75, power = 0.910, delta = 0.9)
power.mmrm.ar1(rho=0.6, ra=ra, sigmaa=1, rb = rb,
lambda = 1, power = 0.903, delta = 0.9)
power.mmrm.ar1(rho=0.6, ra=ra, sigmaa=1, rb = rb,
lambda = 2, power = 0.904, delta = 0.9)
power.mmrm.ar1(N=81, ra=ra, sigmaa=1, rb = rb,
lambda = sqrt(1.25/1.75), power = 0.904, delta = 0.9)
power.mmrm.ar1(N=87, rho=0.6, ra=ra, sigmaa=1, rb = rb,
lambda = 1.25/1.75, power = 0.910)
power.mmrm.ar1(N=80, rho=0.6, ra=ra, sigmaa=1, rb = rb,
lambda = 1, delta = 0.9)
power.mmrm.ar1(N=84, rho=0.6, ra=ra, sigmaa=1, rb = rb,
lambda = 2, power = 0.904, delta = 0.9, sig.level = NULL)
# Extracting paramaters from gls objects with AR1 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.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.