------------------------------------------------------------------------------------------------------------------------------- name: log: c:\acdbookrevision\stata_final_programs_2013\racd11.txt log type: text opened on: 19 Jan 2013, 11:51:18 . . ********** OVERVIEW OF racd11.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 analyzes Patents data in Chapter 11 . * 11.8.1 PATENTS - TABLE 11.1 only (parametric models) . . * NOTE: The rest of chapter 11 is done using R program racd11.R . . * To run you need files . * racd09data.dta . * and user-written Stata addon . * fmm . * in your directory . . ********** SETUP ********** . . set more off . version 11.2 . clear all . * set linesize 82 . set scheme s1mono /* Graphics scheme */ . . ************ . . * NOTE: This program produces just TABLE 11.1 . * The rest of chapter 11 is done using R program racd09.R . . ********** DATA DESCRIPTION . . * The original data is from . * Bronwyn Hall, Zvi Griliches, and Jerry Hausman (1986), . * "Patents and R&D: Is There a Lag?", . * International Economic Review, 27, 265-283. . * See this article for more detailed discussion . * Also see racd09makedata.do for further details . . * NOTE: Here we use just 1979 data . . ********** 11.8.1 PATENTS - PARAMETRIC MODELS . . use racd09data.dta, clear . . * Create log of total R&D over five years . generate LOGRandD = ln(exp(LOGR)+exp(LOGR1)+exp(LOGR2)+exp(LOGR3)+exp(LOGR4)+exp(LOGR5)) . . * Use only 1979 data . keep if YEAR==5 (1384 observations deleted) . . * Regressor list . global XLIST LOGRandD LOGK SCISECT . * global XLIST LOGR LOGR1 LOGR2 LOGR3 LOGR4 LOGR5 LOGK SCISECT . . * Variable descriptions and summary statistics . describe PAT $XLIST storage display value variable name type format label variable label ------------------------------------------------------------------------------------------------------------------------------- PAT float %9.0g Number of (successful) patents applied for this year LOGRandD float %9.0g LOGK float %9.0g Logarithm of the book value of capital in 1972 SCISECT float %9.0g Equals 1 if firm in the scientific sector . summarize PAT $XLIST Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- PAT | 346 32.10116 66.36197 0 515 LOGRandD | 346 3.071489 1.965798 -1.543733 8.723914 LOGK | 346 3.921063 2.095542 -1.76965 9.66626 SCISECT | 346 .4248555 .4950369 0 1 . . ****** POISSON AND NB2: Model estimates and (grouped) predicted probabilities . . * Data frequencies grouped . generate yle0 = PAT <= 0 . generate yle1 = PAT <= 1 . generate yle2 = PAT <= 2 . generate yle5 = PAT <= 5 . generate yle10 = PAT <= 10 . generate yle20 = PAT <= 20 . generate yle50 = PAT <= 50 . generate yle100 = PAT <= 100 . generate yle300 = PAT <= 300 . generate yle500 = PAT <= 500 . . * Poisson regression . poisson PAT $XLIST, vce(robust) Iteration 0: log pseudolikelihood = -3366.1483 Iteration 1: log pseudolikelihood = -3365.8303 Iteration 2: log pseudolikelihood = -3365.8303 Poisson regression Number of obs = 346 Wald chi2(3) = 641.60 Prob > chi2 = 0.0000 Log pseudolikelihood = -3365.8303 Pseudo R2 = 0.7632 ------------------------------------------------------------------------------ | Robust PAT | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- LOGRandD | .4753631 .064298 7.39 0.000 .3493413 .601385 LOGK | .2708473 .0552599 4.90 0.000 .1625398 .3791547 SCISECT | .4672292 .1549081 3.02 0.003 .163615 .7708435 _cons | -.3402801 .1900477 -1.79 0.073 -.7127667 .0322066 ------------------------------------------------------------------------------ . estimates store POISSON . predict ple0, pr(0) . predict ple1, pr(0,1) . predict ple2, pr(0,2) . predict ple5, pr(0,5) . predict ple10, pr(0,10) . predict ple20, pr(0,20) . predict ple50, pr(0,50) . predict ple100, pr(0,100) . predict ple300, pr(0,300) . predict ple500, pr(0,1000) . . * NB2 intercept-only . nbreg PAT, vce(robust) Fitting Poisson model: Iteration 0: log pseudolikelihood = -14212.49 Iteration 1: log pseudolikelihood = -14212.49 Fitting constant-only model: Iteration 0: log pseudolikelihood = -1551.5708 Iteration 1: log pseudolikelihood = -1347.2714 Iteration 2: log pseudolikelihood = -1346.8824 Iteration 3: log pseudolikelihood = -1346.8822 Iteration 4: log pseudolikelihood = -1346.8822 Fitting full model: Iteration 0: log pseudolikelihood = -1346.8822 Iteration 1: log pseudolikelihood = -1346.8822 Negative binomial regression Number of obs = 346 Dispersion = mean Wald chi2(0) = . Log pseudolikelihood = -1346.8822 Prob > chi2 = . ------------------------------------------------------------------------------ | Robust PAT | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- _cons | 3.468892 .1111375 31.21 0.000 3.251067 3.686718 -------------+---------------------------------------------------------------- /lnalpha | 1.263334 .0654783 1.134999 1.391669 -------------+---------------------------------------------------------------- alpha | 3.537195 .2316097 3.11117 4.021557 ------------------------------------------------------------------------------ . predict nbintle0, pr(0) . predict nbintle1, pr(0,1) . predict nbintle2, pr(0,2) . predict nbintle5, pr(0,5) . predict nbintle10, pr(0,10) . predict nbintle20, pr(0,20) . predict nbintle50, pr(0,50) . predict nbintle100, pr(0,100) . predict nbintle300, pr(0,300) . predict nbintle500, pr(0,1000) . . * NB2 regression . nbreg PAT $XLIST, vce(robust) Fitting Poisson model: Iteration 0: log pseudolikelihood = -3366.1483 Iteration 1: log pseudolikelihood = -3365.8303 Iteration 2: log pseudolikelihood = -3365.8303 Fitting constant-only model: Iteration 0: log pseudolikelihood = -1551.5708 Iteration 1: log pseudolikelihood = -1347.2714 Iteration 2: log pseudolikelihood = -1346.8824 Iteration 3: log pseudolikelihood = -1346.8822 Iteration 4: log pseudolikelihood = -1346.8822 Fitting full model: Iteration 0: log pseudolikelihood = -1277.5529 (not concave) Iteration 1: log pseudolikelihood = -1196.6667 Iteration 2: log pseudolikelihood = -1137.5362 Iteration 3: log pseudolikelihood = -1127.8388 Iteration 4: log pseudolikelihood = -1127.5664 Iteration 5: log pseudolikelihood = -1127.5662 Iteration 6: log pseudolikelihood = -1127.5662 Negative binomial regression Number of obs = 346 Dispersion = mean Wald chi2(3) = 616.24 Log pseudolikelihood = -1127.5662 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust PAT | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- LOGRandD | .8159469 .0827754 9.86 0.000 .6537101 .9781837 LOGK | .114636 .0687745 1.67 0.096 -.0201595 .2494315 SCISECT | -.0897712 .1381794 -0.65 0.516 -.3605979 .1810555 _cons | -.8106295 .180078 -4.50 0.000 -1.163576 -.4576832 -------------+---------------------------------------------------------------- /lnalpha | -.1224513 .1134868 -.3448814 .0999788 -------------+---------------------------------------------------------------- alpha | .884749 .1004074 .7083043 1.105148 ------------------------------------------------------------------------------ . estimates store NB2 . predict nble0, pr(0) . predict nble1, pr(0,1) . predict nble2, pr(0,2) . predict nble5, pr(0,5) . predict nble10, pr(0,10) . predict nble20, pr(0,20) . predict nble50, pr(0,50) . predict nble100, pr(0,100) . predict nble300, pr(0,300) . predict nble500, pr(0,1000) . . * And model fit for NB2 . predict yhat, n . summarize yhat PAT Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- yhat | 346 41.60105 130.0129 .1298683 1532.695 PAT | 346 32.10116 66.36197 0 515 . correlate yhat PAT (obs=346) | yhat PAT -------------+------------------ yhat | 1.0000 PAT | 0.6861 1.0000 . display "Squared correlation of y and yhat = " r(rho)^2 Squared correlation of y and yhat = .47070323 . . * Aside: Calculate Li and Racine R2 measure . quietly sum PAT . scalar ymean = r(mean) . generate ycross = (PAT-ymean)*(yhat-ymean) . generate ylessybarsq = (PAT-ymean)^2 . generate yhatlessybarsq = (yhat-ymean)^2 . quietly summarize ycross . scalar Rnumerator = r(sum) . quietly summarize ylessybarsq . scalar Rdenominator1 = r(sum) . quietly summarize yhatlessybarsq . scalar Rdenominator2 = r(sum) . display "Li and Racine R-squared = " Rnumerator^2 / (Rdenominator1*Rdenominator2) Li and Racine R-squared = .46819624 . . * Find the mode probability for each observation (i.e. k than maximizes Pr[y = k] . predict nbfit0, pr(0) . generate mode = 0 . forvalues i = 1/600 { 2. predict nbfit1, pr(`i') 3. quietly replace mode = `i' if nbfit1> nbfit0 4. quietly replace nbfit0 = nbfit1 5. drop nbfit1 6. } . . * Compare the actual count to the predicted mode . recode mode (0=0) (1=1) (2/5=2) (6/20=6) (21/50=31) (51/100=51) /// > (101/200=101) (201/300=201) (301/600=301), gen(modegrouped) (90 differences between mode and modegrouped) . recode PAT (0=0) (1=1) (2/5=2) (6/20=6) (21/50=31) (51/100=51) /// > (101/200=101) (201/300=201) (301/600=301), gen(PATgrouped) (194 differences between PAT and PATgrouped) . tabulate PATgrouped modegrouped RECODE of | PAT | (Number of | (successfu | l) patents | applied | for this | RECODE of mode year) | 0 1 2 6 31 101 | Total -----------+------------------------------------------------------------------+---------- 0 | 73 3 0 0 0 0 | 76 1 | 34 5 2 0 0 0 | 41 2 | 59 6 8 1 0 0 | 74 6 | 18 18 13 2 0 0 | 51 31 | 3 6 23 10 1 0 | 43 51 | 1 2 8 19 0 0 | 30 101 | 0 0 2 12 4 1 | 19 201 | 0 0 0 3 5 0 | 8 301 | 0 0 0 1 1 2 | 4 -----------+------------------------------------------------------------------+---------- Total | 188 40 56 48 11 3 | 346 . tabulate mode mode | Freq. Percent Cum. ------------+----------------------------------- 0 | 188 54.34 54.34 1 | 40 11.56 65.90 2 | 22 6.36 72.25 3 | 13 3.76 76.01 4 | 15 4.34 80.35 5 | 6 1.73 82.08 6 | 5 1.45 83.53 7 | 8 2.31 85.84 8 | 6 1.73 87.57 9 | 5 1.45 89.02 10 | 4 1.16 90.17 11 | 3 0.87 91.04 12 | 2 0.58 91.62 13 | 2 0.58 92.20 15 | 3 0.87 93.06 16 | 4 1.16 94.22 17 | 1 0.29 94.51 18 | 3 0.87 95.38 20 | 2 0.58 95.95 22 | 3 0.87 96.82 23 | 1 0.29 97.11 27 | 1 0.29 97.40 31 | 1 0.29 97.69 32 | 1 0.29 97.98 40 | 1 0.29 98.27 42 | 1 0.29 98.55 43 | 2 0.58 99.13 124 | 1 0.29 99.42 136 | 1 0.29 99.71 176 | 1 0.29 100.00 ------------+----------------------------------- Total | 346 100.00 . count if PAT == mode 80 . count if PATgrouped == modegrouped 90 . . sum yle0 yle1 yle2 yle5 yle10 yle20 yle50 yle100 yle300 yle500 Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- yle0 | 346 .2196532 .4146112 0 1 yle1 | 346 .3381503 .4737651 0 1 yle2 | 346 .4104046 .4926196 0 1 yle5 | 346 .5520231 .4980064 0 1 yle10 | 346 .6184971 .486459 0 1 -------------+-------------------------------------------------------- yle20 | 346 .699422 .4591734 0 1 yle50 | 346 .8236994 .3816276 0 1 yle100 | 346 .9104046 .2860148 0 1 yle300 | 346 .9884393 .1070522 0 1 yle500 | 346 .9971098 .0537603 0 1 . sum ple0 ple1 ple2 ple5 ple10 ple20 ple50 ple100 ple300 ple500 Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- ple0 | 346 .0486209 .1152675 0 .6936235 ple1 | 346 .1176931 .2228108 0 .947369 ple2 | 346 .1886973 .2958484 0 .9937823 ple5 | 346 .3720066 .4080619 0 .9999976 ple10 | 346 .5343497 .4621704 0 1 -------------+-------------------------------------------------------- ple20 | 346 .6622245 .4507554 0 1 ple50 | 346 .8251793 .362435 0 1 ple100 | 346 .9194856 .2596965 0 1 ple300 | 346 .9907562 .0933535 0 1 ple500 | 346 1 0 1 1 . sum nbintle0 nbintle1 nbintle2 nbintle5 nbintle10 nbintle20 nbintle50 /// > nbintle100 nbintle300 nbintle500 Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- nbintle0 | 346 .2617603 0 .2617603 .2617603 nbintle1 | 346 .3351165 0 .3351165 .3351165 nbintle2 | 346 .3817531 0 .3817531 .3817531 nbintle5 | 346 .4698427 0 .4698427 .4698427 nbintle10 | 346 .5567626 0 .5567626 .5567626 -------------+-------------------------------------------------------- nbintle20 | 346 .6591401 0 .6591401 .6591401 nbintle50 | 346 .8063407 0 .8063407 .8063407 nbintle100 | 346 .9063662 0 .9063662 .9063662 nbintle300 | 346 .9907655 0 .9907655 .9907655 nbintle500 | 346 .9999905 0 .9999905 .9999905 . sum nble0 nble1 nble2 nble5 nble10 nble20 nble50 nble100 nble300 nble500 Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- nble0 | 346 .1902577 .2071034 .000288 .8843223 nble1 | 346 .3082096 .2870145 .0006132 .9873318 nble2 | 346 .3886912 .3250769 .0009594 .9986393 nble5 | 346 .5295132 .3637599 .0020745 .9999984 nble10 | 346 .6401065 .3642567 .0040861 1 -------------+-------------------------------------------------------- nble20 | 346 .7386319 .3325746 .0084271 1 nble50 | 346 .8491667 .2573956 .0226604 1 nble100 | 346 .9142927 .1897329 .0480778 1 nble300 | 346 .9751445 .0987155 .1529473 1 nble500 | 346 .995048 .0436341 .4613153 1 . . ***** NB2 MODEL WITH FIRST ORDER POLYNOMIAL - Estimates and (grouped) predicted probabilities . . program lfnb2p1 1. version 11 2. args lnf theta1 alpha a // theta1=x'b, alpha=alpha, a=a 3. tempvar mu 4. local y $ML_y1 // Define y so program more readable 5. generate double `mu' = exp(`theta1') 6. quietly replace `lnf' = lngamma(`y'+(1/`alpha')) - lngamma((1/`alpha')) /// > - lnfactorial(`y') - (`y'+(1/`alpha'))*ln(1+`alpha'*`mu') /// > + `y'*ln(`alpha') + `y'*ln(`mu') /// > + ln((1+`a'*`y')^2) - ln(1+2*`a'*`mu'+(`a'^2)*(`mu'+(1+`alpha')*(`mu'^2))) 7. end . ml model lf lfnb2p1 (PAT = $XLIST) () (), vce(robust) . ml search initial: log pseudolikelihood = - (could not be evaluated) feasible: log pseudolikelihood = -7906.086 rescale: log pseudolikelihood = -2231.5024 rescale eq: log pseudolikelihood = -1394.2889 . * ml init .7 .3 .0 1.0 .0 -.1 .0 -.6 .9 -.5, copy . ml maximize initial: log pseudolikelihood = -1394.2889 rescale: log pseudolikelihood = -1394.2889 rescale eq: log pseudolikelihood = -1394.2889 Iteration 0: log pseudolikelihood = -1394.2889 (not concave) Iteration 1: log pseudolikelihood = -1287.3643 (not concave) Iteration 2: log pseudolikelihood = -1255.5891 Iteration 3: log pseudolikelihood = -1179.8953 (backed up) Iteration 4: log pseudolikelihood = -1121.1922 Iteration 5: log pseudolikelihood = -1118.1782 Iteration 6: log pseudolikelihood = -1116.8733 Iteration 7: log pseudolikelihood = -1116.8135 Iteration 8: log pseudolikelihood = -1116.8123 Iteration 9: log pseudolikelihood = -1116.8123 Number of obs = 346 Wald chi2(3) = 477.31 Log pseudolikelihood = -1116.8123 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust PAT | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- eq1 | LOGRandD | .6113362 .0658558 9.28 0.000 .4822612 .7404112 LOGK | .1104351 .0503529 2.19 0.028 .0117452 .209125 SCISECT | .0064389 .1071066 0.06 0.952 -.2034862 .2163639 _cons | -.9351948 .1991611 -4.70 0.000 -1.325543 -.5448463 -------------+---------------------------------------------------------------- eq2 | _cons | 1.50628 .2368957 6.36 0.000 1.041973 1.970587 -------------+---------------------------------------------------------------- eq3 | _cons | .0817346 .0258581 3.16 0.002 .0310536 .1324156 ------------------------------------------------------------------------------ . estimates store NB2P1 . estat ic ----------------------------------------------------------------------------- Model | Obs ll(null) ll(model) df AIC BIC -------------+--------------------------------------------------------------- NB2P1 | 346 . -1116.812 6 2245.625 2268.703 ----------------------------------------------------------------------------- Note: N=Obs used in calculating BIC; see [R] BIC note . . * Use the results to get predicted probabilities . matrix b = e(b) . scalar a1 = b[1,e(k)] . scalar alpha = b[1,e(k)-1] . matrix b = b[1,1..e(k)-2] . generate mu = 0 . generate one = 1 . local i 1 . foreach var of varlist $XLIST one { 2. scalar beta = b[1,`i'] 3. quietly replace mu = mu + beta*`var' 4. local i = `i' + 1 5. } . replace mu = exp(mu) (346 real changes made) . scalar ainv = 1/alpha . forvalues i = 0/500 { 2. generate pfitnb2p1`i' = lngamma(`i'+ainv) - lngamma(ainv) - lnfactorial(`i') + /// > ainv*ln(ainv/(ainv+mu)) + `i'*ln(mu/(ainv+mu)) + ln((1+a1*`i')^2) - ln(1+2*a1*mu+(a1^2)*(mu+(1+alpha)*(mu^2))) 3. quietly replace pfitnb2p1`i' = exp(pfitnb2p1`i') 4. } . generate nb2p1le0 = pfitnb2p10 . generate nb2p1le1 = pfitnb2p10 + pfitnb2p11 . generate nb2p1le2 = pfitnb2p10 + pfitnb2p11 + pfitnb2p12 . generate nb2p1le5=pfitnb2p10+pfitnb2p11+pfitnb2p12+pfitnb2p13+pfitnb2p13+pfitnb2p15 . generate nb2p1le10 = 0 . generate nb2p1le20 = 0 . generate nb2p1le50 = 0 . generate nb2p1le100 = 0 . generate nb2p1le300 = 0 . generate nb2p1le500 = 0 . forvalues i = 0/500 { 2. quietly replace nb2p1le10 = nb2p1le10 + pfitnb2p1`i' if `i'<=10 3. quietly replace nb2p1le20 = nb2p1le20 + pfitnb2p1`i' if `i'<=20 4. quietly replace nb2p1le50 = nb2p1le50 + pfitnb2p1`i' if `i'<=50 5. quietly replace nb2p1le100 = nb2p1le100 + pfitnb2p1`i' if `i'<=100 6. quietly replace nb2p1le300 = nb2p1le300 + pfitnb2p1`i' if `i'<=300 7. quietly replace nb2p1le500 = nb2p1le500 + pfitnb2p1`i' if `i'<=500 8. } . . * The following gives NB2P1 parametric results corresponding to Table 11.2 . summarize nb2p1le* Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- nb2p1le0 | 346 .2242861 .2117458 .0000254 .8455491 nb2p1le1 | 346 .3320578 .2863826 .000045 .9712303 nb2p1le2 | 346 .4016375 .3239715 .0000638 .9943737 nb2p1le5 | 346 .5332051 .3782282 .0001228 1.01478 nb2p1le10 | 346 .618834 .3800833 .0002459 1 -------------+-------------------------------------------------------- nb2p1le20 | 346 .7114312 .3614461 .0006089 1 nb2p1le50 | 346 .8293393 .2899445 .002928 1 nb2p1le100 | 346 .9087492 .2087086 .0118645 1 nb2p1le300 | 346 .9811061 .0934584 .1129382 1 nb2p1le500 | 346 .9920137 .0624729 .278686 1 . . ****** FINITE MIXTURE MODEL NB2-2 component - Estimates . . * Following does not converge . fmm PAT $XLIST, mixtureof(negbin2) components(2) vce(robust) Fitting Negative Binomial-2 model: Iteration 0: log likelihood = -3366.1483 Iteration 1: log likelihood = -3365.8303 Iteration 2: log likelihood = -3365.8303 Iteration 0: log likelihood = -1551.5708 Iteration 1: log likelihood = -1347.2714 Iteration 2: log likelihood = -1346.8824 Iteration 3: log likelihood = -1346.8822 Iteration 4: log likelihood = -1346.8822 Iteration 0: log likelihood = -1277.5529 (not concave) Iteration 1: log likelihood = -1196.6667 Iteration 2: log likelihood = -1137.5362 Iteration 3: log likelihood = -1127.8388 Iteration 4: log likelihood = -1127.5664 Iteration 5: log likelihood = -1127.5662 Iteration 6: log likelihood = -1127.5662 Fitting 2 component Negative Binomial-2 model: Iteration 0: log pseudolikelihood = -1127.5701 (not concave) Iteration 1: log pseudolikelihood = -1127.5046 (not concave) Iteration 2: log pseudolikelihood = -1118.0772 (not concave) Iteration 3: log pseudolikelihood = -1115.7068 Iteration 4: log pseudolikelihood = -1115.4488 (not concave) Iteration 5: log pseudolikelihood = -1109.4996 (not concave) Iteration 6: log pseudolikelihood = -1107.1628 Iteration 7: log pseudolikelihood = -1106.4365 Iteration 8: log pseudolikelihood = -1106.1363 Iteration 9: log pseudolikelihood = -1106.1329 Iteration 10: log pseudolikelihood = -1106.1329 2 component Negative Binomial-2 regression Number of obs = 346 Wald chi2(6) = 1496.35 Log pseudolikelihood = -1106.1329 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust PAT | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- component1 | LOGRandD | .9181075 .0634763 14.46 0.000 .7936963 1.042519 LOGK | .0653257 .056885 1.15 0.251 -.0461668 .1768182 SCISECT | .1028518 .1416338 0.73 0.468 -.1747454 .380449 _cons | -1.273755 .1607671 -7.92 0.000 -1.588853 -.9586572 -------------+---------------------------------------------------------------- component2 | LOGRandD | .3519998 .0896255 3.93 0.000 .1763371 .5276625 LOGK | .2771243 .0617372 4.49 0.000 .1561216 .398127 SCISECT | -.5731597 .3459274 -1.66 0.098 -1.251165 .1048455 _cons | 1.516435 .1911198 7.93 0.000 1.141847 1.891023 -------------+---------------------------------------------------------------- /imlogitpi1 | 2.571207 .3377545 7.61 0.000 1.90922 3.233194 /lnalpha1 | -.4287745 .1222261 -3.51 0.000 -.6683332 -.1892157 /lnalpha2 | -2.77605 .7533299 -3.69 0.000 -4.252549 -1.29955 ------------------------------------------------------------------------------ alpha1 | .6513068 .0796067 .5125622 .827608 alpha2 | .0622841 .0469204 .0142279 .2726544 pi1 | .9289854 .0222822 .8709315 .9620645 pi2 | .0710146 .0222822 .0379355 .1290685 ------------------------------------------------------------------------------ . estimates store FM2NB2 . . ****** POISSON MODEL WITH POLYNOMIALS HAS LOW LOG-LIKELIHOOD .... . . *** POISSON MODEL WITH FIRST ORDER POLYNOMIAL . . program lfpp1 1. version 11 2. args lnf theta1 a1 // theta1=x'b, alpha=alpha, a=a 3. tempvar mu m1 m2 4. local y $ML_y1 // Define y so program more readable 5. generate double `mu' = exp(`theta1') 6. generate double `m1' = `mu' 7. generate double `m2' = `mu'+`mu'^2 8. quietly replace `lnf' = -`mu' + `y'*`theta1' - lnfactorial(`y') /// > + ln((1+`a1'*`y')^2) - ln(1 + 2*`a1'*`m1' + (`a1'^2)*`m2' ) 9. end . ml model lf lfpp1 (PAT = $XLIST) (), vce(robust) . ml search initial: log pseudolikelihood = -41980.474 alternative: log pseudolikelihood = -35993.508 rescale: log pseudolikelihood = -17754.154 rescale eq: log pseudolikelihood = -16081.453 . * ml init .5 .5 -.04 .5 .03 -.06 .03 -2.3 5.3 .8, copy . ml maximize, iter(20) initial: log pseudolikelihood = -16081.453 rescale: log pseudolikelihood = -16081.453 rescale eq: log pseudolikelihood = -16081.453 Iteration 0: log pseudolikelihood = -16081.453 Iteration 1: log pseudolikelihood = -6497.119 Iteration 2: log pseudolikelihood = -3435.8267 Iteration 3: log pseudolikelihood = -3373.8173 Iteration 4: log pseudolikelihood = -3372.346 Iteration 5: log pseudolikelihood = -3371.9521 Iteration 6: log pseudolikelihood = -3371.9431 Iteration 7: log pseudolikelihood = -3371.9431 Number of obs = 346 Wald chi2(3) = 684.54 Log pseudolikelihood = -3371.9431 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust PAT | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- eq1 | LOGRandD | .4728655 .0594023 7.96 0.000 .3564391 .5892919 LOGK | .2715676 .0541807 5.01 0.000 .1653754 .3777597 SCISECT | .4700587 .1436189 3.27 0.001 .1885708 .7515465 _cons | -.3292864 .1932391 -1.70 0.088 -.7080281 .0494553 -------------+---------------------------------------------------------------- eq2 | _cons | -.0020155 .0000812 -24.84 0.000 -.0021745 -.0018564 ------------------------------------------------------------------------------ . estat ic ----------------------------------------------------------------------------- Model | Obs ll(null) ll(model) df AIC BIC -------------+--------------------------------------------------------------- . | 346 . -3371.943 5 6753.886 6773.118 ----------------------------------------------------------------------------- Note: N=Obs used in calculating BIC; see [R] BIC note . estimates store PP1 . . *** POISSON MODEL WITH SECOND ORDER POLYNOMIAL . . program lfpp2 1. version 11 2. args lnf theta1 a1 a2 // theta1=x'b, alpha=alpha, a=a 3. tempvar mu m1 m2 m3 m4 4. local y $ML_y1 // Define y so program more readable 5. generate double `mu' = exp(`theta1') 6. generate double `m1' = `mu' 7. generate double `m2' = `mu'+`mu'^2 8. generate double `m3' = `mu'+3*`mu'^2+`mu'^3 9. generate double `m4' = `mu'+7*`mu'^2+6*`mu'^3+`mu'^4 10. quietly replace `lnf' = -`mu' + `y'*`theta1' - lnfactorial(`y') /// > + ln((1+`a1'*`y'+`a2'*`y'^2)^2) /// > - ln(1 + 2*`a1'*`m1' + (`a1'^2+2*`a2')*`m2' + 2*`a1'*`a2'*`m3' + `a2'^2*`m4') 11. end . ml model lf lfpp2 (PAT = $XLIST) () (), vce(robust) . ml search initial: log pseudolikelihood = -41980.474 alternative: log pseudolikelihood = -35399.008 rescale: log pseudolikelihood = -19171.476 rescale eq: log pseudolikelihood = -17437.722 . * ml init .5 .5 -.04 .5 .03 -.06 .03 -2.3 5.3 .8, copy . ml maximize, iter(20) initial: log pseudolikelihood = -17437.722 rescale: log pseudolikelihood = -17437.722 rescale eq: log pseudolikelihood = -17437.722 Iteration 0: log pseudolikelihood = -17437.722 (not concave) Iteration 1: log pseudolikelihood = -8088.9579 (not concave) Iteration 2: log pseudolikelihood = -7470.3323 (not concave) Iteration 3: log pseudolikelihood = -7045.7402 (not concave) Iteration 4: log pseudolikelihood = -6060.8943 (not concave) Iteration 5: log pseudolikelihood = -5520.7862 (not concave) Iteration 6: log pseudolikelihood = -5139.7123 (not concave) Iteration 7: log pseudolikelihood = -4883.5333 (not concave) Iteration 8: log pseudolikelihood = -4757.1676 (not concave) Iteration 9: log pseudolikelihood = -4558.9023 (not concave) Iteration 10: log pseudolikelihood = -4437.0907 (not concave) Iteration 11: log pseudolikelihood = -4387.8102 (not concave) Iteration 12: log pseudolikelihood = -4348.9356 (not concave) Iteration 13: log pseudolikelihood = -4317.0448 (not concave) Iteration 14: log pseudolikelihood = -4291.5092 (not concave) Iteration 15: log pseudolikelihood = -4270.4285 (not concave) Iteration 16: log pseudolikelihood = -4253.3678 (not concave) Iteration 17: log pseudolikelihood = -4239.2011 (not concave) Iteration 18: log pseudolikelihood = -4227.6118 (not concave) Iteration 19: log pseudolikelihood = -4217.9416 (not concave) Iteration 20: log pseudolikelihood = -4209.9458 (not concave) convergence not achieved Number of obs = 346 Wald chi2(3) = 609.27 Log pseudolikelihood = -4209.9458 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust PAT | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- eq1 | LOGRandD | .5265669 .0692926 7.60 0.000 .3907559 .662378 LOGK | .2172299 .0602154 3.61 0.000 .0992099 .3352499 SCISECT | .4074079 .1625726 2.51 0.012 .0887715 .7260444 _cons | -.3115985 .1942118 -1.60 0.109 -.6922466 .0690496 -------------+---------------------------------------------------------------- eq2 | _cons | -20.43389 . . . . . -------------+---------------------------------------------------------------- eq3 | _cons | .2091537 .0005961 350.89 0.000 .2079854 .210322 ------------------------------------------------------------------------------ Warning: convergence not achieved . estat ic ----------------------------------------------------------------------------- Model | Obs ll(null) ll(model) df AIC BIC -------------+--------------------------------------------------------------- . | 346 . -4209.946 5 8429.892 8449.124 ----------------------------------------------------------------------------- Note: N=Obs used in calculating BIC; see [R] BIC note . estimates store PP2 . . *** POISSON MODEL WITH THIRD ORDER POLYNOMIAL . . program lfpp3 1. version 11 2. args lnf theta1 a1 a2 a3 // theta1=x'b, alpha=alpha, a=a 3. tempvar mu m1 m2 m3 m4 m5 m6 4. local y $ML_y1 // Define y so program more readable 5. generate double `mu' = exp(`theta1') 6. generate double `m1' = `mu' 7. generate double `m2' = `mu'+`mu'^2 8. generate double `m3' = `mu'+3*`mu'^2+`mu'^3 9. generate double `m4' = `mu'+7*`mu'^2+6*`mu'^3+`mu'^4 10. generate double `m5' = `mu'+15*`mu'^2+25*`mu'^3+10*`mu'^4+`mu'^5 11. generate double `m6' = `mu'+31*`mu'^2+90*`mu'^3+65*`mu'^4+15+`mu'^5+`mu'^6 12. quietly replace `lnf' = -`mu' + `y'*`theta1' - lnfactorial(`y') /// > + ln((1+`a1'*`y'+`a2'*`y'^2+`a3'*`y'^3)^2) /// > - ln(1 + 2*`a1'*`m1' + (`a1'^2+2*`a2')*`m2' + (2*`a1'*`a2'+2*`a3')*`m3' /// > + (2*`a1'*`a3'+`a2'^2)*`m4' + 2*`a2'*`a3'*`m5' + `a3'^2*`m6') 13. end . ml model lf lfpp3 (PAT = $XLIST) () () (), vce(robust) . ml search initial: log pseudolikelihood = -41980.474 alternative: log pseudolikelihood = -35031.215 rescale: log pseudolikelihood = -20579.661 rescale eq: log pseudolikelihood = -17373.859 . * ml init .72 .36 -.01 1.1 .02 -.06 .04 -.66 .94 -.52, copy . ml maximize, iter(20) initial: log pseudolikelihood = -17373.859 rescale: log pseudolikelihood = -17373.859 rescale eq: log pseudolikelihood = -17373.859 Iteration 0: log pseudolikelihood = -17373.859 (not concave) Iteration 1: log pseudolikelihood = -10136.582 (not concave) Iteration 2: log pseudolikelihood = -6559.9486 (not concave) Iteration 3: log pseudolikelihood = -5970.06 (not concave) Iteration 4: log pseudolikelihood = -5671.4683 (not concave) Iteration 5: log pseudolikelihood = -5429.3822 (not concave) Iteration 6: log pseudolikelihood = -5206.126 (not concave) Iteration 7: log pseudolikelihood = -5090.7878 (not concave) Iteration 8: log pseudolikelihood = -4987.0244 (not concave) Iteration 9: log pseudolikelihood = -4847.1549 (not concave) Iteration 10: log pseudolikelihood = -4730.1409 (not concave) Iteration 11: log pseudolikelihood = -4546.1279 (not concave) Iteration 12: log pseudolikelihood = -4483.7424 (not concave) Iteration 13: log pseudolikelihood = -4482.0319 (not concave) Iteration 14: log pseudolikelihood = -4480.4065 (not concave) Iteration 15: log pseudolikelihood = -4474.3321 (not concave) Iteration 16: log pseudolikelihood = -4470.7284 (not concave) Iteration 17: log pseudolikelihood = -4469.6513 (not concave) Iteration 18: log pseudolikelihood = -4468.5954 (not concave) Iteration 19: log pseudolikelihood = -4466.7069 (not concave) Iteration 20: log pseudolikelihood = -4465.1042 (not concave) convergence not achieved Number of obs = 346 Wald chi2(3) = 498.23 Log pseudolikelihood = -4465.1042 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust PAT | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- eq1 | LOGRandD | .5344883 .0687703 7.77 0.000 .399701 .6692756 LOGK | .2566202 .0573377 4.48 0.000 .1442404 .3690001 SCISECT | .4742346 .1714748 2.77 0.006 .1381501 .8103191 _cons | -.6767714 .2280623 -2.97 0.003 -1.123765 -.2297775 -------------+---------------------------------------------------------------- eq2 | _cons | 230.9255 . . . . . -------------+---------------------------------------------------------------- eq3 | _cons | -6.381776 . . . . . -------------+---------------------------------------------------------------- eq4 | _cons | -.0518383 .0017536 -29.56 0.000 -.0552753 -.0484012 ------------------------------------------------------------------------------ Warning: convergence not achieved . estat ic ----------------------------------------------------------------------------- Model | Obs ll(null) ll(model) df AIC BIC -------------+--------------------------------------------------------------- . | 346 . -4465.104 5 8940.208 8959.441 ----------------------------------------------------------------------------- Note: N=Obs used in calculating BIC; see [R] BIC note . estimates store PP3 . . ****** RESULTS . . *** TABLE 11.1: PARAMETRIC MODEL ESTIMATES . . estimates table POISSON NB2 FM2NB2 NB2P1 PP1, /// > b(%7.4f) se(%7.3f) stats(N ll) stfmt(%9.1f) -------------------------------------------------------------------------- Variable | POISSON NB2 FM2NB2 NB2P1 PP1 -------------+------------------------------------------------------------ PAT | LOGRandD | 0.4754 0.8159 | 0.064 0.083 LOGK | 0.2708 0.1146 | 0.055 0.069 SCISECT | 0.4672 -0.0898 | 0.155 0.138 _cons | -0.3403 -0.8106 | 0.190 0.180 -------------+------------------------------------------------------------ lnalpha | _cons | -0.1225 | 0.113 -------------+------------------------------------------------------------ component1 | LOGRandD | 0.9181 | 0.063 LOGK | 0.0653 | 0.057 SCISECT | 0.1029 | 0.142 _cons | -1.2738 | 0.161 -------------+------------------------------------------------------------ component2 | LOGRandD | 0.3520 | 0.090 LOGK | 0.2771 | 0.062 SCISECT | -0.5732 | 0.346 _cons | 1.5164 | 0.191 -------------+------------------------------------------------------------ imlogitpi1 | _cons | 2.5712 | 0.338 -------------+------------------------------------------------------------ lnalpha1 | _cons | -0.4288 | 0.122 -------------+------------------------------------------------------------ lnalpha2 | _cons | -2.7760 | 0.753 -------------+------------------------------------------------------------ eq1 | LOGRandD | 0.6113 0.4729 | 0.066 0.059 LOGK | 0.1104 0.2716 | 0.050 0.054 SCISECT | 0.0064 0.4701 | 0.107 0.144 _cons | -0.9352 -0.3293 | 0.199 0.193 -------------+------------------------------------------------------------ eq2 | _cons | 1.5063 -0.0020 | 0.237 0.000 -------------+------------------------------------------------------------ eq3 | _cons | 0.0817 | 0.026 -------------+------------------------------------------------------------ Statistics | N | 346 346 346 346 346 ll | -3365.8 -1127.6 -1106.1 -1116.8 -3371.9 -------------------------------------------------------------------------- legend: b/se . . ********** CLOSE OUTPUT . . * log close . * exit . * clear . . end of do-file . exit, clear