How do implement gsdesign usermethod with simulations?

[Andrew Xavier]

Dear all,

Please could you help with gsdesign usermethod?

The manual says that "Your program can be as complicated as you would like; you can even use simulations to compute your results."

How can this be done? It is not clear to me. Below is some code I found in the manual to demonstrate the power command with simulations. But then I am unclear how to then implement this with gsdesign. The simulations are based by defining first the sample size so the power can be computed. However the premise of gsdesign is to specify a power so the sample size for a group sequential design is defined given the looks and alpha.

Can you help me how to implement gsdesign usermethod with simulations as noted in the manual?

capture program drop simttest
program simttest, rclass
    version 18
    syntax, n(integer)          ///  Sample size
          [ alpha(real 0.05)    ///  Alpha level
            m0(real 0)          ///  Mean under the null
            ma(real 1)          ///  Mean under the alternative
            sd(real 1) ]        //   Standard deviation
    drawnorm y, n(`n') means(`ma') sds(`sd') clear
    ttest y = `m0'
    return scalar reject = (r(p)<`alpha') 
end

capture program drop power_cmd_simttest
program power_cmd_simttest, rclass
    version 18
    // DEFINE THE INPUT PARAMETERS AND THEIR DEFAULT VALUES
    syntax, n(integer)          ///  Sample size
          [ alpha(real 0.05)    ///  Alpha level
            m0(real 0)          ///  Mean under the null
            ma(real 1)          ///  Mean under the alternative
            sd(real 1)          ///  Standard deviation
            reps(integer 100)]  //   Number of repetitions
                           
    // GENERATE THE RANDOM DATA AND TEST THE NULL HYPOTHESIS
    quietly simulate reject=r(reject), reps(`reps'):    ///
                     simttest, n(`n') m0(`m0') ma(`ma') sd(`sd') alpha(`alpha')
    quietly summarize reject
         
    // RETURN RESULTS
    return scalar power = r(mean)
    return scalar N = `n'
    return scalar alpha = `alpha'
    return scalar m0 = `m0'
    return scalar ma = `ma'
    return scalar sd = `sd'
end

power simttest, n(100) m0(70) ma(75) sd(15)

gsdesign simttest, n(100)

Thank you

Andrew

[Andrew Xavier]

I found this solution provided by stata support:


Simulation-based sample size calculation for group sequential designs is an advanced use that requires some additional considerations. First a little background: internally, the -gsdesign- command runs the
-gsbounds- command to calculate stopping boundaries and the -power- command to calculate the required sample size for an equivalently powered fixed study design (see
https://eur01.safelinks.protection.o...%3D&reserved=0 for details).

The -power- command can be a user-written sample-size evaluator, but it
*must* calculate sample size, not power. However, simulation is used to estimate the power of a test at a given sample size -- simulation cannot directly estimate the required sample size to attain a desired power.
This introduces another step when you want to use simulation to calculate sample size.

There are two ways you can address this issue. The first option is to guess-and-check various sample sizes until you find one that yields the desired power -- this gives you the fixed-study sample size. Then use
-gsbounds- to calculate stopping boundaries, and calculate sample sizes at interim looks yourself.

The second option is to write a program that automates the guess-and-check procedure. This user-written sample-size evaluator will accept the desired power as input and return the required sample size as output. This takes more programming, but the resulting sample-size evaluator can be used directly with -gsdesign-.

********************************************************************************
//program simtest, //unchanged from your 18 Jan email 
capture program drop simttest 
program simttest, rclass
    version 18
    syntax, n(integer)          ///  Sample size
          [ alpha(real 0.05)    ///  Alpha level
            m0(real 0)          ///  Mean under the null
            ma(real 1)          ///  Mean under the alternative
            sd(real 1) ]        //   Standard deviation
    drawnorm y, n(`n') means(`ma') sds(`sd') clear
    ttest y = `m0'
    return scalar reject = (r(p)<`alpha') 
    
end

    
//program power_cmd_simttest, //unchanged from your 18 Jan email 
capture program drop power_cmd_simttest 
program power_cmd_simttest, rclass
    version 18
    // DEFINE THE INPUT PARAMETERS AND THEIR DEFAULT VALUES
    syntax, n(integer)          ///  Sample size
          [ alpha(real 0.05)    ///  Alpha level
            m0(real 0)          ///  Mean under the null
            ma(real 1)          ///  Mean under the alternative
            sd(real 1)          ///  Standard deviation
            reps(integer 100)]  //   Number of repetitions

    // GENERATE THE RANDOM DATA AND TEST THE NULL HYPOTHESIS
    quietly simulate reject=r(reject), reps(`reps'):    ///
                     simttest, n(`n') m0(`m0') ma(`ma') sd(`sd') alpha(`alpha')
    quietly summarize reject

    // RETURN RESULTS
    return scalar power = r(mean)
    return scalar N = `n'
    return scalar alpha = `alpha'
    return scalar m0 = `m0'
    return scalar ma = `ma'
    return scalar sd = `sd'
end

************************************************************
// Option 1: guess-and-check method (easier to implement)
************************************************************

// Say we want 90% power. We will guess-and-check various sample sizes until // we find one that attains about 90% power. Then we will use -gsbounds- // to calculate stopping boundaries. For more information, see // https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.stata.com%2Fmanuals%2Fadaptgsbounds.pdf&data=05%7C02%7Candre.lopes%40ucl.ac.uk%7C9f115a63c9dd4326b56a08dc1ce9b347%7C1faf88fea9984c5b93c9210a11d9a5c2%7C0%7C0%7C638417037583591953%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=JxT8%2BAsjA79JZlY1sVHbcKZq2wNtHfuR0R%2Bp6neDuKc%3D&reserved=0

// set the random seed for reproducibility 
set seed 123

// start with N=100. Use 1000 reps for better accuracy 

power simttest, n(100) m0(70) ma(75) sd(15) reps(1000)

// We attained a power of 91.8%, so we decrease the sample size and try again 

power simttest, n(95) m0(70) ma(75) sd(15) reps(1000)

// This yields 90.4% power, which is close enough for me. Knowing that the 
// sample required by an equivalent fixed design is 95 subjects, we can 
// use -gsbounds- to calculate stopping boundaries. Say we want 4 evenly-spaced 
// looks with a two-sided O'Brien-Fleming-style error-spending efficacy bound.
gsbounds, nlooks(4) efficacy(errobfleming) graphbounds

// The maximum sample required by the group sequential trial is equal to the 
// fixed-study sample size multiplied by the information ratio, which is 
// stored as r(info_ratio) 
scalar Nmax = 95 * r(info_ratio)

// Calculate N at each look using stored matrix r(info_frac) 
matrix sampsize = Nmax * r(info_frac) 
matrix rownames sampsize = "Sample:size"

// Transpose the vector of per-look sample sizes and concatenate it with the 
// matrix of stopping boundaries, which is stored as r(bounds) 
matrix boundmat = (r(bounds), sampsize') 
matlist boundmat

************************************************************
// Option 2: Add program -power_cmd_guessncheck- to automate the 
// guess-and-check process. This is more work to implement, but allows 
// you to use -gsdesign- directly. Here we will use a bisection algorithm 
// to determine the necessary sample size.
************************************************************

// -power_cmd_guessncheck-, a simulation-based sample-size evaluator
capture program drop power_cmd_guessncheck
program power_cmd_guessncheck, rclass
    version 18
    syntax, power(real)           /// Requested power
          [ alpha(real 0.05)      /// Alpha level
            m0(real 0)            /// Mean under the null
            ma(real 1)            /// Mean under the alternative
            sd(real 1)            /// Standard deviation
            reps(integer 100)     /// Number of repetitions for power calculation
            startval(integer 100) /// Starting value for sample size
            maxval(integer 1000)  /// Maximum value for sample size
            minval(integer 1)     /// Minimum value for sample size
            powtol(real 0.01)     /// Tolerance for power search (default 1%)
            biter(integer 100)    /// Max iterations for bisection algorithm
            seed(string)          /// Random seed (no default seed)
            nfractional ]         //  Unused, but required by -gsdesign-

    // Set the random seed (if requested)
    capture confirm number `seed'
    if (!_rc) {
        set seed `seed'
    }

    // Test -powtol()- to ensure it's not set too low
    assert (`powtol' > 0)
    // Test -maxval()- to ensure it's not set too low 
    power_cmd_simttest, n(`maxval') m0(`m0') ma(`ma') sd(`sd') reps(`reps')
    if (r(power) < `power') {
        di as err "error: {bf:maxval()} too small"
        di as error "{p 4 4 2}With `maxval' subjects, the attained power is " ///
                    r(power) ", which is less than the requested power of "   ///
                    "`power'.{p_end}"
        exit 480
    }
    // Test -minval()- to ensure it's not set too high
    power_cmd_simttest, n(`minval') m0(`m0') ma(`ma') sd(`sd') reps(`reps')
    if (r(power) > `power') {
        di as err "error: {bf:minval()} too large"
        di as error "{p 4 4 2}With `minval' subjects, the attained power is " ///
                    r(power) ", which is more than the requested power of "   ///
                    "`power'.{p_end}"
        exit 480
    }

    // Compare requested power to power attained via simulation 
    // at the starting value
    local thisN = `startval'
    local lastN = .
    power_cmd_simttest, n(`thisN') m0(`m0') ma(`ma') sd(`sd') reps(`reps')
    local pdiff = r(power) - `power'

    // Repeat guess-n-check process using bisection algorithm
    local i = 1
    while (abs(`pdiff') > `powtol') {
        local ++i
        if (`i' > `biter') {
            di as err "convergence not achieved"
            di as err "{p 4 4 2}To continue iterating, re-run with a "   ///
                      "larger value of {bf:biter()}.{p_end}"
            exit 430
        }

        // Bisection algorithm
        local lastN = `thisN'
        if (`pdiff' > 0) { // Overpowered: decrease N
            local maxval = `thisN'
        }
        else {             // Underpowered: increase N
            local minval = `thisN'
        }
        local thisN = ceil((`maxval' + `minval') / 2)
        if (`thisN' == `lastN') { // Already converged
            continue, break
        }
        power_cmd_simttest, n(`thisN') m0(`m0') ma(`ma') sd(`sd') reps(`reps')
        local pdiff = r(power) - `power'
    }

    // Return results
    return scalar power = r(power)
    return scalar N = `thisN'
    return scalar alpha = `alpha'
    return scalar m0 = `m0'
    return scalar ma = `ma'
    return scalar sd = `sd'
end

// Try it as a -power- command
power guessncheck, power(0.9) m0(70) ma(75) sd(15) startval(150) reps(500) seed(9)

// Now with -gsdesign-
gsdesign guessncheck, power(0.9) m0(70) ma(75) sd(15) startval(150) reps(500) ///
         seed(9) nlooks(4) efficacy(errobfleming) graphbounds

********************************************************************************

​​​​​​​