------------------------------------------------------------------------------------------------------------------------------- name: log: c:\Imbook\courses\2013_bgpe_germany\day4asymptotic_nonlinear\ct_nonlinear.txt log type: text opened on: 7 Jul 2013, 10:48:40 . . * STATA Program . * For A. Colin Cameron and Pravin Trivedi "Lectures in Microeconometrics" . * Nonlinear regression (M-estimation) . . * To run you need file . * mus10data.dta . * in your directory . . ********** SETUP ********** . . set more off . version 12.0 . set mem 10m set memory ignored. Memory no longer needs to be set in modern Statas; memory adjustments are performed on the fly automatically. . set scheme s1mono /* Graphics scheme */ . . ********** DATA DESCRIPTION ********** . . * 2002 Medical Expenditure Panel Survey (MEPS) . * U.S. individuals aged 25-64 years working in private sector . * but not self-employed . * and not receiving public insurance (Medicare and Medicaid) . * Data due to Deb, Munkin and Trivedi (2006) . . ******* NONLINEAR REGRESSION - POISSON EXAMPLE . . * Read in data, select 2002 data, describe and summarize key variables . use mus10data.dta, clear . quietly keep if year02==1 . describe docvis private chronic female income storage display value variable name type format label variable label ------------------------------------------------------------------------------------------------------------------------------- docvis int %8.0g number of doctor visits private byte %8.0g = 1 if private insurance chronic byte %8.0g = 1 if a chronic condition female byte %8.0g = 1 if female income float %9.0g Income in $ / 1000 . summarize docvis private chronic female income Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- docvis | 4412 3.957389 7.947601 0 134 private | 4412 .7853581 .4106202 0 1 chronic | 4412 .3263826 .4689423 0 1 female | 4412 .4718948 .4992661 0 1 income | 4412 34.34018 29.03987 -49.999 280.777 . histogram docvis if docvis < 40, width(1) (bin=39, start=0, width=1) . . * Poisson regression with default standard errors . poisson docvis private chronic female income Iteration 0: log likelihood = -18504.413 Iteration 1: log likelihood = -18503.549 Iteration 2: log likelihood = -18503.549 Poisson regression Number of obs = 4412 LR chi2(4) = 8852.71 Prob > chi2 = 0.0000 Log likelihood = -18503.549 Pseudo R2 = 0.1930 ------------------------------------------------------------------------------ docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .7986652 .027719 28.81 0.000 .744337 .8529934 chronic | 1.091865 .0157985 69.11 0.000 1.060901 1.12283 female | .4925481 .0160073 30.77 0.000 .4611744 .5239218 income | .003557 .0002412 14.75 0.000 .0030844 .0040297 _cons | -.2297262 .0287022 -8.00 0.000 -.2859814 -.173471 ------------------------------------------------------------------------------ . . * Poisson regression with robust standard errors . poisson docvis private chronic female income, vce(robust) Iteration 0: log pseudolikelihood = -18504.413 Iteration 1: log pseudolikelihood = -18503.549 Iteration 2: log pseudolikelihood = -18503.549 Poisson regression Number of obs = 4412 Wald chi2(4) = 594.72 Prob > chi2 = 0.0000 Log pseudolikelihood = -18503.549 Pseudo R2 = 0.1930 ------------------------------------------------------------------------------ | Robust docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .7986652 .1090014 7.33 0.000 .5850263 1.012304 chronic | 1.091865 .0559951 19.50 0.000 .9821167 1.201614 female | .4925481 .0585365 8.41 0.000 .3778187 .6072774 income | .003557 .0010825 3.29 0.001 .0014354 .0056787 _cons | -.2297262 .1108732 -2.07 0.038 -.4470338 -.0124186 ------------------------------------------------------------------------------ . . * Comparison of standard errors . quietly poisson docvis private chronic female income . estimates store DEFAULT . quietly poisson docvis private chronic female income, vce(robust) . estimates store ROBUST . estimates table DEFAULT ROBUST, b(%9.4f) se(%9.3f) stats(N r2 F) -------------------------------------- Variable | DEFAULT ROBUST -------------+------------------------ private | 0.7987 0.7987 | 0.028 0.109 chronic | 1.0919 1.0919 | 0.016 0.056 female | 0.4925 0.4925 | 0.016 0.059 income | 0.0036 0.0036 | 0.000 0.001 _cons | -0.2297 -0.2297 | 0.029 0.111 -------------+------------------------ N | 4412 4412 r2 | F | -------------------------------------- legend: b/se . . * Marginal effects AME . margins, dydx(*) Average marginal effects Number of obs = 4412 Model VCE : Robust Expression : Predicted number of events, predict() dy/dx w.r.t. : private chronic female income ------------------------------------------------------------------------------ | Delta-method | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | 3.160629 .4352572 7.26 0.000 2.307541 4.013717 chronic | 4.320935 .2757872 15.67 0.000 3.780402 4.861468 female | 1.949204 .2313726 8.42 0.000 1.495722 2.402686 income | .0140765 .0043457 3.24 0.001 .0055591 .0225939 ------------------------------------------------------------------------------ . * Marginal effects MEM . margins, dydx(*) atmean Conditional marginal effects Number of obs = 4412 Model VCE : Robust Expression : Predicted number of events, predict() dy/dx w.r.t. : private chronic female income at : private = .7853581 (mean) chronic = .3263826 (mean) female = .4718948 (mean) income = 34.34018 (mean) ------------------------------------------------------------------------------ | Delta-method | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | 2.4197 .3057397 7.91 0.000 1.820462 3.018939 chronic | 3.308002 .1794551 18.43 0.000 2.956277 3.659728 female | 1.492263 .1675949 8.90 0.000 1.163783 1.820743 income | .0107766 .0033149 3.25 0.001 .0042796 .0172737 ------------------------------------------------------------------------------ . . * Old commands superceded by margins . * mfx . * margeff . . * Marginal effects using finite difference for binary regressors . quietly poisson docvis i.private i.chronic i.female income, vce(robust) . margins, dydx(*) Average marginal effects Number of obs = 4412 Model VCE : Robust Expression : Predicted number of events, predict() dy/dx w.r.t. : 1.private 1.chronic 1.female income ------------------------------------------------------------------------------ | Delta-method | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- 1.private | 2.404721 .2438573 9.86 0.000 1.926769 2.882672 1.chronic | 4.599174 .2886176 15.94 0.000 4.033494 5.164854 1.female | 1.900212 .2156694 8.81 0.000 1.477508 2.322917 income | .0140765 .0043457 3.24 0.001 .0055591 .0225939 ------------------------------------------------------------------------------ Note: dy/dx for factor levels is the discrete change from the base level. . . * Wald test of nonlinear hypothesis . quietly poisson docvis private chronic female income, vce(robust) . testnl _b[female]/_b[private]=1 (1) _b[female]/_b[private] = 1 chi2(1) = 13.51 Prob > chi2 = 0.0002 . . * Delta method confidence interval . nlcom _b[female]/_b[private] - 1 _nl_1: _b[female]/_b[private] - 1 ------------------------------------------------------------------------------ docvis | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- _nl_1 | -.383286 .1042734 -3.68 0.000 -.587658 -.1789139 ------------------------------------------------------------------------------ . . * Nonlinear least-squares regression (command nl) . generate one = 1 . nl (docvis = exp({xb: private chronic female income one})), vce(robust) nolog (obs = 4412) Nonlinear regression Number of obs = 4412 R-squared = 0.3046 Adj R-squared = 0.3038 Root MSE = 7.407479 Res. dev. = 30185.68 ------------------------------------------------------------------------------ | Robust docvis | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- /xb_private | .7105104 .1086194 6.54 0.000 .4975618 .9234589 /xb_chronic | 1.057318 .0558352 18.94 0.000 .947853 1.166783 /xb_female | .4320225 .0694662 6.22 0.000 .2958338 .5682112 /xb_income | .002558 .0012544 2.04 0.041 .0000988 .0050173 /xb_one | -.040563 .1126216 -0.36 0.719 -.2613578 .1802319 ------------------------------------------------------------------------------ . . * Newton-Raphson in Mata for Poisson MLE . * Set up data and local macros for dependent variable and regressors . generate cons = 1 . local y docvis . local xlist private chronic female income cons . * Mata commands for Poisson MLE NR iterations . mata: ------------------------------------------------- mata (type end to exit) ----------------------------------------------------- : st_view(y=., ., "`y'") // read in stata data to y and X : st_view(X=., ., tokens("`xlist'")) : b = J(cols(X),1,0) // compute starting values : n = rows(X) : iter = 1 // initialize number of iterations : cha = 1 // initialize stopping criterion : do { > mu = exp(X*b) > grad = X'(y-mu) // k x 1 gradient vector > hes = cross(X, mu, X) // negative of the k x k Hessian matrix > bold = b > b = bold + cholinv(hes)*grad > cha = (bold-b)'(bold-b)/(bold'bold) > iter = iter + 1 > } while (cha > 1e-16) // end of iteration loops : mu = exp(X*b) : hes = cross(X, mu, X) : vgrad = cross(X, (y-mu):^2, X) : vb = cholinv(hes)*vgrad*cholinv(hes)*n/(n-cols(X)) : iter // num iterations 13 : cha // stopping criterion 1.11467e-24 : st_matrix("b",b') // pass results from Mata to Stata : st_matrix("V",vb) // pass results from Mata to Stata : end ------------------------------------------------------------------------------------------------------------------------------- . * Present results, nicely formatted using Stata command ereturn . matrix colnames b = `xlist' . matrix colnames V = `xlist' . matrix rownames V = `xlist' . ereturn post b V . ereturn display ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .7986654 .1090509 7.32 0.000 .5849295 1.012401 chronic | 1.091865 .0560205 19.49 0.000 .9820669 1.201663 female | .4925481 .058563 8.41 0.000 .3777666 .6073295 income | .003557 .001083 3.28 0.001 .0014344 .0056796 cons | -.2297263 .1109236 -2.07 0.038 -.4471325 -.0123202 ------------------------------------------------------------------------------ . . ********** CLOSE OUTPUT . * log close . * clear . * exit . . . end of do-file . exit, clear