C-------------------------------------------------------- C COMPUTE PARTIAL CORRELATION COEFFICIENT AND C SIGNIFICANCE FOR CENSORED DATA C C AUGUST 1995 C C THE CODE IS BASED ON THE METHODOLOGY PRESENTED IN C 'A test for partial correlation with censored C astronomical data' C BY C M.G.Akritas and J.Siebert C Monthly Notices of the Royal Astronomical Society C 278, 919-924, 1996 C------------------------------------------------------- program partial_tau common /data/ dat(500,3),idat(500,3) common ntot common /kx/ k1,k2,k3 C------------------------------------------------------------- C INPUT DATA FILE CALLED 'DATA'. C CURRENTLY THE DATA FORMAT IS FIXED TO 3(f10.4,1x,i2,1x). C 1ST, 3RD AND 5TH COLUMN ARE INDEPENDENT, DEPENDENT AND TEST C VARIABLE, RESPECTIVELY. 2ND, 4TH AND 6TH COLUMN DENOTE C CENSORING WITH 1 = DETECTION, 0= UPPER LIMIT C EXAMPLE: C ' 26.9800 1 44.4340 0 -1.0714 1 ' C------------------------------------------------------------- open(10,file='data',status='old') C------------------------------------------------- C READ IN DATA: C DAT(I,K) = MEASUREMENT I OF VARIABLE K C IDAT(I,K) = CENSORING INDICATOR FOR DATA POINT (I,K) C DETECTION --> IDAT(I,K)=1 C UPPER LIMIT --> IDAT(I,K)=0 C-------------------------------------------------- i=1 1 read(10,110,end=99) dat(i,1),idat(i,1),dat(i,2), # idat(i,2),dat(i,3),idat(i,3) 110 format(3(f10.4,1x,i2,1x)) dat(i,1)=-dat(i,1) ! CHANGE TO RIGHT CENSORING dat(i,2)=-dat(i,2) ! CHANGE TO RIGHT CENSORING dat(i,3)=-dat(i,3) ! CHANGE TO RIGHT CENSORING i=i+1 goto 1 99 ntot=i-1 k1=1 ! INDEPENDENT VARIABLE = 1.COL OF DAT k2=2 ! DEPENDENT VARIABLE = 2.COL OF DAT k3=3 ! THIRD VARIABLE = 3.COL OF DAT call tau123(res) ! COMPUTE PARTIAL KENDALLS TAU write(6,*) 'Tau(1,2):',tau(k1,k2) write(6,*) 'Tau(1,3):',tau(k1,k3) write(6,*) 'Tau(2,3):',tau(k2,k3) write(6,*) '--> Partial Kendalls tau:', res write(6,*) ' ' write(6,*) 'Calculating variance...this takes some time....' write(6,*) ' ' call sigma(sig) ! COMPUTE VARIANCE write(6,*) 'Square root of variance (sigma):',sig write(6,*) ' ' if(abs(res/sig).gt.1.96) then write(6,*) 'Zero partial correlation rejected at level 0.05' else write(6,*) 'Null hypothesis cannot be rejected!' write(6,*) '(--> No correlation present, if influence of #third variable is excluded)' endif stop end C------------------------------------------------------ C-------- SUBROUTINES AND FUNCTIONS ------------------- C------------------------------------------------------ C------------ TAU123 --------------------------------- C-------- PARTIAL KENDALLS TAU ----------------------- subroutine tau123(res) common /kx/ k1,k2,k3 res= (tau(k1,k2)-tau(k1,k3)*tau(k2,k3))/ # sqrt((1.-tau(k1,k3)**2)*(1.-tau(k2,k3)**2)) end C------------ SIGMA ----------------------------------- C-------- VARIANCE OF STATISTIC ----------------------- subroutine sigma(sigres) common ntot common /kx/ k1,k2,k3 sig2=an( )/(ntot*(1.-tau(k1,k3)**2)*(1.-tau(k2,k3)**2)) sigres=sqrt(sig2) end C------------ AN --------------------------------------- C------- COMPUTES VALUE FOR A_N ------------------------- function an( ) double precision aasum(500) common ntot common /data/ dat(500,3),idat(500,3) common /kx/ k1,k2,k3 c1=16./(float(ntot)-1.) c2=6./((float(ntot)-1.)*(float(ntot)-2.)*(float(ntot)-3.)) asum=0.0 ave = 0.0 do 5 i=1,ntot aasum(i)=0.0 5 continue do 10 i1=1,ntot ! OUTER SUMMATION (I1) write(6,*) i1 do 11 j1=1,ntot-2 ! INNER SUMMATION WITH if(j1.eq.i1) goto 11 ! J1