/* mmrmT module for SAS/IML An implementation of the method of Lu, et al. (2008) in SAS/IML. Lu K, Luo X, Chen PY. Sample size estimation for repeated measures analysis in randomized clinical trials with missing data. International Journal of Biostatistics 2008; 4(1): Article 9. Copyright (c) 2014 Kengo NAGASHIMA This software is released under the MIT License, see . */ proc iml; start mmrmT( Ra, pa, sa, Rb, pb, sb, lambda, delta, side, alpha, power, test, print ); if Rb = . then Rb = Ra; if pb = . then pb = pa; if sb = . then sb = sa; if lambda = . then lambda = 1; if side = . then side = 2; if test = . then test = 1; if print = . then print = 1; Ia = J(nrow(Ra), nrow(Ra), 0); do j = 1 to nrow(Ra) - 1; Ia = Ia + (pa[j] - pa[j + 1])* block(inv(Ra[1:j, 1:j]), J(nrow(Ra) - j, nrow(Ra) - j, 0)); end; Ia = Ia + pa[nrow(Ra)]*inv(Ra); Ib = J(nrow(Rb), nrow(Rb), 0); do j = 1 to nrow(Rb) - 1; Ib = Ib + (pb[j] - pb[j + 1])* block(inv(Rb[1:j, 1:j]), J(nrow(Rb) - j, nrow(Rb) - j, 0)); end; Ib = Ib + pb[nrow(Rb)]*inv(Rb); phia = inv(Ia)[nrow(Ra), nrow(Ra)]; phib = inv(Ib)[nrow(Rb), nrow(Rb)]; sigma = (sa + sb)*0.5; Nas = (1 + lambda * phib/phia)* (quantile("Normal", alpha/side) + quantile("Normal", 1-power))**2* sigma**2/delta**2 - 1; if test = 1 then do; print "Sample Size for MMRM (t-test)"; tpower = 0; do while(tpower < power); Nas = Nas + 1; Nbs = Nas*phia/(lambda*phib); Nu = Nas + Nbs - 2; SE = sqrt(sigma**2 * (1/Nas + 1/Nbs)); t = quantile("t", 1 - alpha/side, Nu); if side = 1 then tpower = 1 - cdf("t", t, Nu, delta/SE); else tpower = 1 - cdf("t", t, Nu, delta/SE) + cdf("t", -t, Nu, delta/SE); end; Nah = Nas; Nal = Nas - 1; do i = 1 to 50; Nas = (Nal + Nah)*0.5; Nbs = Nas*phia/(lambda*phib); Nu = Nas + Nbs - 2; SE = sqrt(sigma**2 * (1/Nas + 1/Nbs)); t = quantile("t", 1 - alpha/side, Nu); if side = 1 then tpower = 1 - cdf("t", t, Nu, delta/SE); else tpower = 1 - cdf("t", t, Nu, delta/SE) + cdf("t", -t, Nu, delta/SE); if tpower - power < 0 then Nal = Nas; else Nah = Nas; end; power = tpower; end; else do; Nas = Nas + 1; Nbs = Nas*phia/(lambda*phib); print "Sample Size for MMRM (z-test)"; end; Na = phia*Nas; Nb = phib*Nbs; Results = Na//Nb//lambda//alpha//power//delta//sigma//side; if print = 1 then print Results[rowname={ "Na =" "Nb =" "lambda =" "alpha =" "power =" "delta =" "sigma =" "side =" "test ="} label = ""]; return(Results); finish; store module = mmrmT; quit; /* Example */ proc iml; load module = mmrmT; /* mmrmT(Ra, pa, sa, Rb, pb, sb, lambda, delta, side, alpha, power, test, print); Ra: Correlation matrix for group A pa: Retention in group A sa: Standard deviation of group A Rb: Correlation matrix for group B (if not specified, Rb = Ra) pb: Retention in group B (if not specified, pb = pa) sb: Standard deviation of group B (if not specified, sb = sa) delta: Effect size lambda: Allcation ratio (default: 1) side: 1 = one-sided test, 2 = two-sided test (default) alpha: Type I error rate power: Power test: 0 = z-test, 1 = t-test (default) print: 0 = not print result, 1 = print result (default) */ Ra = { 1.00 0.25 0.25 0.25, 0.25 1.00 0.25 0.25, 0.25 0.25 1.00 0.25, 0.25 0.25 0.25 1.00 }; pa = { 1, 0.90, 0.80, 0.70 }; sa = 1; delta = 0.5; lambda = 1; N = mmrmT(Ra, pa, sa, ., ., ., lambda, delta, 2, 0.05, 0.8, 0, 1); N = mmrmT(Ra, pa, sa, ., ., ., lambda, delta, 2, 0.05, 0.8, 1, 1); lambda = 2; N = mmrmT(Ra, pa, sa, ., ., ., lambda, delta, 2, 0.05, 0.8, 1, 1); quit;