tikistitch

proc iml;
** create beta and the covariance matrix for beta;
use _parmest_;
read all var _num_ into beta where(_type_="PARMS");
read all var _num_ into covbeta where(_type_="COV");
close _parmest_;
beta=t(beta);

** create design matrix with intercept and covariates ;
use _design_ nobs _rcnobs_;
read all into _x_;
close _design_;

** compute predictions and variance estimates
for when the risk factor is present;
_x1_ = _x_ || j(_rcnobs_,1,1);
_p1_ = 1
/ (1 + exp(- _x1_*beta) );
avgpred1 = sum(_p1_)/_rcnobs_;
deriv1 = ( _p1_ # ( 1 - _p1_ ) # _x1_ / _rcnobs_ );
var_avg1 = sum( deriv1 * covbeta * t(deriv1) );
free _x1_ _p1_ ;

** compute predictions and variance estimates
for when the risk factor is absent;
_x0_ = _x_ || j(_rcnobs_,1,0);
_p0_ = 1
/ (1 + exp(- _x0_*beta) );
avgpred0 = sum(_p0_)/_rcnobs_;
deriv0 = ( _p0_ # ( 1 - _p0_ ) # _x0_ / _rcnobs_ );
var_avg0 = sum( deriv0 * covbeta * t(deriv0) );
free _x0_ _p0_;

_rcycdif_ = avgpred1 - avgpred0 ;
deriv = deriv1 - deriv0 ;
var_dif = sum( deriv * covbeta * t(deriv) );

se_avg1 = sqrt( var_avg1 );
se_avg0 = sqrt( var_avg0 );
se_dif = sqrt( var_dif );

** Report covariate adjusted estimates
for when risk factor is present and
for when it is absent. ;
print avgpred1 se_avg1;
print avgpred0 se_avg0;
** Now present the difference between these means and a standard error for the difference.;
print _rcycdif_ se_dif;
quit;
run;