%let maximum = 1000;

%macro parameter_estimate(percent,index);

ods listing close;

proc nlmixed data = mcar&percent&index cov;

parms b0 = -0.0645 b_group = -0.1433 b_diabbase = -0.04 b_hdbase = 0.1224 b_age = -0.0066

b_base_bpcontrolled = 1.1487 b_sex = 0.0873 s2u = 0.5;

eta = b0 + b_group*group + b_diabbase*diabbase + b_hdbase*hdbase + b_age*age

+ b_base_bpcontrolled*base_bpcontrolled + b_sex*sex + u;

expeta = exp(eta);

p = expeta/(1+expeta);

model outcome ~ binary(p);

random u ~ normal(0,s2u) subject = assfpid;

ods output ParameterEstimates = parameter&percent&index

CovMatParmEst = covariance&percent&index;

run;

data parameter&percent&index;

set parameter&percent&index;

keep estimate;

run;

data covariance&percent&index;

set covariance&percent&index;

drop row parameter;

run;

%mend parameter_estimate;

%macro mvn(percent, index, n);

/* arguments for the macro:

1. varcov: data set for variance-covariance matrix

2. means: data set for mean vector

3. n: sample size

4. myMVN: output data set name */

proc iml;

use covariance&percent&index;/* read in data for variance-covariance matrix */

read all into sigma;

use parameter&percent&index;/* read in data for means */

read all into mu;

p = nrow(sigma);/* calculate number of variables */

n = &n;

l = t(half(sigma));/* calculate cholesky root of cov matrix */

z = normal(j(p,&n,1234));/* generate nvars*samplesize normals */

y = l*z;/* premultiply by cholesky root */

yall = t(repeat(mu,1,&n)+y);/* add in the means */

varnames = { b0 b_group b_diabbase b_hdbase b_age b_base_bpcontrolled b_sex s2u};

create myMVN&percent&index from yall (|colname = varnames|);

append from yall;

quit;

%mend mvn;

%macro mi_random_effect(percent, index);

%parameter_estimate(&percent, &index);

%mvn(&percent, &index, 5);

proc iml symsize = 512;

use mymvn&percent&index;

read all into mvndata;

use mcar&percent&index;

read all var {ptid DIABBASE HDBASE base_bpcontrolled last_bpimproved sex age assfpid_num group missing outcome} into temp_data;

log_icca_cov = j(7700,12,0);

do i = 0 to 4;

do j = 1 to 1540;

do k = 1 to 11;

log_icca_cov[i*1540+j,k] = temp_data[j,k];

end;

log_icca_cov[i*1540+j,12] = i+1;

end;

end;

do i = 1 to 7700;

if log_icca_cov[i, 11] = . then do;

num = log_icca_cov[i, 12];

logit_p = mvndata[num, 1] + mvndata[num, 2]*log_icca_cov[i, 9]

+ mvndata[num, 3]*log_icca_cov[i, 2] + mvndata[num, 4]*log_icca_cov[i,3]

+ mvndata[num, 5]*log_icca_cov[i, 7] + mvndata[num, 6]*log_icca_cov[i,4]

+ mvndata[num, 7]*log_icca_cov[i, 6] + rand('NORMAL', 0, sqrt(mvndata[num, 8]));

log_icca_cov[i, 11] = rand('BERNOULLI', exp(logit_p)/(1+exp(logit_p)));

end;

end;

varnames = {ptid DIABBASE HDBASE base_bpcontrolled last_bpimproved sex age assfpid_num group missing outcome _imputation_};

create log_icca_cov&percent&index from log_icca_cov (|colname = varnames|);

append from log_icca_cov;

quit;

%mend mi_random_effect;

%macro mi_icca_log(percent, index);

ods listing close;

%mi_random_effect(&percent, &index);

data log_icca_cov&percent&index;

set log_icca_cov&percent&index;

if outcome > = 1 then outcome = 1;

else if outcome < 1 then outcome = 0;

run;

proc freq data = log_icca_cov&percent&index;

table last_bpimproved*outcome/kappa;

ods output SimpleKappa = log_icca_kappapool&percent&index;

run;

data log_icca_kappapool&percent&index;

set log_icca_kappapool&percent&index;

if Label1 = 'Kappa';

run;

proc sort data = log_icca_cov&percent&index;

by _imputation_;

run;

proc genmod data = log_icca_cov&percent&index;

class outcome assfpid_num;

model outcome = group diabbase hdbase age base_bpcontrolled sex/D = B link = logit;

repeated subject = assfpid_num/type = exch covb;

by _imputation_;

ods output GEEEmpPEst = log_icca_geepar&percent&index

GEERCov = log_icca_geecov&percent&index;

run;

data log_icca_geepar&percent&index;

set log_icca_geepar&percent&index;

if Parameter~ = 'Scale';

if Parm = 'Prm' then Parm = 'Prm1';

else if Parm = 'GROUP' then Parm = 'Prm2';

else if Parm = 'DIABBASE' then Parm = 'Prm3';

else if Parm = 'HDBASE' then Parm = 'Prm4';

else if Parm = 'AGE' then Parm = 'Prm5';

else if Parm = 'BASE_BPCONTROLLED' then Parm = 'Prm6';

else if Parm = 'SEX' then Parm = 'Prm7';

run;

proc mianalyze parms = log_icca_geepar&percent&index covb = log_icca_geecov&percent&index;

modeleffects Prm2;

ods output ParameterEstimates = pool_log_icca_gee&percent&index;

run;

proc nlmixed data = log_icca_cov&percent&index cov;

by _imputation_;

parms b0 = -0.0645 b_group = -0.1433 b_diabbase = -0.04 b_hdbase = 0.1224 b_age = -0.0066

b_base_bpcontrolled = 1.1487 b_sex = 0.0873 s2u = 0.5;

eta = b0 + b_group*group + b_diabbase*diabbase + b_hdbase*hdbase + b_age*age

+ b_base_bpcontrolled*base_bpcontrolled + b_sex*sex + u;

expeta = exp(eta);

p = expeta/(1+expeta);

model outcome ~ binary(p);

random u ~ normal(0,s2u) subject = assfpid_num;

ods output ParameterEstimates = log_icca_repar&percent&index

CovMatParmEst = log_icca_recov&percent&index;

run;

proc mianalyze parms = log_icca_repar&percent&index covb = log_icca_recov&percent&index;

modeleffects b_group;

ods output ParameterEstimates = pool_log_icca_re&percent&index;

run;

ods listing;

%mend mi_icca_log;

%macro append_log_icca(percent);

%do index = 1%to &maximum;

%if &index = 1%then%do;

data pool_log_icca_re&percent;

set pool_log_icca_re&percent&index;

run;

data pool_log_icca_gee&percent;

set pool_log_icca_gee&percent&index;

run;

data log_icca_kappa&percent;

set log_icca_kappapool&percent&index;

run;

%end;

%else%do;

proc append base = pool_log_icca_re&percent data = pool_log_icca_re&percent&index;

run;

proc append base = pool_log_icca_gee&percent data = pool_log_icca_gee&percent&index;

run;

proc append base = log_icca_kappa&percent data = log_icca_kappapool&percent&index;

run;

%end;

%end;

%mend append_log_icca;

%macro collect_result_log_icca(percent);

%do index = 1%to &maximum;

%mi_icca_log(&percent,&index);

%end;

%append_log_icca(&percent);

proc univariate data = log_icca_kappa&percent;

var nValue1;

run;

proc univariate data = pool_log_icca_gee&percent;

var Estimate StdErr;

run;

proc univariate data = pool_log_icca_re&percent;

var Estimate StdErr;

run;

%mend collect_result_log_icca;

filename junk dummy;

proc printto log = junk;run;

%collect_result_log_icca(05);

%collect_result_log_icca(10);

%collect_result_log_icca(15);

%collect_result_log_icca(30);

%collect_result_log_icca(50);

proc printto; run;