|
|
沙发
楼主 |
发表于 2005-7-5 10:01:54
|
只看该作者
Re: SAS怎么做RIDIT分析啊？
[b:88781]one example[/b:88781]
[quote:88781]%let group=3; /* 三种药物 */
%let rank=4; /* 每种药物的疗效有四个等级 */
%let col_row=%eval(&group*&rank);
DATA _NULL_;
row=&rank;
col=&group;
ARRAY num(&rank,&group);
ARRAY cpd(&group);
ARRAY ridit(&group);
ARRAY trow(&rank);
ARRAY tcol(&group);
ARRAY y(&rank);
ARRAY ucpd(&group);
ARRAY uridit(&group);
ARRAY t(&group);
FILE PRINT;
n=0; buf1=0;
buf4=0; buf5=0;
INPUT num1-num&col_row;
DO j=1 TO col;
tcol(j)=0;
DO i=1 TO row;
tcol(j)=tcol(j)+num(i,j);
END;
n=n+tcol(j);
END;
buf3=n;
DO i=1 TO row;
trow(i)=0;
DO j=1 TO col;
trow(i)=trow(i)+num(i,j);
END;
buf4=buf4+trow(i)*trow(i)*trow(i);
buf2=trow(i);
y(i)=buf3-buf1-buf2;
buf3=y(i);
buf1=trow(i);
END;
DO j=1 TO col;
cpd(j)=0;
DO i=1 TO row;
cpd(j)=cpd(j)+num(i,j)*y(i);
END;
ridit(j)=0.5-cpd(j)/(2*tcol(j)*n);
t(j)=tcol(j)*(n+1)/2-cpd(j)/2;
ucpd(j)=cpd(j)/sqrt(tcol(j)*(n-tcol(j))*(n*n*n-buf4)/(3*n*(n-1)));
uridit(j)=(0.5-ridit(j))/sqrt((n*n*n-buf4)/(12*n*n*(n-1)*tcol(j)));
buf5=buf5+t(j)*t(j)/tcol(j); END;
Hc=(12/(n*(n+1))*buf5-3*(n+1))/(1-(buf4-n)/(n*(n*n-1)));
prob_Hc=1-PROBCHI(Hc,col-1);
PUT @5 "秩和检验结果"
#3 @5 '秩和检验统计量值' @30 '秩和检验P值'
#4 @5 Hc 10.4 @30 prob_Hc 6.4;
PUT #7 @5 'Ridit分析结果'
#9 @5 'Ridit值' @15 "检验Ridit所对应的U统计量的值";
DO j=1 TO col;
PUT #(j+9) @5 ridit(j) 7.4 @15 uridit(j) 8.4;
END;
buf1=probit(0.05/2/col); buf2=-buf1;
buf3=probit(0.05/2); buf4=-buf3;
PUT #(col+12) @5 "U 统计量值的95％置信区间"
@30 buf3 7.4 @40 buf4 7.4;
buf1=probit(0.01/2/col); buf2=-buf1;
buf3=probit(0.01/2); buf4=-buf3;
PUT #(col+14) @5 "U 统计量值的99％置信区间"
@30 buf3 7.4 @40 buf4 7.4;
CARDS;
15 4 1
49 9 15
31 50 45
5 22 24
;
RUN;[/quote:88781] |
|
