/*

   Marginal effects to interpret models 
   M Perraillon
   October 2020

*/



// ---- Analytical standard errors
* Get data 
use https://www.stata-press.com/data/r16/cattaneo2, clear

* Generate interactions so e(V) more compact
gen sm_age = mbsmoke * mage
qui reg bweight mbsmoke mage sm_age
matrix list e(V)

* var(beta2)
di  e(V)[2,2] 
* var (beta3)
di  e(V)[3,3]
* covariance b2,b3
di e(V)[3,2]

* whem smoke is 0 variance of marginal effect is  
di sqrt(e(V)[2,2])
* when smoke is 1
di sqrt(e(V)[2,2] + e(V)[3,3] + 2*e(V)[3,2])


qui reg bweight i.mbsmoke##c.mage
margins, dydx(mage) at(mbsmoke=(0 1)) vsquish
matrix list r(Jacobian)
matrix list r(V)

reg bweight i.mbsmoke##c.mage,  nofvlabel
matrix list e(V)

// ---- Obtain the delta method variance matrix
qui reg bweight i.mbsmoke##c.mage
margins mbsmoke, nofvlabel
* Stored Jacobian from margins
matrix list r(Jacobian)
matrix J = r(Jacobian)
* Stored variance-covariance matrix from reg model 
matrix list e(V)
* Replicate delta method variance (of margins)
matrix Vrep = J*e(V)*J'
matrix list Vrep
* match margins output SE for nonsmoker
di sqrt(Vrep[1,1])
* matches stored variance 
matrix list r(V)


// --- Replicating marginal effects 
bcuse bwght, clear
desc 

* Define low weight (not based on any clinical criteria)
gen lw = 0
replace lw = 1 if bwght < 100 & bwght ~= .

* Logit model 
logit lw cigs faminc motheduc, nolog

* odds ratios
logit, or

* LPM 
reg lw cigs faminc motheduc, robust

* Smooth with lpoly -- fairly linear
lpoly lw cigs, bw(5) deg(1)  
graph export l.png, replace


// --- Average marginal effects (they are all averages, by the way)

* Compute the initial "small change" h
qui sum cigs
scalar h = (abs(r(mean))+.0001)*.0001
di h  

preserve 
   qui logit lw cigs faminc motheduc, nolog
   * as is
   predict double lw_0 if e(sample)   
   * Change cigs by a bit
   replace cigs = cigs + scalar(h)
   predict lw_1 if e(sample)   
   * For each ob calculate dydx
   gen double dydx = (lw_1-lw_0)/scalar(h)  
  * Average
   sum dydx
restore 

* Same with margins
qui logit lw cigs faminc motheduc, nolog
margins, dydx(cigs)

* The Stata way, two-sided derivative 
preserve 
   qui logit lw cigs faminc motheduc
   * Define small change for cigs 
   qui sum cigs
   scalar h = (abs(r(mean))+0.0001)*0.0001

   * Duplicte variable 
   clonevar cigs_c = cigs
   
   * Small negative change 	
   replace cigs = cigs_c - scalar(h) 
   predict double lw_0 if e(sample)
   * Small positive change
   replace cigs = cigs_c + scalar(h)
   predict double lw_1 if e(sample)
   
   gen double dydx = (lw_1-lw_0)/(2*scalar(h))
   sum dydx
restore 

// --- Incremental marginal effects
gen     smoked = 0
replace smoked = 1 if cigs > 0 & cigs ~=.

* Incremental change 
preserve 
   qui logit lw smoked faminc motheduc
   * Nobody smokes
   replace smoked = 0 
   predict double lw_0 if e(sample)
   
   * Everybody smokes
   replace smoked = 1 
   predict double lw_1 if e(sample)
   
   gen double dydx = (lw_1-lw_0)
   sum dydx
restore 

* Wrong syntax
qui logit lw smoked faminc motheduc, nolog
margins, dydx(smoked)

* Right
qui logit lw i.smoked faminc motheduc, nolog
margins, dydx(smoked)


// --- Predicted margins or adjusted means 
qui logit lw i.smoked faminc motheduc, nolog
margins smoked, post 

ereturn list
matrix list e(b)

* marginal effects 
di e(b)[1,2] - e(b)[1,1]

* marginal effects 
di e(b)[1,2] - e(b)[1,1]

* relative risks 
di e(b)[1,2] / e(b)[1,1]

* odds ratios 
di (e(b)[1,2]/(1-e(b)[1,2]) )/ (e(b)[1,1]/(1-e(b)[1,1]))

* Model with no covariates (unadjusted)
logit lw i.smoked, nolog or
margins smoked, post 

// --- Marginal effects at means 
qui sum faminc
scalar fincbar = r(mean)
qui sum motheduc
scalar medbar = r(mean)

preserve
   qui sum cigs
   scalar h = (abs(r(mean))+0.0001)*0.0001
   
   qui logit lw cigs faminc motheduc, nolog
   
   clonevar cigs_c = cigs
   
   * Replace with means 
   replace faminc = scalar(fincbar)
   replace motheduc = scalar(medbar)
   
   * Small negative change
   replace cigs = cigs_c - scalar(h)
   predict double lw_0 if e(sample)

   * Small positive change change
   replace cigs = cigs_c + scalar(h)
   predict double lw_1 if e(sample)

   gen double dydx = (lw_1-lw_0)/(2*scalar(h))
   sum dydx
restore

qui logit lw cigs faminc motheduc, nolog

* Careful with syntax
margins, dydx(cigs)  at((mean) faminc motheduc) 

* Not same as (here cigs is change to its mean before changing by h)
margins, dydx(cigs)  atmeans


// --- Marginsplot
reg bwght i.motheduc faminc parity

margins motheduc, at((mean) faminc parity)
marginsplot
graph export marplot1.png, replace


marginsplot, recast(line) recastci(rarea) 
graph export marplot2.png, replace

* other neat things you can do 
sum bwght
qui reg bwght i.motheduc faminc parity
margins motheduc, at((mean) faminc parity)
marginsplot, horizontal recast(scatter) xline(119, lcolor(red)) ///
   xscale(range()) yscale(reverse)
graph export mhor.png, replace


// --- Interactions 
gen     inc = 0
replace inc = 1 if faminc > 40 & faminc~= .

logit lw c.cigs##i.inc, nolog
margins, dydx(*)

margins, dydx(cigs) at(inc=(0 1)) vsquish
marginsplot, recastci(rarea)
graph export inter.png, replace

* Test the interaction in the probability scale 
* using delta-method SEs
logit lw c.cigs##i.inc, nolog
margins, dydx(cigs) at(inc=(0 1)) contrast(at)




// --- replicating Karaka-Mandic paper
preserve
   webuse margex, clear
   gen female = (sex==1)

   logit outcome c.age##i.female, nolog
   margins, dydx(*)
   margins, dydx(age) at(female=(0 1)) post
   
   
   logit outcome c.age##i.female, nolog
   margins, dydx(age) at(female=(0 1))
   marginsplot 
   
   margins, dydx(age) at(female=(0 1)) contrast(at)
restore


// -- GLM 
qui glm lw c.cigs##i.inc, family(binomial) link(logit)
margins, dydx(*)
margins, dydx(cigs) at(inc=(0 1)) vsquish

// -- Two-part models

* Get data from example in twopm command 
webuse womenwk, clear

replace wage = 0 if wage==.

* First part logit, second part GLM with Gamma family
qui twopm wage i.married children, firstpart(logit) secondpart(glm, family(gamma) link(log))
margins, dydx(married)

gen     nozeroc = 0 
replace nozeroc = 1 if wage >0 
preserve  
   clonevar marr  = married
   * not married
   qui logit nozeroc i.married children
   replace married = 0
   predict double fp0 
   replace married = marr 
   qui glm wage i.married children if wage > 0, f(gamma) l(log)
   replace married = 0
   predict double c0
   gen chat0 = fp0*c0
     
   * married
   replace married = marr
   qui logit nozeroc i.married children
   replace married = 1
   predict double fp1
   replace married = marr 
   qui glm wage i.married children if wage > 0, f(gamma) l(log)
   replace married = 1
   predict double c1
   gen chat1 = fp1*c1
   
   gen dif = chat1 - chat0
   sum dif 
restore 


// --- Subpopulations (the "if" version in margins)
bcuse bwght, clear
gen     smoked = 0
replace smoked = 1 if cigs > 0 & cigs ~=.

qui reg bwght i.smoked##(c.fatheduc c.faminc), robust
*ATET
margins r.smoked, subpop(smoked) vce(unconditional)
teffects ra (bwght fatheduc faminc)(smoked), atet