Bookmark and Share

Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

st: Re: Dunnett's P-value and negative binomial regression


From   "Joseph Coveney" <[email protected]>
To   <[email protected]>
Subject   st: Re: Dunnett's P-value and negative binomial regression
Date   Fri, 17 Jan 2014 00:34:51 +0900

Garry Anderson wrote:

The command -pwcompare trt,mcomp(dunnett) effect- is not allowed after -nbreg-
or -stcox- because of the dunnett method of comparison.  dunnett is only 
allowed with results from anova, manova, regress and mvreg.

The R package multcomp allows the dunnett option for generalized linear models 
and Cox models, however I would prefer to use Stata.

Is it possible to use Stata to estimate dunnett P-values and confidence 
intervals after a negative binomial regression?

The dunnett method compares multiple treatments with a control treatment.

--------------------------------------------------------------------------------

You could use the -dunnettprob()- function to knit your own Dunnett's tests, as
shown below.  It should work with linear predictors of generalized linear 
models.  You'll need to decide what to do about denominator degrees of freedom 
and any imbalances between groups.

Joseph Coveney

. clear *

. set more off

. set seed `=date("2014-01-16", "YMD")'

. set obs 30
obs was 0, now 30

. 
. * linear model
. generate byte group = mod(_n, 3) + 1

. generate double response = rnormal()

. quietly regress response i.group

. scalar define k = e(df_m)

. scalar define df = e(df_r)

. quietly pwcompare i.group, mcompare(dunnett)

. matrix list r(table_vs)

r(table_vs)[9,2]
              2vs1.       3vs1.
             group       group
     b    .2322887   .08838348
    se   .42676226   .42676226
     t   .54430469   .20710239
pvalue   .80959746   .96911369
    ll  -.76352298   -.9074282
    ul   1.2281004   1.0841952
    df          27          27
  crit   2.3334108   2.3334108
 eform           0           0

. forvalues i = 2/3 {
  2.     quietly lincom `i'.group
  3.     display in smcl as result 1 - ///
>       dunnettprob(k, df, abs(r(estimate) / r(se)))
  4. }
.80959746
.96911369

. 
. * negative binomial model
. webuse rod93, clear

. quietly regress deaths i.cohort

. scalar define k = e(df_m)

. scalar define df = e(df_r)

. nbreg deaths i.cohort, exposure(exp) nolog

Negative binomial regression                      Number of obs   =         21
                                                  LR chi2(2)      =       0.40
Dispersion     = mean                             Prob > chi2     =     0.8171
Log likelihood = -131.3799                        Pseudo R2       =     0.0015

------------------------------------------------------------------------------
      deaths |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      cohort |
          2  |  -.2676187   .7237203    -0.37   0.712    -1.686085    1.150847
          3  |  -.4573957   .7236651    -0.63   0.527    -1.875753    .9609618
             |
       _cons |  -2.086731   .5118559    -4.08   0.000     -3.08995   -1.083511
ln(exposure) |          1  (exposure)
-------------+----------------------------------------------------------------
    /lnalpha |   .5939963   .2583615                       .087617    1.100376
-------------+----------------------------------------------------------------
       alpha |   1.811212   .4679475                       1.09157    3.005294
------------------------------------------------------------------------------
Likelihood-ratio test of alpha=0:  chibar2(01) = 4056.27 Prob>=chibar2 = 0.000

. forvalues i = 2/3 {
  2.     quietly lincom `i'.cohort
  3.     display in smcl as text _newline(1) ///
>       "with ersatz df: " _continue
  4.     display in smcl as result 1 - ///
>       dunnettprob(k, df, abs(r(estimate) / r(se)))
  5.     display in smcl as text "asymptotic: " _continue
  6.     display in smcl as result 1 - ///
>       dunnettprob(k, 1e6, abs(r(estimate) / r(se)))
  7. }

with ersatz df: .90589551
asymptotic: .90531712

with ersatz df: .75545379
asymptotic: .75174037

. 
. 
end of do-file

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/


© Copyright 1996–2018 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   Site index