Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • power for Dunnett many to one comparisons?

    I was trying to figure out the power for a Dunnett many to one comparisons design, not just the omnibus F-test There is not a noncentral distribution understanding for the Dunnett t distribution that I know of, which is a multivariate normal type. So I did something that matched my simulation below, using Stata's power oneway command, using the contrast option, with an alpha determined from a critical Dunnett t statistic using the ttail function. The example gets the adjusted two tailed probability for an ANOVA design with k=7 treatment groups and 1 control group with n=4 subjects per group. I'm not aware of others approaching power for Dunnett's comparisons in this way, but it seems to have worked. I was going to translate this to a full power command, but wanted to know others opinions on how I've approached this.

    Code:
    local k = 7
    local n = 4
    local v = (`k'+1)*(`n'-1)
    local dunnett_t = invdunnettprob(`k', `v', 0.95)
    scalar adj_p = 2*(1-t(`v', `dunnett_t'))
    display adj_p
    Here is the full simulation and matching power oneway calculations:

    Code:
    local k = 7
    local n = 4
    local v = (`k'+1)*(`n'-1)
    local dunnett_t = invdunnettprob(`k', `v', 0.95)
    scalar adj_p = 2*(1-t(`v', `dunnett_t'))
    display adj_p
    
    local kd "40 35 70 20 30 50 55"
    local t1 = 5 // control group dct
    local i = 2
    local dct = "5"
    foreach l of local kd {
        local t`i' = `t1' + ln((100-`l')/100)/ln(2)*-1
        display `t`i''
        local dct = "`dct'" + " `t`i''"
        local ++i
    }
    numlist "`dct'"
    local mymeans `r(numlist)'
    display "`mymeans'"
    power oneway `mymeans', ///
        varerror(0.25) npergroup(4) // power .97
    power oneway `mymeans', ///
        varerror(0.25) contrast(-1 1 0 0 0 0 0 0) ///
        npergroup(4) alpha(`=adj_p') // power .26
    power oneway `mymeans', ///
        varerror(0.25) contrast(-1 0 1 0 0 0 0 0) ///
        npergroup(4) alpha(`=adj_p') // power .17
    power oneway `mymeans', ///
        varerror(0.25) contrast(-1 0 0 1 0 0 0 0) ///
        npergroup(4) alpha(`=adj_p') // power .98
    power oneway `mymeans', ///
        varerror(0.25) contrast(-1 0 0 0 1 0 0 0) ///
        npergroup(4) alpha(`=adj_p') // power .04
    power oneway `mymeans', ///
        varerror(0.25) contrast(-1 0 0 0 0 1 0 0) ///
        npergroup(4) alpha(`=adj_p') // power .11
    power oneway `mymeans', ///
        varerror(0.25) contrast(-1 0 0 0 0 0 1 0) ///
        npergroup(4) alpha(`=adj_p') // power .52
    power oneway `mymeans', ///
        varerror(0.25) contrast(-1 0 0 0 0 0 0 1) ///
        npergroup(4) alpha(`=adj_p') // power .67
    
    // compare power of F-test and all Dunnett p-value
    frame reset
    frame create results double(ftest dunnett2vs1 dunnett3vs1 dunnett4vs1 ///
        dunnett5vs1 dunnett6vs1 dunnett7vs1 dunnett8vs1)
    quietly {
        forvalues i = 1(1)5000 {
            frame create simdata
            frame change simdata
            set obs 4
            local kd "40 35 70 20 30 50 55"
            local t1 = 5 // control group dct
            local j = 2
            foreach l of local kd {
                local t`j' = `t1' + ln((100-`l')/100)/ln(2)*-1
                local ++j
            }
            local ve = 0.5
            local ntx :word count `kd'
            local groups = `ntx' + 1
            generate y1 = rnormal(`t1', `ve')
            forvalues j = 2/`groups' {
                generate y`j' = rnormal(`t`j'', `ve')
            }
            generate row = _n
            reshape long y, i(row) j(group)
            regress y i.group
            scalar ftest = Ftail(e(df_m), e(df_r), e(F))
            pwcompare group, effects mcompare(dunnett)
            forvalues j = 1/`ntx' {
                local g = `j' + 1
                scalar dunnett`g'vs1 = el(r(table_vs),4,`j') // *1/2 if one-sided
            }
            frame post results ///
                (scalar(ftest)) ///
                (scalar(dunnett2vs1)) ///
                (scalar(dunnett3vs1)) ///
                (scalar(dunnett4vs1)) ///
                (scalar(dunnett5vs1)) ///
                (scalar(dunnett6vs1)) ///
                (scalar(dunnett7vs1)) ///
                (scalar(dunnett8vs1))
            frame change default
            frame drop simdata
            if mod(`i',100) == 0 noisily display `i'
        }
    }
    frame change results
    save results_power_2.dta, replace
    generate sig_ftest = ftest < 0.05 // power = .97, above = .97
    generate sig_2vs1 = dunnett2vs1 < 0.05 // power = .26, above = .26
    generate sig_3vs1 = dunnett3vs1 < 0.05 // power = .18, above = .17
    generate sig_4vs1 = dunnett4vs1 < 0.05 // power = .97, above = .98
    generate sig_5vs1 = dunnett5vs1 < 0.05 // power = .04, above = .04
    generate sig_6vs1 = dunnett6vs1 < 0.05 // power = .12, above = .11
    generate sig_7vs1 = dunnett7vs1 < 0.05 // power = .52, above = .52
    generate sig_8vs1 = dunnett8vs1 < 0.05 // power = .67, above = .67
    summarize sig*
    frame change default
    frame drop results

  • #2
    Hi Dave Airey. I don't have time right now (and possibly not the expertise either!) to review what you have done. But your post prompted me to check on whether PASS has a procedure for Dunnett's test--and it does. I thought, therefore, you might find something useful in its documentation, which you can view here: Cheers,
    Bruce
    --
    Bruce Weaver
    Email: [email protected]
    Web: http://sites.google.com/a/lakeheadu.ca/bweaver/
    Version: Stata/MP 18.0 (Windows)

    Comment


    • #3
      Thank you Bruce Weaver. I read the NCSS PASS help page, and they do the sample size or power estimation by simulation. I generally do power and sample size by simulation unless there is a noncentral distribution solution. Here I think I found a short cut for power estimation for a many to one design. I am suggesting doing power estimation by changing the problem to use a noncentral distribution for a given contrast. That's represented in the first snippet of code. That's the crux of my question, whether you buy this approach or not.

      The rest of the code was really just to prove to myself that I could get the same power for a given contrast using this approach vs using a simulation and collecting the p-values from Stata's pwcompare...mcompare(dunnett) command. Stata's procedure of getting the Dunnett twosided p-value appears to be to rely on the normal t statistic and use it's internal function dunnettprob() which is based on a multivariate normal estimation of the Dunnett t distribution. So I figured if I could go from the Student's t to a Dunnett's probability, then I could go from Dunnett's t to a Student's probability and then the power calculation (for a given Dunnett's comparison) could used a noncentral distribution (either t of F which is known).

      All that said, I did learn from your link how they speeded up their simulation code by up to 70% and that I'm calculating "any pair" power rather than "all pairs" power. I will also check what I did vs their in help example.

      Below is some output showing how the Dunnett p-value and critical values can be calculated by hand to familiarize you with these invdunnettprob and dunnettprob functions.

      Code:
      . regress y i.group
      
            Source |       SS           df       MS      Number of obs   =        24
      -------------+----------------------------------   F(3, 20)        =      1.24
             Model |  1.65438846         3  .551462821   Prob > F        =    0.3229
          Residual |  8.92361586        20  .446180793   R-squared       =    0.1564
      -------------+----------------------------------   Adj R-squared   =    0.0299
             Total |  10.5780043        23  .459913231   Root MSE        =    .66797
      
      ------------------------------------------------------------------------------
                 y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
      -------------+----------------------------------------------------------------
             group |
                2  |   .5787694   .3856513     1.50   0.149    -.2256851    1.383224
                3  |   .4849888   .3856513     1.26   0.223    -.3194657    1.289443
                4  |   .6867247   .3856513     1.78   0.090    -.1177298    1.491179
                   |
             _cons |   5.148036   .2726967    18.88   0.000       4.5792    5.716871
      ------------------------------------------------------------------------------
      
      . matrix list r(table)
      
      r(table)[9,5]
                      1b.          2.          3.          4.            
                   group       group       group       group       _cons
           b           0   .57876945   .48498885   .68672474   5.1480355
          se           .   .38565131   .38565131   .38565131   .27269665
           t           .   1.5007584   1.2575838   1.7806882    18.87825
      pvalue           .   .14904077   .22302887   .09015986   3.241e-14
          ll           .  -.22568509  -.31946568  -.11772979   4.5792003
          ul           .    1.383224   1.2894434   1.4911793   5.7168708
          df          20          20          20          20          20
        crit   2.0859634   2.0859634   2.0859634   2.0859634   2.0859634
       eform           0           0           0           0           0
      
      . display invt(20, 0.975) // matches crit
      2.0859634
      
      . display 2*(1-t(20,el(r(table),3,2))) // matches pvalue
      .14904077
      
      . pwcompare group, effects mcompare(dunnett)
      
      Pairwise comparisons of marginal linear predictions
      
      Margins: asbalanced
      
      ---------------------------
                   |    Number of
                   |  comparisons
      -------------+-------------
             group |            3
      ---------------------------
      
      ------------------------------------------------------------------------------
                   |                             Dunnett              Dunnett
                   |   Contrast   Std. err.      t    P>|t|     [95% conf. interval]
      -------------+----------------------------------------------------------------
             group |
           2 vs 1  |   .5787694   .3856513     1.50   0.330    -.4009193    1.558458
           3 vs 1  |   .4849888   .3856513     1.26   0.465    -.4946999    1.464678
           4 vs 1  |   .6867247   .3856513     1.78   0.211     -.292964    1.666413
      ------------------------------------------------------------------------------
      
      . matrix list r(table_vs)
      
      r(table_vs)[9,3]
                    2vs1.       3vs1.       4vs1.
                   group       group       group
           b   .57876945   .48498885   .68672474
          se   .38565131   .38565131   .38565131
           t   1.5007584   1.2575838   1.7806882
      pvalue   .32999048   .46542638   .21050642
          ll  -.40091928  -.49469988  -.29296398
          ul   1.5584582   1.4646776   1.6664135
          df          20          20          20
        crit   2.5403485   2.5403485   2.5403485
       eform           0           0           0
      
      . display invdunnettprob(3, 20, 0.95) // matches crit
      2.5403485
      
      . display invdunnettprob(3, 20, 0.90) // onesided crit
      2.1917274
      
      . display 1-dunnettprob(3, 20, abs(el(r(table_vs),3,1))) // matches pvalue
      .32999048

      Comment


      • #4
        I got positive feedback from another statistician and will role this up into a power "dunnett_anova" command. Thanks for the feedback and link to NCSS PASS Bruce Weaver.

        Comment


        • #5
          Just a follow up to note that Hsu (1996) and Pan and Kupper (1999) each developed analytical solutions to this problem, and both are implemented in NCSS PASS, not just simulation. NCSS PASS provides both simulation and analytic solutions for the common many-to-one control ANOVA design. I was under the mistaken impression that only simulation solutions were possible. Tsk.

          Hsu JC. 1996. Multiple Comparisons: Theory and Methods. London: Chapman & Hall
          Pan Z, Kupper LL. 1999. Sample size determination for multiple comparison studies treating confidence interval width as random. Stat. Med. 18:1475–88

          Comment

          Working...
          X