Yet Another Blog in Statistical Computing

I can calculate the motion of heavenly bodies but not the madness of people. -Isaac Newton

Calculating Marginal Effects in Zero-Inflated Beta Model

libname data 'c:\projects\sgf14\data';

ods output parameterestimates = _parms;
proc nlmixed data = data.full tech = trureg alpha = 0.01;
  parms a0 = 0  a1 = 0  a2 = 0  a3 = 0  a4 = 0  a5 = 0  a6 = 0  a7 = 0
        b0 = 0  b1 = 0  b2 = 0  b3 = 0  b4 = 0  b5 = 0  b6 = 0  b7 = 0 
        c0 = 1  c1 = 0  c2 = 0  c3 = 0  c4 = 0  c5 = 0  c6 = 0  c7 = 0;
  xa = a0 + a1 * x1 + a2 * x2 + a3 * x3 + a4 * x4 + a5 * x5 + a6 * x6 + a7 * x7;
  xb = b0 + b1 * x1 + b2 * x2 + b3 * x3 + b4 * x4 + b5 * x5 + b6 * x6 + b7 * x7;
  xc = c0 + c1 * x1 + c2 * x2 + c3 * x3 + c4 * x4 + c5 * x5 + c6 * x6 + c7 * x7;
  mu_xa = 1 / (1 + exp(-xa));
  mu_xb = 1 / (1 + exp(-xb));
  phi = exp(xc);
  w = mu_xb * phi;
  t = (1 - mu_xb) * phi;
  if y = 0 then lh = 1 - mu_xa;
  else lh = mu_xa * (gamma(w + t) / (gamma(w) * gamma(t)) * (y ** (w - 1)) * ((1 - y) ** (t - 1)));
  ll = log(lh);
  model y ~ general(ll);
  *** calculate components for marginal effects ***;
  _mfx_a = exp(xa) / ((1 + exp(xa)) ** 2);
  _mfx_b = exp(xb) / ((1 + exp(xb)) ** 2);
  _p_a = 1 / (1 + exp(-xa));
  _p_b = 1 / (1 + exp(-xb));
  predict _mfx_a out = _marg1 (rename = (pred = _mfx_a) keep = id pred y);
  predict _p_a   out = _marg2 (rename = (pred = _p_a)   keep = id pred);
  predict _mfx_b out = _marg3 (rename = (pred = _mfx_b) keep = id pred);
  predict _p_b   out = _marg4 (rename = (pred = _p_b)   keep = id pred);
run;

data _null_;
  set _parms;
  call symput(parameter, put(estimate, 20.15));
run;

data _marg5;
  merge _marg1 _marg2 _marg3 _marg4;
  by id;

  _marg_x1 = _mfx_a * &a1 * _p_b + _mfx_b * &b1 * _p_a;
  _marg_x2 = _mfx_a * &a2 * _p_b + _mfx_b * &b2 * _p_a;
  _marg_x3 = _mfx_a * &a3 * _p_b + _mfx_b * &b3 * _p_a;
  _marg_x4 = _mfx_a * &a4 * _p_b + _mfx_b * &b4 * _p_a;
  _marg_x5 = _mfx_a * &a5 * _p_b + _mfx_b * &b5 * _p_a;
  _marg_x6 = _mfx_a * &a6 * _p_b + _mfx_b * &b6 * _p_a;
  _marg_x7 = _mfx_a * &a7 * _p_b + _mfx_b * &b7 * _p_a;
run;

proc means data = _marg5 mean;
  var _marg_x:;
run;
/*
Variable            Mean
------------------------
_marg_x1      -0.0037445
_marg_x2       0.0783118
_marg_x3       0.0261884
_marg_x4      -0.3105482
_marg_x5     0.000156693
_marg_x6    -0.000430756
_marg_x7      -0.0977589
------------------------
*/
Advertisements

Written by statcompute

October 12, 2013 at 4:17 pm

%d bloggers like this: