/* Methods I, MLE October 2020 */ cd "H:\Teaching\Methods 2020\lectures\Week 6 marginal effects mle\MLE\code" clear set seed 1234567 set obs 100 // --- Logistic regression * Simulate a Bernoulli gen bernie = uniform()<0.4 *list bernie, c sum bernie scalar p = r(mean) * Calculate highest log-likelihood di 100*0.46*ln(0.46) + (100-100*0.46)*ln(1-0.46) * same as di 100*p*ln(p) + (100-100*p)*ln(1-p) * Verify with logit logit bernie * get in robability scale di 1/(1+exp(.1603427 )) * Plot -log-likelihood twoway function y= -(100*x*ln(x) + (100-100*x)*ln(1-x)), range(0 1) /// xtitle("p") ytitle("-Ln P") saving(l100.gph, replace) graph export l100.png, replace // --- Logistic response function twoway function y=exp(x) / (1+ exp(x)), range(-10 10) saving(l1.gph, replace) twoway function y=exp(-x) / (1+ exp(-x)), range(-10 10) saving(l2.gph, replace) graph combine l1.gph l2.gph, xsize(20) ysize(10) graph export lboth.png, replace // --- Linear model, normal regression with MLE * Normal clear set seed 1234567 set obs 100 gen ynorm = rnormal(100, 10) sum ynorm reg ynorm sysuse auto, clear qui reg price weight mpg ereturn list qui reg price weight mpg * Save sample size and SSE local N = e(N) local rss = e(rss) * Use formula local ll = -0.5*`N'*(ln(2*_pi)+ln(`rss'/`N')+1) display %20.6f `ll' display %20.6f e(ll) // --- The mlexp command mlexp (ln(normalden(price, {xb: weight mpg _cons}, {sigma}))) // --- MLE example use "H:\Teaching\Methods 2018\DNM book\heus\heus_mepssample.dta", clear gen lexp = log(exp_tot +1) reg lexp age female /* mlexp (ln(normalden(lexp,{sigma}))) mlexp (ln(normalden(lexp, {b0=-2000}+{b1=120}*age,{sigma=9000}))), diff */ mlexp (ln(normalden(lexp, {xb: age female _cons} , {sigma}))) * Use a small program capture program drop lfols program lfols args lnf xb lnsigma local y "\$ML_y1" quietly replace `lnf' = ln(normalden(`y', `xb',exp(`lnsigma'))) end ml model lf lfols (xb: lexp = age female) (lnsigma:) ml maximize display exp([lnsigma]_cons) // --- Likelihood function negative twoway function y =log(x), range(-2 2) xline(0 1) yline(0) /// color(red) title("y = log(x)") // --- Likelihood ratio test use "H:\Teaching\Methods 2017\data\GPA1.dta", clear rename colGPA colgpa rename hsGPA hsgpa quietly { reg colgpa est sto m1 reg colgpa hsgpa est sto m2 reg colgpa hsgpa skipped est sto m3 } est table m1 m2 m3, star stat(r2 r2_a ll bic aic) b(%7.3f) * LRT lrtest m3 m2 lrtest m3 m1 * LRT by hand qui reg colgpa est sto m0 scalar ll0 = e(ll) reg colgpa male campus est sto m1 scalar ll1 = e(ll) lrtest m0 m1 * By hand di -2*[ll0 - ll1] // --- Logistic vs normal clear set seed 123456 set obs 5000 gen u = uniform() * Simulate logistic distribution gen l = -ln((1 - u)/u) sum l * Simulated normal with same parameters gen n = rnormal(r(mean), r(sd)) * Plot kdensity l, bw(0.3) gen(xl dl) kdensity n, bw(0.3) gen(xn dn) line dl xl, sort color(red) || line dn xn, sort /// title("Logistic (red) vs normal distribution") ytitle("Density") /// xtitle("x") legend(off) graph export logvsnorm.png, replace // *** New dataset bcuse mroz, clear lowess inlf nwifeinc, gen(lflow) nograph scatter inlf nwifeinc, jitter(5) msize(small) || line lflow nwifeinc, sort /// legend(off) saving(lblow.gph, replace) graph export lblow.png, replace logit inlf nwifeinc, nolog * LRT qui logit inlf nwifeinc, nolog est sto full qui logit inlf, nolog est sto redu lrtest full redu qui logit inlf nwifeinc, nolog scalar ll_cm = e(ll) qui logit inlf, nolog scalar ll_n = e(ll) di 1 - (ll_cm/ll_n) di "cm: " ll_cm " " "null: " ll_n " " "(ll_cm/ll_n): " (ll_cm/ll_n) gen hsp = 0 replace hsp = 1 if educ > 12 & educ ~= . // -- Predictions qui logit inlf hsp, nolog * Predictions for logit manually gen phat_manu = 1/(1+exp(-(_b[_cons] +_b[hsp]*hsp))) *Same as using the inverse logit function gen phat_invl = invlogit(_b[_cons] +_b[hsp]*hsp) * Same as default of predict command predict phat_pred sum phat_manu phat_invl phat_pred *probit qui probit inlf hsp, nolog * use inverse normal gen phat1_norm = normal(_b[_cons] + _b[hsp]*hsp) predict phat1_predprob predict zscore, xb sum phat1* zscore di normal(0.1760512)