version 8.2
#delimit ;
clear, replace;
program drop _all;
set mem 100M;

infix str10  county    1- 9
      str2   state    10-11
      long   pop      12-25
      float  pct6574  26-31
      float  pct75_   33-40
      int    hosps    42-47
      long   beds     48-55
      double income   56-63
      using "c:\misc\tmp\classes\countyhosp.prn";

gen anyhosp = hosps>0 & hosps<.;
gen logpop = log(pop);
gen loginc = log(income);

des;
sum;

logit anyhosp logpop pct75_ loginc;
mfx compute;
predict phat, p;
matrix Var = e(V);
matrix list Var;

gen me_inc = phat*(1-phat)*_b[loginc];

sum me_inc;

gen dmedb1   = (1-2*phat)*phat*(1-phat)*_b[loginc];
gen dmedbinc = dmedb1*loginc+phat*(1-phat);
gen dmedbpct = dmedb1*pct75_;
gen dmedbpop = dmedb1*logpop;

sum dmedb1 dmedbpop dmedbpct dmedbinc;

matrix grad = J(1,4,0);
matrix list grad;
quietly sum dmedb1;   matrix grad[1,1]=r(mean);
quietly sum dmedbpop; matrix grad[1,2]=r(mean);
quietly sum dmedbpct; matrix grad[1,3]=r(mean);
quietly sum dmedbinc; matrix grad[1,4]=r(mean);

matrix list grad;

matrix v_me_inc = grad*Var*(grad');

matrix list v_me_inc;

display sqrt(v_me_inc[1,1]);

bootstrap "logit anyhosp logpop pct75_ loginc" _b, reps(100) bca;

program me_income, rclass;
  version 8.2;
  logit anyhosp logpop pct75_ loginc;
  predict phat_temp, p;
  gen me_inc_temp = phat_temp*(1-phat_temp)*_b[loginc];
  quietly sum me_inc_temp;
  local meinc = r(mean);
  drop phat_temp me_inc_temp;
  return scalar marg_eff = `meinc';
end;

bootstrap "me_income" r(marg_eff), reps(100);