1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
|
C Output from Public domain Ratfor, version 1.01
C------------------------------------------------------------------------
C Compute c-index (c) and Brown-Hollander-Krowar-Goodman-Kruskal-Somer
C rank correlation (gamma) between X and Y with censoring indicator E
C Also returns number of relevant, concordant, and uncertain pairs
C (nrel, nconc, nuncert) and estimated s.d. of gamma (sd) using
C Quade formula (see SAS PROC MATPAR). Pairs with tied x are
C excluded if outx=.TRUE.
C
C F. Harrell 27Nov90
C Modification of SAS Procedure KGKC (1980)
C-------------------------------------------------------------------------
subroutine cidxcn(x,y,e,n,nrel,nconc,nuncert,c,gamma,sd,outx)
implicit double precision (a-h,o-z)
double precision x(n),y(n),dx,dy
logical e(n),outx
double precision nrel,nuncert,nconc
nconc=0d0
nrel=0d0
nuncert=0d0
sumr=0d0
sumr2=0d0
sumw=0d0
sumw2=0d0
sumrw=0d0
do23000 i=1,n
wi=0d0
ri=0d0
do23002 j=1,n
if(j.ne.i)then
dx=x(i)-x(j)
dy=y(i)-y(j)
if(dx.ne.0. .or. .not.outx)then
if((e(i).and.dy.lt.0.).or.(e(i).and..not.e(j).and.dy.eq.0.))then
if(dx.lt.0.)then
nconc=nconc+1d0
wi=wi+1d0
else
if(dx.eq.0.)then
nconc=nconc+.5d0
else
wi=wi-1d0
endif
endif
nrel=nrel+1d0
ri=ri+1d0
else
if((e(j).and.dy.gt.0.).or.(e(j).and..not.e(i).and.dy.eq.0.))then
if(dx.gt.0.)then
nconc=nconc+1d0
wi=wi+1d0
else
if(dx.eq.0.)then
nconc=nconc+.5d0
else
wi=wi-1d0
endif
endif
nrel=nrel+1d0
ri=ri+1d0
else
if(.not.(e(i).and.e(j)))then
nuncert=nuncert+1d0
endif
endif
endif
endif
endif
23002 continue
23003 continue
sumr=sumr+ri
sumr2=sumr2+ri*ri
sumw=sumw+wi
sumw2=sumw2+wi*wi
sumrw=sumrw+ri*wi
23000 continue
23001 continue
c=nconc/nrel
gamma=2.*(c-.5)
Ccall dblepr('sumr',4,sumr,1)
Ccall dblepr('sumw',4,sumw,1)
Ccall dblepr('sumr2',5,sumr2,1)
Ccall dblepr('sumw2',5,sumw2,1)
Ccall dblepr('sumrw',5,sumrw,1)
sd=sumr2*sumw**2-2d0*sumr*sumw*sumrw+sumw2*sumr**2
sd=2.*dsqrt(sd)/sumr/sumr
return
end
|