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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
C A copy from Hmisc:
C Package: Hmisc
C Version: 2.2-3
C Date: 2004-09-12
C Title: Harrell Miscellaneous
C Author: Frank E Harrell Jr <f.harrell@vanderbilt.edu>, with
C contributions from many other users.
C Maintainer: Frank E Harrell Jr <f.harrell@vanderbilt.edu>
C Depends: R (>= 1.9), grid, lattice, acepack
C Description: The Hmisc library contains many functions useful for data
C analysis, high-level graphics, utility operations, functions for
C computing sample size and power, importing datasets,
C imputing missing values, advanced table making, variable clustering,
C character string manipulation, conversion of S objects to LaTeX code,
C and recoding variables.
C License: GPL version 2 or newer
C URL: http://biostat.mc.vanderbilt.edu/s/Hmisc,
C http://biostat.mc.vanderbilt.edu/twiki/pub/Main/RS/sintro.pdf,
C http://biostat.mc.vanderbilt.edu/twiki/pub/Main/StatReport/summary.pdf
C Packaged: Sun Sep 12 18:06:07 2004; hornik
subroutine rcorr(xx, n, p, itype, dmat, npair, x, y, rx, ry, work,
& iwork)
implicit real*8 (a-h,o-z)
integer p, npair(p,p)
real*8 xx(n,p), dmat(p,p), x(1), y(1), rx(1), ry(1), work(1)
integer*4 iwork(1)
real*8 sumx,sumx2,sumy,sumy2,sumxy,z,a,b
do 23000 i=1, p
np=0
do 23002 k=1, n
if(.not.(xx(k,i).lt.1d30))goto 23004
np=np+1
23004 continue
23002 continue
npair(i,i)=np
do 23006 j=(i+1),p
m=0
if(.not.(itype.eq.1))goto 23008
sumx=0d0
sumy=0d0
sumx2=0d0
sumy2=0d0
sumxy=0d0
23008 continue
do 23010 k=1,n
xk=xx(k,i)
yk=xx(k,j)
if(.not.(xk.lt.1d30 .and. yk.lt.1d30))goto 23012
m=m+1
if(.not.(itype.eq.1))goto 23014
a=xk
b=yk
sumx=sumx+a
sumx2=sumx2+a*a
sumy=sumy+b
sumy2=sumy2+b*b
sumxy=sumxy+a*b
goto 23015
23014 continue
x(m)=xk
y(m)=yk
23015 continue
23012 continue
23010 continue
npair(i,j)=m
if(.not.(m.gt.1))goto 23016
if(.not.(itype.eq.1))goto 23018
z=m
d=(sumxy-sumx*sumy/z)/dsqrt((sumx2-sumx*sumx/z)*(sumy2-sumy*sumy/
&z))
goto 23019
23018 continue
call docorr(x, y, m, d, rx, ry, work, iwork)
23019 continue
dmat(i,j)=d
goto 23017
23016 continue
dmat(i,j)=1d30
23017 continue
23006 continue
23000 continue
do 23020 i=1,p
dmat(i,i)=1d0
do 23022 j=(i+1),p
dmat(j,i)=dmat(i,j)
npair(j,i)=npair(i,j)
23022 continue
23020 continue
return
end
subroutine docorr(x, y, n, d, rx, ry, work, iwork)
implicit real*8 (a-h,o-z)
real*8 x(1), y(1), rx(1), ry(1)
integer*4 iwork(1)
real*8 sumx,sumx2,sumy,sumy2,sumxy,a,b,z
call rank(n, x, work, iwork, rx)
call rank(n, y, work, iwork, ry)
sumx=0d0
sumx2=0d0
sumy=0d0
sumy2=0d0
sumxy=0d0
do 23024 i=1,n
a=rx(i)
b=ry(i)
sumx=sumx+a
sumx2=sumx2+a*a
sumy=sumy+b
sumy2=sumy2+b*b
sumxy=sumxy+a*b
23024 continue
z=n
d=(sumxy-sumx*sumy/z)/dsqrt((sumx2-sumx*sumx/z)*(sumy2-sumy*sumy/
&z))
return
end
C ------------------------------------------------------------------------------
|
