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

[Garry Anderson <g.anderson@unimelb.edu.au>]

Many thanks Joseph for your most helpful response to my posting.

I was considering how to obtain the Dunnett 95% confidence intervals from the negative binomial example.
In the example is  -lincom `i'.cohort-
Since this is the unadjusted P-value and confidence interval, is it possible to obtain only Dunnett 95% confidence intervals for the incidence rate ratio by using 
-lincom `i'.cohort, eform level(xx.xx)

However, do you know please what I should use for the value of xx.xx ?

I have 4 comparisons and 737 residual degrees of freedom.

Kind regards, Garry


________________________________________
From: Joseph Coveney [stajc2@gmail.com]
Sent: Friday, 17 January 2014 2:34 AM
To: statalist@hsphsun2.harvard.edu
Cc: Garry Anderson
Subject: Re: Dunnett's P-value and negative binomial regression

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/