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/