/*Setup the global parameters*/ %LET nimpute = 100; *number of imputations; %LET nboot = 1000; *number of bootstraps; %LET alpha = 0.95; *confidence level; %LET seed = 2010; *random number seed; %LET seed2 = 2011; %LET filename="data.txt"; %LET varname=x m y; %LET missing=99999; /*Read data into sas*/ data dset; infile &filename; input &varname; array nvarlist &varname; do over nvarlist; if nvarlist = &missing then nvarlist = .; end; run; /*Initial imputation to obtain model parameter estimates*/ /*Impute data*/ PROC MI data=dset seed=&seed nimpute=&nimpute out=imputed NOPRINT; var &varname; run; quit; /*Estimate model parameters*/ proc reg data=imputed OUTEST= est NOPRINT; model y = x m; model m = x; by _Imputation_; run; quit; data temp; set est; id =INT((_N_-.1)/2)+1; modelnum = MOD(_N_+1, 2)+1; run; /*Orgnize the model parameters*/ DATA temp1; SET temp; array int[2] iY iM; array xpar[2] c a; array mpar[2] b tmp1; array sigma[2] sy sm; retain a b c iY iM sy sm; BY id; IF FIRST.id THEN DO I = 1 to 2; int[I] = .; xpar[I] = .; mpar[I]=.; sigma[I]=.; END; int[modelnum] = intercept; xpar[modelnum] = x; mpar[modelnum] = m; sigma[modelnum] = _RMSE_; IF LAST.id THEN OUTPUT; KEEP _imputation_ a b c iY iM sy sm; RUN; data temp2; set temp1; ab=a*b; run; /*Combine multiple imputation results*/ proc means data=temp2 NOPRINT; var a b c ab iY iM sy sm; output out=iniest mean(a b c ab iY iM sy sm)=a b c ab iY iM sy sm; run; /*bootstrap data*/ data bootsamp0; do sampnum = 1 to 500; /* To create 1000 bootstrap replications */ do i = 1 to nobs; ran = round(ranuni(&seed) * nobs); set dset nobs = nobs point = ran; output; end; end; stop; run; quit; /*Impute data*/ PROC MI data=bootsamp0 seed=&seed nimpute=&nimpute out=imputed0 NOPRINT; EM maxiter = 500; var &varname; by sampnum; run; quit; /*bootstrap data*/ data bootsamp1; do sampnum = 501 to 1000; /* To create 1000 bootstrap replications */ do i = 1 to nobs; ran = round(ranuni(&seed2) * nobs); set dset nobs = nobs point = ran; output; end; end; stop; run; quit; /*Impute data*/ PROC MI data=bootsamp1 seed=&seed2 nimpute=&nimpute out=imputed1 NOPRINT; EM maxiter = 500; var &varname; by sampnum; run; quit; /*Combine two set of imputations*/ data imputed; set imputed0 imputed1; run; /*Estimate model parameters*/ proc reg data=imputed OUTEST= est NOPRINT; model y = x m; model m = x; by sampnum _Imputation_; run; quit; data temp; set est; id =INT((_N_-.1)/2)+1; modelnum = MOD(_N_+1, 2)+1; run; /*Orgnize the model parameters*/ DATA temp1; SET temp; array int[2] iY iM; array xpar[2] c a; array mpar[2] b tmp1; array sigma[2] sy sm; retain a b c iY iM sy sm; BY id; IF FIRST.id THEN DO I = 1 to 2; int[I] = .; xpar[I] = .; mpar[I]=.; sigma[I]=.; END; int[modelnum] = intercept; xpar[modelnum] = x; mpar[modelnum] = m; sigma[modelnum] = _RMSE_; IF LAST.id THEN OUTPUT; KEEP sampnum _imputation_ a b c ab iY iM sy sm; RUN; data temp2; set temp1; ab=a*b; run; /*Combine multiple imputation results*/ proc means data=temp2 NOPRINT; by sampnum; var a b c ab iY iM sy sm; output out=bootest mean(a b c ab iY iM sy sm)=a b c ab iY iM sy sm; run; /*Calculate the BC intervals*/ proc iml; start main; use bootest; read all into X; use iniest; read all into Y; n=nrow(X); m=ncol(X); *print m,n; bc_lo=J(1,m-3,0); bc_up=J(1,m-3,0); alphas=1-(1-&alpha)/2; zcrit = probit(alphas); do j=1 to m-3; number=0; do i=1 to n; if X[i,j+3]n then up=n; bc_lo[j]=vec[low]; bc_up[j]=vec[up]; end; *print bc_lo, bc_up; res=Y || bc_lo || bc_up; *print res; create bcci from res; append from res; finish; run main; quit; data _NULL_; set bcci; file "bcci.txt"; put col1-col26; run; quit;