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

["Joseph Coveney" <stajc2@gmail.com>]

Garry Anderson wrote:

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.

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

You'd use the inverse function, -invdunnettprob()-.  See below.

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

. scalar define Dunnett = invdunnettprob(k, df, 0.95)

. forvalues i = 2/3 {
  2.         display in smcl as text _newline(1) "Group `i' versus Group 1"
  3.     quietly lincom `i'.group
  4.         scalar define HalfCI = r(se) * Dunnett
  5.     display in smcl as result r(estimate) - HalfCI
  6.         display in smcl as result r(estimate) + HalfCI
  7. }

Group 2 versus Group 1
-.76352298
1.2281004

Group 3 versus Group 1
-.9074282
1.0841952

. 
. 
. * negative binomial model
. webuse rod93, clear

. nbreg deaths i.cohort, exposure(exp) irr 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 |        IRR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      cohort |
          2  |   .7651995   .5537904    -0.37   0.712     .1852434    3.160869
          3  |   .6329298   .4580292    -0.63   0.527     .1532395    2.614209
             |
       _cons |   .1240922   .0635173    -4.08   0.000     .0455042    .3384052
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

. scalar define Dunnett = invdunnettprob(e(df_m), 1e6, 0.95)

. forvalues i = 2/3 {
  2.         display in smcl as text _newline(1) "Group `i' versus Group 1"
  3.     quietly lincom `i'.cohort
  4.         scalar define HalfCI = r(se) * Dunnett
  5.         display in smcl as text "Estimate = " exp(r(estimate))
  6.     display in smcl as text "LCL = " as result exp(r(estimate) - HalfCI)
  7.         display in smcl as text "UCL = " as result exp(r(estimate) +
HalfCI)
  8. }

Group 2 versus Group 1
Estimate = .76519946
LCL = .1543422
UCL = 3.7937143

Group 3 versus Group 1
Estimate = .63292982
LCL = .12767877
UCL = 3.1375629

. 
. exit

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/