帮忙看看这个程序怎么啦！-经管之家官网！

以下程序在运行的时候总是出现：
NOTE: IML Ready
ERROR: (execution) Invalid subscript or subscript out of range.
operation : [ at line 14952 column 114
operands: Z_VALUE1, J, *LIT1067
Z_VALUE1 1 row 3 cols (numeric)
1.1337026 1.27154 1.1668986
J 1 row 1 col (numeric)
2
*LIT1067 1 row 1 col (numeric)
1
statement : ASSIGN at line 14952 column 61
ERROR: (execution) Matrices do not conform to the operation.
有高手帮忙看看啦！
%macro stratCI(data=, treatment=, strata=, response=, level=);
%if %sysfunc(exist(&data))=0 %then
%do;
%put " ERROR: Input dataset not specified"; %goto exit;
%end;
%if %length(&treatment) = 0 %then
%do;
%put " ERROR: Treatment group must be specified";
%goto exit;
%end;
%if %length(&strata) = 0 %then
%do;
%put " ERROR: Strata must be specified";
%goto exit;
%end;
%if %length(&response) = 0 %then
%do; %put " ERROR: Response must be specified";
%goto exit;
%end;
proc sql noprint;
create table temp
as select
&treatment as treatment,
&strata as strata,
&response as response
from &data;
quit;
%let alpha= 1-&level/2;
proc sort data=temp;
by strata treatment;
run;
proc freq data=temp noprint ;
by strata;
tables treatment/riskdiff nopercent nocol norow out=freq_strata;
run;
proc freq data=temp;
by strata;
tables treatment*response/riskdiff nopercent nocol out=freq;
ods output RiskDiffCol2 = riskstrata;
run;
* Total Frequency;
proc transpose data=freq_strata out=freq_strata1 (drop=_NAME_ _LABEL_) prefix=n_trt;
by strata;
var count;
run;
* Proportion;
proc transpose data=riskstrata out=riskstrata1 (drop=_NAME_) prefix=prop;
where Row in ('第 1 行','第 2 行');
by strata;
var Risk;
run;
* Risk difference;
data riskstrata2;
set riskstrata(rename=(risk=r_diff));
if Row = '差值';
keep r_diff ASE strata;
run;
* Merge;
data ready;
merge freq_strata1 riskstrata1 riskstrata2;
by strata;
run;
*STRATA-ADJUSTED RISK DIFFERENCE;
data ready1;
set ready;
*calculate weight;
CMH_part = (n_trt1*n_trt2)/(n_trt1+n_trt2);/*Cochran-Mantel-Haenszel weight*/
invar_part =1/ASE**2; /*inverse of variance weight*/
min_part1= (n_trt1+n_trt2);
min_part2= min_part1*r_diff;
min_part3= invar_part*r_diff;
run;
proc sql;
create table ready2 as select strata,
sum(CMH_part) as CMH_sum,
sum(invar_part) as invar_sum,
sum(min_part1) as min_sum1,
sum(min_part2) as min_sum2,
sum(min_part3) as min_sum3
from ready1;
quit;
data ready3;
merge ready1 ready2;
by strata;
min_alpha = r_diff*invar_sum - min_sum3;
min_beta = invar_part * (1 + min_alpha*min_sum2/min_sum1);
min_part4=r_diff*min_beta;
min_part5 = r_diff*min_alpha*invar_part; run;
proc sql;
create table ready3_1 as select *,
sum(min_part4) as min_sum4,
sum(min_part5) as min_sum5 from ready3;
quit;
data weights;
set ready3_1;
CMH_weight = CMH_part/CMH_sum;
INVAR_weight = invar_part/invar_sum;
MIN_weight = min_beta/invar_sum-(min_alpha*invar_part/(invar_sum + min_sum5))*(min_sum4/invar_sum);
run;
/* all weights are calculated. */
proc iml;
use weights;
read all;
read all var {CMH_weight invar_weight min_weight} into w;
n=n_trt1||n_trt2;
n_strata=nrow(n); /*表示组别的数*/
totaln=sum(n);
p= prop1||prop2;
/* unadjusted CI :Pool */
p1_1 = n_trt1 # prop1;
p2_1 = n_trt2 # prop2;
n1=n_trt1[+,];
n2=n_trt2[+,];
p1=p1_1[+,]/n1;
p2=p2_1[+,]/n2;
Crude_Lower = (p1-p2) - probit(&alpha)# sqrt(p1 # (1-p1)/n1 + p2 # (1-p2)/n2);
Crude_Upper = (p1-p2) + probit(&alpha)# sqrt(p1 # (1-p1)/n1 + p2 # (1-p2)/n2);
/* start to calculate Wald CI */
adj_rdiff= w # r_diff;
wald_part1= adj_rdiff[+,];
wald_var = 1/invar_part;
sq_w = w ## 2;
adj_var = sq_w # wald_var;
wald_part2 = adj_var[+,];
Wald_L = wald_part1 - probit(&alpha)*sqrt(wald_part2);
Wald_U = wald_part1 + probit(&alpha) *sqrt(wald_part2);
Wald_Lower=Wald_L`;
Wald_Upper=Wald_U`;
/* end of Wald CI */
* 1. calculate z_gamma;/*表示两组，三个分层因素*/
/*partz_a1 = J(n_strata,3,1);
partz_a2 = J(n_strata,3,1);
partz_b1 = J(n_strata,3,1);
partz_b2 = J(n_strata,3,1); */
partz_a1 = J(n_strata,3,1); /*表示两组，三个分层因素*/
partz_a2 = J(n_strata,3,1);
partz_b1 = J(n_strata,3,1);
partz_b2 = J(n_strata,3,1);
do i=1 to n_strata;
do j=1 to 3;
partz_a1[i,j] = (sq_w[i,j]/n_trt1[i,1]) # (p[i,1] # (1-p[i,1]));
partz_a2[i,j] = (sq_w[i,j]/n_trt2[i,1]) # (p[i,2] # (1-p[i,2]));
partz_b1[i,j] = w[i,j] # sqrt(p[i,1] # (1-p[i,1])/n_trt1[i,1]);
partz_b2[i,j] = w[i,j] # sqrt(p[i,2] # (1-p[i,2])/n_trt2[i,1]);
end;
end;
z_value1=(probit(&alpha)*sqrt(partz_a1[+,])/partz_b1[+,]);
z_value2=(probit(&alpha)*sqrt(partz_a2[+,])/partz_b2[+,]);
PRINT z_value1 z_value2;
* 2. calculate L1 L2 U1 U2 based on stratified Wilson CI;
/*wilson_part1a= J(n_strata,3,1);
wilson_part1b= J(n_strata,3,1);
wilson_part2a= J(n_strata,3,1);
wilson_part2b= J(n_strata,3,1);
*/
wilson_part1a= J(n_strata,3,1);
wilson_part1b= J(n_strata,3,1);
wilson_part2a= J(n_strata,3,1);
wilson_part2b= J(n_strata,3,1);
do i=1 to n_strata;
do j=1 to 3; /*表示有多少分层*/
wilson_part1a[i,j]=(n[i,1] # p[i,1] + 0.5 # (Z_value1[j,1] ## 2))/(n[i,1] + Z_value1[j,1] ## 2);
wilson_part1b[i,j]=((n[i,1] # Z_value1[j,1])/(n[i,1]+Z_value1[j,1]##2)) # (sqrt(4 # n[i,1] # p[i,1] # (1-p[i,1]) + Z_value1[j,1] ##2)/(2 # n[i,1]) );
wilson_part2a[i,j]=(n[i,2] # p[i,2] + 0.5 # (Z_value2[j,1] ## 2))/(n[i,2] + Z_value2[j,1] ## 2);
wilson_part2b[i,j]=((n[i,2] # Z_value2[j,1])/(n[i,2]+Z_value2[j,1]##2)) # (sqrt(4 # n[i,2] # p[i,2] # (1-p[i,2]) + Z_value2[j,1] ##2)/(2 # n[i,2]) );
end;
end;
wilson_L1 = w #( wilson_part1a - wilson_part1b); wilson_U1 = w #( wilson_part1a + wilson_part1b);
wilson_L2 = w #( wilson_part2a - wilson_part2b); wilson_U2 = w #( wilson_part2a + wilson_part2b);
newcombe_L1 = wilson_L1[+,];
newcombe_U1 = wilson_U1[+,];
newcombe_L2 =wilson_L2[+,];
newcombe_U2 =wilson_U2[+,];
* 3. Calculate lambda1 and 2;
lam1= sq_w/n_trt1;
lam2 = sq_w/n_trt2;
lambda1 = lam1[+,];
lambda2 = lam2[+,];
* 4. calculate lower/upper limit for newcombe CI;
NEWCOMBE_L = wald_part1 - probit(&alpha)#sqrt(lambda1 # newcombe_L1 #(1- newcombe_L1)+lambda2#newcombe_U2#(1-newcombe_U2));
NEWCOMBE_U = wald_part1 + probit(&alpha)#sqrt(lambda2 # newcombe_L2 #(1-newcombe_L2)+lambda1#newcombe_U1#(1-newcombe_U1));
Newcombe_Lower=newcombe_L;
Newcombe_Upper=newcombe_U;
/* end of Newcombe CI */
Stratum_nm = {'Strata','Sample size in Trt1','Sample size in Trt2','Proportion in Trt1', 'Proportion in Trt2','Proportion difference (Trt1 - Trt2)','CMH weight', 'Inverse of Variance weight', 'Minimum risk weight' };
Adj_nm = {"Adjusted proportion difference", "Wald CI (lower)","Wald CI (upper)", "Newcombe CI (lower)","Newcombe CI (upper)"} ;
Weights= { "CMH", "Inverse of Variance", "Minimum risk"};
CI = {"Crude CI (lower)","Crude CI (upper)"};
Strata_sum = (1:n_strata)`||n_trt1|| n_trt2 || prop1 || prop2 || r_diff || CMH_weight || invar_weight || min_weight;
Adj_sum = wald_part1`||wald_Lower || Wald_Upper ||newcombe_Lower || newcombe_Upper ; Crude_CI = crude_Lower` || crude_Upper`;
file print; put @35 "Strata-adjusted confidence interval"//;
put @35 "Summary for dataset "/; put @30 "Sample size " @75 totaln / @30 "Number of strata " @75 n_strata;
print crude_CI[colname=CI label= "Unadjusted summary" format=6.4]; print Adj_sum [rowname=weights colname=Adj_nm label= "Adjusted summary" format=6.4];
quit;
%exit:
%mend stratCI;
DATA TEMP;
INPUT ID treatment strata response @@;
cards;
1 1 2 0
2 2 2 1
3 1 2 0
4 1 2 1
5 2 3 0
6 1 1 1
6 2 1 1
6 2 1 1
7 1 2 1
8 1 3 1
9 1 1 1
10 1 1 1
11 1 2 1
12 1 3 0
1 1 2 0
2 2 2 1
3 1 2 1
4 1 2 0
5 2 3 1
6 1 1 0
7 1 2 1
8 1 3 1
9 1 1 0
10 1 1 0
11 1 2 0
12 1 3 1
;
run;
%stratCI(data=TEMP, treatment=treatment, strata=strata, response=response, level=0.05);