Gradient field of predicted levels with data overlay

gradient-field.do

// Figure 3. Gradient
use "${git}/constructed/mental_c.dta", clear

  replace hh_logconscap = 8 if hh_logconscap < 8 // For scaling
  replace hh_logconscap = 13 if hh_logconscap >13 &!missing(hh_logconscap) // For scaling

  local fault_controls "hh_epidist hh_slope hh_district_1 hh_district_2 hh_district_3"
  local indiv_controls "indiv_age indiv_male indiv_married indiv_edu_binary"

    reg mh_irt ///
      c.hh_faultdist##c.hh_logconscap  `indiv_controls' `fault_controls' ///
    , cl(village_code)

    margins, at(hh_faultdist=(0(2)60) hh_logconscap=(8(0.2)13)) vce(unconditional) ///
      saving("${git}/outputs/predictions.dta", replace)
    use "${git}/outputs/predictions.dta", clear
      append using "$git/data/analysis_mental.dta"

        replace hh_logconscap = 8 if hh_logconscap < 8 // For scaling
        replace hh_logconscap = 13 if hh_logconscap >13 &!missing(hh_logconscap) // For scaling

    lab var _margin "Predicted Mental Health (Higher = Better)"

    twoway ///
      (contour _margin _at1 _at2    if _at < . ///
        , ccuts(-0.5(0.1)0.5) crule(intensity) ecolor(dkgreen)) ///
      (scatter hh_faultdist hh_logconscap   , ///
        m(O) mc(black) mlw(vvthin) mfc(none) msize(tiny) ///
        jitter(1) jitterseed(123456))  ///
      if hh_faultdist< 60 | hh_faultdist == . ///
    , legend(on order(1 "Sample") ring(0) pos(5)) yscale(reverse) ///
      xtit("(log) Household Consumption per Capita") ///
      ytit("Distance to Activated Fault Line") ///
      ylab(0 "Fault Line" 20 "20km" 40 "40km" 60 "60km") ///
      yscale(noline) xscale(noline)

      graph export "${git}/outputs/F_gradient.png", replace width(4000)