*! version 1.0)

help trialcounts

Title

trialcounts -- trialcounts is a program for calculating the number of subjects recruited to,
and the number having reached subsequent scheduled visits, in a clinical trial based on supplied linear piecewise
recruitment functions, whilst additionally accounting for potential dropout.
From these are derived the 'cohort weights'
which can be used in the mixedpower program in order to calculate power or sample size for linear mixed models,
accounting for the recruitment and dropout functions established in trialcounts at a given point in trial time.

Syntax

trialcounts
,
schedule(numlist)
ends(numlist)
rates(numlist)
[options]

* schedule(numlist) list of visit times for the proposed study

* ends(varname|numlist) list of timepoints ends of each piecewise (linear) recruitment function

* rates(varname numlist|numlist) list of recruitment rates for each piecewise (linear) function

timelist(numlist) single time value to evaluate recruitment function; or timelist to search over when using search

quadratic set the first piecewise recruitment function to be quadratic

maxn(#) maximum number to be recruited to the trial, overrriding recruitment function

search(# #) search for the timepoint that reaches target #_1 at visit #_2

dropouts(numlist) list of proportions of those who only reached visit k, and no further, due to dropout

drfunction(drf_type, p(#) [l(#) k(#) s(#)]) define dropout function using a parametric survival model

rgraph request a plot of the recruitment function

dgraph request a plot each of the parametric survival and hazard dropout functions

nolines suppress red count lines in recruitment function

rgopts(twoway_options) recruitment function plot options, except graph name

dhgopts(twoway_options) dropout hazard function plot options, except graph name

dsgopts(twoway_options) dropout survival function plot options, except graph name

rename(name_option) recruitment function plot name option

dhname(name_option) dropout hazard function plot name option

dsname(name_option) dropout survival function plot name option

compact display table of results in compact form

disp2 include the rescaled probability weights in the displayed table of results

*this option is required.

drf_type

weibull weibull dropout model

lognormal log normal dropout model

gompertz gompertz dropout model

loglogistic log logistic dropout model

ggamma generalized gamma dropout model

Description

trialcounts calculates the number of subjects recruited to a trial, as well as the number having reached subsequent
scheduled visits, based on a user-supplied recruitment function. The reported numbers are given as both the total number
having reached each visit, and the unique number who have only reached each visit but no further.

The user defines the recruitment function in
a list of piecewise linear functions in two required options; ends which defines that endpoint t for
each consecutive piecewise function and rates which defines the recruitment rates for each of those
functions i.e number recruited per unit time, regardless of length of the function (difference between ends).
The recruitment function can therefore be arbitrarily complex. The ends and rates options
can be supplied either directly as number lists or as variable names containing the successive values.

The other required option is schedule which informs trialcounts the visit timepoints where the
counts need to be evaluated at. Another crucial option, though not required, timelist indicates at what point in
'trial time' the counts are to be calculated for. Other useful options are maxn which will override the
defined recruitment function in the sense that you can set a limit to the number of total recruited to the trial.
You can also use the search option to establish at what point in time do a target number of subjects reach
a specific visit, which is when timelist needs to be an actual list to search over. The quadratic option
can be used to easily replicate an increasing recruitment rate in its early phase.

Based on the list of counts, trialcounts also produces 'probability weights' representing the
proportion who have reached the various visit schedule points,
either out of all those who will be recruited according to the trialcounts recruitment function
('final weights' in the output) or out of all those currently recruited according to the time value used
to evaluate the recruitment function ('final weights rescaled' in the output). Either of the weight lists
may be used, according to required needs, in the mixedpower, strecruitment option to calculate power
or sample size accounting for partial subject follow-up at interim study timepoints.

trialcounts also incorporates dropout rates and integrates the dropout and recruitment functions to produce
the necessary counts and probability weights, accounting for both factors. Dropout probabilities can be supplied
as a list within the dropouts option or using the drfunction which calculates dropout probabilities
for the user based on common parametric survival models, selected according to preferential hazard shapes. Plots
can also be produced by trialcounts to visualise both dropout and recruitment function selections.

Options

Required option

schedule(numlist) specifies the visit times for the proposed trial. The supplied number list must always be increasing and
not negative. The length of the list is limited to 101 visits including baseline, though the table of counts output will not
be printed with visit length>20. Even so, the output will not be presented prettily before that number is reached. All important
information will in supplied in the returned list regardless of visit length.

ends(varname|numlist) list of timepoints that relate to the 'end' of each piecewise (linear) recruitment function. The list can be directly supplied as a number list,
allowing the usual shorthand conventions, or as a (numeric) variable name that contain the numbers therein, with relevant rows
as indexed by the rates option. However supplied, the number list can be just one value and a maximum of 120 values
but must always be increasing
and not negative. The number of values must match the number supplied in rates (guaranteed if using the varname
method for both ends and rates)
, but not necessarily that in schedule.
The first piecewise function then, spans time zero (baseline) to the first value in ends,
and the second piecewise from first value in ends to the second, and so on.

rates(varname numlist|numlist) list of recruitment rates that relate to each piecewise (linear) recruitment function.
The list can be directly supplied as a number list,
allowing the usual shorthand conventions, or as a (numeric) variable name that contain the numbers therein, with relevant rows
to include additionally indexed by a number list. The number list used for either method can be just one value and a maximum of
120 values
and must not be negative. However, only for the variable name method (where its role is as a row index) must the number list
be always increasing.
The number of values must match the number supplied in ends (guaranteed if using the varname
method for both ends and rates), but not necessarily that in schedule. The rate values
themselves should represent the number recruited per unit time (hence rate) in that piecewise function span,
regardless of the length of the span. One could even supply a long list of values if one was attempting to
replicate the actual recruitment rates of a trial through time (e.g. one rate per month so far, and perhaps supplemented with predicted
rates into the future), and this can be easily performed with the variable name method if such rates exist in a data file.

If using the variable name input method for ends then one must do so also for rates; if using the variable name
input method for rates then one may use either the variable name or number list input method for ends, assuming the length matches for both.
With the variable name method one may select non-consecutive rows to input (e.g. rates(1/5 8/10 14 20)), if
consecutive rows have the same rate. This may speed-up computation if there are many rows (say, recruitment rates for every
month of a trial) and the search option is being utilised. Also be aware that including an empty row will cause the
program to exit.

Other options

timelist(# [# ...]) despite the option name one will typically supply just one value that indicates the point in trial
time where the recruitment function will be evaluated and counts of subjects at each visit made.
However, you need not request any time value, or alternatively you may specify a needlessly large value far beyond
the length of the trial. In which case the functions
are evaluated at 'the end of the trial' which is defined as the time of the last value in ends plus the last value
in schedule i.e. time at the last visit for the last subject recruited. In which case, a warning message is given telling
the user no time value has been supplied.

When using the search option then the timelist option must be now be a list of at least 2 increasing
numbers for the program to search over. You will probably want to use a shorthand convention such as timelist(4(0.01)5).

quadratic this is to specify that the first (and only first) piecewise recruitment function is to be a quadratic, rather
than linear, function. This is to enable a convenient shortcut to approximating a smooth and increasing recruitment rate in the
early phase of the trial, rather than trying to replicate through a series of linear piecewise functions, each increasing in rates
slightly. The first value in the rates list will be be interpreted as if the coefficient of a quadratic function, and so
will typically require a smaller value than for the linear rates.

maxn(#) to specify a maximum size of the trial. This will override the recruitment function as supplied in ends
and rates. One may be able to estimate a recruitment rate stretching into the future but be unsure for how long one needs
to indicate this function i.e. what should the last value from ends be? If you just know that the trial only needs to
recruit a certain number of subjects then use maxn(#). Once the supplied recruitment function hits the maxn(#) limit,
that point becomes the time value of the last recruited subjects (and hence the last value of ends
in this restricted recruitment function). Counts and recruitment plots (including the limit of the trial) will reflect this.

search(# #) allows one to search for the point in trial time where a target number (the first #) of subjects reaches a
particular visit number (the second #). Hence the second # must not be greater than the length of schedule. If the search
option is requested then a list of increasing time values must be supplied to timelist for the program to search over.
If the time list is first specified in large intervals you could rerun the command with a more focused, granular time list to obtain an
accurate estimate of the objective time value. That being said, the trialcounts, search
option is typically fast even when searching over many values, although there is a limit of 2500 for the number of elements Stata allows in a
number list. If the supplied list is entirely below the objective time value then
trialcounts will produce a warning message to that effect and report the counts and weights at the last supplied time value.
If the supplied list is entirely above the objective time value then trialcounts will produce a warning message
to that effect and report the counts and weights at the first supplied time value.

dropouts(numlist) specifies the estimated proportion of dropouts you
expect at 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. Therefore, the sum of all proportions should not exceed one. There is a built in tolerance of 0.000001
to deal with numerical inaccuracies within Stata.

drfunction(drf_type, p(#) [l(#) k(#) s(#)]) allows the user to define in a more convenient,
and perhaps plausible, manner the dropout function compared to the direct entry of dropouts. Instead a parametric survival model
is chosen with a hazard function that reflects anticipated dropout rates over the trial history (for example increasing, decreasing,
non-monotonic). The function choices are weibull, log normal, gompertz, log logistic and generalized gamma.

All models except the generalized gamma are specified using 2 of 3
possible parameters p(#), l(#) and s(#). The letters p and l are taken from the Weibull model relating
to the shape parameter p and the rate parameter lambda. The other three models have a 'shape' or 'scale' parameter,
and a 'rate' parameter
of a sort where the linear predictor would enter in a survival regression model (although typically known with differing Greek letters
but have been given the same indications here for convenience).
The third parameter s is an additional parameter where the user can alternatively specify the
survival function probability at the last vist in the schedule list, instead of l(#) which is rarely intuitive. The
value of l is then back-calculated in order to generate the dropout probabilities. Hence s(#) must be between 0 and 1.
For the generalized gamma model, there is no s(#) as it is not possible to back-calculate the other parameters from a survival
probability. Instead the generalized gamma has an additional parameter kappa k(#)

weibull where p is the shape parameter p (>0) and l the rate parameter lambda (>0). p<1
results in an decreasing monotonic hazard function; p=1 results in a constant hazard function
i.e. exponential survival model at value of l(#); p=2 results in an positive linear hazard function;
p>1 results in an increasing monotonic hazard function

lognormal where p is the scale parameter sigma (>0) and l the 'rate' parameter mu. p>0
results in an non-monotonic hazard function that increases and then decreases; p dictates the magnitude of the
rise and fall; depending on the parameter values the decrease in hazard might not be observed for the range of the trial schedule

gompertz where p is the shape parameter gamma and l the rate parameter lambda (>0). p<0
results in an decreasing monotonic hazard function starting from value of l(#); p=0 results in a constant
hazard function i.e. exponential survival model at value of l(#); p>0 results in an increasing monotonic
hazard function starting from value of l(#)

loglogistic where p is the scale parameter gamma (>0) and l the rate parameter lambda (>0).
p>0 & <1 results in an non-monotonic hazard function that increases and then decreases; depending on the parameter
values the decrease in hazard might not be oberved for the range of the trial schedule; p>=1 results in a monotonic
decreasing hazard function

ggmamma where p is the shape parameter sigma (>0), k is another ancillary parameter
and l the 'rate' parameter mu (>0). ggmama is a flexible model, with the Weibull (k=1), exponential
(k=1 and p=1) and log-normal k=0 as special cases.

rgraph produces a graph of the recruitment function up to the final timepoint of the trial, defined
as time of last entry in ends list plus time of last entry in schedule list. The furthest right vertical
red line is for the first time value in schedule
plotted at the specified time value, and meets the line for number recruited.
Each horizontal line to the left, in turn, is distanced equal
to the time interval between visits and so corresponds to the horizontal for the counts at each visit,
evaluated at the supplied time. Red horizontal
lines hence reflect, in order from highest line to lowest, the number of subjects who have reached each visit in the
schedule list, in ascending time order, at the time value supplied in timelist.
Think of the red vertical lines moving from left to right as trial time progresses
with distances between them kept fixed, and where they meet the recruitment function reflecting the trialcounts
sums at a given timelist(#).
These red lines can be excluded with option nolines. If no one has reached a given visit that
line will not be present. If the search function is called the lines
are drawn for the time value in timelist that meets the target of the search. Additional
information regarding the recruitment function is also presented. Note, that if dropout is indicated then this will
not impact this graph which exclusively concerns recruitment, unlike the table result of counts which integrates
in the dropout function or list.

dgraph produces two graphs of the dropout survival and hazard functions specified in drfunction with
additional information concerning parameter values.

nolines suppresses the red vertical and horizontal lines depicting the visit counts, as described above in rgraph

rgopts(twoway_options) supply (and possibly override) any options allowed for twoway plots for the recruitment graph,
except for name_option in which case use rename.

dhgopts(twoway_options) supply (and possibly override) any options allowed for twoway plots for the dropout hazard graph,
except for name_option in which case use dhname.

dsgopts(twoway_options) supply (and possibly override) any options allowed for twoway plots for the dropout survival graph,
except for name_option in which case use dsname.

rename(name_option) override the default name given to the recruitment function graph.

dhname(name_option) override the default name given to the dropout hazard function graph.

dsname(name_option) override the default name given to the dropout survival function graph.

compact display the table of results in compact form. If the table is not displaying nicely into your output window
because it is too wide, this option may allow the table to fit. Vist index and time value labelling is abbreviated and the results
are presented to fewer decimal places.

disp2 additionally print in the output the list of final cohort probability weights that have been rescaled so as to sum to one.
Note, the unscaled 'final weights' may sum<1 if not every subject has been recruited at time(#).

Output

The output first provides a summary with the total number who are planned to be recruited to the trial based purely on the recruitment
function, and limited by the maxn option if necessary. Also given is the number who have have been recruited given the recruitment
function (i.e. have had at least a
baseline visit, so that dropout is not a factor) at the timepoint given by the option time or evaluated by search.

Following this is a table of the cohort probability weights and counts for each visit. The first line gives the final cohort weights, integrated
over dropout and recruitment, which may be used for power and sample size calculation in mixedpower. They reflect the probabilities
of having reached visit k, and no further, when considering an overall final sample size i.e. the first number in the summary output. The second line
(if requested with disp2) are the final cohort weights rescaled so they sum to one, and are appropriate for use in mixedpower if considering
an intermediate sample size i.e. the second number in the summary output at time(#).

Following the probability weights are the counts per visit. 'Unique counts' indicates the number having reached visit k, and no further, at
time(#), whilst 'sum counts' are the total number who have reached at least visit k at time(#), both taking dropout into account.

Examples

Specify a recruitment function with 5 piecewise linear functions and request recruitment graph. No time value specified

`". trialcounts, sched(0 0.5 1 1.5 2 2.5 3) ends(0.5 1 1.5 3 5) rates(5 7 12 16 25) rgr"'

Same recruitment function but supply a time to evaluate function. Also use disp2to the rescaled weights in table

`". trialcounts, sched(0 0.5 1 1.5 2 2.5 3) ends(0.5 1 1.5 3 5) rates(5 7 12 16 25) time(3.5) rgr disp2 "'

Use maxn to place a limit to recruitment of n=60 and integrate in a weibull dropout function with high early hazard (p=0.7)
and a 'survival' probability of 0.75 at the end of the trial. Also request dropout graphs

`". trialcounts, sched(0 0.5 1 1.5 2 2.5 3) ends(0.5 1 1.5 3 5) rates(5 7 12 16 25) time(5) rgr disp2 drfunction(weibull, p(0.7) s(0.75)) maxn(60) dgr"'

A small trial with slow exponential growth in recruitment

`". trialcounts,sched(0(1)10) rates(5(0.5)15) ends(0.5(0.5)10.5) rgr time(9)"'

`". trialcounts,sched(0(1)10) rates(5(0.5)15) ends(0.5(0.5)10.5) rgr time(9) maxn(80)"'

Use quad to indicate first recruitment piece is interpreted as a quadratic function and recreate increasing early recruitment,
followed by a static recruitment rate and giving the plot a title. Note, only 2 piecewise functions required

`". trialcounts, sched(0 0.5 1 1.5 2 2.5 3) ends(3 5) rates(40 250) time(3.5) rgr quad rgopts(title(recruitment function for new trial))"'

Use search function to find when the 5th visit is reached by 80 subjects. We try searching over years 3 to 4 in 0.01 unit
increments. At least 2 values are required to search over. Also include a lognormal dropout function

`". trialcounts, sched(0 0.5 1 1.5 2 2.5 3) ends(3 5) rates(40 250) time(3(0.01)4) quad drfunction(lognormal, p(1.3) s(0.6)) search(80 5) dgr"'

The same set-up (plus a gompertz dropout function) but we search over a range too low to reach the target.
Note the warning message

`". trialcounts, sched(0 0.5 1 1.5 2 2.5 3) ends(3 5) rates(40 250) time(3(0.01)3.5) quad drfunction(gompertz, p(-0.3) l(0.2)) search(80 5) dgr"'

The same set-up (plus a log-logistic dropout function) but we search over a range too high and surpass the target with the first time value. Note the warning message

`". trialcounts, sched(0 0.5 1 1.5 2 2.5 3) ends(3 5) rates(40 250) time(4(0.01)5) quad drfunction(loglogistic, p(0.5) s(0.75)) search(80 5) dgr"'

COMBINING MIXEDPOWER WITH TRIALCOUNTS

An advanced example looking at combining trialcounts with mixedpower, or rather in this example,
mvmixedpower and using the variable name entry method.
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.

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

`". est use mvmixed_ms.ster"'

`". mvmixedpower, trtspec(slope) sched(0 0.5 1 1.5 2 2.5 3) diff(0.05) auto multi(3) weighting(1 1 1) n(686) alpha(0.3) 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

trialcounts stores the following in r():

Scalars

r(max_time) maximum time length of the trial

r(N_limit) maximum sample size

r(N_current) number recruited evaluated at the supplied time value

r(time_value) supplied (or evaluated if using search) time value

r(fwgtsum) sum of the recruitment weights based on r(N_limit)

r(drate) dropout function rate parameter value

r(dshape) dropout function shape parameter value

r(dSurv_max) dropout function survival probability at last visit

Matrices

r(schedlist) supplied schedule list

r(origends) supplied ends list

r(origrates) supplied rates list

r(ends) effective ends list, different if maxn used

r(rates) effective rates list, different if maxn used

r(dropwgts) supplied dropout weights

r(recweights) recruitment weights based on recruitment function

r(rescalefwgts) rescaled final weights so sum to 1

r(uniqcts) the unique count making each visit

r(cumcts) the cumulative count making each visit

trialcounts also returns all the above as macros, as well as the command line,
the full timelist if using search and the dropout function name in r():

Author

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