------------------------------------------------------------------------------------------------------------------------------- name: log: c:\acdbookrevision\stata_final_programs_2013\racd08.txt log type: text opened on: 19 Jan 2013, 11:50:03 . . ********** OVERVIEW OF racd08.do ********** . . * STATA Program . * copyright C 2013 by A. Colin Cameron and Pravin K. Trivedi . * used for "Regression Analyis of Count Data" SECOND EDITION . * by A. Colin Cameron and Pravin K. Trivedi (2013) . * Cambridge University Press . . * This STATA program does examples for chapter 8.8 and 8.9 . * 8.5 COPULAS: CLAYTON AND GUMBEL . * 8.9 EXAMPLE: BIVARIATE COUNT ANALYSIS . * (1) INDEPENDENCE TESTS BASED ON ORTHOGONAL POLYNOMIALS . * (2) NLSUR: NUNLINEAR SEEMINGLY UNRELATED REGRESSION ESTIMATOR . * (3) MULTIVARIATE NEGATIVE BINOMIAL ESTIMATION ESTIMATED BY ML . . * To run you need file . * racd06data1healthcare.dta . * in your directory . . ********** SETUP ********** . . set more off . version 12 . clear all . * set linesize 82 . set scheme s1mono /* Graphics scheme */ . . ********** DATA DESCRIPTION for CHAPTER 8.9 . . * The data are extracted from the 1987-88 National Medical Expenditure Survey (NMES). . * The extract and analysis are in P. Deb and P.K. Trivedi (1997), . * Demand for Medical Care by the Elderly: A Finite Mixture Approach" . * Journal of Applied Econometrics, 12, 313-326. . * See this article for more detailed discussion . * Also see racd06makedata1healthcare.do for further details . . * This STATA program does the analysis for chapter 9 . * 8.5 COPULA . * 8.9 EMPIRICAL EXAMPLE (EMR and HOSP) . . * To run you need file . * racd06data1healthcare.dta . * in your directory . . ********** DATA DESCRIPTION . . * The data are extracted from the 1987-88 National Medical Expenditure Survey (NMES). . * The extract and analysis are in P. Deb and P.K. Trivedi (1997), . * Demand for Medical Care by the Elderly: A Finite Mixture Approach" . * Journal of Applied Econometrics, 12, 313-326. . * See this article for more detailed discussion . . * The simulation to show copula generated data . * is based on code from P.K. Trivedi and D.M. Zimmer (2005) . * Copula Modeling: An Introduction for Practitioners . * Foundations and Trends in Econometrics Vol. 1, No 1 1-111. . . ********** 8.5 COPULAS: CLAYTON AND GUMBEL . . * The data are generated data . * y1 ~ Poisson(10) . * y2 ~ Poisson(10) . * Copula is . * Clayton theta = 2 so Kendall's tau = 2/(2+2) = 0.5 . * Gumbel theta = 2 so Kendall's tau = (2-1)/2 = 0.5 . . *** CLAYTON COPULA . . * Code from David Zimmer . clear . set obs 1000 obs was 0, now 1000 . set seed 10101 . gen y1=0 . gen y2=0 . mata: ------------------------------------------------- mata (type end to exit) ----------------------------------------------------- : obs = 1000 : mean1 = 10 : mean2 = 10 : theta = 2 : y1 = J(obs,1,-999) : y2 = J(obs,1,-999) : v1 = uniform(obs,1) : v2 = uniform(obs,1) : u1 = v1 : u2 = ( (v1:^(-theta)) :* (v2:^(-theta/(theta+1)) :- 1) :+ 1 ) :^(-1/theta) : for(j=1; j<=obs; j++) { > p10 = exp(-mean1) > p20 = exp(-mean2) > s1 = p10 > s2 = p20 > for(i=0; i<50; i++) { > if (u1[j,] < s1) y1[j,]=i > if (u1[j,] < s1) i=50 > if (u1[j,] > s1) p10 = mean1*p10/(i+1) > if (u1[j,] > s1) s1 = s1 + p10 > if (u1[j,] > s1) y1[j,]=i+1 > } > for(i=0; i<50; i++) { > if (u2[j,] < s2) y2[j,]=i > if (u2[j,] < s2) i=50 > if (u2[j,] > s2) p20 = mean2*p20/(i+1) > if (u2[j,] > s2) s2 = s2 + p20 > if (u2[j,] > s2) y2[j,]=i+1 > } > } : st_store(., "y1", y1) : st_store(., "y2", y2) : end ------------------------------------------------------------------------------------------------------------------------------- . summarize Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- y1 | 1000 10.15 3.279368 1 22 y2 | 1000 10.174 3.222438 1 23 . ktau y1 y2 Number of obs = 1000 Kendall's tau-a = 0.5114 Kendall's tau-b = 0.5601 Kendall's score = 255464 SE of score = 10453.860 (corrected for ties) Test of Ho: y1 and y2 are independent Prob > |z| = 0.0000 (continuity corrected) . graph twoway (scatter y1 y2) (lfit y1 y2, lwidth(medthick)), legend(off) /// > title("Sample from Clayton Copula") ytitle("y1") scale(1.2) saving(clayton, replace) (file clayton.gph saved) . . *** GUMBEL COPULA . . * Code from David Zimmer . clear . set obs 1000 obs was 0, now 1000 . set seed 10101 . gen y1=0 . gen y2=0 . mata: ------------------------------------------------- mata (type end to exit) ----------------------------------------------------- : obs = 1000 : mean1 = 10 : mean2 = 10 : theta = 4 : y1 = J(obs,1,-999) : y2 = J(obs,1,-999) : thet = uniform(2000,1) :* 3.1415 : w = uniform(2000,1) : ww = -1:*ln(w) : alph = 1/theta : z1 = sin((1-alph):*thet) :* (sin(alph:*thet)):^(alph/(1-alph)) : z2 = (sin(thet)):^(1/(1-alph)) : z = z1:/z2 : xx = (z:/ww):^((1-alph)/alph) : v1 = uniform(2000,1) : v2 = uniform(2000,1) : u1 = exp( -1:*((-1:*ln(v1):/xx):^(1/theta)) ) : u2 = exp( -1:*((-1:*ln(v2):/xx):^(1/theta)) ) : for(j=1; j<=obs; j++) { > p10 = exp(-mean1) > p20 = exp(-mean2) > s1 = p10 > s2 = p20 > for(i=0; i<50; i++) { > if (u1[j,] < s1) y1[j,]=i > if (u1[j,] < s1) i=50 > if (u1[j,] > s1) p10 = mean1*p10/(i+1) > if (u1[j,] > s1) s1 = s1 + p10 > if (u1[j,] > s1) y1[j,]=i+1 > } > for(i=0; i<50; i++) { > if (u2[j,] < s2) y2[j,]=i > if (u2[j,] < s2) i=50 > if (u2[j,] > s2) p20 = mean2*p20/(i+1) > if (u2[j,] > s2) s2 = s2 + p20 > if (u2[j,] > s2) y2[j,]=i+1 > } > } : st_store(., "y1", y1) : st_store(., "y2", y2) : end ------------------------------------------------------------------------------------------------------------------------------- . sum Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- y1 | 1000 10.038 3.158882 2 21 y2 | 1000 10.005 3.150063 2 22 . ktau y1 y2 Number of obs = 1000 Kendall's tau-a = 0.7191 Kendall's tau-b = 0.7897 Kendall's score = 359182 SE of score = 10448.973 (corrected for ties) Test of Ho: y1 and y2 are independent Prob > |z| = 0.0000 (continuity corrected) . graph twoway (scatter y1 y2) (lfit y1 y2, lwidth(medthick)), legend(off) /// > title("Sample from Gumbel Copula") ytitle("y1") scale(1.2) saving(gumbel, replace) (file gumbel.gph saved) . . *** FIGURE 8.1: CLAYTON AND GUMBEL COPULAS EXAMPLE . . graph combine clayton.gph gumbel.gph, ycommon xcommon ysize(3) xsize(6) . graph export racd08fig1.wmf, replace (file c:\acdbookrevision\stata_final_programs_2013\racd08fig1.wmf written in Windows Metafile format) . graph export racd08fig1.eps, replace (file racd08fig1.eps written in EPS format) . . ********** 8.9 EMPIRICAL EXAMPLE . . ****** READ IN DATAA AND SUMMARIZE . . use racd06data1healthcare.dta, clear . . * Variable descriptions and summary statistics . describe Contains data from racd06data1healthcare.dta obs: 4,406 vars: 22 7 Jun 2011 10:39 size: 387,728 ------------------------------------------------------------------------------------------------------------------------------- storage display value variable name type format label variable label ------------------------------------------------------------------------------------------------------------------------------- OFP float %9.0g Number of physician office visits OFNP float %9.0g Number of non-physician office visits OPP float %9.0g Number of physician outpatient visits OPNP float %9.0g Number of non-physician outpatient visits EMR float %9.0g Number of emergency room visits HOSP float %9.0g Number hospitalizations EXCLHLTH float %9.0g Equals 1 if self perceived health is excellent POORHLTH float %9.0g Equals 1 if self perceived health is poor NUMCHRON float %9.0g Number of chronic conditions ADLDIFF float %9.0g Equals 1 if the person has a condition that limits activities of daily living NOREAST float %9.0g Equals 1 if the person lives in northeastern U.S. MIDWEST float %9.0g Equals 1 if the person lives in the midwestern U.S. WEST float %9.0g Equals 1 if the person lives in the western U.S. AGE float %9.0g Age in years (divided by 10) BLACK float %9.0g Equals 1 if the person is African American MALE float %9.0g Equals 1 if the person is male MARRIED float %9.0g Equals 1 if the person is married SCHOOL float %9.0g Number of years of education FAMINC float %9.0g Family income in $10,000 EMPLOYED float %9.0g Equals 1 if the person is employed PRIVINS float %9.0g Equals 1 if the person is covered by private health insurance MEDICAID float %9.0g Equals 1 if the person is covered by Medicaid ------------------------------------------------------------------------------------------------------------------------------- Sorted by: . summarize Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- OFP | 4406 5.774399 6.759225 0 89 OFNP | 4406 1.618021 5.317056 0 104 OPP | 4406 .7507944 3.652759 0 141 OPNP | 4406 .5360872 3.879506 0 155 EMR | 4406 .2635043 .7036586 0 12 -------------+-------------------------------------------------------- HOSP | 4406 .2959601 .7463978 0 8 EXCLHLTH | 4406 .0778484 .2679633 0 1 POORHLTH | 4406 .1257376 .3315911 0 1 NUMCHRON | 4406 1.541988 1.349632 0 8 ADLDIFF | 4406 .2040399 .4030441 0 1 -------------+-------------------------------------------------------- NOREAST | 4406 .1899682 .3923203 0 1 MIDWEST | 4406 .2625965 .4400949 0 1 WEST | 4406 .1811167 .3851585 0 1 AGE | 4406 7.402406 .6334051 6.6 10.9 BLACK | 4406 .117113 .3215914 0 1 -------------+-------------------------------------------------------- MALE | 4406 .4035406 .4906631 0 1 MARRIED | 4406 .5460735 .4979292 0 1 SCHOOL | 4406 10.29029 3.738736 0 18 FAMINC | 4406 2.527132 2.924648 -1.0125 54.8351 EMPLOYED | 4406 .1032683 .3043435 0 1 -------------+-------------------------------------------------------- PRIVINS | 4406 .7764412 .4166769 0 1 MEDICAID | 4406 .0912392 .2879817 0 1 . . summarize EMR HOSP, detail Number of emergency room visits ------------------------------------------------------------- Percentiles Smallest 1% 0 0 5% 0 0 10% 0 0 Obs 4406 25% 0 0 Sum of Wgt. 4406 50% 0 Mean .2635043 Largest Std. Dev. .7036586 75% 0 8 90% 1 8 Variance .4951354 95% 1 11 Skewness 5.066106 99% 3 12 Kurtosis 49.81375 Number hospitalizations ------------------------------------------------------------- Percentiles Smallest 1% 0 0 5% 0 0 10% 0 0 Obs 4406 25% 0 0 Sum of Wgt. 4406 50% 0 Mean .2959601 Largest Std. Dev. .7463978 75% 0 8 90% 1 8 Variance .5571097 95% 2 8 Skewness 3.964427 99% 3 8 Kurtosis 25.80247 . summarize EMR HOSP Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- EMR | 4406 .2635043 .7036586 0 12 HOSP | 4406 .2959601 .7463978 0 8 . correlate EMR HOSP (obs=4406) | EMR HOSP -------------+------------------ EMR | 1.0000 HOSP | 0.4761 1.0000 . correlate EMR HOSP if EMR > 0 & HOSP > 0 (obs=454) | EMR HOSP -------------+------------------ EMR | 1.0000 HOSP | 0.3987 1.0000 . tabulate EMR Number of | emergency | room visits | Freq. Percent Cum. ------------+----------------------------------- 0 | 3,602 81.75 81.75 1 | 588 13.35 95.10 2 | 137 3.11 98.21 3 | 54 1.23 99.43 4 | 11 0.25 99.68 5 | 7 0.16 99.84 6 | 2 0.05 99.89 7 | 1 0.02 99.91 8 | 2 0.05 99.95 11 | 1 0.02 99.98 12 | 1 0.02 100.00 ------------+----------------------------------- Total | 4,406 100.00 . tabulate HOSP Number | hospitaliza | tions | Freq. Percent Cum. ------------+----------------------------------- 0 | 3,541 80.37 80.37 1 | 599 13.60 93.96 2 | 176 3.99 97.96 3 | 48 1.09 99.05 4 | 20 0.45 99.50 5 | 12 0.27 99.77 6 | 5 0.11 99.89 7 | 1 0.02 99.91 8 | 4 0.09 100.00 ------------+----------------------------------- Total | 4,406 100.00 . . ****** INDEPENDENCE TESTS BASED ON ORTHOGONAL POLYNOMIALS . . global XLIST EXCLHLTH POORHLTH NUMCHRON ADLDIFF NOREAST MIDWEST WEST AGE /// > BLACK MALE MARRIED SCHOOL FAMINC EMPLOYED PRIVINS MEDICAID . global y1 EMR . global y2 HOSP . . * Following is for NB2 model . * Generate first two orthogonal polynomials for y1 . nbreg $y1 $XLIST Fitting Poisson model: Iteration 0: log likelihood = -2811.3847 Iteration 1: log likelihood = -2810.6856 Iteration 2: log likelihood = -2810.6853 Iteration 3: log likelihood = -2810.6853 Fitting constant-only model: Iteration 0: log likelihood = -2850.4693 Iteration 1: log likelihood = -2811.7131 Iteration 2: log likelihood = -2802.979 Iteration 3: log likelihood = -2802.9781 Iteration 4: log likelihood = -2802.9781 Fitting full model: Iteration 0: log likelihood = -2685.7858 Iteration 1: log likelihood = -2662.9502 Iteration 2: log likelihood = -2662.2099 Iteration 3: log likelihood = -2662.2093 Iteration 4: log likelihood = -2662.2093 Negative binomial regression Number of obs = 4406 LR chi2(16) = 281.54 Dispersion = mean Prob > chi2 = 0.0000 Log likelihood = -2662.2093 Pseudo R2 = 0.0502 ------------------------------------------------------------------------------ EMR | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- EXCLHLTH | -.6370736 .1957708 -3.25 0.001 -1.020777 -.2533699 POORHLTH | .4939609 .102036 4.84 0.000 .293974 .6939478 NUMCHRON | .2212737 .0267322 8.28 0.000 .1688795 .2736678 ADLDIFF | .4036297 .0925641 4.36 0.000 .2222074 .5850521 NOREAST | .046639 .10587 0.44 0.660 -.1608624 .2541403 MIDWEST | .0210809 .0975606 0.22 0.829 -.1701343 .2122961 WEST | .1770833 .1070172 1.65 0.098 -.0326667 .3868332 AGE | .0724237 .0606626 1.19 0.233 -.0464728 .1913203 BLACK | .174874 .1163038 1.50 0.133 -.0530773 .4028253 MALE | .068292 .0841144 0.81 0.417 -.0965693 .2331532 MARRIED | -.1076015 .0876875 -1.23 0.220 -.2794658 .0642628 SCHOOL | -.017987 .0110336 -1.63 0.103 -.0396125 .0036384 FAMINC | .0038962 .0135328 0.29 0.773 -.0226275 .03042 EMPLOYED | .1722555 .1319031 1.31 0.192 -.0862699 .4307809 PRIVINS | .059815 .1052137 0.57 0.570 -.1464 .26603 MEDICAID | .1810393 .135732 1.33 0.182 -.0849906 .4470692 _cons | -2.430219 .4916736 -4.94 0.000 -3.393881 -1.466556 -------------+---------------------------------------------------------------- /lnalpha | .4992164 .1016006 .3000828 .6983499 -------------+---------------------------------------------------------------- alpha | 1.64743 .1673799 1.349971 2.010433 ------------------------------------------------------------------------------ Likelihood-ratio test of alpha=0: chibar2(01) = 296.95 Prob>=chibar2 = 0.000 . predict mu1, n . scalar alpha1 = exp([lnalpha]_cons) . generate Q1y1 = $y1 - mu1 . generate Q2y1 = ($y1-mu1)^2 - (1+2*alpha1*mu1)*($y1-mu1) - (1+alpha1*mu1)*mu1 . * Generate first two orthogonal polynomials for y2 . nbreg $y2 $XLIST Fitting Poisson model: Iteration 0: log likelihood = -3028.5825 Iteration 1: log likelihood = -3027.7741 Iteration 2: log likelihood = -3027.7737 Iteration 3: log likelihood = -3027.7737 Fitting constant-only model: Iteration 0: log likelihood = -3067.9878 Iteration 1: log likelihood = -3020.9814 Iteration 2: log likelihood = -3009.627 Iteration 3: log likelihood = -3009.6246 Iteration 4: log likelihood = -3009.6246 Fitting full model: Iteration 0: log likelihood = -2874.2435 Iteration 1: log likelihood = -2847.0582 Iteration 2: log likelihood = -2846.414 Iteration 3: log likelihood = -2846.4137 Iteration 4: log likelihood = -2846.4137 Negative binomial regression Number of obs = 4406 LR chi2(16) = 326.42 Dispersion = mean Prob > chi2 = 0.0000 Log likelihood = -2846.4137 Pseudo R2 = 0.0542 ------------------------------------------------------------------------------ HOSP | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- EXCLHLTH | -.6894085 .1939777 -3.55 0.000 -1.069598 -.3092192 POORHLTH | .511255 .0992636 5.15 0.000 .3167019 .7058081 NUMCHRON | .2764486 .0267385 10.34 0.000 .224042 .3288552 ADLDIFF | .3377104 .0909624 3.71 0.000 .1594273 .5159935 NOREAST | -.0002047 .1040003 -0.00 0.998 -.2040415 .2036321 MIDWEST | .1344601 .0931499 1.44 0.149 -.0481103 .3170305 WEST | .1348615 .1045674 1.29 0.197 -.0700868 .3398098 AGE | .1720692 .0595326 2.89 0.004 .0553874 .288751 BLACK | .1028964 .1187966 0.87 0.386 -.1299406 .3357334 MALE | .2131734 .081686 2.61 0.009 .0530718 .373275 MARRIED | -.0342766 .085295 -0.40 0.688 -.2014518 .1328986 SCHOOL | .0011724 .0106513 0.11 0.912 -.0197038 .0220486 FAMINC | .0007634 .0133114 0.06 0.954 -.0253265 .0268533 EMPLOYED | .0435324 .1306011 0.33 0.739 -.212441 .2995058 PRIVINS | .1842365 .1052341 1.75 0.080 -.0220185 .3904916 MEDICAID | .1411514 .1388956 1.02 0.310 -.1310789 .4133817 _cons | -3.511593 .485745 -7.23 0.000 -4.463636 -2.55955 -------------+---------------------------------------------------------------- /lnalpha | .537769 .0920903 .3572754 .7182626 -------------+---------------------------------------------------------------- alpha | 1.712183 .1576754 1.42943 2.050867 ------------------------------------------------------------------------------ Likelihood-ratio test of alpha=0: chibar2(01) = 362.72 Prob>=chibar2 = 0.000 . predict mu2, n . scalar alpha2 = exp([lnalpha]_cons) . generate Q1y2 = $y2 - mu2 . generate Q2y2 = ($y2-mu2)^2 - (1+2*alpha2*mu2)*($y2-mu2) - (1+alpha2*mu2)*mu2 . . /* > * Following is for NB1 model - not reported in book > * Note that here we use Var = alpha*mu whereas Stata sets Var = mu + delta*mu > * Generate first two orthogonal polynomials for y1 > nbreg $y1 $XLIST, dispersion(constant) > predict mu1, n > scalar delta1 = exp([lndelta]_cons) > scalar phi1 = phi1 + 1 > generate Q1y1 = $y1 - mu1 > generate Q2y1 = ($y1-mu1)^2 - (2*phi1)*($y1-mu1) - phi1*mu1 > * Generate first two orthogonal polynomials for y2 > nbreg $y2 $XLIST, dispersion(constant) > predict mu2, n > scalar delta2 = exp([lndelta]_cons) > scalar phi2 = delta2 + 1 > generate Q1y2 = $y2 - mu2 > generate Q2y2 = ($y2-mu2)^2 - (2*phi2)*($y2-mu1) - phi2*mu1 > */ . . /* > * Following is for Poisson model - not reported in book > * Generate first two orthogonal polynomials for y1 > poisson $y1 $XLIST > predict mu1, n > generate Q1y1 = $y1 - mu1 > generate Q2y1 = ($y1-mu1)^2 - $y1 > * Generate first two orthogonal polynomials for y2 > poisson $y2 $XLIST > predict mu2, n > generate Q1y2 = $y2 - mu2 > generate Q2y2 = ($y2-mu2)^2 - $y2 > */ . . * Now perform tests based on crossproducts of Q1 and Q2 . generate one = 1 . generate Q1y1Q1y2 = Q1y1*Q1y2 . quietly regress one Q1y1Q1y2, noconstant . scalar test11 = e(N)*e(r2) . generate Q2y1Q2y2 = Q2y1*Q2y2 . quietly regress one Q2y1Q2y2, noconstant . scalar test22 = e(N)*e(r2) . generate Q1y1Q2y2 = Q1y1*Q2y2 . quietly regress one Q1y1Q2y2, noconstant . scalar test12 = e(N)*e(r2) . generate Q2y1Q1y2 = Q2y1*Q1y2 . quietly regress one Q2y1Q1y2, noconstant . scalar test21 = e(N)*e(r2) . . *** RESULTS GIVEN IN TEXT: Display the four test statisics . . display "Test based on Q1y1 x Q1y2 = " test11 " and p = " chi2tail(1,test11) Test based on Q1y1 x Q1y2 = 88.44908 and p = 5.216e-21 . display "Test based on Q1y1 x Q2y2 = " test12 " and p = " chi2tail(1,test12) Test based on Q1y1 x Q2y2 = 18.049421 and p = .00002152 . display "Test based on Q2y1 x Q1y2 = " test21 " and p = " chi2tail(1,test21) Test based on Q2y1 x Q1y2 = .01426886 and p = .904917 . display "Test based on Q2y1 x Q2y2 = " test22 " and p = " chi2tail(1,test22) Test based on Q2y1 x Q2y2 = 3.9304097 and p = .04742039 . . sum Q* Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- Q1y1 | 4406 .0002089 .675437 -1.737044 11.23508 Q2y1 | 4406 .0141149 1.862419 -18.2109 84.94746 Q1y2 | 4406 -.0041016 .7178018 -2.185327 7.487564 Q2y2 | 4406 -.0024391 1.800491 -31.8428 34.47506 Q1y1Q1y2 | 4406 .2096681 1.465054 -4.303709 61.43302 -------------+-------------------------------------------------------- Q2y1Q2y2 | 4406 .4279603 14.32393 -615.7653 395.9419 Q1y1Q2y2 | 4406 -.227011 3.539936 -126.2112 57.97785 Q2y1Q1y2 | 4406 -.0104336 5.79841 -34.7777 299.7224 . correlate Q* (obs=4406) | Q1y1 Q2y1 Q1y2 Q2y2 Q1y1Q1y2 Q2y1Q2y2 Q1y1Q2y2 Q2y1Q1y2 -------------+------------------------------------------------------------------------ Q1y1 | 1.0000 Q2y1 | 0.1316 1.0000 Q1y2 | 0.4326 -0.0078 1.0000 Q2y2 | -0.1867 0.1277 -0.0723 1.0000 Q1y1Q1y2 | 0.5449 0.5107 0.4189 0.0133 1.0000 Q2y1Q2y2 | -0.1777 -0.6290 -0.0334 -0.0608 -0.4003 1.0000 Q1y1Q2y2 | -0.2877 -0.2603 0.0132 0.3471 -0.3388 0.2951 1.0000 Q2y1Q1y2 | 0.3565 0.7573 0.1258 -0.0331 0.6904 -0.7480 -0.4050 1.0000 . . /* To apply to original Cameron and Trivedi (1993) example use > use racd03data.dta, clear > global XLIST SEX AGE AGESQ INCOME LEVYPLUS FREEPOOR FREEREPA ILLNESS /// > ACTDAYS HSCORE CHCOND1 CHCOND2 > global y1 HOSPADMI > global y2 HOSPDAYS > */ . . ****** NLSUR: NUNLINEAR SEEMINGLY UNRELATED REGRESSION ESTIMATOR . . * Global for the regressors . global XLIST EXCLHLTH POORHLTH NUMCHRON ADLDIFF NOREAST MIDWEST WEST AGE /// > BLACK MALE MARRIED SCHOOL FAMINC EMPLOYED PRIVINS MEDICAID . . generate CONSTANT = 1 . . *** TABLE 8.3: NLSUR RESULTS (second half of table) . . nlsur (EMR = exp({xb1: $XLIST CONSTANT})) /// > (HOSP = exp({xb2: $XLIST CONSTANT})), vce(robust) nolog (obs = 4406) Calculating NLS estimates... Calculating FGNLS estimates... FGNLS regression --------------------------------------------------------------------- Equation | Obs Parms RMSE R-sq Constant ----------------+---------------------------------------------------- 1 EMR | 4406 17 .672962 0.1977* (none) 2 HOSP | 4406 17 .7131008 0.2111* (none) --------------------------------------------------------------------- * Uncentered R-sq ------------------------------------------------------------------------------ | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- /xb1_EXCLH~H | -.426552 .2535779 -1.68 0.093 -.9235557 .0704516 /xb1_POORH~H | .6386311 .1107226 5.77 0.000 .4216189 .8556434 /xb1_NUMCH~N | .2670969 .043998 6.07 0.000 .1808623 .3533314 /xb1_ADLDIFF | .4295973 .1657307 2.59 0.010 .1047711 .7544235 /xb1_NOREAST | -.1194534 .1735414 -0.69 0.491 -.4595883 .2206814 /xb1_MIDWEST | -.0339255 .1707602 -0.20 0.843 -.3686095 .3007584 /xb1_WEST | -.0778962 .1833873 -0.42 0.671 -.4373286 .2815362 /xb1_AGE | -.1484804 .0896563 -1.66 0.098 -.3242035 .0272428 /xb1_BLACK | .2788804 .2128257 1.31 0.190 -.1382503 .6960111 /xb1_MALE | .0787842 .1583821 0.50 0.619 -.2316389 .3892073 /xb1_MARRIED | -.1596034 .1510262 -1.06 0.291 -.4556093 .1364024 /xb1_SCHOOL | -.0012683 .0191043 -0.07 0.947 -.038712 .0361754 /xb1_FAMINC | .0085027 .0213362 0.40 0.690 -.0333155 .050321 /xb1_EMPLO~D | .3199558 .204672 1.56 0.118 -.0811939 .7211056 /xb1_PRIVINS | .0019576 .154603 0.01 0.990 -.3010587 .3049738 /xb1_MEDIC~D | .0508229 .2049068 0.25 0.804 -.3507871 .452433 /xb1_CONST~T | -1.026209 .6258008 -1.64 0.101 -2.252756 .2003385 /xb2_EXCLH~H | -.6465456 .2060276 -3.14 0.002 -1.050352 -.2427389 /xb2_POORH~H | .6365103 .0950259 6.70 0.000 .4502629 .8227577 /xb2_NUMCH~N | .2416503 .0300529 8.04 0.000 .1827478 .3005529 /xb2_ADLDIFF | .4033449 .1084431 3.72 0.000 .1908003 .6158894 /xb2_NOREAST | -.0715548 .1333755 -0.54 0.592 -.332966 .1898564 /xb2_MIDWEST | .0823602 .1401352 0.59 0.557 -.1922998 .3570201 /xb2_WEST | -.0543309 .1331784 -0.41 0.683 -.3153558 .206694 /xb2_AGE | -.0568123 .0759019 -0.75 0.454 -.2055773 .0919527 /xb2_BLACK | .1295524 .1374455 0.94 0.346 -.1398358 .3989405 /xb2_MALE | .0475138 .1204375 0.39 0.693 -.1885393 .2835669 /xb2_MARRIED | .0245081 .1082608 0.23 0.821 -.1876792 .2366953 /xb2_SCHOOL | .0087595 .016803 0.52 0.602 -.0241739 .0416928 /xb2_FAMINC | .0147487 .0125029 1.18 0.238 -.0097565 .039254 /xb2_EMPLO~D | .0991653 .1527081 0.65 0.516 -.2001372 .3984678 /xb2_PRIVINS | .3213019 .1253396 2.56 0.010 .0756408 .5669629 /xb2_MEDIC~D | .2585363 .1421873 1.82 0.069 -.0201457 .5372184 /xb2_CONST~T | -1.951981 .6011933 -3.25 0.001 -3.130298 -.7736634 ------------------------------------------------------------------------------ . estimates store NLSUR . . * Save the parameter estimates to use as starting values later for Bivariate Ne . matrix bnlsur = e(b) . . * Calculate the error correlation in two ways . matrix Sigma = e(Sigma) . matrix list Sigma symmetric Sigma[2,2] EMR HOSP EMR .45284894 HOSP .20755806 .50839241 . scalar rho = Sigma[1,2] / sqrt(Sigma[1,1]*Sigma[2,2]) . display rho .43257697 . predict u1hat, equation(#1) residuals . predict u2hat, equation(#2) residuals . correlate u1hat u2hat (obs=4406) | u1hat u2hat -------------+------------------ u1hat | 1.0000 u2hat | 0.4331 1.0000 . . ****** MULTIVARIATE NEGATIVE BINOMIAL ESTIMATION ESTIMATED BY ML . . * Bivariate Negbin ML program lfnbbivariate to be called by command ml method lf . program lfnbbivariate 1. version 10.1 2. args lnf theta1 theta2 a // theta1=x1'b1, theta2=x2'b2 a=alpha, lnf=lnf(y) 3. tempvar mu1 mu2 4. local y1 $ML_y1 // Define y1 so program more readable 5. local y2 $ML_y2 // Define y2 so program more readable 6. generate double `mu1' = exp(`theta1') 7. generate double `mu2' = exp(`theta2') 8. quietly replace `lnf' = lngamma(`y1'+`y2'+(1/`a')) - lngamma((1/`a')) /// > - lnfactorial(`y1') - lnfactorial(`y2') /// > + `y1'*ln(`mu1') + `y2'*ln(`mu2') /// > - (`y1'+`y2'+(1/`a'))*ln(1+`mu1'+`mu2') 9. end . . * Initial values - for betas use initial values from preceding NLSUR . * For alpha use alpha = 2 . ml model lf lfnbbivariate (EMR = $XLIST) (HOSP = $XLIST) (), vce(robust) . ml init bnlsur 2.0, copy . . *** TABLE 8.3: ML RESULTS . . ml maximize initial: log pseudolikelihood = -5545.08 rescale: log pseudolikelihood = -5545.08 rescale eq: log pseudolikelihood = -5294.3623 Iteration 0: log pseudolikelihood = -5294.3623 Iteration 1: log pseudolikelihood = -5225.7803 Iteration 2: log pseudolikelihood = -5222.9883 Iteration 3: log pseudolikelihood = -5222.9437 Iteration 4: log pseudolikelihood = -5222.9437 Number of obs = 4406 Wald chi2(16) = 263.35 Log pseudolikelihood = -5222.9437 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- eq1 | EXCLHLTH | -.618348 .2127487 -2.91 0.004 -1.035328 -.2013683 POORHLTH | .488027 .0973097 5.02 0.000 .2973034 .6787505 NUMCHRON | .2405193 .0281647 8.54 0.000 .1853176 .2957211 ADLDIFF | .404453 .0916514 4.41 0.000 .2248196 .5840865 NOREAST | .0435052 .1048306 0.42 0.678 -.1619589 .2489693 MIDWEST | .0276668 .0975625 0.28 0.777 -.1635522 .2188857 WEST | .1817455 .1146418 1.59 0.113 -.0429483 .4064394 AGE | .0988923 .0613506 1.61 0.107 -.0213526 .2191372 BLACK | .1737047 .1171529 1.48 0.138 -.0559107 .4033201 MALE | .1106915 .0876456 1.26 0.207 -.0610908 .2824738 MARRIED | -.1115078 .0900464 -1.24 0.216 -.2879954 .0649799 SCHOOL | -.0174469 .0111293 -1.57 0.117 -.0392598 .0043661 FAMINC | -.0008557 .0145005 -0.06 0.953 -.0292761 .0275647 EMPLOYED | .185813 .1260944 1.47 0.141 -.0613275 .4329535 PRIVINS | .0294672 .1018669 0.29 0.772 -.1701883 .2291227 MEDICAID | .1456729 .1346714 1.08 0.279 -.1182782 .409624 _cons | -1.875552 .5066264 -3.70 0.000 -2.868522 -.8825828 -------------+---------------------------------------------------------------- eq2 | EXCLHLTH | -.7012096 .1875118 -3.74 0.000 -1.068726 -.3336933 POORHLTH | .5005783 .0968527 5.17 0.000 .3107506 .690406 NUMCHRON | .2716642 .0242844 11.19 0.000 .2240676 .3192607 ADLDIFF | .3366176 .0921516 3.65 0.000 .1560038 .5172314 NOREAST | .0242009 .1037591 0.23 0.816 -.1791631 .227565 MIDWEST | .1585836 .0974262 1.63 0.104 -.0323683 .3495355 WEST | .1612162 .1075588 1.50 0.134 -.0495952 .3720277 AGE | .1814116 .0587687 3.09 0.002 .066227 .2965962 BLACK | .095689 .1136181 0.84 0.400 -.1269983 .3183764 MALE | .2062592 .0819623 2.52 0.012 .0456161 .3669023 MARRIED | -.0145049 .0880294 -0.16 0.869 -.1870393 .1580296 SCHOOL | .0015326 .0112992 0.14 0.892 -.0206135 .0236786 FAMINC | -.001753 .0112715 -0.16 0.876 -.0238448 .0203387 EMPLOYED | .0508569 .1315916 0.39 0.699 -.2070579 .3087718 PRIVINS | .1885835 .1099582 1.72 0.086 -.0269307 .4040976 MEDICAID | .1821086 .1361372 1.34 0.181 -.0847153 .4489326 _cons | -2.824792 .4910653 -5.75 0.000 -3.787262 -1.862321 -------------+---------------------------------------------------------------- eq3 | _cons | 2.165871 .1316891 16.45 0.000 1.907765 2.423977 ------------------------------------------------------------------------------ . estimates store MLBVNB . . * Aside: Confidence interval for 1/alpha . nlcom 1/[eq3]_cons _nl_1: 1/[eq3]_cons ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- _nl_1 | .4617079 .0280727 16.45 0.000 .4066864 .5167295 ------------------------------------------------------------------------------ . . *** ALTERNATIVE SINGLE EQUATION ESTIMATORS (Not given in book) . . nl (EMR = exp({xb3: $XLIST CONSTANT})), vce(robust) nolog (obs = 4406) Nonlinear regression Number of obs = 4406 R-squared = 0.1977 Adj R-squared = 0.1946 Root MSE = .6742425 Res. dev. = 9013.268 ------------------------------------------------------------------------------ | Robust EMR | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- /xb3_EXCLH~H | -.4506454 .2488543 -1.81 0.070 -.9385255 .0372346 /xb3_POORH~H | .6232095 .1104892 5.64 0.000 .4065949 .8398242 /xb3_NUMCH~N | .2626472 .0438513 5.99 0.000 .1766766 .3486178 /xb3_ADLDIFF | .4321203 .1604543 2.69 0.007 .1175488 .7466917 /xb3_NOREAST | -.1044931 .1714978 -0.61 0.542 -.4407152 .2317291 /xb3_MIDWEST | -.0513287 .1691808 -0.30 0.762 -.3830085 .2803511 /xb3_WEST | -.0710776 .1835512 -0.39 0.699 -.4309306 .2887754 /xb3_AGE | -.1397973 .0902243 -1.55 0.121 -.3166824 .0370878 /xb3_BLACK | .2694216 .2088309 1.29 0.197 -.1399924 .6788357 /xb3_MALE | .0712789 .1573248 0.45 0.651 -.237157 .3797148 /xb3_MARRIED | -.1524561 .149551 -1.02 0.308 -.4456515 .1407393 /xb3_SCHOOL | -.0010404 .0191956 -0.05 0.957 -.0386735 .0365926 /xb3_FAMINC | .0092552 .0213445 0.43 0.665 -.0325908 .0511012 /xb3_EMPLO~D | .2935965 .1976169 1.49 0.137 -.0938323 .6810254 /xb3_PRIVINS | .0132496 .1556414 0.09 0.932 -.2918861 .3183852 /xb3_MEDIC~D | .062663 .2043911 0.31 0.759 -.3380468 .4633727 /xb3_CONST~T | -1.07612 .6324334 -1.70 0.089 -2.316009 .1637682 ------------------------------------------------------------------------------ . nl (HOSP = exp({xb4: $XLIST CONSTANT})), vce(robust) nolog (obs = 4406) Nonlinear regression Number of obs = 4406 R-squared = 0.2113 Adj R-squared = 0.2082 Root MSE = .7143959 Res. dev. = 9523.02 ------------------------------------------------------------------------------ | Robust HOSP | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- /xb4_EXCLH~H | -.6971222 .2095509 -3.33 0.001 -1.107948 -.2862967 /xb4_POORH~H | .6142087 .0935843 6.56 0.000 .4307361 .7976812 /xb4_NUMCH~N | .2320049 .0290481 7.99 0.000 .175056 .2889538 /xb4_ADLDIFF | .4160149 .104057 4.00 0.000 .2120107 .6200191 /xb4_NOREAST | -.0559222 .1328301 -0.42 0.674 -.3163362 .2044919 /xb4_MIDWEST | .0641432 .1404636 0.46 0.648 -.2112364 .3395228 /xb4_WEST | -.0325414 .1329101 -0.24 0.807 -.2931123 .2280294 /xb4_AGE | -.0369077 .0751163 -0.49 0.623 -.1841735 .1103582 /xb4_BLACK | .1053548 .1328203 0.79 0.428 -.1550401 .3657496 /xb4_MALE | .0422865 .1178446 0.36 0.720 -.1887484 .2733213 /xb4_MARRIED | .0323321 .1055873 0.31 0.759 -.1746723 .2393364 /xb4_SCHOOL | .0053263 .0164276 0.32 0.746 -.02688 .0375327 /xb4_FAMINC | .0160215 .0120372 1.33 0.183 -.0075774 .0396204 /xb4_EMPLO~D | .059557 .1532933 0.39 0.698 -.2409753 .3600892 /xb4_PRIVINS | .3228933 .1221102 2.64 0.008 .0834957 .5622908 /xb4_MEDIC~D | .2513341 .1364241 1.84 0.065 -.016126 .5187943 /xb4_CONST~T | -2.028423 .5948647 -3.41 0.001 -3.194658 -.8621875 ------------------------------------------------------------------------------ . quietly poisson EMR $XLIST, vce(robust) nolog . estimates store Poi_EMR . quietly poisson HOSP $XLIST, vce(robust) nolog . estimates store Poi_HOSP . quietly nbreg EMR $XLIST, vce(robust) nolog . estimates store NB_EMR . quietly nbreg HOSP $XLIST, vce(robust) nolog . estimates store NB_HOSP . estimates table Poi_EMR Poi_HOSP NB_EMR NB_HOSP, b(%9.3f) se(%9.2f) eq(1) -------------------------------------------------------------- Variable | Poi_EMR Poi_HOSP NB_EMR NB_HOSP -------------+------------------------------------------------ #1 | EXCLHLTH | -0.625 -0.710 -0.637 -0.689 | 0.21 0.19 0.21 0.19 POORHLTH | 0.516 0.535 0.494 0.511 | 0.10 0.09 0.10 0.10 NUMCHRON | 0.223 0.251 0.221 0.276 | 0.03 0.02 0.03 0.02 ADLDIFF | 0.419 0.345 0.404 0.338 | 0.09 0.09 0.09 0.09 NOREAST | 0.025 -0.009 0.047 -0.000 | 0.11 0.10 0.10 0.10 MIDWEST | -0.001 0.112 0.021 0.134 | 0.10 0.10 0.10 0.10 WEST | 0.142 0.099 0.177 0.135 | 0.12 0.11 0.11 0.11 AGE | 0.036 0.117 0.072 0.172 | 0.06 0.06 0.06 0.06 BLACK | 0.171 0.096 0.175 0.103 | 0.12 0.11 0.12 0.12 MALE | 0.060 0.154 0.068 0.213 | 0.09 0.08 0.09 0.08 MARRIED | -0.126 -0.027 -0.108 -0.034 | 0.09 0.08 0.09 0.09 SCHOOL | -0.018 0.002 -0.018 0.001 | 0.01 0.01 0.01 0.01 FAMINC | 0.008 0.007 0.004 0.001 | 0.01 0.01 0.01 0.01 EMPLOYED | 0.176 0.039 0.172 0.044 | 0.13 0.13 0.13 0.13 PRIVINS | 0.058 0.215 0.060 0.184 | 0.10 0.10 0.10 0.11 MEDICAID | 0.144 0.180 0.181 0.141 | 0.13 0.12 0.13 0.14 _cons | -2.148 -3.080 -2.430 -3.512 | 0.49 0.47 0.49 0.49 -------------+------------------------------------------------ lnalpha | _cons | 0.499 0.538 | 0.11 0.10 -------------------------------------------------------------- legend: b/se . . ******* TABLE in the BOOK . . *** TABLE 8.3: ML and NLSUR bivariate estimates . * Note: Following gives default se's for NB1FE and not jackknife se's (given above) . estimates table MLBVNB NLSUR, b(%7.4f) se(%7.3f) stats(N ll) stfmt(%9.1f) modelwidth(9) -------------------------------------- Variable | MLBVNB NLSUR -------------+------------------------ eq1 | EXCLHLTH | -0.6183 | 0.213 POORHLTH | 0.4880 | 0.097 NUMCHRON | 0.2405 | 0.028 ADLDIFF | 0.4045 | 0.092 NOREAST | 0.0435 | 0.105 MIDWEST | 0.0277 | 0.098 WEST | 0.1817 | 0.115 AGE | 0.0989 | 0.061 BLACK | 0.1737 | 0.117 MALE | 0.1107 | 0.088 MARRIED | -0.1115 | 0.090 SCHOOL | -0.0174 | 0.011 FAMINC | -0.0009 | 0.015 EMPLOYED | 0.1858 | 0.126 PRIVINS | 0.0295 | 0.102 MEDICAID | 0.1457 | 0.135 _cons | -1.8756 | 0.507 -------------+------------------------ eq2 | EXCLHLTH | -0.7012 | 0.188 POORHLTH | 0.5006 | 0.097 NUMCHRON | 0.2717 | 0.024 ADLDIFF | 0.3366 | 0.092 NOREAST | 0.0242 | 0.104 MIDWEST | 0.1586 | 0.097 WEST | 0.1612 | 0.108 AGE | 0.1814 | 0.059 BLACK | 0.0957 | 0.114 MALE | 0.2063 | 0.082 MARRIED | -0.0145 | 0.088 SCHOOL | 0.0015 | 0.011 FAMINC | -0.0018 | 0.011 EMPLOYED | 0.0509 | 0.132 PRIVINS | 0.1886 | 0.110 MEDICAID | 0.1821 | 0.136 _cons | -2.8248 | 0.491 -------------+------------------------ eq3 | _cons | 2.1659 | 0.132 -------------+------------------------ xb1_EXCLHLTH | _cons | -0.4266 | 0.254 -------------+------------------------ xb1_POORHLTH | _cons | 0.6386 | 0.111 -------------+------------------------ xb1_NUMCHRON | _cons | 0.2671 | 0.044 -------------+------------------------ xb1_ADLDIFF | _cons | 0.4296 | 0.166 -------------+------------------------ xb1_NOREAST | _cons | -0.1195 | 0.174 -------------+------------------------ xb1_MIDWEST | _cons | -0.0339 | 0.171 -------------+------------------------ xb1_WEST | _cons | -0.0779 | 0.183 -------------+------------------------ xb1_AGE | _cons | -0.1485 | 0.090 -------------+------------------------ xb1_BLACK | _cons | 0.2789 | 0.213 -------------+------------------------ xb1_MALE | _cons | 0.0788 | 0.158 -------------+------------------------ xb1_MARRIED | _cons | -0.1596 | 0.151 -------------+------------------------ xb1_SCHOOL | _cons | -0.0013 | 0.019 -------------+------------------------ xb1_FAMINC | _cons | 0.0085 | 0.021 -------------+------------------------ xb1_EMPLOYED | _cons | 0.3200 | 0.205 -------------+------------------------ xb1_PRIVINS | _cons | 0.0020 | 0.155 -------------+------------------------ xb1_MEDICAID | _cons | 0.0508 | 0.205 -------------+------------------------ xb1_CONSTANT | _cons | -1.0262 | 0.626 -------------+------------------------ xb2_EXCLHLTH | _cons | -0.6465 | 0.206 -------------+------------------------ xb2_POORHLTH | _cons | 0.6365 | 0.095 -------------+------------------------ xb2_NUMCHRON | _cons | 0.2417 | 0.030 -------------+------------------------ xb2_ADLDIFF | _cons | 0.4033 | 0.108 -------------+------------------------ xb2_NOREAST | _cons | -0.0716 | 0.133 -------------+------------------------ xb2_MIDWEST | _cons | 0.0824 | 0.140 -------------+------------------------ xb2_WEST | _cons | -0.0543 | 0.133 -------------+------------------------ xb2_AGE | _cons | -0.0568 | 0.076 -------------+------------------------ xb2_BLACK | _cons | 0.1296 | 0.137 -------------+------------------------ xb2_MALE | _cons | 0.0475 | 0.120 -------------+------------------------ xb2_MARRIED | _cons | 0.0245 | 0.108 -------------+------------------------ xb2_SCHOOL | _cons | 0.0088 | 0.017 -------------+------------------------ xb2_FAMINC | _cons | 0.0147 | 0.013 -------------+------------------------ xb2_EMPLOYED | _cons | 0.0992 | 0.153 -------------+------------------------ xb2_PRIVINS | _cons | 0.3213 | 0.125 -------------+------------------------ xb2_MEDICAID | _cons | 0.2585 | 0.142 -------------+------------------------ xb2_CONSTANT | _cons | -1.9520 | 0.601 -------------+------------------------ Statistics | N | 4406 4406 ll | -5222.9 -8811.1 -------------------------------------- legend: b/se . . ********** CLOSE OUTPUT . . * log close . * clear . * exit . end of do-file . exit, clear