--------------------------------------------------------------------------------------------------------------- log: c:\Imbook\transparencies\acdcount\trcount2009.txt log type: text opened on: 27 Mar 2009, 15:01:14 . . ********** OVERVIEW OF trcount2009.do ********** . . * STATA Program to demonstrate Cound Data Regression Models . * Based on mus17p1cnt.do . * A. Colin Cameron and Pravin K. Trivedi (2009) . * "Microeconometrics using Stata", Stata Press . . * 1: COUNT REGRESSION: BASIC CROSS-SECTION . * ANALYSIS WITHOUT REGRESSORS . * BASIC POISSON REGRESSION . * GLM REGRESSION . * NONLINEAR LEAST SQUARES . * DIAGNOSTICS . * NEGATIVE BINOMIAL REGRESSION . * 2: COUNT REGRESSION: ADDITIONAL CROSS-SECTION . * CENSORED AND TRUNCATED . * HURDLE OR TWO-PART . * ZERO-INFLATED . * MODEL COMPARISON . * FINITE MIXTURE . * ENDOGENOUS REGRESSORS . * 3. COUNT REGRESSION: BEYOND CROSS-SECTION . * TIME SERIES - NO EXAMPLE . * PANEL POISSON: POOLED AND PA . * PANEL POISSON: RANDOM EFFECTS . * PANEL POISSON: FIXED EFFECTS . * PANEL NEGATIVE BINOMIAL: POOLED AND PA . * 4: MULTIVARIATE, MAXIMUM SIMULATED LIKELIHOOD, BAYESIAN . * NO DATA EXAMPLES . . * To run you need file mus17data.dta . * in your directory . * Stata user-written commands . * spost9_ado // For prvalue, prcount, listcoef in MUS section 17.3, 17.4 . * hplogit // For hurdle model in MUS section 17.3 . * hnblogit // For hurdle model in MUS section 17.3 . * fmm // For finite mixture models in MUS section 17.3 . * ivpois // For endogenous regressors . * are used . . ********** SETUP ********** . . set more off . version 10.1 . clear all . set mem 10m Current memory allocation current memory usage settable value description (1M = 1024k) -------------------------------------------------------------------- set maxvar 5000 max. variables allowed 1.909M set memory 10M max. data space 10.000M set matsize 400 max. RHS vars in models 1.254M ----------- 13.163M . set scheme s1mono /* Graphics scheme */ . . *********** 1: COUNT REGRESSION: BASICS (POISSON, GLM, NLS, DIAGNOSTICS, NB) . . *** 1. BASICS: ANALYSIS WITHOUT REGRESSORS . . * Summarize doctor visits data (MEPS 2003 data 65-90 years old: see MUS p.557) . use mus17data.dta . histogram docvis if docvis < 40, discrete /// > xtitle("# Doctor visits (for < 40 visits)") (start=0, width=1) . graph export countfig01histraw.wmf, replace (file c:\Imbook\transparencies\acdcount\countfig01histraw.wmf written in Windows Metafile format) . . * Tabulate docvis after recoding values > 10 to ranges 11-40 and 41-60 . generate dvrange = docvis . recode dvrange (11/40 = 40) (41/60 = 60) (dvrange: 784 changes made) . label define dvcounts 40 "11-40" 60 "41-60" . label values dvrange dvcounts . tabulate dvrange dvrange | Freq. Percent Cum. ------------+----------------------------------- 0 | 401 10.91 10.91 1 | 314 8.54 19.45 2 | 358 9.74 29.18 3 | 334 9.08 38.26 4 | 339 9.22 47.48 5 | 266 7.23 54.72 6 | 231 6.28 61.00 7 | 202 5.49 66.49 8 | 179 4.87 71.36 9 | 154 4.19 75.55 10 | 108 2.94 78.49 11-40 | 774 21.05 99.54 41-60 | 14 0.38 99.92 73 | 1 0.03 99.95 106 | 1 0.03 99.97 144 | 1 0.03 100.00 ------------+----------------------------------- Total | 3,677 100.00 . . * Generate Poisson data with sample mean of docvis . quietly summarize docvis . scalar dvmean = r(mean) . set seed 10101 // set the seed ! . generate ypoisson= rpoisson(dvmean) // Draw from Poisson(mu=dvmean) . summarize docvis ypoisson Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- docvis | 3677 6.822682 7.394937 0 144 ypoisson | 3677 6.836008 2.631633 0 18 . tabulate ypoisson ypoisson | Freq. Percent Cum. ------------+----------------------------------- 0 | 4 0.11 0.11 1 | 33 0.90 1.01 2 | 89 2.42 3.43 3 | 214 5.82 9.25 4 | 375 10.20 19.45 5 | 475 12.92 32.36 6 | 546 14.85 47.21 7 | 533 14.50 61.71 8 | 471 12.81 74.52 9 | 374 10.17 84.69 10 | 252 6.85 91.54 11 | 137 3.73 95.27 12 | 92 2.50 97.77 13 | 45 1.22 98.99 14 | 19 0.52 99.51 15 | 10 0.27 99.78 16 | 3 0.08 99.86 17 | 3 0.08 99.95 18 | 2 0.05 100.00 ------------+----------------------------------- Total | 3,677 100.00 . histogram docvis if docvis < 40, discrete addplot(hist ypoisson, fcolor(white) discrete) /// > xtitle("# Doctor visits versus Poisson") legend(off) (start=0, width=1) . graph export countfig02histpoisson.wmf, replace (file c:\Imbook\transparencies\acdcount\countfig02histpoisson.wmf written in Windows Metafile format) . . * Generate negative binomial (mu=1 var=2) generated data . quietly nbreg docvis . scalar alpha = e(alpha) . display alpha .84081223 . set seed 10101 // set the seed ! . generate mugamma = rgamma(1/alpha,alpha) . generate ynegbin = rpoisson(dvmean*mugamma) // NB generated as a Poisson-gamma mixture . summarize docvis ynegbin mugamma Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- docvis | 3677 6.822682 7.394937 0 144 ynegbin | 3677 6.996736 6.999999 0 66 mugamma | 3677 1.027634 .9443734 .0001458 8.468604 . tabulate ynegbin ynegbin | Freq. Percent Cum. ------------+----------------------------------- 0 | 375 10.20 10.20 1 | 411 11.18 21.38 2 | 349 9.49 30.87 3 | 322 8.76 39.62 4 | 248 6.74 46.37 5 | 232 6.31 52.68 6 | 241 6.55 59.23 7 | 185 5.03 64.26 8 | 179 4.87 69.13 9 | 171 4.65 73.78 10 | 98 2.67 76.45 11 | 104 2.83 79.28 12 | 120 3.26 82.54 13 | 90 2.45 84.99 14 | 93 2.53 87.52 15 | 71 1.93 89.45 16 | 52 1.41 90.86 17 | 48 1.31 92.17 18 | 37 1.01 93.17 19 | 35 0.95 94.13 20 | 36 0.98 95.10 21 | 21 0.57 95.68 22 | 22 0.60 96.27 23 | 17 0.46 96.74 24 | 16 0.44 97.17 25 | 13 0.35 97.53 26 | 11 0.30 97.82 27 | 9 0.24 98.07 28 | 10 0.27 98.34 29 | 11 0.30 98.64 30 | 5 0.14 98.78 31 | 4 0.11 98.88 32 | 6 0.16 99.05 33 | 8 0.22 99.27 34 | 4 0.11 99.37 35 | 2 0.05 99.43 36 | 1 0.03 99.46 37 | 4 0.11 99.56 38 | 5 0.14 99.70 41 | 3 0.08 99.78 42 | 1 0.03 99.81 43 | 1 0.03 99.84 45 | 1 0.03 99.86 48 | 2 0.05 99.92 51 | 1 0.03 99.95 59 | 1 0.03 99.97 66 | 1 0.03 100.00 ------------+----------------------------------- Total | 3,677 100.00 . histogram docvis if docvis < 40, discrete addplot(hist ynegbin if ynegbin < 40, fcolor(white) discrete) /// > xtitle("# Doctor visits versus negative binomial") legend(off) (start=0, width=1) . graph export countfig03histnegbin.wmf, replace (file c:\Imbook\transparencies\acdcount\countfig03histnegbin.wmf written in Windows Metafile format) . . *** 1. BASICS: POISSON REGRESSION . . drop _all . . * Summary statistics for doctor visits data . use mus17data.dta . global xlist private medicaid age age2 educyr actlim totchr . describe docvis $xlist storage display value variable name type format label variable label --------------------------------------------------------------------------------------------------------------- docvis float %9.0g # doctor visits private byte %8.0g =1 if has private supplementary insurance medicaid byte %8.0g =1 if has Medicaid public insurance age byte %8.0g Age age2 float %9.0g Age-squared educyr byte %8.0g Years of education actlim byte %8.0g =1 if activity limitation totchr byte %8.0g # chronic conditions . summarize docvis $xlist, sep(10) Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- docvis | 3677 6.822682 7.394937 0 144 private | 3677 .4966005 .5000564 0 1 medicaid | 3677 .166712 .3727692 0 1 age | 3677 74.24476 6.376638 65 90 age2 | 3677 5552.936 958.9996 4225 8100 educyr | 3677 11.18031 3.827676 0 17 actlim | 3677 .333152 .4714045 0 1 totchr | 3677 1.843351 1.350026 0 8 . . * Poisson with preffered robust standard errors . poisson docvis $xlist, vce(robust) nolog // Poisson robust SEs Poisson regression Number of obs = 3677 Wald chi2(7) = 720.43 Prob > chi2 = 0.0000 Log pseudolikelihood = -15019.64 Pseudo R2 = 0.1297 ------------------------------------------------------------------------------ | Robust docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .1422324 .036356 3.91 0.000 .070976 .2134889 medicaid | .0970005 .0568264 1.71 0.088 -.0143773 .2083783 age | .2936722 .0629776 4.66 0.000 .1702383 .4171061 age2 | -.0019311 .0004166 -4.64 0.000 -.0027475 -.0011147 educyr | .0295562 .0048454 6.10 0.000 .0200594 .039053 actlim | .1864213 .0396569 4.70 0.000 .1086953 .2641474 totchr | .2483898 .0125786 19.75 0.000 .2237361 .2730435 _cons | -10.18221 2.369212 -4.30 0.000 -14.82578 -5.538638 ------------------------------------------------------------------------------ . . * Poisson with default ML standard errors . poisson docvis $xlist // Poisson default ML standard errors Iteration 0: log likelihood = -15019.656 Iteration 1: log likelihood = -15019.64 Iteration 2: log likelihood = -15019.64 Poisson regression Number of obs = 3677 LR chi2(7) = 4477.98 Prob > chi2 = 0.0000 Log likelihood = -15019.64 Pseudo R2 = 0.1297 ------------------------------------------------------------------------------ docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .1422324 .0143311 9.92 0.000 .114144 .1703208 medicaid | .0970005 .0189307 5.12 0.000 .0598969 .134104 age | .2936722 .0259563 11.31 0.000 .2427988 .3445457 age2 | -.0019311 .0001724 -11.20 0.000 -.0022691 -.0015931 educyr | .0295562 .001882 15.70 0.000 .0258676 .0332449 actlim | .1864213 .014566 12.80 0.000 .1578726 .2149701 totchr | .2483898 .0046447 53.48 0.000 .2392864 .2574933 _cons | -10.18221 .9720115 -10.48 0.000 -12.08732 -8.277101 ------------------------------------------------------------------------------ . . * Poisson coefficients as incidence ratios . poisson docvis $xlist, irr vce(robust) nolog // Incidence ratios Poisson regression Number of obs = 3677 Wald chi2(7) = 720.43 Prob > chi2 = 0.0000 Log pseudolikelihood = -15019.64 Pseudo R2 = 0.1297 ------------------------------------------------------------------------------ | Robust docvis | IRR Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | 1.152845 .0419128 3.91 0.000 1.073555 1.23799 medicaid | 1.101861 .0626148 1.71 0.088 .9857255 1.231679 age | 1.341344 .0844747 4.66 0.000 1.185587 1.517564 age2 | .9980708 .0004158 -4.64 0.000 .9972562 .998886 educyr | 1.029997 .0049907 6.10 0.000 1.020262 1.039826 actlim | 1.20493 .0477838 4.70 0.000 1.114823 1.30232 totchr | 1.28196 .0161253 19.75 0.000 1.250741 1.313957 ------------------------------------------------------------------------------ . . * MEM and AME for Poisson . quietly poisson docvis $xlist, vce(robust) . mfx // MEM: Marginal effect for Poisson evaluated at average of x Marginal effects after poisson y = predicted number of events (predict) = 6.2674204 ------------------------------------------------------------------------------ variable | dy/dx Std. Err. z P>|z| [ 95% C.I. ] X ---------+-------------------------------------------------------------------- private*| .8926136 .22788 3.92 0.000 .445984 1.33924 .4966 medicaid*| .6281643 .38182 1.65 0.100 -.120185 1.37651 .166712 age | 1.840567 .39247 4.69 0.000 1.07134 2.60979 74.2448 age2 | -.012103 .0026 -4.66 0.000 -.017194 -.007012 5552.94 educyr | .1852412 .03067 6.04 0.000 .125127 .245355 11.1803 actlim*| 1.207039 .26864 4.49 0.000 .680506 1.73357 .333152 totchr | 1.556764 .07602 20.48 0.000 1.40777 1.70575 1.84335 ------------------------------------------------------------------------------ (*) dy/dx is for discrete change of dummy variable from 0 to 1 . margeff // AME: Average marginal effect for Poisson Average partial effects after poisson y = E(docvis) (expected number of counts) ------------------------------------------------------------------------------ variable | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .9701906 .2645543 3.67 0.000 .4516738 1.488707 medicaid | .6830664 .4219606 1.62 0.105 -.1439612 1.510094 age | 2.032557 .4489563 4.53 0.000 1.152618 2.912495 age2 | -.0131753 .0028473 -4.63 0.000 -.0187559 -.0075947 educyr | .201682 .033795 5.97 0.000 .135445 .2679189 actlim | 1.295942 .304792 4.25 0.000 .6985604 1.893323 totchr | 1.712165 .0935177 18.31 0.000 1.528874 1.895456 ------------------------------------------------------------------------------ . . *** 1. BASICS: GLM REGRESSION . . * GLM for Poisson with robust standard errors . glm docvis $xlist, family(poisson) link(log) vce(robust) nolog Generalized linear models No. of obs = 3677 Optimization : ML Residual df = 3669 Scale parameter = 1 Deviance = 18395.14033 (1/df) Deviance = 5.013666 Pearson = 23147.37781 (1/df) Pearson = 6.308906 Variance function: V(u) = u [Poisson] Link function : g(u) = ln(u) [Log] AIC = 8.173859 Log pseudolikelihood = -15019.6398 BIC = -11726.81 ------------------------------------------------------------------------------ | Robust docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .1422324 .036356 3.91 0.000 .070976 .2134889 medicaid | .0970005 .0568264 1.71 0.088 -.0143773 .2083783 age | .2936722 .0629776 4.66 0.000 .1702383 .4171061 age2 | -.0019311 .0004166 -4.64 0.000 -.0027475 -.0011147 educyr | .0295562 .0048454 6.10 0.000 .0200594 .039053 actlim | .1864213 .0396569 4.70 0.000 .1086953 .2641474 totchr | .2483898 .0125786 19.75 0.000 .2237361 .2730435 _cons | -10.18221 2.369212 -4.30 0.000 -14.82578 -5.538638 ------------------------------------------------------------------------------ . . * GLM for Poisson with GLM corrected standard errors . * These are default ML standard errors multiplied by square root Pearson statistic . glm docvis $xlist, family(poisson) link(log) scale(x2) nolog Generalized linear models No. of obs = 3677 Optimization : ML Residual df = 3669 Scale parameter = 1 Deviance = 18395.14033 (1/df) Deviance = 5.013666 Pearson = 23147.37781 (1/df) Pearson = 6.308906 Variance function: V(u) = u [Poisson] Link function : g(u) = ln(u) [Log] AIC = 8.173859 Log likelihood = -15019.6398 BIC = -11726.81 ------------------------------------------------------------------------------ | OIM docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .1422324 .0359961 3.95 0.000 .0716813 .2127836 medicaid | .0970005 .0475493 2.04 0.041 .0038055 .1901954 age | .2936722 .0651959 4.50 0.000 .1658906 .4214538 age2 | -.0019311 .0004331 -4.46 0.000 -.00278 -.0010822 educyr | .0295562 .0047271 6.25 0.000 .0202912 .0388212 actlim | .1864213 .0365861 5.10 0.000 .1147139 .2581288 totchr | .2483898 .0116663 21.29 0.000 .2255243 .2712554 _cons | -10.18221 2.441453 -4.17 0.000 -14.96737 -5.397048 ------------------------------------------------------------------------------ (Standard errors scaled using square root of Pearson X2-based dispersion) . . *** 1. BASICS: NONLINEAR LEAST SQUARES . . * Nonlinear least squares with exponential conditional mean . capture drop one . generate one = 1 . nl (docvis = exp({xb: $xlist one})), vce(robust) nolog (obs = 3677) Nonlinear regression Number of obs = 3677 R-squared = 0.5436 Adj R-squared = 0.5426 Root MSE = 6.804007 Res. dev. = 24528.25 ------------------------------------------------------------------------------ | Robust docvis | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- /xb_private | .1235144 .0395179 3.13 0.002 .0460351 .2009937 /xb_medicaid | .0856747 .0649936 1.32 0.188 -.0417525 .2131018 /xb_age | .2951153 .0720509 4.10 0.000 .1538516 .4363789 /xb_age2 | -.0019481 .0004771 -4.08 0.000 -.0028836 -.0010127 /xb_educyr | .0309924 .0051192 6.05 0.000 .0209557 .0410291 /xb_actlim | .1916735 .0413705 4.63 0.000 .110562 .2727851 /xb_totchr | .2191967 .0151021 14.51 0.000 .1895874 .248806 /xb_one | -10.12438 2.713159 -3.73 0.000 -15.44383 -4.804931 ------------------------------------------------------------------------------ . . * Poisson: Sample vs avg predicted probabilities of y = 0, 1, ..., 10 . countfit docvis $xlist, maxcount(10) prm nograph noestimates nofit Comparison of Mean Observed and Predicted Count Maximum At Mean Model Difference Value |Diff| --------------------------------------------- PRM 0.102 0 0.040 PRM: Predicted and actual probabilities Count Actual Predicted |Diff| Pearson ------------------------------------------------ 0 0.109 0.007 0.102 5168.233 1 0.085 0.030 0.056 387.868 2 0.097 0.063 0.034 69.000 3 0.091 0.095 0.005 0.789 4 0.092 0.116 0.024 17.861 5 0.072 0.121 0.049 72.441 6 0.063 0.114 0.051 84.355 7 0.055 0.099 0.044 73.016 8 0.049 0.082 0.033 50.128 9 0.042 0.065 0.024 31.225 10 0.029 0.051 0.021 33.402 ------------------------------------------------ Sum 0.785 0.844 0.443 5988.318 . . *** 1. BASICS: DIAGNOSTICS . . * Residuals . quietly glm docvis $xlist, family(poisson) link(log) . predict rraw, response . predict rpearson, pearson . predict rdeviance, deviance . predict ranscombe, anscombe . predict hat, hat . predict cooksd, cooksd . summarize rraw rpearson rdeviance ranscombe hat cooksd, sep(10) Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- rraw | 3677 -7.08e-10 6.808178 -29.12996 136.9007 rpearson | 3677 -.0060737 2.509354 -4.914746 51.38051 rdeviance | 3677 -.3438453 2.210397 -6.087067 24.35213 ranscombe | 3677 -.3634087 2.254119 -6.153762 25.72892 hat | 3677 .0021757 .0015322 .0007023 .027556 cooksd | 3677 .0019357 .0177516 6.54e-11 .964198 . correlate rraw rpearson rdeviance ranscombe (obs=3677) | rraw rpearson rdevia~e ransco~e -------------+------------------------------------ rraw | 1.0000 rpearson | 0.9792 1.0000 rdeviance | 0.9454 0.9669 1.0000 ranscombe | 0.9435 0.9661 0.9998 1.0000 . . * An R-squared measure . capture drop yphat . quietly poisson docvis $xlist, vce(robust) . predict yphat, n . quietly correlate docvis yphat . display "Squared correlation between y and yhat = " r(rho)^2 Squared correlation between y and yhat = .1530784 . . * Overdispersion test against V[y|x] = E[y|x] + a*(E[y|x]^2) . quietly poisson docvis $xlist, vce(robust) . predict muhat, n . quietly generate ystar = ((docvis-muhat)^2 - docvis)/muhat . regress ystar muhat, noconstant noheader ------------------------------------------------------------------------------ ystar | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- muhat | .7047319 .1035926 6.80 0.000 .5016273 .9078365 ------------------------------------------------------------------------------ . . *** 1. BASICS: NEGATIVE BINOMIAL REGRESSION . . * Standard negative binomial (NB2) with default SEs . nbreg docvis $xlist, nolog Negative binomial regression Number of obs = 3677 LR chi2(7) = 773.44 Dispersion = mean Prob > chi2 = 0.0000 Log likelihood = -10589.339 Pseudo R2 = 0.0352 ------------------------------------------------------------------------------ docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .1640928 .0332186 4.94 0.000 .0989856 .2292001 medicaid | .100337 .0454209 2.21 0.027 .0113137 .1893603 age | .2941294 .0601588 4.89 0.000 .1762203 .4120384 age2 | -.0019282 .0004004 -4.82 0.000 -.0027129 -.0011434 educyr | .0286947 .0042241 6.79 0.000 .0204157 .0369737 actlim | .1895376 .0347601 5.45 0.000 .121409 .2576662 totchr | .2776441 .0121463 22.86 0.000 .2538378 .3014505 _cons | -10.29749 2.247436 -4.58 0.000 -14.70238 -5.892595 -------------+---------------------------------------------------------------- /lnalpha | -.4452773 .0306758 -.5054007 -.3851539 -------------+---------------------------------------------------------------- alpha | .6406466 .0196523 .6032638 .6803459 ------------------------------------------------------------------------------ Likelihood-ratio test of alpha=0: chibar2(01) = 8860.60 Prob>=chibar2 = 0.000 . . * NB2: Sample vs average predicted probabilities of y = 0, 1, ..., 10 . countfit docvis $xlist, maxcount(10) nbreg nograph noestimates nofit Comparison of Mean Observed and Predicted Count Maximum At Mean Model Difference Value |Diff| --------------------------------------------- NBRM -0.023 1 0.007 NBRM: Predicted and actual probabilities Count Actual Predicted |Diff| Pearson ------------------------------------------------ 0 0.109 0.091 0.018 12.708 1 0.085 0.108 0.023 17.288 2 0.097 0.105 0.008 2.270 3 0.091 0.096 0.005 1.086 4 0.092 0.085 0.007 2.333 5 0.072 0.074 0.001 0.072 6 0.063 0.063 0.000 0.004 7 0.055 0.054 0.001 0.088 8 0.049 0.046 0.003 0.694 9 0.042 0.039 0.003 0.873 10 0.029 0.033 0.004 1.456 ------------------------------------------------ Sum 0.785 0.794 0.073 38.872 . . * Generalized negative binomial with alpha parameterized . gnbreg docvis $xlist, lnalpha($xlist) nolog Generalized negative binomial regression Number of obs = 3677 LR chi2(7) = 733.91 Prob > chi2 = 0.0000 Log likelihood = -10518.541 Pseudo R2 = 0.0337 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- docvis | private | .1421503 .0330792 4.30 0.000 .0773163 .2069843 medicaid | .0896531 .0459235 1.95 0.051 -.0003554 .1796616 age | .2996728 .0597171 5.02 0.000 .1826295 .4167161 age2 | -.0019674 .0003968 -4.96 0.000 -.0027452 -.0011896 educyr | .0286626 .0042215 6.79 0.000 .0203887 .0369366 actlim | .1764802 .0337961 5.22 0.000 .1102412 .2427193 totchr | .2523805 .0113765 22.18 0.000 .2300829 .2746781 _cons | -10.42024 2.234642 -4.66 0.000 -14.80006 -6.040424 -------------+---------------------------------------------------------------- lnalpha | private | -.2143547 .0666325 -3.22 0.001 -.3449519 -.0837574 medicaid | .0561093 .0897534 0.63 0.532 -.1198041 .2320227 age | -.1606495 .1221281 -1.32 0.188 -.4000162 .0787172 age2 | .0010047 .0008126 1.24 0.216 -.000588 .0025974 educyr | -.0085659 .0086239 -0.99 0.321 -.0254685 .0083366 actlim | .096942 .0698175 1.39 0.165 -.0398978 .2337818 totchr | -.2584128 .0235344 -10.98 0.000 -.3045393 -.2122862 _cons | 6.587102 4.563726 1.44 0.149 -2.357637 15.53184 ------------------------------------------------------------------------------ . . * Comparison of Poisson and NB and various standard error estimates . quietly poisson docvis $xlist . estimates store PDEFAULT . quietly poisson docvis $xlist, vce(robust) . estimates store PROBUST . quietly glm docvis $xlist, family(poisson) link(log) scale(x2) . estimates store PPEARSON . quietly nbreg docvis $xlist . estimates store NBDEFAULT . quietly nbreg docvis $xlist, vce(robust) . estimates store NBROBUST . esttab PDEFAULT PROBUST PPEARSON NBDEFAULT NBROBUST, b(%10.4f) se mtitles pr2 -------------------------------------------------------------------------------------------- (1) (2) (3) (4) (5) PDEFAULT PROBUST PPEARSON NBDEFAULT NBROBUST -------------------------------------------------------------------------------------------- docvis private 0.1422*** 0.1422*** 0.1422*** 0.1641*** 0.1641*** (0.0143) (0.0364) (0.0360) (0.0332) (0.0369) medicaid 0.0970*** 0.0970 0.0970* 0.1003* 0.1003 (0.0189) (0.0568) (0.0475) (0.0454) (0.0567) age 0.2937*** 0.2937*** 0.2937*** 0.2941*** 0.2941*** (0.0260) (0.0630) (0.0652) (0.0602) (0.0646) age2 -0.0019*** -0.0019*** -0.0019*** -0.0019*** -0.0019*** (0.0002) (0.0004) (0.0004) (0.0004) (0.0004) educyr 0.0296*** 0.0296*** 0.0296*** 0.0287*** 0.0287*** (0.0019) (0.0048) (0.0047) (0.0042) (0.0049) actlim 0.1864*** 0.1864*** 0.1864*** 0.1895*** 0.1895*** (0.0146) (0.0397) (0.0366) (0.0348) (0.0394) totchr 0.2484*** 0.2484*** 0.2484*** 0.2776*** 0.2776*** (0.0046) (0.0126) (0.0117) (0.0121) (0.0132) _cons -10.1822*** -10.1822*** -10.1822*** -10.2975*** -10.2975*** (0.9720) (2.3692) (2.4415) (2.2474) (2.4241) -------------------------------------------------------------------------------------------- lnalpha _cons -0.4453*** -0.4453*** (0.0307) (0.0378) -------------------------------------------------------------------------------------------- N 3677 3677 3677 3677 3677 pseudo R-sq 0.130 0.130 0.035 0.035 -------------------------------------------------------------------------------------------- Standard errors in parentheses * p<0.05, ** p<0.01, *** p<0.001 . . *********** 2: COUNT REGRESSION: ADDITIONAL (TRUNCATED/CENSORED, HURDLE, ZERO-INFLATED, . **** FINITE MIXTURE, ENDOGENOUS) . . *** 2. ADDITONAL: TRUNCATED AND CENSORED . . use mus17data.dta, clear . global xlist private medicaid age age2 educyr actlim totchr . . * Zero-truncated negative binomial . ztnb docvis $xlist if docvis>0, nolog Zero-truncated negative binomial regression Number of obs = 3276 LR chi2(7) = 509.10 Dispersion = mean Prob > chi2 = 0.0000 Log likelihood = -9452.899 Pseudo R2 = 0.0262 ------------------------------------------------------------------------------ docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .1095567 .0345239 3.17 0.002 .0418911 .1772223 medicaid | .0972309 .0470358 2.07 0.039 .0050425 .1894193 age | .2719032 .0625359 4.35 0.000 .1493352 .3944712 age2 | -.0017959 .000416 -4.32 0.000 -.0026113 -.0009805 educyr | .0265974 .0043937 6.05 0.000 .0179859 .035209 actlim | .1955384 .0355161 5.51 0.000 .1259281 .2651488 totchr | .2226969 .0124128 17.94 0.000 .1983683 .2470254 _cons | -9.19017 2.337591 -3.93 0.000 -13.77176 -4.608576 -------------+---------------------------------------------------------------- /lnalpha | -.5259629 .0418671 -.6080209 -.443905 -------------+---------------------------------------------------------------- alpha | .590986 .0247429 .5444273 .6415264 ------------------------------------------------------------------------------ Likelihood-ratio test of alpha=0: chibar2(01) = 7089.37 Prob>=chibar2 = 0.000 . estimates store ZTNB . . * Censored negative binomial at y <= 10 . generate docviscens10 = docvis . replace docviscens10 = 10 if docvis > 10 (791 real changes made) . . * Censored Negbin ML program lfnbcens10 with censoring from right at 10 . program lfnbcens10 1. version 10.1 2. args lnf theta1 a // theta1=x'b, a=alpha, lnf=lnf(y) 3. tempvar mu sumgammaratios probyle10 cens 4. local y $ML_y1 // Define y so program more readable 5. generate double `mu' = exp(`theta1') 6. generate double `sumgammaratios' = exp(lngamma(0+(1/`a')))/1 /// > + exp(lngamma(1+(1/`a')))/1 + exp(lngamma(2+(1/`a')))/2 /// > + exp(lngamma(3+(1/`a')))/6 + exp(lngamma(4+(1/`a')))/24 /// > + exp(lngamma(5+(1/`a')))/120 + exp(lngamma(6+(1/`a')))/720 /// > + exp(lngamma(7+(1/`a')))/5040 + exp(lngamma(8+(1/`a')))/40320 /// > + exp(lngamma(9+(1/`a')))/362880 + exp(lngamma(10+(1/`a')))/3628800 7. generate double `probyle10' = (`sumgammaratios'/exp(lngamma((1/`a')))) /// > *((1/`a')^(1/`a'))*(`mu'^`y')/((1/`a'+`mu')^(1/`a'+`y')) 8. generate double `cens' = `y' >= 10 9. quietly replace `lnf' = (1-`cens') * (lngamma(`y'+(1/`a')) - lngamma(1/`a') /// > - lnfactorial(`y') - (`y'+(1/`a'))*ln(1+`a'*`mu') /// > + `y'*ln(`a') + `y'*ln(`mu') ) + `cens'*ln(`probyle10') 10. end . * Command lfnbcens10 implemented for negative binomial MLE . ml model lf lfnbcens10 (docviscens10 = $xlist) () . ml maximize, nolog (3677 missing values generated) (3677 missing values generated) initial: log likelihood = - (could not be evaluated) (3677 missing values generated) (3677 missing values generated) feasible: log likelihood = -12293.433 rescale: log likelihood = -8141.2566 rescale eq: log likelihood = -8118.0172 Number of obs = 3677 Wald chi2(7) = 293.34 Log likelihood = -7796.8328 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ docviscens10 | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- eq1 | private | .1258389 .0408548 3.08 0.002 .045765 .2059128 medicaid | .0652582 .0558857 1.17 0.243 -.0442758 .1747922 age | .2275648 .0735028 3.10 0.002 .0835019 .3716277 age2 | -.0014782 .0004892 -3.02 0.003 -.002437 -.0005194 educyr | .0196804 .0052061 3.78 0.000 .0094767 .0298842 actlim | .097664 .0429777 2.27 0.023 .0134293 .1818988 totchr | .2147348 .0150827 14.24 0.000 .1851732 .2442964 _cons | -7.8 2.746158 -2.84 0.005 -13.18237 -2.41763 -------------+---------------------------------------------------------------- eq2 | _cons | 1.01341 .0378353 26.78 0.000 .9392546 1.087566 ------------------------------------------------------------------------------ . estimates store NBCENS10 . . /* Following did not work > generate double `sumgammaratios' = 0 > forvalues i = 0/10 { > quietly replace `sumgammaratios' = `sumgammaratios' + exp(lngamma(`i'+(1/`a')))/exp(lngamma(`i'+1)) > } > */ . . * Comparison with regular negative binomial . quietly nbreg docvis $xlist, nolog . estimates store NBREG . estimates table NBREG ZTNB NBCENS10, equation(1) b(%10.4f) se stats(N ll) ----------------------------------------------------- Variable | NBREG ZTNB NBCENS10 -------------+--------------------------------------- #1 | private | 0.1641 0.1096 0.1258 | 0.0332 0.0345 0.0409 medicaid | 0.1003 0.0972 0.0653 | 0.0454 0.0470 0.0559 age | 0.2941 0.2719 0.2276 | 0.0602 0.0625 0.0735 age2 | -0.0019 -0.0018 -0.0015 | 0.0004 0.0004 0.0005 educyr | 0.0287 0.0266 0.0197 | 0.0042 0.0044 0.0052 actlim | 0.1895 0.1955 0.0977 | 0.0348 0.0355 0.0430 totchr | 0.2776 0.2227 0.2147 | 0.0121 0.0124 0.0151 _cons | -10.2975 -9.1902 -7.8000 | 2.2474 2.3376 2.7462 -------------+--------------------------------------- lnalpha | _cons | -0.4453 -0.5260 | 0.0307 0.0419 -------------+--------------------------------------- eq2 | _cons | 1.0134 | 0.0378 -------------+--------------------------------------- Statistics | N | 3677 3276 3677 ll | -1.059e+04 -9452.8990 -7796.8328 ----------------------------------------------------- legend: b/se . . * Ordered logit . * Do for docviscens10 which has 10 categories . ologit docviscens10 $xlist, nolog Ordered logistic regression Number of obs = 3677 LR chi2(7) = 882.51 Prob > chi2 = 0.0000 Log likelihood = -7883.0551 Pseudo R2 = 0.0530 ------------------------------------------------------------------------------ docviscens10 | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .3343672 .06573 5.09 0.000 .2055389 .4631956 medicaid | .1576941 .0910149 1.73 0.083 -.0206919 .33608 age | .5883239 .1163094 5.06 0.000 .3603616 .8162862 age2 | -.0038291 .0007736 -4.95 0.000 -.0053453 -.0023128 educyr | .0541312 .0084618 6.40 0.000 .0375464 .0707159 actlim | .2771719 .0707753 3.92 0.000 .1384549 .415889 totchr | .6188858 .0258162 23.97 0.000 .5682869 .6694846 -------------+---------------------------------------------------------------- /cut1 | 22.06108 4.348472 13.53823 30.58393 /cut2 | 22.81988 4.348907 14.29618 31.34358 /cut3 | 23.42779 4.349421 14.90308 31.9525 /cut4 | 23.89882 4.349988 15.373 32.42464 /cut5 | 24.34237 4.350657 15.81524 32.8695 /cut6 | 24.68387 4.351224 16.15562 33.21211 /cut7 | 24.98738 4.351741 16.45813 33.51664 /cut8 | 25.26841 4.352181 16.73829 33.79853 /cut9 | 25.53461 4.352499 17.00387 34.06535 /cut10 | 25.78217 4.352753 17.25093 34.31341 ------------------------------------------------------------------------------ . . *** 2. ADDITONAL: HURDLE MODEL . . * Hurdle logit-nb model manually . logit docvis $xlist, nolog Logistic regression Number of obs = 3677 LR chi2(7) = 453.08 Prob > chi2 = 0.0000 Log likelihood = -1040.3258 Pseudo R2 = 0.1788 ------------------------------------------------------------------------------ docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .6586978 .1264608 5.21 0.000 .4108393 .9065563 medicaid | .0554225 .1726693 0.32 0.748 -.2830032 .3938482 age | .5428779 .2238845 2.42 0.015 .1040724 .9816834 age2 | -.0034989 .0014957 -2.34 0.019 -.0064304 -.0005673 educyr | .047035 .0155706 3.02 0.003 .0165171 .0775529 actlim | .1623927 .1523743 1.07 0.287 -.1362553 .4610408 totchr | 1.050562 .0671922 15.64 0.000 .9188676 1.182256 _cons | -20.94163 8.335137 -2.51 0.012 -37.2782 -4.605058 ------------------------------------------------------------------------------ . predict dvplogit, p // Logit Pr[y=1] . * Second step uses positives only . ztnb docvis $xlist if docvis>0, nolog Zero-truncated negative binomial regression Number of obs = 3276 LR chi2(7) = 509.10 Dispersion = mean Prob > chi2 = 0.0000 Log likelihood = -9452.899 Pseudo R2 = 0.0262 ------------------------------------------------------------------------------ docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .1095567 .0345239 3.17 0.002 .0418911 .1772223 medicaid | .0972309 .0470358 2.07 0.039 .0050425 .1894193 age | .2719032 .0625359 4.35 0.000 .1493352 .3944712 age2 | -.0017959 .000416 -4.32 0.000 -.0026113 -.0009805 educyr | .0265974 .0043937 6.05 0.000 .0179859 .035209 actlim | .1955384 .0355161 5.51 0.000 .1259281 .2651488 totchr | .2226969 .0124128 17.94 0.000 .1983683 .2470254 _cons | -9.19017 2.337591 -3.93 0.000 -13.77176 -4.608576 -------------+---------------------------------------------------------------- /lnalpha | -.5259629 .0418671 -.6080209 -.443905 -------------+---------------------------------------------------------------- alpha | .590986 .0247429 .5444273 .6415264 ------------------------------------------------------------------------------ Likelihood-ratio test of alpha=0: chibar2(01) = 7089.37 Prob>=chibar2 = 0.000 . predict dvztnbcm, cm // ztnb E[y|y>0] = E[y](1-/Pr[|y=0] . /* Following computes this manually > scalar alpha = e(alpha) > predict dvztnb, n > generate pryeq0 = (1+alpha*dvztnb)^(-1/alpha) > generate dvztnbcmman = dvztnb/(1-pryeq0) > */ . generate dvhurdle = dvplogit*dvztnbcm . . * Same hurdle logit-nb model using the user-written hnblogit command . hnblogit docvis $xlist, nolog Negative Binomial-Logit Hurdle Regression Number of obs = 3677 Wald chi2(7) = 309.90 Log likelihood = -10493.225 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- logit | private | .6586978 .1264608 5.21 0.000 .4108393 .9065563 medicaid | .0554225 .1726694 0.32 0.748 -.2830032 .3938483 age | .542878 .2238845 2.42 0.015 .1040724 .9816835 age2 | -.0034989 .0014957 -2.34 0.019 -.0064304 -.0005673 educyr | .047035 .0155706 3.02 0.003 .0165171 .0775529 actlim | .1623927 .1523743 1.07 0.287 -.1362554 .4610408 totchr | 1.050562 .0671922 15.64 0.000 .9188676 1.182256 _cons | -20.94163 8.335138 -2.51 0.012 -37.2782 -4.605058 -------------+---------------------------------------------------------------- negbinomial | private | .1095566 .0345239 3.17 0.002 .041891 .1772222 medicaid | .0972308 .0470358 2.07 0.039 .0050423 .1894193 age | .2719031 .0625359 4.35 0.000 .149335 .3944712 age2 | -.0017959 .000416 -4.32 0.000 -.0026113 -.0009805 educyr | .0265974 .0043937 6.05 0.000 .0179859 .035209 actlim | .1955384 .0355161 5.51 0.000 .125928 .2651487 totchr | .2226967 .0124128 17.94 0.000 .1983681 .2470252 _cons | -9.190165 2.337592 -3.93 0.000 -13.77176 -4.608569 -------------+---------------------------------------------------------------- /lnalpha | -.525962 .0418671 -12.56 0.000 -.60802 -.443904 ------------------------------------------------------------------------------ AIC Statistic = 5.712 . estimates store HURDLENB . . *** 2. ADDITONAL: ZERO-INFLATED . . * Zero-inflated negative binomial . zinb docvis $xlist, inflate($xlist) vuong nolog Zero-inflated negative binomial regression Number of obs = 3677 Nonzero obs = 3276 Zero obs = 401 Inflation model = logit LR chi2(7) = 588.93 Log likelihood = -10492.88 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- docvis | private | .1289797 .032987 3.91 0.000 .0643264 .193633 medicaid | .1091956 .044511 2.45 0.014 .0219556 .1964356 age | .2847325 .0589577 4.83 0.000 .1691776 .4002874 age2 | -.0018781 .0003922 -4.79 0.000 -.0026469 -.0011093 educyr | .0253991 .0041432 6.13 0.000 .0172786 .0335196 actlim | .1737716 .0336464 5.16 0.000 .1078258 .2397173 totchr | .229991 .0120795 19.04 0.000 .2063156 .2536663 _cons | -9.680235 2.204161 -4.39 0.000 -14.00031 -5.36016 -------------+---------------------------------------------------------------- inflate | private | -.9152675 .2758402 -3.32 0.001 -1.455904 -.3746307 medicaid | .3487142 .3372848 1.03 0.301 -.3123519 1.00978 age | -.4357439 .5156094 -0.85 0.398 -1.44632 .5748319 age2 | .002805 .0034886 0.80 0.421 -.0040326 .0096426 educyr | -.08423 .0339273 -2.48 0.013 -.1507263 -.0177336 actlim | -.8241735 .4825621 -1.71 0.088 -1.769978 .1216309 totchr | -2.985208 .6860952 -4.35 0.000 -4.32993 -1.640486 _cons | 17.09618 18.97318 0.90 0.368 -20.09057 54.28294 -------------+---------------------------------------------------------------- /lnalpha | -.5848279 .0349792 -16.72 0.000 -.6533859 -.5162699 -------------+---------------------------------------------------------------- alpha | .5572017 .0194905 .5202812 .5967423 ------------------------------------------------------------------------------ Vuong test of zinb vs. standard negative binomial: z = 6.48 Pr>z = 0.0000 . estimates store ZINB . predict dvzinb, n . . quietly nbreg docvis $xlist . estimates store NBREG . predict dvnbreg, n . . *** 2. ADDITONAL: MODEL COMPARISON . . * Comparison of NB and ZINB using countfit . countfit docvis $xlist, nbreg zinb nograph noestimates Comparison of Mean Observed and Predicted Count Maximum At Mean Model Difference Value |Diff| --------------------------------------------- NBRM -0.023 1 0.007 ZINB 0.009 4 0.003 NBRM: Predicted and actual probabilities Count Actual Predicted |Diff| Pearson ------------------------------------------------ 0 0.109 0.091 0.018 12.708 1 0.085 0.108 0.023 17.288 2 0.097 0.105 0.008 2.270 3 0.091 0.096 0.005 1.086 4 0.092 0.085 0.007 2.333 5 0.072 0.074 0.001 0.072 6 0.063 0.063 0.000 0.004 7 0.055 0.054 0.001 0.088 8 0.049 0.046 0.003 0.694 9 0.042 0.039 0.003 0.873 ------------------------------------------------ Sum 0.756 0.761 0.070 37.416 ZINB: Predicted and actual probabilities Count Actual Predicted |Diff| Pearson ------------------------------------------------ 0 0.109 0.113 0.004 0.627 1 0.085 0.087 0.002 0.170 2 0.097 0.093 0.004 0.782 3 0.091 0.090 0.001 0.038 4 0.092 0.083 0.009 3.931 5 0.072 0.074 0.002 0.155 6 0.063 0.065 0.002 0.304 7 0.055 0.057 0.002 0.176 8 0.049 0.049 0.000 0.000 9 0.042 0.042 0.000 0.001 ------------------------------------------------ Sum 0.756 0.753 0.027 6.184 Tests and Fit Statistics ------------------------------------------------------------------------- NBRM BIC= -8935.061 AIC= 5.765 Prefer Over Evidence ------------------------------------------------------------------------- vs ZINB BIC= -9062.301 dif= 127.240 ZINB NBRM Very strong AIC= 5.717 dif= 0.048 ZINB NBRM Vuong= 6.484 prob= 0.000 ZINB NBRM p=0.000 . . * Compare LL, AIC, BIC across models . estimates table NBREG HURDLENB ZINB, equation(1) b(%12.4f) se /// > stats(N ll aic bic) stfmt(%12.1f) newpanel ----------------------------------------------------------- Variable | NBREG HURDLENB ZINB -------------+--------------------------------------------- #1 | private | 0.1641 0.6587 0.1290 | 0.0332 0.1265 0.0330 medicaid | 0.1003 0.0554 0.1092 | 0.0454 0.1727 0.0445 age | 0.2941 0.5429 0.2847 | 0.0602 0.2239 0.0590 age2 | -0.0019 -0.0035 -0.0019 | 0.0004 0.0015 0.0004 educyr | 0.0287 0.0470 0.0254 | 0.0042 0.0156 0.0041 actlim | 0.1895 0.1624 0.1738 | 0.0348 0.1524 0.0336 totchr | 0.2776 1.0506 0.2300 | 0.0121 0.0672 0.0121 _cons | -10.2975 -20.9416 -9.6802 | 2.2474 8.3351 2.2042 -------------+--------------------------------------------- lnalpha | _cons | -0.4453 -0.5260 -0.5848 | 0.0307 0.0419 0.0350 -------------+--------------------------------------------- negbinomial | private | 0.1096 | 0.0345 medicaid | 0.0972 | 0.0470 age | 0.2719 | 0.0625 age2 | -0.0018 | 0.0004 educyr | 0.0266 | 0.0044 actlim | 0.1955 | 0.0355 totchr | 0.2227 | 0.0124 _cons | -9.1902 | 2.3376 -------------+--------------------------------------------- inflate | private | -0.9153 | 0.2758 medicaid | 0.3487 | 0.3373 age | -0.4357 | 0.5156 age2 | 0.0028 | 0.0035 educyr | -0.0842 | 0.0339 actlim | -0.8242 | 0.4826 totchr | -2.9852 | 0.6861 _cons | 17.0962 | 18.9732 ----------------------------------------------------------- ----------------------------------------------------------- Statistics | NBREG HURDLENB ZINB -------------+--------------------------------------------- N | 3677 3677 3677 ll | -10589.3 -10493.2 -10492.9 aic | 21196.7 21020.4 21019.8 bic | 21252.6 21126.0 21125.3 ----------------------------------------------------------- legend: b/se . . * Compare predicted means across models . summarize docvis dvnbreg dvhurdle dvzinb Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- docvis | 3677 6.822682 7.394937 0 144 dvnbreg | 3677 6.890034 3.486562 2.078925 41.31503 dvhurdle | 3677 6.840676 3.134925 1.35431 31.86874 dvzinb | 3677 6.838704 3.135122 .9473827 32.98153 . correlate docvis dvnbreg dvhurdle dvzinb (obs=3677) | docvis dvnbreg dvhurdle dvzinb -------------+------------------------------------ docvis | 1.0000 dvnbreg | 0.3870 1.0000 dvhurdle | 0.3990 0.9894 1.0000 dvzinb | 0.3983 0.9882 0.9982 1.0000 . . *** 2. ADDITONAL: FINITE MIXTURE MODEL . . * Finite-mixture model using fmm command with constant probabilities . use mus17data.dta, clear . fmm docvis $xlist, vce(robust) components(2) mixtureof(poisson) Fitting Poisson model: Iteration 0: log likelihood = -15019.656 Iteration 1: log likelihood = -15019.64 Iteration 2: log likelihood = -15019.64 Fitting 2 component Poisson model: Iteration 0: log pseudolikelihood = -14985.068 (not concave) Iteration 1: log pseudolikelihood = -12233.072 (not concave) Iteration 2: log pseudolikelihood = -11752.598 Iteration 3: log pseudolikelihood = -11518.01 Iteration 4: log pseudolikelihood = -11502.758 Iteration 5: log pseudolikelihood = -11502.686 Iteration 6: log pseudolikelihood = -11502.686 2 component Poisson regression Number of obs = 3677 Wald chi2(14) = 576.86 Log pseudolikelihood = -11502.686 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- component1 | private | .2077415 .0560256 3.71 0.000 .0979333 .3175497 medicaid | .1071618 .0964233 1.11 0.266 -.0818245 .2961481 age | .3798087 .100821 3.77 0.000 .1822032 .5774143 age2 | -.0024869 .0006711 -3.71 0.000 -.0038022 -.0011717 educyr | .029099 .0067908 4.29 0.000 .0157893 .0424087 actlim | .1244235 .0558883 2.23 0.026 .0148844 .2339625 totchr | .3191166 .0184744 17.27 0.000 .2829074 .3553259 _cons | -14.25713 3.759845 -3.79 0.000 -21.62629 -6.887972 -------------+---------------------------------------------------------------- component2 | private | .138229 .0614901 2.25 0.025 .0177106 .2587474 medicaid | .1269723 .1329626 0.95 0.340 -.1336297 .3875742 age | .2628874 .1140355 2.31 0.021 .0393819 .486393 age2 | -.0017418 .0007542 -2.31 0.021 -.00322 -.0002636 educyr | .0241679 .0076208 3.17 0.002 .0092314 .0391045 actlim | .1831598 .0622267 2.94 0.003 .0611977 .3051218 totchr | .1970511 .0263763 7.47 0.000 .1453545 .2487477 _cons | -8.051256 4.28211 -1.88 0.060 -16.44404 .3415266 -------------+---------------------------------------------------------------- /imlogitpi1 | .877227 .0952018 9.21 0.000 .690635 1.063819 ------------------------------------------------------------------------------ pi1 | .7062473 .0197508 .6661082 .7434197 pi2 | .2937527 .0197508 .2565803 .3338918 ------------------------------------------------------------------------------ . scalar pi = e(pi1_est) . . * Predict y for two components . quietly fmm docvis $xlist, vce(robust) components(2) mixtureof(poisson) . quietly predict dvfit1, equation(component1) . quietly predict dvfit2, equation(component2) . quietly predict dvcombinedfit . summarize dvfit1 dvfit2 dvcombinedfit docvis Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- dvfit1 | 3677 3.801692 2.176922 .9815563 27.28715 dvfit2 | 3677 13.95943 5.077463 5.615584 55.13366 dvcombined~t | 3677 6.785555 3.013985 2.342815 35.46714 docvis | 3677 6.822682 7.394937 0 144 . . * Create histograms of fitted values . quietly histogram dvfit1, name(comp_1, replace) . quietly histogram dvfit2, name(comp_2, replace) . quietly graph combine comp_1 comp_2, ysize(3) xsize(6) ycommon xcommon . quietly graph export countfig21finitemixture.wmf, replace . . * 2-component mixture of NB1 . fmm docvis $xlist, vce(robust) components(2) mixtureof(negbin1) Fitting Negative Binomial-1 model: Iteration 0: log likelihood = -15019.656 Iteration 1: log likelihood = -15019.64 Iteration 2: log likelihood = -15019.64 Iteration 0: log likelihood = -12739.566 Iteration 1: log likelihood = -11125.786 Iteration 2: log likelihood = -10976.314 Iteration 3: log likelihood = -10976.058 Iteration 4: log likelihood = -10976.058 Iteration 0: log likelihood = -10976.058 Iteration 1: log likelihood = -10566.829 Iteration 2: log likelihood = -10531.205 Iteration 3: log likelihood = -10531.054 Iteration 4: log likelihood = -10531.054 Fitting 2 component Negative Binomial-1 model: Iteration 0: log pseudolikelihood = -10531.611 (not concave) Iteration 1: log pseudolikelihood = -10529.012 (not concave) Iteration 2: log pseudolikelihood = -10515.85 (not concave) Iteration 3: log pseudolikelihood = -10500.668 (not concave) Iteration 4: log pseudolikelihood = -10495.501 (not concave) Iteration 5: log pseudolikelihood = -10494.709 Iteration 6: log pseudolikelihood = -10493.449 Iteration 7: log pseudolikelihood = -10493.333 Iteration 8: log pseudolikelihood = -10493.324 Iteration 9: log pseudolikelihood = -10493.324 2 component Negative Binomial-1 regression Number of obs = 3677 Wald chi2(14) = 560.31 Log pseudolikelihood = -10493.324 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- component1 | private | .137827 .0610423 2.26 0.024 .0181863 .2574676 medicaid | .0379753 .0628139 0.60 0.545 -.0851377 .1610883 age | .253357 .0633567 4.00 0.000 .12918 .3775339 age2 | -.0016569 .0004261 -3.89 0.000 -.002492 -.0008218 educyr | .0228524 .0055063 4.15 0.000 .0120602 .0336446 actlim | .1060655 .0514198 2.06 0.039 .0052845 .2068464 totchr | .2434641 .0294843 8.26 0.000 .1856759 .3012523 _cons | -8.645394 2.352187 -3.68 0.000 -13.2556 -4.035192 -------------+---------------------------------------------------------------- component2 | private | .372013 .5124233 0.73 0.468 -.6323182 1.376344 medicaid | .3344168 .856897 0.39 0.696 -1.34507 2.013904 age | .5260549 .6902627 0.76 0.446 -.8268352 1.878945 age2 | -.0034424 .0047508 -0.72 0.469 -.0127539 .005869 educyr | .0457671 .0499026 0.92 0.359 -.0520402 .1435743 actlim | .3599301 .3852059 0.93 0.350 -.3950595 1.11492 totchr | .4150389 .1332826 3.11 0.002 .1538097 .6762681 _cons | -19.3304 25.16197 -0.77 0.442 -68.64696 29.98615 -------------+---------------------------------------------------------------- /imlogitpi1 | 2.382195 2.159316 1.10 0.270 -1.849987 6.614377 /lndelta1 | 1.210492 .2343343 5.17 0.000 .7512047 1.669778 /lndelta2 | 2.484476 .7928709 3.13 0.002 .9304772 4.038474 ------------------------------------------------------------------------------ delta1 | 3.355133 .7862229 2.119552 5.310991 delta2 | 11.99483 9.510352 2.535719 56.7397 pi1 | .9154595 .1671169 .1358744 .9986608 pi2 | .0845405 .1671169 .0013392 .8641256 ------------------------------------------------------------------------------ . . *** 2. ADDITONAL: ENDOGENOUS REGRESSORS . . use mus17data.dta, clear . global xlist private medicaid age age2 educyr actlim totchr . . * Nonlinear 2SLS IV estimator for Poisson: computation using command optimize . generate cons = 1 . local y docvis . local xlist private medicaid age age2 educyr actlim totchr cons . local zlist income ssiratio medicaid age age2 educyr actlim totchr cons . mata ------------------------------------------------- mata (type end to exit) ------------------------------------- : void pgmm(todo, b, y, X, Z, Qb, g, H) > { > Xb = X*b' > mu = exp(Xb) > h = Z'(y-mu) > W = cholinv(cross(Z,Z)) > Qb = h'W*h > if (todo == 0) return > G = -(mu:*Z)'X > g = (G'W*h)' > if (todo == 1) return > H = G'W*G > _makesymmetric(H) > } : st_view(y=., ., "`y'") : st_view(X=., ., tokens("`xlist'")) : st_view(Z=., ., tokens("`zlist'")) : S = optimize_init() : optimize_init_which(S,"min") : optimize_init_evaluator(S, &pgmm()) : optimize_init_evaluatortype(S, "d2") : optimize_init_argument(S, 1, y) : optimize_init_argument(S, 2, X) : optimize_init_argument(S, 3, Z) : optimize_init_params(S, J(1,cols(X),0)) : optimize_init_technique(S,"nr") : b = optimize(S) Iteration 0: f(p) = 156836.04 Iteration 1: f(p) = 21765.741 Iteration 2: f(p) = 2087.4467 Iteration 3: f(p) = 186.55764 Iteration 4: f(p) = 182.298 Iteration 5: f(p) = 182.29546 Iteration 6: f(p) = 182.29541 Iteration 7: f(p) = 182.29541 : // Compute robust estimate of VCE : Xb = X*b' : mu = exp(Xb) : h = Z'(y-mu) : W = cholinv(cross(Z,Z)) : G = -(mu:*Z)'X : Shat = ((y-mu):*Z)'((y-mu):*Z)*rows(X)/(rows(X)-cols(X)) : Vb = luinv(G'W*G)*G'W*Shat*W*G*luinv(G'W*G) : st_matrix("b",b) : st_matrix("Vb",Vb) : end --------------------------------------------------------------------------------------------------------------- . . * Nonlinear IV estimator for Poisson: formatted results . matrix colnames b = `xlist' . matrix colnames Vb = `xlist' . matrix rownames Vb = `xlist' . ereturn post b Vb . ereturn display ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .5920658 .3401151 1.74 0.082 -.0745475 1.258679 medicaid | .3186961 .1912099 1.67 0.096 -.0560685 .6934607 age | .3323219 .0706128 4.71 0.000 .1939233 .4707205 age2 | -.002176 .0004648 -4.68 0.000 -.003087 -.001265 educyr | .0190875 .0092318 2.07 0.039 .0009935 .0371815 actlim | .2084997 .0434233 4.80 0.000 .1233916 .2936079 totchr | .2418424 .013001 18.60 0.000 .2163608 .267324 cons | -11.86341 2.735737 -4.34 0.000 -17.22535 -6.50146 ------------------------------------------------------------------------------ . . * ivpois is also NL2SLS but works with multiplicative form of Mullahy's condition . ivpois docvis private $xlist2, exog(income ssiratio) endog(private) ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- docvis | private | .0412665 .1334695 0.31 0.757 -.2203288 .3028619 _cons | 1.898924 .0708494 26.80 0.000 1.760062 2.037786 ------------------------------------------------------------------------------ . . * Control function approach to endogeneity . * First stage is reduced form to get predicted residuals . global xlist2 medicaid age age2 educyr actlim totchr . regress private $xlist2 income ssiratio, vce(robust) Linear regression Number of obs = 3677 F( 8, 3668) = 249.61 Prob > F = 0.0000 R-squared = 0.2108 Root MSE = .44472 ------------------------------------------------------------------------------ | Robust private | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- medicaid | -.3934477 .0173623 -22.66 0.000 -.4274884 -.3594071 age | -.0831201 .0293734 -2.83 0.005 -.1407098 -.0255303 age2 | .0005257 .0001959 2.68 0.007 .0001417 .0009098 educyr | .0212523 .0020492 10.37 0.000 .0172345 .02527 actlim | -.0300936 .0176874 -1.70 0.089 -.0647718 .0045845 totchr | .0185063 .005743 3.22 0.001 .0072465 .0297662 income | .0027416 .0004736 5.79 0.000 .0018131 .0036702 ssiratio | -.0647637 .0211178 -3.07 0.002 -.1061675 -.0233599 _cons | 3.531058 1.09581 3.22 0.001 1.3826 5.679516 ------------------------------------------------------------------------------ . predict lpuhat, residual . . * Second-stage Poisson with robust SEs . poisson docvis private $xlist2 lpuhat, vce(robust) nolog Poisson regression Number of obs = 3677 Wald chi2(8) = 718.87 Prob > chi2 = 0.0000 Log pseudolikelihood = -15010.614 Pseudo R2 = 0.1303 ------------------------------------------------------------------------------ | Robust docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .5505541 .2453175 2.24 0.025 .0697407 1.031368 medicaid | .2628822 .1197162 2.20 0.028 .0282428 .4975217 age | .3350604 .0696064 4.81 0.000 .1986344 .4714865 age2 | -.0021923 .0004576 -4.79 0.000 -.0030893 -.0012954 educyr | .018606 .0080461 2.31 0.021 .002836 .034376 actlim | .2053417 .0414248 4.96 0.000 .1241505 .286533 totchr | .24147 .0129175 18.69 0.000 .2161523 .2667878 lpuhat | -.4166838 .249347 -1.67 0.095 -.9053949 .0720272 _cons | -11.90647 2.661445 -4.47 0.000 -17.1228 -6.69013 ------------------------------------------------------------------------------ . . * Program and bootstrap for Poisson two-step estimator . program endogtwostep, eclass 1. version 10.1 2. tempname b 3. capture drop lpuhat2 4. regress private $xlist2 income ssiratio 5. predict lpuhat2, residual 6. poisson docvis private $xlist2 lpuhat2 7. matrix `b' = e(b) 8. ereturn post `b' 9. end . bootstrap _b, reps(400) seed(10101) nodots nowarn: endogtwostep Bootstrap results Number of obs = 3677 Replications = 400 ------------------------------------------------------------------------------ | Observed Bootstrap Normal-based | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .5505541 .2567815 2.14 0.032 .0472716 1.053837 medicaid | .2628822 .1205813 2.18 0.029 .0265473 .4992172 age | .3350604 .0707275 4.74 0.000 .1964371 .4736838 age2 | -.0021923 .0004667 -4.70 0.000 -.0031071 -.0012776 educyr | .018606 .0083042 2.24 0.025 .0023301 .034882 actlim | .2053417 .0412756 4.97 0.000 .124443 .2862405 totchr | .24147 .0134522 17.95 0.000 .2151042 .2678359 lpuhat2 | -.4166838 .2617964 -1.59 0.111 -.9297953 .0964276 _cons | -11.90647 2.698704 -4.41 0.000 -17.19583 -6.617104 ------------------------------------------------------------------------------ . . *********** 3: COUNT REGRESSION: BEYOND CROSS-SECTION . . *** 3. BEYOND CROSS-SECTION: PANEL . . * Describe dependent variables and regressors . use mus18data.dta, clear . describe mdu lcoins ndisease female age lfam child id year storage display value variable name type format label variable label --------------------------------------------------------------------------------------------------------------- mdu float %9.0g number face-to-fact md visits lcoins float %9.0g log(coinsurance+1) ndisease float %9.0g count of chronic diseases -- ba female float %9.0g female age float %9.0g age that year lfam float %9.0g log of family size child float %9.0g child id float %9.0g person id, leading digit is sit year float %9.0g study year . . * Summarize dependent variables and regressors . summarize mdu lcoins ndisease female age lfam child id year Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- mdu | 20186 2.860696 4.504765 0 77 lcoins | 20186 2.383588 2.041713 0 4.564348 ndisease | 20186 11.2445 6.741647 0 58.6 female | 20186 .5169424 .4997252 0 1 age | 20186 25.71844 16.76759 0 64.27515 -------------+-------------------------------------------------------- lfam | 20186 1.248404 .5390681 0 2.639057 child | 20186 .4014168 .4901972 0 1 id | 20186 357971.2 180885.6 125024 632167 year | 20186 2.420044 1.217237 1 5 . . * Panel description of dataset . xtset id year panel variable: id (unbalanced) time variable: year, 1 to 5, but with gaps delta: 1 unit . xtdescribe id: 125024, 125025, ..., 632167 n = 5908 year: 1, 2, ..., 5 T = 5 Delta(year) = 1 unit Span(year) = 5 periods (id*year uniquely identifies each observation) Distribution of T_i: min 5% 25% 50% 75% 95% max 1 2 3 3 5 5 5 Freq. Percent Cum. | Pattern ---------------------------+--------- 3710 62.80 62.80 | 111.. 1584 26.81 89.61 | 11111 156 2.64 92.25 | 1.... 147 2.49 94.74 | 11... 79 1.34 96.07 | ..1.. 66 1.12 97.19 | .11.. 33 0.56 97.75 | ..111 33 0.56 98.31 | .1111 29 0.49 98.80 | ...11 71 1.20 100.00 | (other patterns) ---------------------------+--------- 5908 100.00 | XXXXX . . * Panel summary of dependent variable . xtsum mdu Variable | Mean Std. Dev. Min Max | Observations -----------------+--------------------------------------------+---------------- mdu overall | 2.860696 4.504765 0 77 | N = 20186 between | 3.785971 0 63.33333 | n = 5908 within | 2.575881 -34.47264 40.0607 | T-bar = 3.41672 . . * Year-to-year transitions in doctor visits . generate mdushort = mdu . replace mdushort = 4 if mdu >= 4 (4039 real changes made) . xttrans mdushort | mdushort mdushort | 0 1 2 3 4 | Total -----------+-------------------------------------------------------+---------- 0 | 58.87 19.61 9.21 4.88 7.42 | 100.00 1 | 33.16 24.95 17.58 10.14 14.16 | 100.00 2 | 23.55 24.26 17.90 12.10 22.19 | 100.00 3 | 17.80 20.74 18.55 12.14 30.77 | 100.00 4 | 8.79 11.72 12.32 11.93 55.23 | 100.00 -----------+-------------------------------------------------------+---------- Total | 31.81 19.27 13.73 9.46 25.73 | 100.00 . corr mdu L.mdu (obs=14266) | L. | mdu mdu -------------+------------------ mdu | --. | 1.0000 L1. | 0.6184 1.0000 . . * Simple time-series plot for first 20 individuals (= first 85 obs) . quietly xtline mdu if _n<=85, overlay legend(off) . graph export countfig31panelplot.wmf, replace (file c:\Imbook\transparencies\acdcount\countfig31panelplot.wmf written in Windows Metafile format) . . *** 3. PANEL POISSON: POOLED AND POPULATION AVERAGED . . * Pooled Poisson estimator with cluster-robust standard errors . poisson mdu lcoins ndisease female age lfam child, vce(cluster id) Iteration 0: log pseudolikelihood = -62580.248 Iteration 1: log pseudolikelihood = -62579.401 Iteration 2: log pseudolikelihood = -62579.401 Poisson regression Number of obs = 20186 Wald chi2(6) = 476.93 Prob > chi2 = 0.0000 Log pseudolikelihood = -62579.401 Pseudo R2 = 0.0609 (Std. Err. adjusted for 5908 clusters in id) ------------------------------------------------------------------------------ | Robust mdu | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- lcoins | -.0808023 .0080013 -10.10 0.000 -.0964846 -.0651199 ndisease | .0339334 .0026024 13.04 0.000 .0288328 .039034 female | .1717862 .0342551 5.01 0.000 .1046473 .2389251 age | .0040585 .0016891 2.40 0.016 .000748 .0073691 lfam | -.1481981 .0323434 -4.58 0.000 -.21159 -.0848062 child | .1030453 .0506901 2.03 0.042 .0036944 .2023961 _cons | .748789 .0785738 9.53 0.000 .5947872 .9027907 ------------------------------------------------------------------------------ . . * Without cluster-robust . poisson mdu lcoins ndisease female age lfam child Iteration 0: log likelihood = -62580.248 Iteration 1: log likelihood = -62579.401 Iteration 2: log likelihood = -62579.401 Poisson regression Number of obs = 20186 LR chi2(6) = 8122.52 Prob > chi2 = 0.0000 Log likelihood = -62579.401 Pseudo R2 = 0.0609 ------------------------------------------------------------------------------ mdu | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- lcoins | -.0808023 .0020358 -39.69 0.000 -.0847924 -.0768121 ndisease | .0339334 .0005581 60.80 0.000 .0328395 .0350273 female | .1717862 .0086984 19.75 0.000 .1547376 .1888348 age | .0040585 .0004128 9.83 0.000 .0032494 .0048677 lfam | -.1481981 .0082043 -18.06 0.000 -.1642783 -.1321179 child | .1030453 .0146897 7.01 0.000 .0742539 .1318366 _cons | .748789 .0208494 35.91 0.000 .7079249 .789653 ------------------------------------------------------------------------------ . . * Poisson PA estimator with unstructured error correlation and robust VCE . xtpoisson mdu lcoins ndisease female age lfam child, pa corr(unstr) vce(robust) Iteration 1: tolerance = .01585489 Iteration 2: tolerance = .00034066 Iteration 3: tolerance = 2.334e-06 Iteration 4: tolerance = 1.939e-08 GEE population-averaged model Number of obs = 20186 Group and time vars: id year Number of groups = 5908 Link: log Obs per group: min = 1 Family: Poisson avg = 3.4 Correlation: unstructured max = 5 Wald chi2(6) = 508.61 Scale parameter: 1 Prob > chi2 = 0.0000 (Std. Err. adjusted for clustering on id) ------------------------------------------------------------------------------ | Semi-robust mdu | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- lcoins | -.0804454 .0077782 -10.34 0.000 -.0956904 -.0652004 ndisease | .0346067 .0024238 14.28 0.000 .0298561 .0393573 female | .1585075 .0334407 4.74 0.000 .0929649 .2240502 age | .0030901 .0015356 2.01 0.044 .0000803 .0060999 lfam | -.1406549 .0293672 -4.79 0.000 -.1982135 -.0830962 child | .1013677 .04301 2.36 0.018 .0170696 .1856658 _cons | .7764626 .0717221 10.83 0.000 .6358897 .9170354 ------------------------------------------------------------------------------ . . *** 3. PANEL POISSON: RANDOM EFFECTS . . * Poisson random-effects estimator with default standard errors . xtpoisson mdu lcoins ndisease female age lfam child, re Fitting Poisson model: Iteration 0: log likelihood = -62580.248 Iteration 1: log likelihood = -62579.401 Iteration 2: log likelihood = -62579.401 Fitting full model: Iteration 0: log likelihood = -43248.161 Iteration 1: log likelihood = -43240.57 Iteration 2: log likelihood = -43240.556 Iteration 3: log likelihood = -43240.556 Random-effects Poisson regression Number of obs = 20186 Group variable: id Number of groups = 5908 Random effects u_i ~ Gamma Obs per group: min = 1 avg = 3.4 max = 5 Wald chi2(6) = 637.49 Log likelihood = -43240.556 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ mdu | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- lcoins | -.0878258 .0068682 -12.79 0.000 -.1012873 -.0743642 ndisease | .0387629 .0022046 17.58 0.000 .034442 .0430839 female | .1667192 .0286298 5.82 0.000 .1106058 .2228325 age | .0019159 .0011134 1.72 0.085 -.0002663 .0040982 lfam | -.1351786 .0260022 -5.20 0.000 -.186142 -.0842152 child | .1082678 .0341477 3.17 0.002 .0413396 .1751961 _cons | .7574177 .0618346 12.25 0.000 .6362241 .8786112 -------------+---------------------------------------------------------------- /lnalpha | .0251256 .0209586 -.0159526 .0662038 -------------+---------------------------------------------------------------- alpha | 1.025444 .0214919 .984174 1.068444 ------------------------------------------------------------------------------ Likelihood-ratio test of alpha=0: chibar2(01) = 3.9e+04 Prob>=chibar2 = 0.000 . . * Poisson random-effects estimator with cluster-robust standard errors . * xtpoisson mdu lcoins ndisease female age lfam child, re vce(boot, reps(400) seed(10101) nodots) . * Previous command commented out as takes a long time . . * Poisson random-effects estimator with normal intercept and default standard errors . xtpoisson mdu lcoins ndisease female age lfam child, re normal Fitting comparison Poisson model: Iteration 0: log likelihood = -62580.248 Iteration 1: log likelihood = -62579.401 Iteration 2: log likelihood = -62579.401 Fitting full model: tau = 0.0 log likelihood = -62579.401 tau = 0.1 log likelihood = -48555.969 tau = 0.2 log likelihood = -45691.914 tau = 0.3 log likelihood = -44608.086 tau = 0.4 log likelihood = -44155.322 tau = 0.5 log likelihood = -44187.327 Iteration 0: log likelihood = -43995.824 Iteration 1: log likelihood = -43258.374 Iteration 2: log likelihood = -43227.009 Iteration 3: log likelihood = -43226.889 Iteration 4: log likelihood = -43226.889 Random-effects Poisson regression Number of obs = 20186 Group variable: id Number of groups = 5908 Random effects u_i ~ Gaussian Obs per group: min = 1 avg = 3.4 max = 5 Wald chi2(6) = 798.96 Log likelihood = -43226.889 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ mdu | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- lcoins | -.1145337 .0072788 -15.74 0.000 -.1287998 -.1002676 ndisease | .0408695 .0022941 17.81 0.000 .0363731 .045366 female | .208436 .0304848 6.84 0.000 .148687 .268185 age | .002674 .0011961 2.24 0.025 .0003297 .0050182 lfam | -.1443327 .0265365 -5.44 0.000 -.1963432 -.0923222 child | .0737146 .0344697 2.14 0.032 .0061553 .141274 _cons | .2872796 .0641625 4.48 0.000 .1615234 .4130359 -------------+---------------------------------------------------------------- /lnsig2u | .0549982 .0254991 2.16 0.031 .005021 .1049755 -------------+---------------------------------------------------------------- sigma_u | 1.027881 .013105 1.002514 1.05389 ------------------------------------------------------------------------------ Likelihood-ratio test of sigma_u=0: chibar2(01) = 3.9e+04 Pr>=chibar2 = 0.000 . . * Poisson random-effects estimator with normal intercept and normal slope for one parameter . * xtmepoisson mdu lcoins ndisease female age lfam child || id: NDISEASE . * Previous command commented out as takes a long time . . *** 3. PANEL POISSON: FIXED EFFECTS . . * Poisson fixed-effects estimator with WRONG default standard errors . xtpoisson mdu lcoins ndisease female age lfam child, fe i(id) note: 265 groups (265 obs) dropped because of only one obs per group note: 666 groups (2130 obs) dropped because of all zero outcomes note: lcoins omitted because it is constant within group note: ndisease omitted because it is constant within group note: female omitted because it is constant within group Iteration 0: log likelihood = -24182.852 Iteration 1: log likelihood = -24173.211 Iteration 2: log likelihood = -24173.211 Conditional fixed-effects Poisson regression Number of obs = 17791 Group variable: id Number of groups = 4977 Obs per group: min = 2 avg = 3.6 max = 5 Wald chi2(3) = 19.24 Log likelihood = -24173.211 Prob > chi2 = 0.0002 ------------------------------------------------------------------------------ mdu | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- age | -.0112009 .0039024 -2.87 0.004 -.0188494 -.0035523 lfam | .0877134 .0554606 1.58 0.114 -.0209874 .1964141 child | .1059867 .0437744 2.42 0.015 .0201905 .1917829 ------------------------------------------------------------------------------ . . * Poisson fixed-effects estimator with cluster-robust standard errors . * xtpoisson mdu lcoins ndisease female age lfam child, fe vce(boot, reps(400) seed(10101) nodots) . * Previous command commented out as takes a long time . . *** 3. PANEL POISSON: ESTIMATOR COMPARISON . . * Comparison of Poisson panel estimators . quietly poisson mdu lcoins ndisease female age lfam child, vce(cluster id) . estimates store POOLED . quietly xtpoisson mdu lcoins ndisease female age lfam child, pa corr(unstr) vce(robust) . estimates store POPAVE . xtpoisson mdu lcoins ndisease female age lfam child, re vce(boot, reps(100) seed(10101)) (running xtpoisson on estimation sample) Bootstrap replications (100) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .................................................. 50 .................................................. 100 Random-effects Poisson regression Number of obs = 20186 Group variable: id Number of groups = 5908 Random effects u_i ~ Gamma Obs per group: min = 1 avg = 3.4 max = 5 Wald chi2(6) = 529.10 Log likelihood = -43240.556 Prob > chi2 = 0.0000 (Replications based on 5908 clusters in id) ------------------------------------------------------------------------------ | Observed Bootstrap Normal-based mdu | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- lcoins | -.0878258 .0086097 -10.20 0.000 -.1047004 -.0709511 ndisease | .0387629 .0026904 14.41 0.000 .0334899 .0440359 female | .1667192 .0379216 4.40 0.000 .0923942 .2410442 age | .0019159 .0016242 1.18 0.238 -.0012675 .0050994 lfam | -.1351786 .0308529 -4.38 0.000 -.1956492 -.0747079 child | .1082678 .0495487 2.19 0.029 .0111541 .2053816 _cons | .7574177 .0754536 10.04 0.000 .6095314 .905304 -------------+---------------------------------------------------------------- /lnalpha | .0251256 .0270297 -.0278516 .0781029 -------------+---------------------------------------------------------------- alpha | 1.025444 .0277175 .9725326 1.081234 ------------------------------------------------------------------------------ Likelihood-ratio test of alpha=0: chibar2(01) = 3.9e+04 Prob>=chibar2 = 0.000 . * quietly xtpoisson mdu lcoins ndisease female age lfam child, re . estimates store RE_GAMMA . quietly xtpoisson mdu lcoins ndisease female age lfam child, re normal . estimates store RE_NORMAL . xtpoisson mdu lcoins ndisease female age lfam child, fe vce(boot, reps(100) seed(10101)) (running xtpoisson on estimation sample) Bootstrap replications (100) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .................................................. 50 .................................................. 100 Conditional fixed-effects Poisson regression Number of obs = 17791 Group variable: id Number of groups = 4977 Obs per group: min = 2 avg = 3.6 max = 5 Wald chi2(3) = 4.64 Log likelihood = -24173.211 Prob > chi2 = 0.2002 (Replications based on 4977 clusters in id) ------------------------------------------------------------------------------ | Observed Bootstrap Normal-based mdu | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- age | -.0112009 .0095077 -1.18 0.239 -.0298356 .0074339 lfam | .0877134 .1125783 0.78 0.436 -.132936 .3083627 child | .1059867 .0738452 1.44 0.151 -.0387472 .2507206 ------------------------------------------------------------------------------ . * quietly xtpoisson mdu lcoins ndisease female age lfam child, fe . estimates store FIXED . estimates table POOLED POPAVE RE_GAMMA RE_NORMAL FIXED, equations(1) b(%8.4f) se stats(N ll) stfmt(%8.0f) --------------------------------------------------------------------- Variable | POOLED POPAVE RE_GAMMA RE_NOR~L FIXED -------------+------------------------------------------------------- #1 | lcoins | -0.0808 -0.0804 -0.0878 -0.1145 | 0.0080 0.0078 0.0086 0.0073 ndisease | 0.0339 0.0346 0.0388 0.0409 | 0.0026 0.0024 0.0027 0.0023 female | 0.1718 0.1585 0.1667 0.2084 | 0.0343 0.0334 0.0379 0.0305 age | 0.0041 0.0031 0.0019 0.0027 -0.0112 | 0.0017 0.0015 0.0016 0.0012 0.0095 lfam | -0.1482 -0.1407 -0.1352 -0.1443 0.0877 | 0.0323 0.0294 0.0309 0.0265 0.1126 child | 0.1030 0.1014 0.1083 0.0737 0.1060 | 0.0507 0.0430 0.0495 0.0345 0.0738 _cons | 0.7488 0.7765 0.7574 0.2873 | 0.0786 0.0717 0.0755 0.0642 -------------+------------------------------------------------------- lnalpha | _cons | 0.0251 | 0.0270 -------------+------------------------------------------------------- lnsig2u | _cons | 0.0550 | 0.0255 -------------+------------------------------------------------------- Statistics | N | 20186 20186 20186 20186 17791 ll | -62579 -43241 -43227 -24173 --------------------------------------------------------------------- legend: b/se . . *** 3. PANEL NEGATIVE BINOMIAL . . /* > * Negative binomial pooled estimator with default standard errors > nbreg mdu lcoins ndisease female age lfam child > * Negative binomial pooled estimator with het robust standard errors > nbreg mdu lcoins ndisease female age lfam child, vce(robust) > * Negative binomial pooled estimator with cluster-robust standard errors > nbreg mdu lcoins ndisease female age lfam child, cluster(id) > * Negative binomial population-averaged estimator with equicorrelated errors > xtnbreg lcoins ndisease female age lfam child, pa corr(exch) vce(robust) > * Negative binomial random-effects estimator with default standard errors > xtnbreg mdu ndisease female age lfam child, re i(id) > * Negative binomial fixed-effects estimator with default standard errors > xtnbreg mdu ndisease female age lfam child, fe i(id) > > * Comparison of negative binomial panel estimators > quietly xtpoisson mdu lcoins ndisease female age lfam child, pa corr(exch) vce(robust) > estimates store PPA_ROB > quietly xtnbreg mdu lcoins ndisease female age lfam child, pa corr(exch) vce(robust) > estimates store NBPA_ROB > quietly xtnbreg mdu lcoins ndisease female age lfam child, re > estimates store NBRE > quietly xtnbreg mdu lcoins ndisease female age lfam child, fe > estimates store NBFE > estimates table PPA_ROB NBPA_ROB NBRE NBFE, equations(1) b(%8.4f) se stats(N ll) stfmt(%8.0f) > */ . . *********** 4: MULTIVARIATE, MAXIMUM SIMULATED LIKELIHOOD, BAYESIAN . . * NO DATA EXAMPLES . . ********** CLOSE OUTPUT . . * log close . * clear . * exit . . end of do-file . exit, clear