Below is the macro, apologies for the unfortunate formatting.
**Data created from Hanley and Shapiro (1994 Journal of Statistical Education 2(1) );
data Flies;
input ID treatment thorax death40;
cards;
1 0 0.64 1
2 0 0.68 1
3 0 0.68 0
4 0 0.72 0
5 0 0.72 0
6 0 0.76 1
7 0 0.76 0
8 0 0.76 0
9 0 0.76 0
10 0 0.76 0
11 0 0.80 0
12 0 0.80 0
13 0 0.80 0
14 0 0.84 0
15 0 0.84 0
16 0 0.84 0
17 0 0.84 0
18 0 0.84 0
19 0 0.84 0
20 0 0.88 0
21 0 0.88 0
22 0 0.92 0
23 0 0.92 0
24 0 0.92 0
25 0 0.94 0
26 1 0.64 1
27 1 0.64 1
28 1 0.68 1
29 1 0.72 1
30 1 0.72 1
31 1 0.74 1
32 1 0.76 1
33 1 0.76 0
34 1 0.76 0
35 1 0.78 1
36 1 0.80 1
37 1 0.80 1
38 1 0.82 1
39 1 0.82 0
40 1 0.84 1
41 1 0.84 1
42 1 0.84 0
43 1 0.84 0
44 1 0.88 0
45 1 0.88 0
46 1 0.88 0
47 1 0.88 0
48 1 0.88 0
49 1 0.88 0
50 1 0.92 0
;
****************************************************************************
* Note: Always remember to change data set name as well as variable names *
****************************************************************************;
%CounterFactual(data =Flies, /*data set to be analyzed*/
yvar = death40, /*outcome variabe*/
trt = treatment, /* Exposure or treatment variable */
xvar = thorax, /* independent variables to be adjusted */
link = logit, /* link function used by PROC GENMOD: logit, probit, cloglog */
alpha=0.05); /* alpha level for confidence interval, default is 0.05*/
%macro CounterFactual(data=, /*data set to be analyzed*/
yvar =, /*outcome variabe*/
trt =, /* Exposure or treatment variable */
xvar =, /* independent variables to be adjusted */
link =, /* link function used in PROC GENMOD: logit, probit, cloglog */
alpha=0.05 /* alpah level for confidence interval, default is 0.05*/
);
%if &link= probit %then %let dist=’normal’;
%else %if &link=logit %then %let dist=’logistic’;
%else %if &link=cloglog %then %let dist=’extreme’;
proc genmod desc data= &data;
model &yvar = &trt &xvar/covb dist=bin link= &link;
ods output ParameterEstimates=est(keep=estimate)
covB=cov(drop=rowname);
proc iml;
start CI4P(Pest, Var, n);
%let z=probit(1-&alpha/2);
*Wald;
W_L = pest - &z*sqrt(var) ;
W_U = pest + &z*sqrt(var) ;
*log, usd by Flanders and Rhodes (1987 J Chron Dis 40: 697-704);
lOG_L = pest*exp(- &z*sqrt(var)/(pest));
log_U = pest*exp(+ &Z*sqrt(var)/(pest));
*logit;
varlgt = var/( pest*(1-pest))**2 ;
l=log(pest/(1-pest)) - &z*sqrt(varlgt);
lgt_L = exp(l)/(1+exp(l));
u=log(pest/(1-pest)) + &z*sqrt(varlgt);
lgt_U = exp(u)/(1+exp(u));
*Wilson, based on Newcombe (2001 Am Stat 55: 200-202);
ll=log(pest/(1-pest)) - 2*arsinh( &z/2 * sqrt(varlgt));
Wln_L = exp(ll)/(1+exp(ll));
uu=log(pest/(1-pest)) + 2*arsinh( &z/2 * sqrt(varlgt));
Wln_U = exp(uu)/(1+exp(uu));
return(pest||W_L||W_U||log_L||log_U||lgt_L||lgt_U||wln_L||wln_U);
finish CI4P;
start MOVER(p1, l1, u1, p2, l2, u2, corr);
**Baded on Zou (2008 Am J Epidemiol 162: 212-224);
point = p1-p2;
L = p1- p2 - sqrt(max(0, (p1-l1)**2 -2*corr*(p1-l1)*(u2-p2) + (u2-p2)**2));
U = p1- p2 + sqrt(max(0, (u1-p1)**2 -2*corr*(u1-p1)*(p2-l2) + (p2-l2)**2));
return(point||L||U);
finish MOVER;
**Bring in data for prediction;
use &data;
read all var{&xvar} into X;
n = nrow(X);
m = ncol(X);
X1 = J(n,2,1)||X;
X0 = J(n,1,1)||J(n,1,0)||X;
use cov;
read all var _num_ into V;
use est;
read all var _num_ into beta;
beta=beta[1:(m+2)];
if &dist = ’extreme’ then do;
p_1 = 1-exp(-exp(X1 * beta));
p_0 = 1-exp(-exp(X0 * beta));
piece1 = exp(X1 * beta -exp(X1 * beta) )# X1;
piece0 = exp(X0 * beta -exp(X0 * beta) )# X0;
end;
else do;
p_1 = cdf(&dist, X1 * beta); **predicted prob, if exposed;
p_0 = cdf(&dist, X0 * beta); **predicted prob, if unexposed;
piece1 = pdf(&dist, X1 * beta)# X1;
piece0 = pdf(&dist, X0 * beta)# X0;
end;
p1est = sum(p_1)/n;
p0est = sum(p_0)/n;
V1=0; V0=0; COV =0;
do i =1 to n;
do j=1 to n;
Xi1 = X1[i,]; Xj1 = X1[j,];
Xi0 = X0[i,]; Xj0 = X0[j,];
if &dist = ’extreme’ then do;
V1 = V1 + 1/N**2 * exp(Xi1 * beta - exp(Xi1 * beta) )*
exp(Xj1 * beta - exp(Xj1 * beta) )* (xi1*V*T(xj1));
V0 = V0 + 1/N**2 * exp(Xi0 * beta - exp(Xi0 * beta) )*
exp(Xj0 * beta - exp(Xj0 * beta) )* (xi0*V*T(xj0));
COV = cov + 1/N**2 * exp(Xi1 * beta - exp(Xi1 * beta) )*
exp(Xj0 * beta - exp(Xj0 * beta) )* (xi1*V*T(xj0));
end;
else do;
V1 = V1 + 1/N**2*pdf(&dist, xi1* beta) *
pdf(&dist, xj1* beta) * (xi1*V*T(xj1));
V0 = V0 + 1/N**2*pdf(&dist, xi0* beta) *
pdf(&dist, xj0* beta) * (xi0*V*T(xj0));
COV = COV + 1/N**2* pdf(&dist, xi1* beta)*
pdf(&dist, xj0* beta) * (xi1*V*T(xj0));
end;
end; *END J;
end; *END I;
ll = p1est-p0est - &z*sqrt(v1+v0-2*cov);
uu = p1est-p0est + &z*sqrt(v1+v0-2*cov);
group1 = CI4P(p1est, V1, n); **CI for exposed risk;
group0 = CI4P(p0est, V0, n); **CI for unexposed risk;
print ’==== CI for risks ===’;
name = {estimate Wald_L wald_U log_L log_U lgt_L lgt_U Wlsn_L Wlsn_U};
print group1[colname=name],
group0[colname=name];
**Confidenc einterval for difference;
rho = COV/sqrt(V1*V0);
dWald = mover(group1[1], group1[2], group1[3], group0[1], group0[2], group0[3], rho);
dlgt = mover(group1[1], group1[6], group1[7], group0[1], group0[6], group0[7], rho);
dWln = mover(group1[1], group1[8], group1[9], group0[1], group0[8], group0[9], rho);
**The following is based on Zou and Donner (2004 Controlled Clin Trials 25: 3-12);
vd = V1 + V0 - 2*cov;
d = p1est - p0est;
F_LL = log( (1+ d)/(1-d) ) - &z *2*sqrt(vd)/(1- d**2);
F_L = (exp(F_LL) - 1)/(exp(F_LL) + 1);
F_Uu = log( (1+ d)/(1-d) ) + &z *2* sqrt(vd)/(1-d**2);
F_U = (exp(F_uu) - 1)/(exp(F_uu) + 1);
Fisher = d||F_L||F_U;
print ’===CI for difference==’;
print dwald, dlgt, dwln, Fisher;
** CI for Ratio;
**log-delta;
RRest = p1est/p0est;
var = V1/p1est**2 + V0/p0est**2 - 2 * cov/(p1est*p0est);
Lower = RRest*exp(-&z*sqrt(var));
Upper = RRest*exp(+&z*sqrt(var));
RlogDelta = RRest||lower||upper;
* Zou and Donner (2008 Stat Med 27: 1693-1702);
Rlgt = exp(mover(log(group1[1]), log(group1[6]), log(group1[7]),
log(group0[1]), log(group0[6]), log(group0[7]), rho));
RWln = exp(mover(log(group1[1]), log(group1[8]), log(group1[9]),
log(group0[1]), log(group0[8]), log(group0[9]), rho));
print ’=== CI for RR ===’;
print RlogDelta, Rlgt, Rwln;
quit;
%mend CounterFactual;
