File: cidxcn.f

package info (click to toggle)
hmisc 4.2-0-1
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 3,332 kB
  • sloc: asm: 27,116; fortran: 606; ansic: 411; xml: 160; makefile: 2
file content (88 lines) | stat: -rw-r--r-- 2,371 bytes parent folder | download | duplicates (9)
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