Linear mixed model sample size calculations.

Description

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.

Usage

1
2
3
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)

Arguments

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, sigmab is assumed same as sigmaa. If not NULL, sigmaa and sigmab are averaged.

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.

Details

See Lu, Luo, & Chen (2008).

Value

The number of subject required per arm to attain the specified power given sig.level and the other parameter estimates.

Author(s)

Michael C. Donohue

References

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)

See Also

power.mmrm.ar1, lmmpower, diggle.linear.power

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
# 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)