*! version 1.0)

help mvmixedpower

Title

mvmixedpower -- mvmixedpower is a program for calculating power or sample size analytically for multivariate
linear mixed-effects
models, typically for the design of a randomised clinical trial where you wish to consider more than one outcome jointly.

See also mixedpower, dmmixedpower, trialcounts

Syntax

mvmixedpower
,
trtspec(trtspec_type)
schedule(numlist)
multi(#)
[options]

* trtspec(trtspec_type) specification of treatment effect type, either slope
or intercept

* schedule(numlist) the time values for the visit schedule

* multi(#) the number of outcomes m (dimensions)

alpha(#) significance level; default is 0.05

twosided request a properly two-sided test

power(#) power; default is 0.8, required to compute sample size

n(#) total sample size; required to compute power

difference(# [# ...]) 1 or m treatment effect size(s) in model parameter terms - alternative to effectiveness.

effectiveness(# [# ...]) 1 or m treatment effect size(s) for the slope, in
proportionate terms relative to the control group slope

contslope(# [# ...]) 1 or m mean control group slope(s) that the effectiveness option is relative to.
Not required if difference used

weighting(# # ...|ivw|nivw) the weighting used in a linear combination test of
each effect size per dimension

covariance(matrix) user input of the random effect covariance matrix

errorvar(matrix) user input of the random error variance

auto alternative specification of covariance and error variance parameters that indicates automatic input of variance estimates from a
multivariate mixed model in memory

scale(#) the ratio of the timescale used for the variance components and the timescale in schedule;
default is 1

aratio(# #) allocation ratio between groups; default is equal allocation (1 1)

dropouts(numlist) the proportion who only reached visit k of schedule, and no further, due to dropout.
Must sum to 1

drop2(numlist) the dropout proportions in the treatment group, if different.
Must sum to 1

strecruitment(numlist) the proportion who only reached visit k of schedule, and no further, due to staggered recruitment. Need not sum to 1

noheader suppress display of mvmixedpower header banner

nosyntax suppress display of the mixed model syntax that mvmixedpower is calculating power/sample size for

notable suppress display of table count per visit given for incomplete follow-up

xmat(#) return the X matrix from 'cohort' number # in the return list

rmat(#) return the R matrix from 'cohort' number # in the return list

zmat(#) return the Z matrix from 'cohort' number # in the return list

gmat(#) return the G matrix from 'cohort' number # in the return list

bvarn(#) return the covariance matrix of fixed effects based on 'cohort' number # in the return list

*this option is required.

trtspec_type

Currently only slope and intercept terms allowed
slope proportionate slope effect i.e. constrained to equal control slope at baseline

intercept intercept effect i.e. parallel shift

Description

mvmixedpower performs analytic sample size or power calculations for a proposed RCT using the asymptotic formula for the
fixed effects variance-covariance matrix of a linear mixed model: Var(B_hat)=(X'(Sigma^-1)X)^-1,
where Sigma=R+ZGZ'. The variance-covariance matrix is calculated for a nominal 2-person trial* and includes the
2-person variance* for the proposed treatment effect parameters per outcome. In mvmixedpower these are then linearly combined
using a weighting scheme of choice, with correct inference based on the relevant variance and covariance terms from Var(B_hat).
The subsequent effect size may then be used in a power/sample size calculation in the conventional
manner, suitably adjusting for the proposed or required sample size.
Further details can be found in the help file for mixedpower.
Integration with trialcounts for dropout and partial follow-up is feasible similarly to mixedpower.

* assuming allocation ratio is (1 1), otherwise covariance matrix is calculated for a N1+N2 person trial, where N1 and N2 are
the integers specified in aratio(N1 N2). If there is assumed incomplete follow-up then the covariance matrix is further
adjusted to reflect the proportions of predicted data patterns

Options

Required options

trtspec(trtspec_type) allows the user to specify different
parameterisations of the treatment effect - currently only a slope or intercept effect allowed.
The parameterisation must be the same for all m
joint outcomes, even if the magnitude of the treatment effects need not be.

slope
is for a proportionate slope effect implying that the control and treatment group slopes have a common overall intercept.

intercept is for
an intercept effect implying a constant parallel difference between control and treatment groups.

schedule(numlist) specifies the visit times for the proposed trial. A baseline
visit at time 0 is not assumed and must be given if required. The supplied number list must always be increasing and
not negative.

multi(#) the number m of joint outcomes (dimensions). Must be of at least dimension 2

Basic options

alpha(#) significance level of the test; default is 0.05. Be aware that when you specify alpha(0.05)
you are in fact obtaining results for a 0.025 one-sided test, unless specifically
requesting twosided.

twosided request that sample size or power is calculated for a properly two-sided test, including power to detect a difference
in the 'opposite' direction to treatment 'improvement'. This will be denoted in the output with an asterisk next to 'alpha'.
The default is a one-sided test with type I error of alpha/2. Usually there will be no substantive difference between the two,
but that difference
will not be trace if alpha is high and/or the effect size is small.

power(#) power; default is 0.8, required to compute sample size. Sample size will be rounded up to nearest even
number and a similar principle is applied when aratio(# #) is used, see option aratio. Note, the exact
'fractional' sample size is also returned as a scalar.

n(#) total sample size; required to compute power. If the provided number is odd, the actual sample size used will be
rounded down to nearest even number. A similar principle is applied when aratio(# #) is used, see option aratio.
Note, the exact power corresponding to the rounded n, is printed as output with target power returned as a scalar as 'fractional' power
- corresponding to 'fractional' sample size.

Effect size and testing options

difference(# [# ...]) the magnitude of the m treatment effects, specified as the values of the relevant regression parameters.
Hence for trtspec(slope) the value supplied is the absolute change of the treatment group relative to control group per unit time per m.
One can either give m treatment differences as a numer list or a single value that applies for all m outcomes. This may be the
a reasonable option when all outcomes have been standardised or are naturally commensurate.

effectiveness(# [# ...]) for the trtspec(slope) selection only, instead of using difference
the user may alternatively specify
the treatment magnitudes in proportionate terms relative to the control group slope.
Like difference, one can provide either m effectiveness
values. 1 per outcome, or a single value applying to all m outcomes. The value(s) supplied should be a real
number >0 and <=1. For example, effectiveness(0.3 0.3 0.3)
for a 30% change towards a null slope for each outcome.
This option also requires use of contslope to represent the control group slopes the
effect is relative to. If you wish to specify an effect that i) results in an effect in the opposite
direction (positive to negative slope or vice-versa) or ii) increase (decreases) further an already positive (negative)
slope then you will need to use difference.

contslope(# [# ...]) this option is to be used in conjunction with effectiveness
and represents the mean control group
slopes that the effectiveness option is relative to, entered as real numbers. One may either
provide m different control slopes per outcome, or a single control slope to be repeated for all m outcomes,
regardless of whether
one or m values were supplied in effectiveness. Not required if difference used.

weighting(# # ...|ivw|nivw) the weighting used in a linear combination test of each slope or intercept effect per dimension.
The number of elements of weighting must be the same as number of dimensions m entered in option multi.
Also allowed are the text entries ivw indicating inverse-variance weighting and nivw indicating 'negative' inverse-variance
weighting. The latter
is necessary if some of the outcomes trend in opposite directions, meaning the effect sizes will typically be given
with opposite signs in difference or assumed with effectiveness. The text selection nivw acknowledges these
differentially signed effects as all indicating improvement whilst respecting the signs of the off-diagonals in the covariance matrix.

Variance options

covariance(matrix) the covariance matrix of the random effects, entered as a 2mx2m symmetrical matrix where
m is the number of joint outcomes. The entries are ordered at the highest level by outcome, with intercept and slope variance per
outcome then entered in 2x2 blocks, with all the other off-diagonal possibly non-zero if desired, to induce correlation between outcomes at the
random effect level.

See option auto for a shortcut to using covariance and errorvar together.

errorvar(matrix) the error variance terms, entered as a mxm matrix,
where m is the number of joint outcomes, and error covariances across outcomes may be non-zero.

See option auto for a shortcut to using covariance and errorvar together.

auto use of this option will automatically transfer the random effect and error (co-)variance values from a suitable
mixed model in memory, and is an auto-input alternative to using covariance and errorvar. It is the user's
responsibility to ensure the model is 'suitable'.

With mvmixedpower the user must be careful how the multivariate mixed model in memory
has been fitted: specifically the m random intercept terms are all named before the m random slope terms
after the first double piping. For example if m=3:

. mixed depvar subscale_1 subscale_2 subscale_3 c.time_1 c.time_2 c.time_3 , nocons || id_level2:
subscale_1 subscale_2 subscale_3 time_1 time_2 time_3, cov(unstr) nocons || visit: , nocons residuals(unstr, t(subscale))

then the three subscale_# terms after the first double piping are the three outcome-specific random intercepts, given before the 3
outcome-specific random coefficients time_#, along with the , nocons option. Note, these variables
have also been given as fixed effects, again with the , nocons option. The covariance structure between outcomes can be
any of those options allowed by mixed and auto will enter the 6x6 (in this example)
covariance matrix of random effects.
At a lower level of hierarchy within visit
the residuals have also been modelled here as unstructured
across the three outcomes and auto will enter the 3x3 error matrix. The auto
option also accepts independent errors but no other residual
structure allowed by the command mixed.

Adjustment options

scale(#) use this option if you wish to specify visit times in a different scale to that relating to the
inputted random effect (co-)variance parameters.
For example, if your variance parameter estimates came from a model where time was specified in years,
but you wished to calculate power for a trial that had a visit schedule specified in months you could use scale(`=1/12').
Scale-dependent effect sizes (slope-based) should be kept unchanged in difference - the fundamental aspects
(schedule and effect size) are specified in the scale you wish, but variance parameters from a source with a different scale
are rescaled by scale.

aratio(# #) specifies relative group sizes (allocation ratio). aratio requires two integers reflecting the control group
to treatment group allocation ratio. The default is aratio(1 1) i.e equal group size, though there is no limit to the
maximum value of either #.
When calculating required sample size,
the smallest n achieving target power is calculated that results in integer group sizes.
Hence actual power may be slightly more than if 'fractional'
subjects were allowed, which could meet target power exactly. When calculating power, the actual sample size used is the largest
n <= n(#) that gives integer group sizes. Hence, actual power may be slightly increased if any 'spare' samples were allocated
to one of the groups. Note, this rule is always followed, so if aratio(3 3) was chosen, then depending on n(#),
power could be slightly less than aratio(1 1) even though the allocation ratio itself is equivalent,
as actual n may need to be smaller to ensure integer group sizes.
Because you can enter any integer value, aratio(# #) can also be used to accurately calculate power after recruitment is finished,
knowing that in fact the allocation was not perfectly balanced. For example, enter aratio(152 147).
Note 'fractional' power and 'fractional' sample size
are also returned in the stored results.

dropouts(numlist) specifies the estimated proportion of dropouts out of the original sample you
expect following each study visit. It must correspond exactly to the length of the schedule list.
Each number in the list should be between 0 and 1 and is the proportion who only reached visit k, and no further,
due to dropout. To be clear, the proportion is not out of those remaining at visit k but of those eligible for full
follow-up. Therefore one does not to consider the impact of strecuitment when specifying the dropout list values.
Furthermore, the sum of all proportions should be equal to one, as it makes no sense to have subjects
withdraw before the trial begins. There is a built-in tolerance of 0.00001
to deal with numerical inaccuracies within Stata.

drop2(numlist) specifies a set of dropout proportions specifically for the treatment group, if it is believed they
will be different from the control. Rules are exactly as for dropouts.

strecruitment(numlist) specifies the proportion of the original sample who have only reached as far as visit k,
and no further,
due to the effect of staggered recruitment to the trial. Each number in the list should be between 0 and 1 and
the sum of all proportions should not exceed one. There is a built-in tolerance of 0.00001
to deal with numerical inaccuracies within Stata. Sums less than one are allowed however, and for staggered recruitment make sense, if one
wishes to calculate power for an interim stage of a trial with an established final sample size, but at this interim point not everyone
has yet been randomised. If both strecruitment and dropouts are specified then the integrated cohort proportions
(weights) are also returned in the stored results. Use of companion program trialcounts
is recommended as an accessible way to generate appropriate strecruitment proportions.

It is also worth mentioning that
strecruitment is more relevant for calculating power rather than sample size, as it reflects recruitment constraints. Used for
a sample size calculation, it implies that needing a bigger sample can be achieved by simply recruiting more in a given timeframe. If
that is not the case, then use trialcounts to work with mvmixedpower to find when a sufficient sample size can be achieved.

Reporting options

noheader prevents display of the mvmixedpower header banner.

nosyntax prevents the default display of the mixed model syntax for the requested model that mvmixedpower has calculated
power or sample size for. All terms in italics are placeholder names for variables
names that the user would supply when fitting that model. The syntax given may not be the only way to fit a particular model, but is
a recommended parameterisation when matching with mvmixedpower and, in particular, for utilising the auto-input methods.

notable prevents display of the table of counts by visit, otherwise produced when follow-up is incomplete.
If either dropouts or
strecuitment is supplied indicating that not all subjects have all been recorded for all outcome measures, then a table is
produced that lists the (rounded) number of subjects who have been observed for each visit by treatment arm (so different if drop2 or aratio meaningully supplied).

xmat(#) returns the X design matrix for a particular 'cohort' identified by visit number # from the schedule list. This is only
possible for cohorts with non-zero probability weightings. This means that if both dropouts and strecruitment are not used then
one can only request the X matrix appropriate for the last visit (number #) in schedule where all individuals have a full complement
of measures. Note r(xmat) is not automatically returned as there is potential for the matrix to be large. When displaying r(xmat) there
are rownames to help identify rows. The tags include 't' for (treatment) group (0 for control, 1 for treatment),
'v' for visit (the visit index number in the schedule list) and 'm' for outcome, in order.
The columns are not labelled other than X1, X2... but terms are order at the highest level by
outcome and then within outcome all control term variables come
before the treatment term variable, and
intercepts come before slopes.
One may use xmat(#), rmat(#), zmat(#), gmat(#) and bvarn(#) together
and # can be different for all options, assuming the selections are individually permissible.

rmat(#) returns the R error matrix for a particular 'cohort' identified by visit number # from the schedule list. All information
given in xmat applies here, except for column identification, where instead column matches row. Displaying r(rmat) may be
particularly useful to ensure the order of error variances and covariances is correct when applying auto.

zmat(#) returns the Z design matrix of random effects for a particular 'cohort' identified by visit number #
from the schedule list. All information
given in xmat applies here, except for column identification, where instead pairs of
columns are represent the random intercept and
slope - in that order - per outcome, firstly for the control and then treatment group.

gmat(#) returns the G covariance matrix of random effects for a particular 'cohort' identified by
visit number # from the schedule list. The matrix r(gmat) has no row- or column-name information other than G1, G2..
and a covariance-block
of random effects is given for the control then treatment group.
Within that covariance-block paired random effects are ordered by outcome. For the first outcome pair
the random intercept is the first variance term on the diagonal followed by the
random slope variance, and similarly for subsequent outcomes. Displaying r(gmat) may be
particularly useful to ensure the order of random effect variances and covariances is correct when applying auto.

bvarn(#) returns the covariance matrix of the estimated fixed parameters calculated for a particular 'cohort'
identified by visit number # from the schedule list. This is only possible for cohorts with non-zero probability weightings.
Note that the option bvarn returns not just matrix r(bvar_n), but also the working correlation matrix r(wcorr_n), and
the working
covariance matrix r(sigma_n). For r(bvar_n) there is
no row- or column-name information, but the order of variables is exactly as described in option xmat(#).

Examples

An example of mvmixedpower for use with multivariate mixed model taken from a real interim analysis plan for a trial in Parkinson's where the 3 joint
outcomes are Parts I, II and III of the MDS-UPDRS, although we have ignored the incomplete follow-up and higer alpha here.
This is hard work with user-entry and easy to make a mistake.
The mvmixedpower command calculates power for a common effectiveness of 30% reduction in slope to all 3 control slopes.
Inverse variance weighting has been used as the linear combination of treatment effects

`". mat err=(6.107172,2.194946,2.105618\2.194946,7.459157,5.858216\2.105618,5.858216,41.01714)"'

`". mat cov=(9.67,0.34,7.01, 0.56,5.23,1.06\0.34,0.63,0.52,0.58,0.3,0.52\7.01,0.52,15.06,0.43,13.62,0.36\0.56,0.58,0.43,1.12,0.16,	1.51\5.23,0.3,13.62,0.16,70.21,-3.01\1.06,0.52,0.36,1.51,-3.01,4.49)"'

`". mvmixedpower , trtspec(slope) sched(0 0.25 0.5(0.5)3) conts(0.865 0.975 1.876) eff(0.3) cov(cov) error(err) m(3) weighting(ivw) n(800)"'

It is much quicker and less error-prone to use auto-entry, as done here. Firstly we load a previously saved estimation results file from a
multivariate mixed model and replay
the results. The slight discrepancy in power compared to the user-entry version is because it was
necessary to reduce the precision of the covariance entries for the input version for it to be displayed properly in the help file

`". estimates use mvmixed_pd.ster"'

`". mixed"'

`". mvmixedpower, trtspec(slope) sched(0 0.25 0.5(0.5)3) conts(0.865 0.975 1.876) eff(0.3) auto m(3) weighting(ivw) n(800)"'

We now look at combining mvmixedpower with trialcounts with an advanced example.
We first load a dataset containing recruitment rates for a real MS trial.
The variable recr2 contains the recruitment rate per month for consecutive 2 month intervals
in the control arm. The first 10 entries are actual monthly rates,
with the subsequent 17 entries being
best estimates into the future. The interim analysis is to be held at 38 months so we really only need to consider the
first 19 entries. However, for exposition we will include all rows in the file, but use the time option to limit the calculations
at month 38. Doing this also allows us to show how one could reduce the number of defined ends and rates by using entering the
number list rates(recr2 1(1)21 27), so excluding rows 22-26 as the rates here are identical to 27. In a different example, this may
reduce the computational burden, especially if combining with search. We also enter the name of the variable containing the ends.

`". use ms_recruitment, clear"'

We will run trialcounts with months as the unit of time and the command
will evalute the supplied recruitment function at 38 months. In addition, we integrate in a Weibull
dropout function that implies 90% will still be in the trial by the last visit (at 36 months). The variance parameters
for a power calculation are loaded from a
saved multivariate mixed model, which we will auto-input into mvmixedpower. We use the rescaled weights
from the trialcounts output as the entries in the strecruitment list because we defined a recruitment function that resulted in
planned sample size greater than the sample size at the time of the interim analysis. By evaluating the local
rescaled weights macro (r(fw_rescale)) returned by trialcounts one can
enable seamless .do file use. We supply the 3 control slopes from the fitted model into the contslope option, together with a relative
effect size of 30% improvement for each scale. That is, the 3 slopes are lessened in gradient by 40% towards the null. We wish to combine the
3 effect sizes with IVW but because the first scale increases with disability progression and the other two decrease, we actually
indicate the option
weighting(nivw) so as to respect the different coefficient signs indicating improvement. As this is an interim analysis looking for
a signal of IMP activity we tolerate a high (one-sided) alpha in order to maintain high power.

`". est use mvmixed_ms.ster"'

`". mixed"'

`". trialcounts, sched(0(6)36) ends(month) rates(recr2 1(1)21 27) time(38) drfunction(weibull, p(0.75) s(0.9)) disp2 rgr dgr"'

`". mvmixedpower, trtspec(slope) sched(0 0.5 1 1.5 2 2.5 3) contslope(.1075532 -.1258883 -.0793978) eff(0.4) auto multi(3) weighting(nivw) n(686) alpha(0.6) dropout(`= r(fw_rescale)')"'

Note, even though use of mvmixedpower will overwrite the r-class returned list from
trialcounts, mixedpower
also returns the (integrated) weightings with the same macro names (either r(final_wgts) or r(fw_rescale)) so one may
continue entering the same mixedpower command when altering some other parameter of the mixed model. Because the mixed model was fitted in years we will return to years as the time unit for the schedule list and
difference effect sizes, which will be linearly combined with equal weighting. Note, we indicate the sample size is 343*2=686 to reflect both groups now.

Stored results

mvmixedpower stores the following in r():

Scalars

r(fractional_ss) 'fractional' sample size without rounding

r(samplesize) calculated or supplied sample size

r(fractional_power) power without rounding of supplied sample size

r(power) calculated or supplied power

r(effectsize) the treatment effect standardised by the variance=r(wgt_trteff)/r(var_trteff)

r(wgt_trteff) the overall (weighted) treatment effect, linearly combined from m outcome effects

r(var_trteff) the variance of r(wgt_trteff) for a (N1+N2)-person trial adjusted for incomplete follow-up

r(se_trial) the standard error of r(trt_eff) for the trial reflecting final sample size

r(scale) scale used

r(alpha) type I error rate (significance level)

Macros

r(cmdline) command line

r(cmd) command

r(schedule) schedule list

r(final_wgts) integrated probability weights of dropout and staggered recruitment

r(final_wgts2) integrated probability weights for treatment group, if different

r(fw_rescale) rescaled integrated weights to sum to 1

r(fw_rescale2) rescaled integrated weights for treatment group, if different

r(cont_spec) model parameterisation of control group

r(trt_spec) model parameterisation of treatment group

r(twosided) if calculation based on properly two-sided test or not

r(aratio) allocation ratio

r(trialcounts_cmd) trialcounts command line, if previous rclass command

r(fr_current) current frame at command

Matrices

r(sched_list) schedule list

r(cont_num) control numbers reaching each visit

r(trt_num) treatment numbers reaching each visit

r(contrd_num) rounded control numbers reaching each visit

r(trtrd_num) rounded treatment numbers reaching each visit

r(m_weights) weights used for linear combination of m treatment effects

r(trtbeta_true) true treatment group beta coefficients

r(beta_var) covariance matrix of fixed effects

r(wcorr_n) working correlation matrix if requested, for a given cohort

r(sigma_n) working covariance matrix if requested, for a given cohort

r(bvar_n) covariance matrix of fixed effects if requested, for a given cohort

r(xmat) X matrix if requested, for a given cohort

r(rmat) R matrix if requested, for a given cohort

r(zmat) Z matrix if requested, for a given cohort

r(gmat) G matrix if requested, for a given cohort

r(cohort_wgts) integrated probability weights of dropout and staggered recruitment

r(cohort_wgts2) integrated probability weights for treatment group, if different

r(st_rec) staggered recruitment probability weights

r(dropouts) dropout probability weights

r(dropouts2) dropout probability weights for treatment group, if different

References

Chris Frost, Michael G. Kenward, Nick C. Fox. Optimizing the design of clinical
trials where the outcome is a rate. Can estimating a baseline rate in a run-in
period increase efficiency? Statist. Med. 2008; 27:3717–3731 doi: 10.1002/sim.3280

Author

Matthew Burnell
MRC Centre of Research Excellence in Clinical Trial Innovation
University College London
London, UK
m.burnell@ucl.ac.uk