NOTE: Some users have encountered difficulty compiling ASURV on Linux and MacOS systems. This problem may be solved by renaming subroutines STAT and UNPACK throughout the code. #! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'README' <<'END_OF_FILE' XREADME for ASURV X XThis shell archive contains source code and documentation for ASURV, XAstronomy SURVival analysis (Rev. 1.3). ASURV implements a suite of Xstatistical methods for the analysis of censored data; i.e. data Xwhich are known to lie above or below some limit. It was written Xspecifically to treat left-censoring arising in observational astronomy Xwhen objects are observed but sometimes not detected due to sensitivity Xlimits. However, the methods can be useful to researchers in other Xdisciplines, as the code includes techniques that are often omitted Xfrom commercial survival analysis packages. X XASURV computes: the maximum-likelihood Kaplan-Meier estimator; several Xunivariate two-sample tests (Gehan, Peto-Peto, Peto-Prentice); three Xbivariate correlation coefficients (Cox regression, generalized Kendall's Xtau and Spearman's rho); and three linear regressions (EM algorithm Xassuming normal residuals, Buckley-James line, Schmitt line). The user Xis strongly encouraged to read the manual and references below if they Xare unfamiliar with these methods. The program is stand-alone and does Xnot call any specialized library. X XThis archive contains: this README file; asurv_code.f with 59 Xsurbroutines in FORTRAN 77; asurv.doc and asurv.tex with the Users' XManual in ASCII and LaTeX respectively; and asurv.etc with test files, Xsubroutine list and Bug Report form. It is distributed via the World XWide Web (ftp://www.astro.psu.edu/users/edf/asurv.shar) or by email Xrequest to code@stat.psu.edu. The archive is unpacked using a UNIX Xcommand like `sh asurv.shar'. The source code can be compiled with Xa statement like `f77 -f -o asurv asurv_code.f'. Operation has been Xverified for a variety of computer platforms including UNIX, VMS, VM Xand DOS. X XASURV was written between 1987 and 1992 by Drs. Takashi Isobe (Center Xfor Space Research, MIT), Michael LaValley (formerly at Dept. of XStatistics, Penn State), and Eric Feigelson. Code development was Xsupported by several NASA grants. Questions and problems should be Xaddressed to: Eric Feigelson, Dept. of Astronomy & Astrophysics, XPennsylvania State University, University Park PA USA, FAX 814-863-3399, XEmail edf@astro.psu.edu, WWW http://www.astro.psu.edu/users/edf). X XIMPORTANT: The authors grant researchers and students permission to Xuse and copy ASURV code and associated material for non-commercial Xpurposes. We request that publications resulting from its use cite Xone of the references below. This software is provided `as is' without Xany expressed or implied warranty. X XReferences: XFeigelson, E. D. and Nelson, P. I. ``Statistical Methods for X Astronomical Data with Upper Limits: I. Univariate Distributions', X Astrophyscal Journal 293, 192-206, 1985. XIsobe, T., Feigelson, E. D., and Nelson, P. I. ``Statistical Methods X for Astronomical Data with Upper Limits: II. Correlation and Regression', X Astrophysical Journal, 306, 490-507, 1986. XLaValley, M., Isobe, T. and Feigelson, E.D. ``ASURV'', Bulletin X Amercan Astronomical Society (Software Reports), 22, 917-918, 1990. X XRevisions: XRev. 0 (1987-1990) Incomplete and obsolete. XRev. 1 (1992-present) Essentially identical versions with minor bugs X corrected. X XREADME for asurv.shar written by Eric feigelson (Sept. 1996) X END_OF_FILE if test 3299 -ne `wc -c <'README'`; then echo shar: \"'README'\" unpacked with wrong size! fi # end of 'README' fi if test -f 'asurv_code.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'asurv_code.f'\" else echo shar: Extracting \"'asurv_code.f'\" \(370641 characters\) sed "s/^X//" >'asurv_code.f' <<'END_OF_FILE' XC XC ********************************************************************** XC *********************** SUBROUTINE AARRAY **************************** XC ********************************************************************** XC X SUBROUTINE AARRAY(Z,IND,ISTA,IS,NTOT,NDAT,MVAR,NG1,NG2,XY, X + ID1,ID2,ISAVE) XC XC******* ISIGN,IFULL IS ADDED ON "COMMON' STATEMENT * XC * * XC * INPUT Z(I,J) : DATA TO BE TESTED * XC * IND(I,J): INDICATOR OF CENSORING * XC * ISTA(I) : INDICATOR OF GROUP * XC * IS : IS-TH SUB-DATA SET * XC * NG1 : INDICATOR OF THE FIRST GROUP * XC * NG2 : INDICATOR OF THE SECOND GROUP * XC * NTOT : TOTAL NUMBER OF DATA POINTS * XC * LL : INDICATOR OF OUTPUT FILE * XC * IPR : INDICATOR FOR PRINTING * XC * * XC * OUTPUT N : NTOT * XC * N1 : NUMBER OF DATA POINTS IN GROUP 1 * XC * N2 : NUMBER OF DATA POINTS IN GROUP 2 * XC * NCEN : NUMBER OF CENSORED DATA POINTS * XC * ISIGN : INDICATOR OF LOWER/UPPER LIMITS * XC * * XC * PUT ALL OBS. IN ARRAY XY AND FORM ARRAYS ID1 AND ID2 * XC * ID1(I)=0 : ITH OBS. IS UNCENSORED * XC * 1 : ITH OBS. IS CENSORED * XC * ID2(I)=J : ITH OBS. IS FROM ITH SAMPLE, J=1,2 * XC * * XC * SUBROUTINES * XC * SORT2 * XC XC******* ALTHOUGH THIS SUBROUTINE HAS THE SAME NAME AS A PROGRAM FROM * XC******* "STATISTICAL METHODS FOR SURVIVAL DATA ANALYSIS" BY ELISA T. * XC******* LEE, 1980, LIFETIME LEARNING PUBLICATIONS (BELMONT:CA), * XC******* IT IS DIFFERENT EXCEPT IN THE GENERAL PURPOSE. * XC******* ID1(I) IS ASSIGNED IN THE OPPOSITE WAY SO THAT THE PPROGRAM * XC******* CAN USE THE DATA SETS WHICH ARE MADE FOR OTHER PROGRAMS. * XC X X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X X DIMENSION Z(MVAR,NDAT),IND(MVAR,NDAT),ISTA(NTOT) X DIMENSION XY(NTOT),ID1(NTOT),ID2(NTOT),ISAVE(NTOT) X COMMON /G/ NCOMP,N1,N2,NCEN,ISIGN,IFULL,LO XC X IU=0 X NCEN=0 X NCOMP=0 X N1=0 X N2=0 X ISIGN=1 XC XC * FIND THE CENSORSHIP OF THE DATA SET. -1 FOR UPPER LIMITS * XC * AND 1 FOR LOWER LIMITS * XC X DO 100 I=1,NTOT X ISAVE(I) = 0 X IF(IND(IS,I) .EQ. 0) GOTO 100 X ISIGN=IND(IS,I)/IABS(IND(IS,I)) X ISAVE(I) = ISIGN X 100 CONTINUE X XC * CHECK WHETHER THE UPPER AND LOWER LIMITS ARE MIXED IN THE SAME * XC * VARIABLE. IF SO, THE PROGRAM IS TERMINATED. * XC * THIS TEST WAS ADDED. * XC X DO 110 I = 1, NTOT X IF(ISAVE(I) .EQ. 0) GOTO 110 X IF(ISAVE(I) .NE. ISIGN) THEN X PRINT * X PRINT *,'YOU CANNOT HAVE BOTH UPPER AND LOWER LIMITS' X PRINT *,'IN ONE VARIABLE AT THE SAME TIME.' X PRINT *,'PLEASE CHECK YOUR DATA.' X PRINT *,'THE PROGRAM HAS BEEN TERMINATED.' X PRINT * X STOP X ENDIF X 110 CONTINUE XC XC * COUNT NUMBER OF DATA POINTS IN THE TWO SUBSAMPLES * XC X X DO 400 I = 1, NTOT X IF((ISTA(I) .EQ. NG1) .OR. (ISTA(I) .EQ. NG2)) THEN X NCOMP = NCOMP + 1 X XY(NCOMP) = ISIGN*Z(IS,I) X IF(ISTA(I) .EQ. NG1) ID2(NCOMP) = 1 X IF(ISTA(I) .EQ. NG2) ID2(NCOMP) = 2 X IF(IABS(IND(IS,I)) .NE. 1) THEN X ID1(NCOMP) = 0 X IU = IU + 1 X IF(ID2(NCOMP) .EQ. 1) N1 = N1 + 1 X IF(ID2(NCOMP) .EQ. 2) N2 = N2 + 1 X ELSE X ID1(NCOMP) = 1 X NCEN = NCEN + 1 X IF(ID2(NCOMP) .EQ. 1) N1 = N1 + 1 X IF(ID2(NCOMP) .EQ. 2) N2 = N2 + 1 X ENDIF X ENDIF X 400 CONTINUE X X CALL SORT2(XY, ID1, ID2, NTOT) X X RETURN X END X X XC XC ********************************************************************** XC ********************* FUNCTION AGAUSS ******************************* XC ********************************************************************** XC X FUNCTION AGAUSS(Z) XC XC * EVALUATES THE INTEGRAL OF THE GAUSSIAN PROBABILITY FUNCTION * XC * OBTAINED FROM PROGRAM 3-5 ON P. 35 OF "DATA REDUCTION AND * XC * ERROR ANALYSIS FOR THE PHYSICAL SCIENCES", P. R. BEVINGTON, * XC * 1969, McGRAW HILL, (NY:NY). * XC * * XC XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) XC XC X Z=DABS(Z) X AGAUSS=1.0 XC XC * IF Z>5.0, USE APPROXIMATION FOR PROB TO AVOID ERROR * XC X IF(Z.LE.5.0) THEN X DENOM=1.0 X IF(Z .GT. 0.0) THEN X TERM=0.7071067812D00*Z X SUM=TERM X Y2=(Z**2)/2.0 X 31 DENOM=DENOM+2.0 X TERM=TERM*(Y2*2.0/DENOM) X SUM=SUM+TERM X IF(TERM/SUM-1.0E-10 .GT. 0.0) THEN X GOTO 31 X ELSE X AGAUSS=1.128379167D00*SUM*DEXP(-Y2) X ENDIF X ELSE X AGAUSS = 0.0 X ENDIF X ENDIF X X RETURN X END X XC XC************************************************************************* XC********************* SUBROUTINE AKRANK ********************************* XC************************************************************************* XC XC X SUBROUTINE AKRANK(IND, X, NTOT, IP, R, MVAR,ZU, ZC, X + PL, F, V, FMASS, ITEMP, PTEMP,Z1, X + WRK1,WRK2,WRK3,DWRK1,IWRK1,SWRK1) XC XC * THIS SUBROUTINE COMPUTES AKRITAS' RANK * XC * * XC * REFERENCE * XC * PENN STATE UNIVERSITY, DEPARTMENT OF STATISTICS, * XC * TECHNICAL REPORTS AND PREPRINTS SERIES, NUMBER 87, * XC * "ALIGNED RANK TESTS FOR REGRESSION WITH CENSORED DATA", * XC * MICHAEL G. AKRITAS, SEPTEMBER 1989 * XC * INPUT * XC * IND : INDICATOR OF CENSORSHIP * XC * X : VARIABLE * XC * NTOT : TOTAL NUMBER OF DATA POINTS * XC * * XC * OUTPUT * XC * R : RANK * XC * PL : PL ESTIMATOR * XC * F : 1.0 - PL (DISTRIBUTION FUNCTION) * XC * * XC * OTHER VARIABLES * XC * IP : INDEX OF VARIABLE BEING RANKED * XC * MVAR : NUMBER OF VARIABLES * XC * ZU : DETECTED DATA * XC * ZC : CENSORED DATA * XC * FMASS : JUMPS IN PL ESTIMATOR * XC * IU : NUMBER OF DETECTIONS * XC * IC : NUMBER OF CENSORED DATA POINTS * XC * PTEMP : TEMPORARY STORAGE OF PL ESTIMATOR * XC * * XC * SUBROUTINES * XC * XVAR, PLESTM, SORT1 * XC XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X X DIMENSION IND(MVAR, NTOT), X(MVAR, NTOT), ZU(NTOT), ZC(NTOT) X DIMENSION PL(NTOT), R(MVAR, NTOT), F(NTOT), V(NTOT), FMASS(NTOT) X DIMENSION Z1(MVAR, NTOT), ITEMP(NTOT), PTEMP(NTOT) X DIMENSION IWRK1(NTOT),DWRK1(MVAR,NTOT),SWRK1(MVAR) X DIMENSION WRK1(NTOT),WRK2(NTOT),WRK3(NTOT) X XC XC * CALL SUBROUTINE XVAR : DISTINGUISH DETECTIONS AND CENSORED * XC * DATA POINTS. * XC X X CALL XVAR(IND,X,IP,NTOT,ISIGN,ZU,ZC,IU,IC,IWRK1,WRK1,WRK2,WRK3, X + DWRK1,SWRK1,LTOT,MVAR,ITEMP) X XC XC * CALL PLESTM : PL ESTIMATOR COMPUTATION * XC XC X DO 5 I = 1, NTOT X ITEMP(I) = 0 X Z1(1, I) = 0.0 X IWRK1(I) = IND(IP,I) X 5 CONTINUE X X IF(IU .EQ. 0) THEN X WRITE(6,3) X 3 FORMAT('NO DETECTIONS: PROGRAM IS TERMINATED') X STOP X ENDIF X X CALL SORT1(IWRK1,Z1,ZU,IU,1,ITEMP,SWRK1,MVAR) XC X IF(IC .NE. 0) CALL SORT1(IWRK1,Z1,ZC,IC,1,ITEMP,SWRK1,MVAR) X XC XC >>>> Bug fixed Sept. 1996. NCH was missing from following line <<<< X X CALL PLESTM(ZU, ZC, IU, IC, PL, V, NTOT,SMEAN,SIGM,ICH,NCH,IWRK1) X XC XC * IF THE DATA CONTAINS CENSORED DATA POINTS, THE PRODUCT LIMIT * XC * ESTIMATOR MUST BE ADJUSTED TO INCLUDE CENSORED DATA POINTS. * XC X IF(IC .NE. 0) THEN X XC * IF THE DATA HAS UPPER LIMITS, FIRST THE PRODUCT LIMIT ESTIMATOR* XC * MUST BE ADJUSTED. * X X IF(ISIGN .LT. 0) THEN X FMASS(1) = 1.0 - PL(1) X DO 10 I = 2, IU X FMASS(I) = PL(I-1)-PL(I) X 10 CONTINUE X X J = IU/2 X DO 20 I = 1, J X FTEMP=FMASS(I) X FMASS(I)=FMASS(IU-I+1) X FMASS(IU-I+1)=FTEMP X 20 CONTINUE X X DO 40 I = 1, IU X PTEMP(I) = 1.0 X DO 30 J = 1, I X PTEMP(I)=PTEMP(I)-FMASS(J) X 30 CONTINUE X 40 CONTINUE X ELSE X DO 50 I = 1, IU X PTEMP(I)=PL(I) X 50 CONTINUE X ENDIF X XC * NOW, PRODUCT LIMIT ESTIMATOR VALUES ARE ASSIGNED TO CENSORED * XC * DATA POINTS. * X X IF(IND(IP,1) .EQ. 0) THEN X PL(1) = PTEMP(1) X J = 1 X ELSE X PL(1) = 1.0 X J = 0 X ENDIF X X DO 60 I = 2, NTOT X IF(IND(IP, I) .EQ. 0) THEN X J = J + 1 X PL(I) = PTEMP(J) X ELSE X PL(I) = PL(I-1) X ENDIF X 60 CONTINUE X X ENDIF X XC * THE PRODUCT LIMIT ESTIMATE IS NOW USED TO ESTIMATE THE * XC * DISTRIBUTION FUNCTION (F) AT ALL POINTS. * X X DO 65 I = 1, NTOT X F(I) = 1.0 - PL(I) X 65 CONTINUE X XC XC * COMPUTE HERE AKRITAS' RANK USING F-VALUES * XC X DO 90 I = 1, NTOT X IF(IND(IP, I) .EQ. 0) THEN X R(IP, I) = REAL(NTOT)*F(I) X ELSEIF(IND(IP, I) .GT. 0) THEN X R(IP, I) = REAL(NTOT)*(0.5 + 0.5*F(I)) X ELSE X R(IP, I) = NTOT*(0.5*F(I)) X ENDIF X 90 CONTINUE X X RETURN X END X XC XC ********************************************************************** XC ******************** SUBROUTINE ARISK ******************************* XC ********************************************************************** XC X SUBROUTINE ARISK(R,XM,X,E1,NG,H,XY,ID1,NTOT) XC XC XC * THIS SUBROUTINE COMPUTES THE FOLLOWING FOUR * XC * ARRAYS FOR SUBROUTINE COX, LRANK, AND PWLCXN. * XC * R(I) : NO. OF OBSERVATIONS IN RISK SET AT THE * XC * I-TH DISTINCT FAILURE TIME. * XC * XM(I) : MULTIPLICITY OF THE I-TH DISTINCT * XC * FAILURE TIME. * XC * E1(I) : XM(I)/R(I) * XC * H(I) : KAPLAN AND MEIER'S ESTIMATES OF THE * XC * SURVIVOR FUNCTION * XC * * XC * X(I) : THE ARRAY OF DISTINCT FAILURE TIMES * XC * NG : NO OF X * XC * THIS SUBROUTINE IS OBTAINED FROM ELISA T. LEE, "STATISTICAL * XC * METHODS FOR SURVIVAL DATA ANALYSIS", 1980, LIFETIME LEARNING * XC * PUBLICATIONS (BELMONT:CA); BUT HAS BEEN SIGNIFICANTLY MODIFIED. * XC * * XC X X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X X DIMENSION R(NTOT),XM(NTOT),X(NTOT),H(NTOT),E1(NTOT) X DIMENSION XY(NTOT),ID1(NTOT) X COMMON /G/ NCOMP,N1,N2,NCEN,ISIGN,IFULL,LO XC X L=1 X I=1 X R(L)=REAL(NCOMP) XC XC * COMPUTE RISK SETS, AND OTHER QUANTITIES * XC X 24 IF(ID1(I).NE.0) THEN X R(L)=R(L)-1.0 X I=I+1 X GOTO 24 X ENDIF X 25 XM(L)=1.0 X XNC=0.0 X TEMP=XY(I) X X(L)=TEMP X X 21 IF(I.NE.NCOMP) THEN X I=I+1 XC X IF(ID1(I).NE.1) THEN X IF(TEMP.NE.XY(I)) GOTO 20 X XM(L)=XM(L)+1.0 X GOTO 21 X ENDIF X X 26 XNC=XNC+1.0 X X(L)=TEMP X GOTO 21 X X 20 L=L+1 X R(L)=R(L-1)-XM(L-1)-XNC X GOTO 25 X ENDIF X X 23 X(L)=TEMP X NG=L XC XC * COMPUTE KM ESTIMATOR * X X DO 30 I=1,NG X E1(I)=XM(I)/R(I) X 30 CONTINUE X X H(1)=1.0 X NG1=NG+1 X X DO 31 I=2,NG1 X H(I)=H(I-1)*(1.0-E1(I-1)) X 31 CONTINUE X X RETURN X END X X XC XC XC * ASURV: SURVIVAL ANALYSIS PACKAGE FOR ASTRONOMERS * XC * * XC * DEVELOPED BY: TAKASHI ISOBE * XC * CENTER FOR SPACE RESEARCH * XC * THE MASSACHUSETTS INSTITUTE OF TECHNOLOGY * XC * * XC * MICHAEL LAVALLEY * XC * DEPARTMENT OF STATISTICS * XC * THE PENSYLVANIA STATE UNIVERSITY * XC * 330A CLASSROOM BUILDING, UNIVERSITY PARK PA 16802 * XC * INTERNET: MLV@STAT.PSU.EDU * XC * * XC * ERIC FEIGELSON * XC * DEPARTMENT OF ASTRONOMY AND ASTROPHYSICS * XC * THE PENSYLVANIA STATE UNIVERSITY * XC * 525 DAVEY LAB. UNIVERSITY PARK PA 16802 * XC * * XC * REV. 1.2 SECOND UPDATE SUMMER 1992 * XC * * XC * THIS PACKAGE IS WRITTEN TO PROVIDE SEVERAL * XC * SURVIVAL ANALYSIS METHODS WHICH ARE USEFUL IN ANALYZING * XC * ASTRONOMICAL DATA. SURVIVAL ANALYSIS IS A GROUP OF STATISTICAL * XC * METHODS WHICH TREAT PROBLEMS WITH CENSORED DATA (UPPER OR LOWER * XC * LIMITS). THIS PACKAGE INCLUDES SOME TECHNIQUES DEVELOPED IN * XC * IN OTHER FIELDS (E.G. ECONOMICS, ACTUARIAL SCIENCE, RELIABILITY * XC * MATHEMATICS), AND A FEW METHODS DEVELOPED BY ASTRONOMERS. * XC * * XC * THE METHODS PROVIDED IN THIS PACKAGE ARE : * XC * * XC * UNIVARIATE DISTRIBUTION : KAPLAN-MEIER ESTIMATOR * XC * TWO-SAMPLE TESTS : GEHAN TEST * XC * LOGRANK TEST * XC * PETO AND PETO TEST * XC * PETO AND PRENTICE TEST * XC * CORRELATION TESTS : COX PROPORTIONAL HAZARDS MODEL * XC * GENERALIZED KENDALL'S TAU (BHK METHOD)* XC * GENERALIZED SPEARMAN'S RHO * XC * (AKRITAS' METHOD) * XC * LINEAR REGRESSIONS : EM ALGORITHM WITH NORMAL DISTRIBUTION * XC * BUCKLEY-JAMES METHOD * XC * TWO-DIMENSIONAL KAPLAN-MEIER * XC * REGRESSION FOR DUAL-CENSORED DATA * XC * * XC * * XC * INPUTS * XC * * XC * IS0 : IF 1 : UNIVARIATE PROBLEM * XC * 2 : CORRELATION/REGRESSION PROBLEM * XC * 3 : EXIT * XC * * XC * SUBROUTINES DATA1, UNIVAR, BIVAR * XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X CHARACTER*1 BLANK XC OPEN(6,CARRIAGECONTROL='LIST',STATUS='OLD') XC XC XC X PRINT * X PRINT *,' ***************************************************' X PRINT *,' * *' X PRINT *,' * WELCOME TO ASURV *' X PRINT *,' * SURVIVAL ANALYSIS PACKAGE *' X PRINT *,' * FOR ASTRONOMERS *' X PRINT *,' * *' X PRINT *,' * DEVELOPED BY: *' X PRINT *,' * TAKASHI ISOBE *' X PRINT *,' * (CENTER FOR SPACE RESEARCH, MIT) *' X PRINT *,' * MICHAEL LAVALLEY *' X PRINT *,' * (DEPT. OF STATISTICS, PENN STATE) *' X PRINT *,' * ERIC FEIGELSON *' X PRINT *,' * (DEPT. OF ASTRONOMY & ASTROPHYSICS, PENN STATE) *' X PRINT *,' * *' X PRINT *,' * *' X PRINT *,' * REV 1.2 SUMMER 1992 *' X PRINT *,' ***************************************************' X PRINT * X PRINT * X PRINT * X PRINT * X PRINT *,' (CARRIAGE RETURN TO CONTINUE) ' X READ(5,50) BLANK X 50 FORMAT(A1) X PRINT * XC XC * START CONVERSATION WITH THE USER * XC X PRINT * X PRINT * X PRINT * X 100 PRINT *,' MENU ' X PRINT * X PRINT * X PRINT *,' UNIVARIATE DATA BIVARIATE DATA ' X PRINT * X PRINT * X PRINT *,' DISTRIBUTION FUNCTION CORRELATION ' X PRINT *,' 1 KAPLAN-MEIER ESTIMATOR 1 COX REGRESSION ' X PRINT *,' 2 GEN. KENDALL TAU' X PRINT *,' 3 GEN. SPEARMAN RHO' X PRINT * X PRINT * X PRINT *,' TWO-SAMPLE TESTS LINEAR REGRESSION ' X PRINT *,' 1 GEHAN TESTS 1 EM ALGORITHM WITH ' X PRINT *,' 2 LOGRANK TEST GAUSSIAN RESIDUALS ' X PRINT *,' 3 PETO AND PETO TEST 2 BUCKLEY-JAMES METHOD ' X PRINT *,' 4 PETO AND PRENTICE TEST WITH KM RESIDUALS ' X PRINT *,' 3 SCHMITT METHOD FOR ' X PRINT *,' DUAL CENSORED DATA ' X PRINT * X PRINT * X PRINT * X PRINT *,' (CARRIAGE RETURN TO CONTINUE) ' X READ(5,50) BLANK XC X PRINT * XC XC * CHOICE : UNIVARIATE PROBLEM OR CORRELATION/REGRESSION PROBLEM * XC X PRINT * X PRINT *,' SELECT DATA TYPE: ' X PRINT *,' 1 UNIVARIATE DATA ' X PRINT *,' 2 BIVARIATE DATA ' X PRINT *,' 3 EXIT ' X 200 WRITE(6,210) X 210 FORMAT(' CHOICE ? ') XC 210 FORMAT(' CHOICE ? ',$) XC X CALL DATA1(IS0) XC X IF((IS0.EQ.1).OR.(IS0.EQ.2).OR.(IS0.EQ.3)) GOTO 300 X PRINT *,'PLEASE TYPE ONCE MORE' X GOTO 200 XC X 300 IBACK=0 X IF(IS0.EQ.1) CALL UNIVAR(IBACK) X IF(IS0.EQ.2) CALL BIVAR(IBACK) X IF(IS0.EQ.3) STOP XC X IF(IBACK.EQ.1) GOTO 100 X STOP X END X XC XC ********************************************************************** XC ********************** SUBROUTINE BHK ******************************* XC ********************************************************************** XC X SUBROUTINE BHK(IND,XX,YY,NTOT,OUTPUT,X,Y,IAA,IBB,IP,MVAR) XC XC * GENERALIZED KENDALL'S TAU CORRELATION COEFFICIENT * XC * FOR CENSORED DATA * XC XC * THIS PROGRAM COMPUTES KENDALL'S TAU FOR BIVARIATE DATA * XC * SETS. THE DATA SETS CAN CONTAIN CENSORED POINTS IN THE * XC * INDEPENDENT VARIABLE AND/OR THE DEPENDENT VARIABLE. * XC * ALTHOUGH THIS PROGRAM GIVES ANSWERS FOR DATA SETS WHICH * XC * CONTAIN TIES, IT MAY NOT BE ACCURATE. * XC * PARAMETERS : * XC * INPUT * XC * NTOT : NUMBER OF OBSERVATIONS * XC * XX(1,I) : INDEPENDENT PARAMETER OF I-TH OBSERVATION * XC * YY(I) : DEPENDENT PARAMETER OF I-TH OBSERVATION * XC * IND(I) : INDICATOR OF CENSORED STATUS * XC * EACH POINT MUST BE SPECIFIED ITS CENSORED STATUS : * XC * FOR THE LOWER LIMITS * XC * 0 : DETECTED POINT * XC * 1 : ONLY DEPENDENT VARIABLE IS LOWER LIMIT * XC * 2 : ONLY INDEPENDENT VARIABLE IS LOWER LIMIT * XC * 3 : BOTH VARIABLES ARE LOWER LIMIT * XC * 4 : INDEPENDENT VARIABLE IS LOWER LIMIT AND * XC * DEPENDENT VARIABLE IS UPPER LIMIT * XC * FOR THE UPPER LIMITS, CHANGE THE SIGN OF ABOVE INDICATORS. * XC * * XC * WORK * XC * X(I) : =XX(1,I) * XC * Y(I) : =YY(I) * XC * IP(I) : =IND(I) * XC * IAA(I) : CONCORDANCE INFORMATION FOR X * XC * IBB(I) : CONCORDANCE INFORMATION FOR Y * XC * OUTPUT * XC * PROB : SIGNIFICANCE LEVEL FOR THE HYPOTHESIS THAT * XC * X AND Y ARE NOT CORRELATED UNDER THE * XC * GAUSSIAN DISTRIBUTION * XC * * XC * SUBROUTINES * XC * CENS, COEFF * XC * * XC * REF. BROWN, HOLLANDER, AND KORWAR 1974, IN RELIABILITY * XC * AND BIOMETRY P.327, EQNS 1 TO 8, PROSCHAN AND * XC * SERFLING EDS (SIAM) * XC XC * NOTE: THIS PROGRAM IS QUITE CPU INTENSIVE FOR LARGE DATA * XC * SETS (MORE THAN A FEW HUNDRED POINTS). * XC * * X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X X DIMENSION XX(MVAR,NTOT),YY(NTOT),IND(NTOT),X(NTOT),Y(NTOT) X DIMENSION IP(NTOT),IAA(NTOT),IBB(NTOT) X CHARACTER*9 OUTPUT XC X SIS =0.0 X ASUM =0.0 X BSUM =0.0 X AASUM=0.0 X BBSUM=0.0 XC XC * SUBSTITUE XX AND YY TO X AND Y SO THAT THE ORIGINAL VALUES * XC * WON'T BE CHANGED. * XC X DO 90 I=1,NTOT X X(I) = XX(1,I) X Y(I) = YY(I) X IP(I) = IND(I) X 90 CONTINUE XC XC XC * THE SUBROUTINE CENS ADDS OR SUBTRACTS A SMALL NUMBER * XC * FROM EACH CENSORED POINT SO THAT NO TIES WITH DETECTED * XC * POINTS OCCUR. * XC XC X CALL CENS(X,Y,IP,NTOT) XC XC XC * START MAKING INFORMATION FOR CONCORDANCE * XC XC X DO 1900 I=1,NTOT XC XC * INFORMATION OF CONCORDANCE FOR THE INDEPENDENT VAR. * XC X IA=2 X IB=3 X IC=4 X ID=-2 X IE=-3 X IG=-4 X IH=1 X IJ=-1 XC XC * SUBROUTINE WHICH FINDS CONCORDANCE INFORMATION * XC X CALL COEFF(I,X,IP,NTOT,IAA,IA,IB,IC,ID,IE,IG,IH,IJ) XC XC * INFORMATION OF CONCORDANCE FOR THE DEPENDENT VAR. * XC X IA=1 X IB=3 X IC=-4 X ID=-1 X IE=-3 X IG=4 X IH=2 X IJ=-2 X X CALL COEFF(I,Y,IP,NTOT,IBB,IA,IB,IC,ID,IE,IG,IH,IJ) XC XC * START COMPUTING QUANTITIES IS, IASUM, IBSUM, * XC * IAASUM, AND IBBSUM. * XC X DO 1800 J=1,NTOT X IF((IAA(J).EQ.0).AND.(IBB(J).EQ.0)) GOTO 1800 X SIS=SIS+IAA(J)*IBB(J) X ASUM=ASUM+IAA(J)**2 X BSUM=BSUM+IBB(J)**2 X X 1650 DO 1700 K=1,NTOT X IF(IAA(J).NE.0) THEN X IF(IAA(K).NE.0) THEN X AASUM=AASUM+IAA(J)*IAA(K) X ENDIF X ENDIF X 1670 IF(IBB(J).NE.0) THEN X IF(IBB(K).NE.0) THEN X BBSUM=BBSUM+IBB(J)*IBB(K) X ENDIF X ENDIF X 1700 CONTINUE X 1800 CONTINUE X 1900 CONTINUE XC XC * NOW COMPUTE THE STATISTIC AND THE PROBABILITY * XC X D1=REAL(NTOT*(NTOT-1)) X D2=REAL(D1*(NTOT-2)) X ALP=2.0*(ASUM*BSUM)/D1 X GAM=4.0*((AASUM-ASUM)*(BBSUM-BSUM))/D2 X VAR=ALP+GAM X SIGMA=DSQRT(VAR) X Z=SIS/SIGMA X PROB=1.0-AGAUSS(Z) XC X IF(OUTPUT.EQ.' ') THEN X WRITE(6,2030) X WRITE(6,2003) X WRITE(6,2030) X WRITE(6,2005) Z X WRITE(6,2007) PROB X WRITE(6,2030) X ELSE X WRITE(60,2030) X WRITE(60,2003) X WRITE(60,2030) X WRITE(60,2005) Z X WRITE(60,2007) PROB X WRITE(60,2030) X ENDIF X 2003 FORMAT(5X,'CORRELATION TEST BY GENERALIZED KENDALL`S TAU') X 2005 FORMAT(7X,'Z-VALUE =',F12.3) X 2007 FORMAT(7X,'PROBABILITY =',F13.4,/, X + ' (PROBABILITY THAT A CORRELATION IS NOT PRESENT)') X 2030 FORMAT(' ') X X RETURN X END X XC XC ********************************************************************** XC ********************** SUBROUTINE BIN ******************************* XC ********************************************************************** XC X SUBROUTINE BIN(NTOT,MX,MY,ISKIP,ICENS,DELX,DELY,XORG,YORG,MM, X + M1,M2,M3,M4,M5,M6,M7,M8,INDEX,LP,XT,YT,Z,SWRK1, X + X,Y,NP,XB,YB,F,N,N1,N2,N3,N4,N5,N6,N7,N8,IB,MVAR) X XC XC XC * * XC * THIS SUBROUTINE DOES BINNING AND CHANGES CENSORED POINTS * XC * WHICH DO NOT HAVE DETECTED POINTS ABOVE (OR BELOW) * XC * TO DETECTED POINTS. * XC * * XC * WARNING WARNING WARNING WARNING * XC * * XC * THE USER SHOULD BE WARNED THAT THIS SUBROUTINE ACTUALLY * XC * CHANGES THE DATA!! FIRST, IT REDEFINES SOME LIMITS TO * XC * DETECTIONS. IF THE BINS ARE CHOSEN TO BE TOO NARROW, THEN * XC * VIRTUALLY ALL LIMITS COULD BE CHANGED. SECOND, IT PUSHES * XC * EACH LIMIT INTO THE ADJACENT BIN. IF THE BINS ARE CHOSEN TO * XC * TO BE TOO WIDE, THIS SUBSTANTIALLY ALTERS THE MEASURED VALUES. * XC * THUS, THE USER MUST TREAD A FINE LINE IN CHOSING BIN SIZES. * XC * * XC * * XC * INPUT * XC * X(I) : INDEPENDENT VARIABLE * XC * Y(I) : DEPENDENT VARIABLE * XC * NP(I) : INDICATOR OF CENSORING * XC * NTOT : TOTAL NUMBER OF DATA * XC * MX : NUMBER OF BINS IN X * XC * MY : NUMBER OF BINS IN Y * XC * ISKIP : INDICATOR OF BINNING PROCESS * XC * ICENS : CENSORING STATUS OF THE DATA SET * XC * IF ISKIP>0, THE NEXT VALUES MUST BE PROVIDED : * XC * DELX : BIN SIZE OF X AXIS * XC * DELY : BIN SIZE OF Y AXIS * XC * XORIG : ORIGIN OF X * XC * YORIG : ORIGIN OF Y * XC * * XC * WORK * XC * YT(I) : COPY OF Y(I) FOR SORTING PROGRAM. * XC * M1 : # OF Y LOWER LIMITS CHANGED TO DETECTIONS * XC * M2 : # OF X LOWER LIMITS CHANGED TO DETECTIONS * XC * M3 : # OF DOUBLE LOWER LIMITS CHANGED TO * XC * DETECTIONS * XC * M4 : # OF Y LOWER , X UPPER LIMITS CHANGED TO * XC * DETECTIONS * XC * M5 : # OF Y UPPER LIMITS CHANGED TO DETECTIONS * XC * M6 : # OF X LOWER LIMITS CHANGED TO DETECTIONS * XC * M7 : # OF DOUBLE UPPER LIMITS CHANGED TO * XC * DETECTIONS * XC * M8 : # OF Y UPPER , X LOWER LIMITS CHANGED TO * XC * DETECTIONS * XC * NC1, NC2,...,NC8 : # OF CENSORED POINTS. SEE THE * XC * MAIN PROGRAM FOR THE DEFINITIONS * XC * IB : DIMENSION SIZE OF BINS * XC * * XC * OUTPUT * XC * F(I,J): INITIAL GUESS OF THE PROBABILITY OF THE * XC * BIN(I,J) * XC * N(I,J): NUMBER OF DETECTED POINTS IN THE BIN(I,J) * XC * N1(I,J): NUMBER OF Y LOWER LIMITS IN THE BIN(I,J) * XC * N2(I,J): NUMBER OF X LOWER LIMITS IN THE BIN(I,J) * XC * N3(I,J): NUMBER OF DOUBLE LOWER LIMITS IN THE BIN(I,J) * XC * N4(I,J): NUMBER OF Y LOWER, X UPPER LIMITS IN THE * XC * BIN(I,J) * XC * N5(I,J): NUMBER OF Y UPPER LIMITS IN THE BIN(I,J) * XC * N6(I,J): NUMBER OF X UPPER LIMITS IN THE BIN(I,J) * XC * N7(I,J): NUMBER OF DOUBLE UPPER LIMITS IN THE BIN(I,J) * XC * N8(I,J): NUMBER OF Y UPPER, X LOWER LIMITS IN THE * XC * BIN(I,J) * XC * XB(I) : COORDINATE OF CENTER OF THE BIN IN X * XC * YB(I) : COORDINATE OF CENTER OF THE BINS IN Y * XC * IF ISKIP=0, THE NEXT VALUES ARE OUTPUTS : * XC * DELX : BIN SIZE OF X AXIS * XC * DELY : BIN SIZE OF Y AXIS * XC * XORIG : ORIGIN OF X * XC * YORIG : ORIGIN OF Y * XC * * XC * SUBROUTINES * XC * SORT1 * XC * * XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X X DIMENSION INDEX(NTOT),LP(NTOT),XT(NTOT),YT(NTOT),Z(MVAR,NTOT) X DIMENSION X(NTOT),Y(NTOT),NP(NTOT),XB(IB),YB(IB),SWRK1(MVAR) X DIMENSION F(IB,IB),N(IB,IB),N1(IB,IB),N2(IB,IB),N3(IB,IB) X DIMENSION N4(IB,IB),N5(IB,IB),N6(IB,IB),N7(IB,IB),N8(IB,IB) X COMMON /C1/NC1,NC2,NC3,NC4,NC5,NC6,NC7,NC8 XC XC XC * SUBSTITUE NP, X, AND Y TO LP, XT, AND YT SO THAT THE ORIGINAL DATA* XC * WON'T BE CHANGED. * XC X DO 100 J=1,NTOT X LP(J)=NP(J) X XT(J)=X(J) X YT(J)=Y(J) X Z(1,J)=1.0 X 100 CONTINUE XC XC * CALL THE SUBROUTINE SORT1, AND FIND MIN. AND MAX. OF X AND Y. * XC * IF ISKIP=0, THE ORIGIN AND BIN SIZES ARE ALREADY GIVEN * XC X IF(ISKIP.EQ.0) THEN XC XC * SORTING X * XC X CALL SORT1(LP,Z,XT,NTOT,1,INDEX,SWRK1,MVAR) XC XC * SORTING Y * XC X CALL SORT1(LP,Z,YT,NTOT,1,INDEX,SWRK1,MVAR) XC XC * FIND THE SIZES OF BINS * XC X DELX=XT(NTOT)-XT(1) X DELY=YT(NTOT)-YT(1) X DELX=DELX/FLOAT(MX-2) X DELY=DELY/FLOAT(MY-2) XC XC * FIND THE ORIGIN OF THE GRID * XC X XORG=XT(1)-1.5*DELX X YORG=YT(1)-1.5*DELY X ENDIF XC XC XC * INITIALIZE N, N1,....,N8, AND F * XC X DO 300 I=1,MX X DO 200 J=1,MY X N(I,J) =0 X N1(I,J)=0 X N2(I,J)=0 X N3(I,J)=0 X N4(I,J)=0 X N5(I,J)=0 X N6(I,J)=0 X N7(I,J)=0 X N8(I,J)=0 X F(I,J)=0.0 X 200 CONTINUE X 300 CONTINUE XC X DO 390 I=1,NTOT XC XC * FIND POSITION OF I-TH DATA POINT IN THE GRID AND COUNT * XC * NUMBERS OF N,N1,N2,.....,N8. * XC X IP=INT((X(I)-XORG)/DELX)+1 X JP=INT((Y(I)-YORG)/DELY)+1 X XC XC * FOR CONVENIENCE CENSORED POINTS ARE ASSIGNED TO THE NEXT BIN * XC XC XC * DETECTIONS * XC X IF(NP(I).EQ.0) THEN X N(IP,JP)=N(IP,JP)+1 X XC XC * Y LOWER LIMITS * XC X ELSEIF(NP(I).EQ.1) THEN X N1(IP,JP+1)=N1(IP,JP+1)+1 X XC XC * X LOWER LIMITS * XC X ELSEIF(NP(I).EQ.2) THEN X N2(IP+1,JP)=N2(IP+1,JP)+1 X XC XC * DOUBLE LOWER LIMITS * XC X ELSEIF(NP(I).EQ.3) THEN X N3(IP+1,JP+1)=N3(IP+1,JP+1)+1 X XC XC * Y LOWER LIMITS, X UPPER LIMITS * XC X ELSEIF(NP(I).EQ.4) THEN X N4(IP+1,JP-1)=N4(IP+1,JP-1)+1 X XC XC * Y UPPER LIMITS * XC X ELSEIF(NP(I).EQ.-1) THEN X N5(IP,JP-1)=N5(IP,JP-1)+1 X XC XC * X UPPER LIMITS * XC X ELSEIF(NP(I).EQ.-2) THEN X N6(IP-1,JP)=N6(IP-1,JP)+1 X XC XC * DOUBLE UPPER LIMITS * XC X ELSEIF(NP(I).EQ.-3) THEN X N7(IP-1,JP-1)=N7(IP-1,JP-1)+1 X XC XC * Y UPPER LIMITS, X LOWER LIMITS * XC X ELSEIF(NP(I).EQ.-4) THEN X N8(IP-1,JP+1)=N8(IP-1,JP+1)+1 X X ELSE X PRINT *,' THE CENSORSHIP INDICATOR IS NOT RECOGNIZED' X RETURN X ENDIF X 390 CONTINUE XC XC * SET THE COORDINATES OF THE EACH BIN * XC X DO 410 I=1,MX X XB(I)=XORG+DELX/2.0+DELX*(I-1) X 410 CONTINUE X X DO 420 I=1,MY X YB(I)=YORG+DELY/2.0+DELY*(I-1) X 420 CONTINUE XC XC * START CHECKING THE RELATION BETWEEN CENSORED POINTS AND * XC * DETECTED POINTS. IF THE CENSORED POINTS ARE LOCATED SO * XC * THAT THEY CANNOT GIVE WEIGHT TO DETECTED POINTS, THE * XC * CENSORED POINTS ARE CHANGED TO DETECTIONS. * XC X M1=0 X M2=0 X M3=0 X M4=0 X M5=0 X M6=0 X M7=0 X M8=0 XC XC XC * Y LOWER LIMITS * XC X X IF(NC1.NE.0) THEN X DO 600 I=1,MX X DO 500 J=1,MY X JJ=MY-J+1 X IF(N1(I,JJ).NE.0) THEN X K=JJ X 450 IF(N(I,K).EQ.0) THEN X K=K+1 X IF(K.LE.MY) GOTO 450 X M1=M1+N1(I,JJ) X N(I,JJ)=N(I,JJ)+N1(I,JJ) X N1(I,JJ)=0 X ENDIF X ENDIF X 500 CONTINUE X 600 CONTINUE X ENDIF X XC XC XC * X LOWER LIMITS * XC X IF(NC2.NE.0) THEN X DO 800 J=1,MY X DO 700 I=1,MX X II=MX-I+1 X IF(N2(II,J).NE.0) THEN X L=II X 650 IF(N(L,J).EQ.0) THEN X L=L+1 X IF(L.LE.MX) GOTO 650 X M2=M2+N2(II,J) X N(II,J)=N(II,J)+N2(II,J) X N2(II,J)=0 X ENDIF X ENDIF X 700 CONTINUE X 800 CONTINUE X ENDIF X XC XC XC * DOUBLE LOWER LIMITS * XC X IF(NC3.NE.0) THEN X DO 1000 I=1,MX X II=MX-I+1 X DO 950 J=1,MY X JJ=MY-J+1 X IF(N3(II,JJ).NE.0) THEN X L=II X 850 K=JJ X 900 IF(N(II,JJ).EQ.0) THEN X K=K+1 X IF(K.LE.MY) GOTO 900 X L=L+1 X IF(L.LE.MX) GOTO 850 X M3=M3+N3(II,JJ) X N(II,JJ)=N(II,JJ)+N3(II,JJ) X N3(II,JJ)=0 X ENDIF X ENDIF X 950 CONTINUE X 1000 CONTINUE X ENDIF X XC XC XC * Y LOWER LIMITS, X UPPER LIMITS * XC X IF(NC4.NE.0) THEN X DO 1300 I=1,MX X II=MX-I+1 X DO 1200 J=1,MY X IF(N4(II,J).NE.0) THEN X L=II X 1050 K=J X 1100 IF(N(L,K).EQ.0) THEN X K=K-1 X IF(K.GE.1) GOTO 1100 X L=L+1 X IF(L.LE.MX) GOTO 1050 X M4=M4+N4(II,J) X N(II,J)=N(II,J)+N4(II,J) X N4(II,J)=0 X ENDIF X ENDIF X 1200 CONTINUE X 1300 CONTINUE X ENDIF X XC XC XC * Y UPPER LIMITS * XC X IF(NC5.NE.0) THEN X DO 1600 I=1,MX X DO 1500 J=1,MY X IF(N5(I,J).NE.0) THEN X K=J X 1450 IF(N(I,K).EQ.0) THEN X K=K-1 X IF(K.GE.1) GOTO 1450 X M5=M5+N5(I,J) X N(I,J) = N(I,J) + N5(I,J) X N5(I,J)=0 X ENDIF X ENDIF X 1500 CONTINUE X 1600 CONTINUE X ENDIF X XC XC XC * X UPPER LIMITS * XC X IF(NC6.NE.0) THEN X DO 1800 J=1,MY X DO 1700 I=1,MX X IF(N6(I,J).NE.0) THEN X L=I X 1650 IF(N(L,J).EQ.0) THEN X L=L-1 X IF(L.GE.1) GOTO 1650 X M6=M6+N6(I,J) X N(I,J)=N(I,J)+N6(I,J) X N6(I,J)=0 X ENDIF X ENDIF X 1700 CONTINUE X 1800 CONTINUE X ENDIF X XC XC XC * DOUBLE UPPER LIMITS * XC X IF(NC7.NE.0) THEN X DO 2000 I=1,MX X DO 1950 J=1,MY X IF(N7(I,J).NE.0) THEN X L=I X 1850 K=J X 1900 IF(N(L,K).EQ.0) THEN X K=K-1 X IF(K.GE.1) GOTO 1900 X L=L-1 X IF(L.GE.1) GOTO 1850 X M7=M7+N7(I,J) X N(I,J)=N(I,J)+N7(I,J) X N7(I,J)=0 X ENDIF X ENDIF X 1950 CONTINUE X 2000 CONTINUE X ENDIF X XC XC XC * Y UPPER LIMITS, X LOWER LIMITS * XC X IF(NC8.NE.0) THEN X DO 2300 I=1,MX X DO 2200 J=1,MY X JJ=MY-J+1 X IF(N8(I,JJ).NE.0) THEN X L=I X 2050 K=JJ X 2100 IF(N(L,K).EQ.0) THEN X K=K+1 X IF(K.LE.MY) GOTO 2100 X L=L-1 X IF(L.GE.1) GOTO 2050 X M8=M8+N8(I,JJ) X N(I,JJ) = N(I,JJ)+N8(I,JJ) X N8(I,JJ)=0 X ENDIF X ENDIF X 2200 CONTINUE X 2300 CONTINUE X ENDIF X X XC X MM=M1+M2+M3+M4 XC XC * INITIAL GUESS OF F * XC X SNT=NTOT X DO 2440 I=1,MX X DO 2430 J=1,MY X IF(N(I,J).NE.0) F(I,J)=FLOAT(N(I,J))/SNT X 2430 CONTINUE X 2440 CONTINUE X X RETURN X END X X X XC XC ********************************************************************** XC ******************** SUBROUTINE BJ ********************************* XC ********************************************************************** XC X SUBROUTINE BJ(IND,X,Y,NTOT,TOL,MAX,NVAR,ND,NC,ICENS,OUTPUT, X + ALPHA,SIGMAA, X + IWRK1,IWRK2,IWRK4,IWRK5,IWRK6,IWRK7, X + IWRK8,WRK1,WRK2,WRK3,WRK4,WRK5,WRK6,WRK7,WRK8, X + SWRK1,SWRK2,SWRK3,SWRK4,SWRK5,SWRK6,SWRK7, X + SWRK8,SWRK9,DWRK1,EWRK1,MVAR) XC XC XC * LINEAR REGRESSION WITH CENSORED DATA : BUCKLEY-JAMES METHOD * XC * * XC * USING A NONPARAMETRIC METHOD, THIS PROGRAM CALCULATES * XC * LINEAR REGRESSION COEFFICIENTS FOR DATA WHICH CONTAINS SOME * XC * CENSORED OBSERVATIONS. * XC * * XC * PARAMETERS : * XC * INPUT * XC * NTOT : NUMBER OF OBSERVATIONS * XC * NVAR : NUMBER OF INDEPENDENT VARIABLE * XC * ALPHA : INITIAL REGRESSION COEFFICIENT ESTIMATES * XC * (PLEASE ALWAYS USE 0.0 IN THIS PROGRAM). * XC * TOL : TOLERANCE FOR CONVERGENCE (E.G. 1.0E-05) * XC * MAX : MAXIMUM ITERATION (E.G. 20) * XC * X(J,I) : THE MATRIX CONTAINS THE COEFF.OF J-TH * XC * LOCATION PARAMETER AND I-TH OBSERVATION * XC * Y(I) : DEPENDENT PARAMETER OF I-TH OBSERVATION * XC * IND(I) : INDICATOR OF CENSORED DATA ... * XC * IF IND(I)= 1 : LOWER LIMIT * XC * = 0 : UNCENSORED POINT * XC * =-1 : UPPER LIMIT * XC * OUTPUT * XC * ALPHA(1) : INTERCEPT COEFFICIENT * XC * ALPHA(J) : J>1, J-TH SLOPE COEFFICIENTS * XC * ALPHA(MPLONE) : STANDARD DEVIATION * XC * SIGMAA(J) : STANDARD DEVIATION OF J-TH COEFFICIENT * XC * * XC * SUBROUTINES * XC * BUCKLY * XC XC XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X X DIMENSION IND(NTOT),X(MVAR,NTOT),Y(NTOT),ALPHA(MVAR),SIGMAA(MVAR) X DIMENSION IWRK1(NTOT),IWRK2(NTOT),IWRK4(NTOT) X DIMENSION IWRK5(NTOT),IWRK6(NTOT),IWRK7(NTOT),IWRK8(NTOT) X DIMENSION WRK1(NTOT),WRK2(NTOT),WRK3(NTOT),WRK4(NTOT) X DIMENSION WRK5(NTOT),WRK6(NTOT),WRK7(NTOT),WRK8(NTOT) X DIMENSION SWRK1(MVAR),SWRK2(MVAR),SWRK3(MVAR),SWRK4(MVAR) X DIMENSION SWRK5(MVAR),SWRK6(MVAR),SWRK7(MVAR),SWRK8(MVAR) X DIMENSION SWRK9(MVAR),DWRK1(MVAR,NTOT),EWRK1(MVAR,MVAR) X CHARACTER*9 OUTPUT XC X NVAR1=NVAR+1 XC XC * IF CENSORING IS DUE TO UPPER LIMITS, CHANGE THE SIGNS OF DATA * XC * X(I) AND Y(I) BECAUSE B-J METHOD ASSUMES LOWER LIMITS. * XC X IF(ICENS .LT. 0) THEN X DO 1222 I=1,NTOT X DO 1111 J=1,NVAR X X(J,I)=-X(J,I) X 1111 CONTINUE X Y(I)=-Y(I) X 1222 CONTINUE X ENDIF XC XC XC * BUCKLY : THE SUBROUTINE WHICH PERFORMS THE BUCKLEY AND * XC * JAMES METHOD. * XC XC X CALL BUCKLY(X,Y,ALPHA,IND,TOL,SIGMAA,NTOT,NVAR,ND,NC,MAX,ITE, X + IWRK1,DWRK1,IWRK2,IWRK4,WRK1,WRK2,IWRK5, X + WRK3,WRK4,WRK5,WRK6,WRK7,SWRK1,SWRK2,SWRK3, X + IWRK6,IWRK7,IWRK8,SWRK4,SWRK5,SWRK6,SWRK7, X + SWRK8,SWRK9,EWRK1,WRK8,MVAR) XC XC XC * CORRECT THE SIGNS OF THE DATA TO THE ORIGINAL ONES, IF THE * XC * CENSORING IS UPPER LIMIT. * XC X IF(ICENS.LT.0) THEN X DO 1223 I=1,NTOT X X DO 1113 J=1,NVAR X X(J,I)=-X(J,I) X 1113 CONTINUE X Y(I)=-Y(I) X 1223 CONTINUE XC X ALPHA(1)=-ALPHA(1) XC X ENDIF X X 320 IF(OUTPUT.EQ.' ') THEN X PRINT 1050 X PRINT 1020 X PRINT 1050 XC X PRINT 1200,ALPHA(1) X X DO 452 J=2,NVAR1 X JI=J-1 X PRINT 1250,JI,ALPHA(J),SIGMAA(J) X 452 CONTINUE X X PRINT 1300,ALPHA(NVAR+2) X PRINT 1350,ITE X PRINT 1050 X ELSE X WRITE(60,1050) X WRITE(60,1020) X WRITE(60,1050) XC X WRITE(60,1200) ALPHA(1) X X DO 450 J=2,NVAR1 X JI=J-1 X WRITE(60,1250) JI,ALPHA(J),SIGMAA(J) X 450 CONTINUE X X WRITE(60,1300) ALPHA(NVAR+2) X WRITE(60,1350) ITE X WRITE(60,1050) X ENDIF XC XC X 1020 FORMAT(T5,'LINEAR REGRESSION BY BUCKLEY-JAMES METHOD' ) X 1050 FORMAT(T5,' ') X 1100 FORMAT(T8,'DATA TITLE :',T25,60A1) X 1200 FORMAT(T8,'INTERCEPT COEFF :',F8.4) X 1250 FORMAT(T8,'SLOPE COEFF ',I1,' :',F8.4,T38,'+/-',T41, X + F8.4) X 1300 FORMAT(T8,'STANDARD DEVIATION :',F8.4) X 1350 FORMAT(T8,'ITERATIONS :',I3) X 4000 RETURN X END X XC XC ********************************************************************** XC *********************** SUBROUTINE BUCKLY *************************** XC ********************************************************************** XC X SUBROUTINE BUCKLY(X,Y,ALPHA,IND,TOL,SIGMAA,NTOT, X + NVAR,NU,NC,MAX,ITE,IND2,XX,IPT,IR,ND,TY, X + T,NO,Z,W,WX,ZY,V,TEST,TEST2,BU, X + IWRK1,IWRK2,SWRK1,SWRK2,SWRK3,SWRK4, X + SWRK5,SWRK6,EWRK1,WRK1,MVAR) XC XC XC * THIS IS A SUBPROGRAM WHICH PERFORMS THE BUCKLEY-JAMES * XC * METHOD. THIS SUBROUTINE WAS ADAPTED FROM CODE BY J. HALPERN * XC * (STANFORD UNIVERSITY SCHOOL OF MEDICINE, DEPARTMENT * XC * OF FAMILY, COMMUNITY AND PREVENTIVE MEDICINE.) * XC * * XC * INPUT * XC * X(J,I) : INDEPENDENT VARIABLES * XC * Y(I) : DEPENDENT VARIABLE * XC * IND(I) : INDICATOR OF CENSORING * XC * TOL : TOLERANCE LEVEL * XC * NTOT : NUMBER OF DATA POINTS * XC * NVAR : NUMBER OF INDEPENDENT VARIABLES * XC * NU : NUMBER OF DETECTED POINTS * XC * NC : NUMBER OF CENSORED POINTS * XC * MAX : MAXIMUM ITERATION * XC * * XC * WORK * XC * V(J) : AVERAGE OF J-TH DETECTED INDEPENDENT VARIABLE * XC * BU(J) : VARIANCE OF J-TH DETECTED INDEPENDENT VARIABLE * XC * TEST(J) : STORE OF THE PREVIOUS STEP ESTIMATIONS OF * XC * ALPHA(J) * XC * IR(I) : SORTING ORDER * XC * Z(I) : RESIDUALS * XC * W(I) : KM ESTIMATOR * XC * WX(I) : WEIGHT * XC * * XC * OUTPUT * XC * ALPHA(J) : REGRESSION COEFFICIENTS * XC * SIGMA(J) : ERROR * XC * ITE : ITERATION NUMBER * XC * * XC * SUBROUTINES * XC * SORT1, REGRES * XC XC XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER(I-N) X X DIMENSION IND(NTOT),IND2(NTOT),XX(MVAR,NTOT),IPT(NTOT) X DIMENSION X(MVAR,NTOT),Y(NTOT),IR(NTOT),ND(NTOT),TY(NTOT) X DIMENSION T(NTOT),NO(NTOT),Z(NTOT),W(NTOT),WX(NTOT),ZY(NTOT) X DIMENSION V(MVAR),ALPHA(MVAR),TEST(MVAR),BU(MVAR),SIGMAA(MVAR) X DIMENSION TEST2(MVAR),IWRK1(NTOT),IWRK2(NTOT) X DIMENSION SWRK1(MVAR),SWRK2(MVAR),SWRK3(MVAR),SWRK4(MVAR) X DIMENSION SWRK5(MVAR),SWRK6(MVAR),EWRK1(MVAR,MVAR),WRK1(NTOT) XC XC * INITIALIZATION * XC X ITE=0 X NB=NVAR+1 X DO 343 J=1,NVAR X V(J) =0.0 X BU(J)=0.0 X 343 CONTINUE X X DO 392 IN=1,NB X TEST(IN)=0.0 X TEST2(IN) =0.0 X 392 CONTINUE XC XC * CALCULATE SOME VALUES FOR THE STANDARD DEVIATION * XC X DO 5 I=1,NTOT X IF(IND(I).EQ.0) THEN X DO 63 J=1,NVAR X V(J)=V(J)+X(J,I) X 63 CONTINUE X ENDIF X 5 CONTINUE X X DO 68 J=1,NVAR X V(J)=V(J)/REAL(NU) X 68 CONTINUE X X DO 51 I=1,NTOT X IF(IND(I).EQ.0) THEN X DO 53 J=1,NVAR X BU(J)=BU(J)+(X(J,I)-V(J))**2 X 53 CONTINUE X ENDIF X 51 CONTINUE XC XC * REGRES : SUBPROGRAM FOR LINEAR REGRESSION WITHOUT * XC * CONSIDERING CENSORING STATUS * XC X CALL REGRES(X,Y,NTOT,NVAR,ALPHA,RMUL,SWRK1,SWRK2,WRK1, X + IWRK1,IWRK2,SWRK3,SWRK4,SWRK5,EWRK1,SWRK6,MVAR) XC XC * GET RESIDUALS Z(I) * XC XC * START ITERATION : 2000 LOOP. * XC XC X 2000 DO 31 I=1,NTOT X T(I)=-400.0 X IND2(I)=IND(I) XC X ZS=0.0 X DO 61 J=1,NVAR X JJ=J+1 X ZS=ZS+ALPHA(JJ)*X(J,I) X XX(J,I)=X(J,I) X 61 CONTINUE XC X Z(I)=Y(I)-ZS X 31 CONTINUE XC XC * SORTING .... INCREASING ORDER * XC X CALL SORT1(IND2,XX,Z,NTOT,NVAR,IR,SWRK1,MVAR) XC X DO 311 I=1,NTOT X TY(I)=Y(IR(I)) X ZY(I)=Z(I) X 311 CONTINUE XC XC * THE LARGEST RESIDUAL MUST BE UNCENSORED. * XC X IND2(NTOT)=0 XC XC * ESTIMATE VALUES FOR CENSORED DATA * XC * * XC * TY(I)=YY(I)*DEL+((ALPHA*X+SUM(WXX(K)*Z(K))/(1-W(I)))*(1-DEL) * XC * WHERE * XC * TY : ESTIMATED DEPENDENT VALUE * XC * DEL : IF THE DATA IS UNCENSORED :DEL=1.0 * XC * IF THE DATA IS CENSORED :DEL=0.0 * XC * SUM : SUM OVER UNCENSORED DATA Z(K)0, THESE VALUES MUST BE PROVIDED BY THE USER * XC * IPIRNT : IF 0, NO TWO DIMENSIONAL K-M ESTIMATOR WILL BE * XC * PRINTED * XC * >0, TWO DIMENSIONAL K-M ESTIMATOR WILL BE * XC * PRINTED * XC * * XC * WORKING VARIABLES AND ARRAYS: * XC * NTOT : NUMBER OF DATA POINTS * XC * ND : NUMBER OF DETECTED POINTS * XC * NC1 : NUMBER OF Y LOWER LIMITS * XC * NC2 : NUMBER OF X LOWER LIMITS * XC * NC3 : NUMBER OF DOUBLE LOWER LIMITS * XC * NC4 : NUMBER OF Y LOWER AND X UPPER LIMITS * XC * NC5 : NUMBER OF Y UPPER LIMITS * XC * NC6 : NUMBER OF X UPPER LIMITS * XC * NC7 : NUMBER OF DOUBLE UPPER LIMITS * XC * NC8 : NUMBER OF Y UPPER AND X LOWER LIMITS * XC * ICENS : IF 0, CENSORING IS MIXED * XC * 1, CENSORING IS Y LOWER LIMITS ONLY * XC * -1, CENSORING IS Y UPPER LIMITS ONLY * XC * NYC : NC1+NC2 * XC * NXC : NC3+NC4 * XC * NBC : NC5+NC6+NC7+NC8 * XC * IDO : NXC+NBC * XC * IMUL : INDICATOR OF MULTIVARIATE PROBLEM * XC * XX(J,I) : =X(ICOL,I), EXCEPT FOR MULTI INDEPENDENT VARIABLE * XC * CASE (J=1,NVAR). * XC * IND2(I) : =IND(1,I) * XC * * XC * OUTPUT * XC * COXREG * XC * CHI : GLOBAL CHI-SQUARE * XC * PROB : PROBABILITY FOR NULL * XC * BHK (GNERALIZED KENDALL'S TAU) * XC * Z : DEVIATION * XC * PROB : PROBABILITY FOR NULL * XC * EM ALGORITHM * XC * ALPHA(K) : LINEAR REGRESSION COEFFICIENTS (K=1,NVAR+1) * XC * ALPHA(K+2): STANDARD DEVIATION * XC * SIGMAA(K) : ERROR * XC * ITE : NUMBER OF ITERATION * XC * BUCKLEY-JAMES * XC * ALPHA(K) : LINEAR REGRESSION COEFFICIENTS (K=1,NVAR+1) * XC * ALPHA(K+2): STANDARD DEVIATION * XC * SIGMAA(K) : ERROR * XC * SCHMITT * XC * ALPHA : INTERCEPT COEFFICIENT * XC * BETA : SLOPE COEFFICIENT * XC * ***** ALL OUTPUTS ARE INSIDE OF EACH SUBROUTINE * XC * * XC * SUBROUTINES * XC * DATA1, DATREG, DATA2, MULVAR * XC * * X X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X XC * THIS PARAMETER STATEMENT AND THE ONE IN UNIVAR.F ARE THE ONLY * XC * STATEMENTS THAT NEED TO BE ADJUSTED IF THE USER WISHES TO * XC * ANALYZE DATA SETS OF MORE THAN 500 OBSERVATIONS OR MORE THAN * XC * VARIABLES. * X XC ************************************************************************** X PARAMETER(MVAR=4, NDAT=500, IBIN=50) XC ************************************************************************** X X CHARACTER*1 CHECK,CHAR(4,10) X CHARACTER*7 BB(10),YY X CHARACTER*9 FILE,COMMAND,OUTPUT,COLM X CHARACTER*9 YNAME,DUMMY1 X CHARACTER*80 TITLE X X DIMENSION IND(MVAR,NDAT),X(MVAR,NDAT),Y(NDAT) X DIMENSION IPROG(6),IIND(NDAT) X X DIMENSION DWRK1(MVAR,NDAT),DWRK2(MVAR,NDAT) X DIMENSION DWRK3(MVAR,NDAT),DWRK4(MVAR,NDAT) X DIMENSION DWRK5(MVAR,NDAT),DWRK6(MVAR,NDAT) X DIMENSION DWRK8(MVAR,NDAT) X X DIMENSION EWRK1(4,4),RWRK1(NDAT,MVAR) X X DIMENSION AWRK(5,IBIN) X DIMENSION WWRK1((MVAR+1)+NDAT) X DIMENSION WWRK2((MVAR+1)+NDAT) X DIMENSION VWRK1((MVAR+1)*NDAT) X DIMENSION VWRK2((MVAR+1)*NDAT) X X DIMENSION WRK1(NDAT),WRK2(NDAT),WRK3(NDAT),WRK4(NDAT) X DIMENSION WRK5(NDAT),WRK6(NDAT),WRK7(NDAT),WRK8(NDAT) X DIMENSION WRK9(NDAT),WRK10(NDAT),WRK11(NDAT) X DIMENSION WRK12(NDAT) X X DIMENSION SWRK1(MVAR),SWRK2(MVAR),SWRK3(MVAR) X DIMENSION SWRK4(MVAR),SWRK5(MVAR),SWRK6(MVAR) X DIMENSION SWRK7(MVAR),SWRK8(MVAR),SWRK9(MVAR) X DIMENSION SWRK10(MVAR),SWRK11(MVAR),SWRK17(MVAR) X X DIMENSION LWRK1(MVAR,NDAT), LWRK2(MVAR,NDAT) X DIMENSION LWRK3(MVAR,NDAT) X DIMENSION IWRK1(NDAT),IWRK2(NDAT),IWRK3(NDAT) X DIMENSION IWRK4(NDAT),IWRK5(NDAT),IWRK6(NDAT) X DIMENSION IWRK7(NDAT),IWRK8(NDAT) X X DIMENSION IWRK9(NDAT),CWRK1(IBIN),CWRK2(IBIN) X X DIMENSION IBWRK1(IBIN,IBIN),IBWRK2(IBIN,IBIN) X DIMENSION IBWRK3(IBIN,IBIN),IBWRK4(IBIN,IBIN) X DIMENSION IBWRK5(IBIN,IBIN),IBWRK6(IBIN,IBIN) X DIMENSION IBWRK7(IBIN,IBIN),IBWRK8(IBIN,IBIN) X DIMENSION IBWRK9(IBIN,IBIN) X DIMENSION BWRK1(IBIN,IBIN),BWRK2(IBIN,IBIN) X X LENG = (MVAR+1)+NDAT X LEGWRK = (MVAR+1)*NDAT X X DO 5000 K=1,10 X BB(K)=' X ' X 5000 CONTINUE X YY=' Y ' XC X 6000 PRINT * X PRINT * X PRINT *,' CORRELATION AND REGRESSION CALCULATIONS' X PRINT * X PRINT *,' CORRELATION OPTIONS LINEAR REGRESSION OPTIONS' X PRINT *,' 1. COX HAZARD MODEL 4. EM ALGORITHM WITH ' X PRINT *,' NORMAL DISTRIBUTION' X PRINT *,' 2. GEN. KENDALL`S TAU 5. BUCKLEY-JAMES METHOD' X PRINT *,' 3. GEN. SPEARMAN`S RHO 6. SCHMITT`S BINNING METHOD' X PRINT * X PRINT *,' DATA SETS WITH CENSORING IN ONLY ONE DIRECTION OF THE' X PRINT *,' DEPENDENT VARIABLE CAN USE ALL METHODS.' X PRINT * X PRINT *,' DATA SETS WITH SEVERAL INDEPENDENT AND ONE DEPENDENT' X PRINT *,' VARIABLE CAN USE ONLY THE COX PROPORTIONAL HAZARD' X PRINT *,' MODEL,EM ALGORITHM, OR BUCKLEY-JAMES METHOD. ONLY' X PRINT *,' ONE TYPE OF CENSORING IN THE DEPENDENT VARIABLE IS' X PRINT *,' ALLOWED.' X PRINT * X PRINT *,' IF YOUR DATA SET HAS CENSORED POINTS IN THE ' X PRINT *,' INDEPENDENT VARIABLE AND/OR DUAL CENSORED POINTS,' X PRINT *,' YOU CAN USE ONLY THE GEN. KENDALL`S TAU OR GEN.' X PRINT *,' SPEARMAN`S RHO CORRELATION COEFFICIENT, OR' X PRINT *,' SCHMITT`S BINNED LINEAR REGRESSION.' X PRINT * X 6010 PRINT * XC XC * CHECK WHETHER THE USER WANTS TO USE COMMAND FILE INPUTS. IF SO, * XC * GO TO 6660 * XC X 50 FORMAT(A1) X 1380 FORMAT(A9) XC X OUTPUT=' ' X ICOMM=0 X ICOL=1 X PRINT *,'DO YOU WANT TO READ ALL INFORMATION' X WRITE(6,6020) X 6020 FORMAT(' FROM A COMMAND FILE (Y/N) ? ') X READ(5,50) CHECK X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') GOTO 6660 X IF(CHECK.EQ.'N'.OR.CHECK.EQ.'n') GOTO 6025 X GOTO 6010 XC XC * READ FROM THE TERMINAL * XC X 6025 PRINT * X PRINT *,' OK, LET US READ FROM THE TERMINAL ' XC XC * READ TITLE * XC X 6030 PRINT * X WRITE(6,6040) X 6040 FORMAT('WHAT IS THE TITLE OF THE PROBLEM ? ') X READ(5,6050) TITLE X 6050 FORMAT(A80) XC XC * READ DATA FILE NAME * XC X 6051 PRINT * X WRITE(6,6052) X 6052 FORMAT('WHAT IS THE DATA FILE NAME ? ') X READ(5,1380) FILE XC XC * READ NUMBER OF INDEPENDENT VARIABLES * XC X 6060 PRINT * X WRITE(6,6070) X 6070 FORMAT('HOW MANY INDEPENDENT VARIABLES DO YOU HAVE ? ') X CALL DATA1(NVAR) X IF((NVAR.GE.1).AND.(NVAR.LE.MVAR-2)) GOTO 6080 X PRINT * X PRINT *,' YOUR CHOICE IS NOT ACCEPTABLE. PLEASE TYPE IT AGAIN' X GOTO 6060 XC XC * CALL SUBROUTINE "DATREG" TO READ DATA * XC X 6080 CALL DATREG(NVAR,IND,X,Y,FILE,NTOT,ND,NC1,NC2,NC3,NC4,NC5, X + NC6,NC7,NC8,ICENS,NYC,NXC,NBC,MVAR,NDAT) X DO 6090 I = 1, NTOT X IIND(I) = IND(1,I) X 6090 CONTINUE XC XC * CHECK WHICH METHODS THE USER CAN USE * XC X IDC=NXC+NBC X IF((NVAR.EQ.1).AND.(IDC.EQ.0)) GOTO 6530 X IF((NVAR.NE.1).AND.(IDC.NE.0)) GOTO 6340 X IF((NVAR.EQ.1).AND.(IDC.NE.0)) GOTO 6400 X 6170 PRINT * X WRITE(6,6180) X 6180 FORMAT('IS THIS A MULTIVARIATE PROBLEM (Y/N) ? ') X READ(5,50) CHECK X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') GOTO 6220 X IF(CHECK.EQ.'N'.OR.CHECK.EQ.'n') GOTO 6340 X GOTO 6170 XC XC * DATA SET WITH MORE THAN ONE INDEPENDENT VARIABLES * XC X 6220 PRINT * X PRINT *,' YOU CAN USE THE NEXT METHODS ' X PRINT * X PRINT *,' 1. COX HAZARD METHOD' X PRINT *,' 4. EM ALGORITHM WITH NORMAL DISTRIBUTION' X PRINT *,' 5. BUCKLEY-JAMES METHOD' XC X ICOL=0 X J=1 X 6230 PRINT * X WRITE(6,6240) X 6240 FORMAT('WHICH METHOD DO YOU WANT TO USE ? ') X CALL DATA1(IPROG(J)) X IF((IPROG(J).EQ.1).OR.(IPROG(J).EQ.4).OR.(IPROG(J).EQ.5)) X + GOTO 6245 X GOTO 6230 X 6245 IF(J.EQ.1) GOTO 6260 X J1=J-1 X DO 6250 K=1,J1 X IF(IPROG(K).NE.IPROG(J)) GOTO 6250 X PRINT * X PRINT *,' YOU ALREADY CHOSE THAT METHOD.' X PRINT *,' PLEASE CHOOSE ANOTHER ONE' X GOTO 6230 X 6250 CONTINUE X 6260 IF(J.GE.3) GOTO 6280 X 6265 PRINT * X WRITE(6,6270) X 6270 FORMAT('DO YOU WANT TO USE ANY OTHER METHOD (Y/N) ? ') X READ(5,50) CHECK X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') J=J+1 X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') GOTO 6230 X IF(CHECK.EQ.'N'.OR.CHECK.EQ.'n') GOTO 6280 X GOTO 6265 XC X 6280 NOTEST=J X GOTO 6604 XC XC * FOR THE CASE THAT THE DATA SET CONTAINS MIXED CENSORING * XC * (THAT IS, UPPER AND LOWER LIMITS SIMULTANEOUSLY AND/OR * XC * CENSORING IN BOTH VARIABLES). * XC X 6340 PRINT * X WRITE(6,6350) X 6350 FORMAT('WHICH INDEPENDENT VARIABLE DO YOU WANT TO USE ? ') X CALL DATA1(ICOL) X IF(ICOL.GT.NVAR) GOTO 6340 X IF(ICOL.LE.0) GOTO 6340 XC X 6400 IF(NBC.EQ.0) GOTO 6530 X J=1 X PRINT * X PRINT *,' YOU CAN USE THE FOLLOWING METHODS' X PRINT *,' 2. GEN. KENDALL`S TAU METHOD' X PRINT *,' 3. GEN. SPEARMAN`S RHO METHOD' X PRINT *,' 6. SCHMITT`S BINNING METHOD' X 6410 PRINT * X WRITE(6,6420) X 6420 FORMAT('WHICH METHOD DO YOU WANT TO USE ? ') X CALL DATA1(IPROG(J)) X IF((IPROG(J).EQ.2).OR.(IPROG(J).EQ.3).OR.(IPROG(J).EQ.6)) X + GOTO 6425 X GOTO 6410 X 6425 IF(J.EQ.1) GOTO 6440 X J1=J-1 X DO 6430 K=1,J1 X IF(IPROG(K).NE.IPROG(J)) GOTO 6430 X PRINT * X PRINT *,' YOU ALREADY CHOSE THAT METHOD.' X PRINT *,' PLEASE CHOOSE THE OTHER ONE' X GOTO 6410 X 6430 CONTINUE X 6440 IF(J.EQ.3) GOTO 6600 X 6450 PRINT * X WRITE(6,6460) X 6460 FORMAT('DO YOU WANT TO USE THE OTHER METHOD, TOO (Y/N) ? ') X READ(5,50) CHECK X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') J=J+1 X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') GOTO 6410 X IF(CHECK.EQ.'N'.OR.CHECK.EQ.'n') GOTO 6600 X GOTO 6450 XC XC * FOR THE CASE THAT THE DATA SET CONTAINS ONE INDEPENDENT AND ONE * XC * DEPENDENT VARIABLES AND ONE KIND OF CENSORING IN THE DEPENDENT * XC * VARIABLE. * XC X 6530 PRINT * X PRINT *,' YOU CAN USE THE FOLLOWING METHODS' X PRINT * X PRINT *,' 1. COX HAZARD MODEL 4. EM ALGORITHM WITH' X PRINT *,' NORMAL DISTRIBUTION' X PRINT *,' 2. KENDALL`S TAU 5. BUCKLEY-JAMES REGRESSION' X PRINT *,' 3. SPEARMAN`S RHO', X + ' 6. SCHMITT`S BINNED REGRESSION' X PRINT * X J=1 X 6540 PRINT * X WRITE(6,6550) X 6550 FORMAT('WHICH METHOD DO YOU WANT TO USE ? ') X CALL DATA1(IPROG(J)) X IF(IPROG(J).LT.1) GOTO 6540 X IF(IPROG(J).GT.6) GOTO 6570 X IF(J.EQ.1) GOTO 6580 X J1=J-1 X DO 6560 K=1,J1 X IF(IPROG(K).NE.IPROG(J)) GOTO 6560 X PRINT * X PRINT *,' YOU ALREADY CHOSE THAT METHOD.' X PRINT *,' PLEASE CHOOSE ANOTHER ONE.' X GOTO 6540 X 6560 CONTINUE X 6570 IF(J.GE.6) GOTO 6600 X 6580 PRINT * X WRITE(6,6590) X 6590 FORMAT('DO YOU WANT TO USE ANOTHER METHOD (Y/N) ? ') X READ(5,50) CHECK X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') J=J+1 X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') GOTO 6540 X IF(CHECK.EQ.'N'.OR.CHECK.EQ.'n') GOTO 6600 X GOTO 6580 XC X 6600 NOTEST=J XC XC * NAME THE VARIABLES * XC X 6601 PRINT * X PRINT *,' PLEASE NAME THE VARIABLES : ' X PRINT * X WRITE(6,6602) X 6602 FORMAT('WHAT IS THE NAME OF THE INDEPENDENT VARIABLE ? ') X READ(5,1380) COLM XC X PRINT * X WRITE(6,6603) X 6603 FORMAT('WHAT IS THE NAME OF THE DEPENDENT VARIABLE ? ') X READ(5,1380) YNAME XC X 6604 PRINT * X WRITE(6,6605) X 6605 FORMAT('DO YOU WANT TO PRINT THE ORIGINAL DATA (Y/N) ? ') X READ(5,50) CHECK X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') IDATA=1 X IF(CHECK.EQ.'N'.OR.CHECK.EQ.'n') IDATA=0 X IF((CHECK.EQ.'Y').OR.(CHECK.EQ.'N')) GOTO 6609 X IF((CHECK.EQ.'y').OR.(CHECK.EQ.'n')) GOTO 6609 X GOTO 6604 XC XC * CHECK WHETHER THE USER WANT TO SAVE THE RESULT IN AN OUTPUT FILE * XC X 6609 PRINT * X PRINT *,'DO YOU WANT TO SAVE THE RESULT ' X WRITE(6,6610) X 6610 FORMAT(' IN AN OUTPUT FILE (Y/N) ? ') X READ(5,50) CHECK X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') GOTO 6620 X IF(CHECK.EQ.'N'.OR.CHECK.EQ.'n') GOTO 7116 X GOTO 6609 X 6620 PRINT * X WRITE(6,6630) X 6630 FORMAT('WHAT IS THE NAME OF THE OUTPUT FILE ? ') X READ(5,1380) OUTPUT X GOTO 7116 XC XC XC XC * USE "COMMAND" FILE FOR INPUTS * XC X 6660 PRINT * X WRITE(6,6670) X 6670 FORMAT('WHAT IS THE NAME OF THE COMMAND FILE ? ') X READ(5,1380) COMMAND XC X 6700 OPEN(UNIT=50, FILE=COMMAND, STATUS='OLD', FORM='FORMATTED') X ICOMM=1 XC XC * READ TITLE OF THE PROBLEM ; NAME OF THE DATA FILE * XC X READ(50,6710) TITLE X 6710 FORMAT(A80) X READ(50,1380) FILE XC XC * READ NUMBER OF VARIABLES; WHICH VARIABLE WILL BE USED; AND HOW * XC * MANY METHODS THE USER WANTS TO USE. * XC X READ(50,6720) ((CHAR(I,J),I=1,4),J=1,3) X 6720 FORMAT(20A1) X CALL DATA2(CHAR,1,3,NVAR,LIND) X IF(LIND.EQ.0) GOTO 6750 X 6730 PRINT * X PRINT *,' TOTAL NUMBER OF INDEPENDENT VARIABLES IS NOT CLEAR.' X STOP X 6750 IF(NVAR.LT.1) GOTO 6730 X IF(NVAR.NE.1) GOTO 6760 X ICOL=1 X GOTO 6915 XC X 6760 CALL DATA2(CHAR,50,3,ICOL,LIND) X IF(LIND.EQ.0) GOTO 6900 X 6860 PRINT * X PRINT *,' THE CHOICE OF THE VARIABLE IS NOT CLEAR' X STOP X 6900 IF(ICOL.LE.0) GOTO 6860 X IF(ICOL.GT.NVAR) GOTO 6860 XC XC * CHOICE OF THE METHODS * XC X 6915 CALL DATA2(CHAR,3,3,NOTEST,LIND) X IF(LIND.EQ.0) GOTO 6950 X 6930 PRINT * X PRINT *,' IT IS NOT CLEAR HOW MANY METHODS YOU WANT TO USE ' X STOP X 6950 IF(NOTEST.LE.0) GOTO 6930 X IF(NOTEST.GT.6) GOTO 6930 XC X READ(50,6960) ((CHAR(I,J),I=1,4),J=1,NOTEST) X 6960 FORMAT(30A1) X DO 7020 I=1,NOTEST X CALL DATA2(CHAR,I,NOTEST,IPROG(I),LIND) X IF(LIND.EQ.0) GOTO 7010 X 6970 PRINT * X IF(I.EQ.1) PRINT *,' FIRST PROGRAM NUMBER IS NOT CLEAR' X IF(I.EQ.2) PRINT *,' SECOND PROGRAM NUMBER IS NOT CLEAR' X IF(I.GE.3) WRITE(6,6780) I X 6780 FORMAT(5X,I4,'-TH PROGRAM NUMBER IS NOT CLEAR') X STOP X 7010 IF(IPROG(I).LE.0) GOTO 6970 X IF(IPROG(I).GT.6) GOTO 6970 X 7020 CONTINUE XC XC * READ NAMES OF THE INDEPENDENT AND DEPENDENT VARIABLES. IF THE * XC * PROBLEM HAS MULTI-INDEPENDENT VARIABLES, THESE NAMES WIL BE * XC * IGNORED. * XC X READ(50,7022) COLM,YNAME X 7022 FORMAT(2A9) XC X CLOSE(UNIT=50,STATUS='KEEP') XC XC * CALL SUBROUTINE "DATREG" TO READ DATA * XC X CALL DATREG(NVAR,IND,X,Y,FILE,NTOT,ND,NC1,NC2,NC3,NC4,NC5, X + NC6,NC7,NC8,ICENS,NYC,NXC,NBC,MVAR,NDAT) XC X OPEN(UNIT=50,FILE=COMMAND,STATUS='OLD',FORM='FORMATTED') XC XC * THE NEXT SEVERAL LINES READ IN DUMMY VALUES TO PREVENT READING * XC * THE COMMANDS A SECOND TIME. * XC X READ(50,1380) DUMMY1 X READ(50,1380) DUMMY1 X READ(50,7029) IDUMMY X READ(50,7029) IDUMMY X READ(50,1380) DUMMY1 X 7029 FORMAT(I4) XC XC * CHECK WHETHER THE ASSIGNED METHODS CAN BE USED FOR THE DATA * XC X IF(NVAR.GE.2) GOTO 7070 X IF((NXC.EQ.0).AND.(NBC.EQ.0)) GOTO 7110 XC XC * THE CASE WITH MIXED CENSORING IN DATA * XC X I=1 X 7030 IF(IPROG(I).NE.1) GOTO 7040 X PRINT * X PRINT *,' YOU CANNOT USE COX HAZARD MODEL FOR THIS DATA SET' X PRINT *,' THIS METHOD WILL BE IGNORED' X IPROG(I)=-9 X 7040 IF(IPROG(I).NE.4) GOTO 7050 X PRINT * X PRINT *,' YOU CANNOT USE EM ALGORITHM FOR THIS DATA SET' X PRINT *,' THIS METHOD WILL BE IGNORED' X IPROG(I)=-9 X 7050 IF(IPROG(I).NE.5) GOTO 7060 X PRINT * X PRINT *,' YOU CANNOT USE BUCKLEY-JAMES METHOD FOR THIS' X PRINT *,' DATA SET. THIS METHOD WILL BE IGNORED' X IPROG(I)=-9 X 7060 IF(I.GE.NOTEST) GOTO 7110 X I=I+1 X GOTO 7030 XC XC * THE CASE WITH MORE THAN ONE INDEPENDENT VARIABLES * XC X 7070 I=1 X 7080 IF(IPROG(I).NE.2) GOTO 7085 X PRINT * X PRINT *,' YOU CANNOT USE THE KENDALL`S TAU METHOD FOR' X PRINT *,' THIS DATA SET' X PRINT *,' THIS METHOD WILL BE IGNORED' X IPROG(I)=-9 X 7085 IF(IPROG(I).NE.3) GOTO 7090 X PRINT * X PRINT *,' YOU CANNOT USE SPEARMAN`S RHO FOR THIS DATA SET' X PRINT *,' THIS METHOD WILL BE IGNORED' X IPROG(I)=-9 X 7090 IF(IPROG(I).NE.6) GOTO 7100 X PRINT * X PRINT *,' YOU CANNOT USE SCHMITT`S BINNED REGRESSION' X PRINT *,' THIS METHOD WILL BE IGNORED' X IPROG(I)=-9 X 7100 IF(I.EQ.NOTEST) GOTO 7110 X I=I+1 X GOTO 7080 XC XC * READ PRINT OUT INDICATOR FOR THE DATA * XC X 7110 READ(50,6960) (CHAR(I,1),I=1,4) X CALL DATA2(CHAR,1,1,IDATA,LIND) X IF(LIND.EQ.0) GOTO 7114 X 7112 PRINT * X PRINT *,' THE VALUE FOR "IDATA" IS NOT ACCEPTABLE' X STOP X 7114 IF((IDATA.EQ.0).OR.(IDATA.EQ.1)) GOTO 7115 X GOTO 7112 XC XC * READ OUTPUT FILE NAME * XC X 7115 READ(50,1380) OUTPUT XC XC * CALL SUBROUTINE "MULVAR" TO COMPUTE CORRELATION/REGRESSION PROBLEMS* XC X X 7116 IF(OUTPUT .NE. ' ') OPEN(UNIT=60,FILE=OUTPUT, X + STATUS='NEW' XC THE FOLLOWING LINE SHOULD BE COMMENTED OUT ON ALL MACHINES OTHER THAN XC VAX/VMS MACHINES. XC + ,CARRIAGECONTROL='LIST' X + ) X X X CALL MULVAR(X,Y,IND,NTOT,ICOL,NVAR,NOTEST,IPROG,ICOMM, X + OUTPUT,COLM,FILE,YNAME,TITLE,ND,NYC,ICENS, X + NC1,NC2,NC3,NC4,NC5,NC6,NC7,NC8,MVAR, X + LENG,LEGWRK,IBIN,DWRK1,IWRK9,SWRK17,DWRK2, X + DWRK3,DWRK4,DWRK5,DWRK6,DWRK8,RWRK1, X + EWRK1,AWRK,WWRK1,WWRK2, X + VWRK1,VWRK2,WRK1,WRK2,WRK3,WRK4,WRK5,WRK6, X + WRK7,WRK8,WRK9,WRK10,WRK11,WRK12, X + SWRK1,SWRK2,SWRK3,SWRK4,SWRK5,SWRK6,SWRK7, X + SWRK8,SWRK9,SWRK10,SWRK11,LWRK1,LWRK2,LWRK3, X + IWRK1,IWRK2,IWRK3,IWRK4,IWRK5,IWRK6,IWRK7, X + IWRK8,CWRK1,CWRK2,IBWRK1,IBWRK2,IBWRK3, X + IBWRK4,IBWRK5,IBWRK6,IBWRK7,IBWRK8,IBWRK9, X + BWRK1,BWRK2) XC X IF(IDATA.EQ.0) GOTO 7219 X IF(OUTPUT.NE.' ') WRITE(60,7140) X IF(OUTPUT.NE.' ') WRITE(60,7117) FILE X IF(OUTPUT.EQ.' ') PRINT 7140 X IF(OUTPUT.EQ.' ') PRINT 7117, FILE X 7117 FORMAT(5X,'INPUT DATA FILE : ',A9) X IF(ICOL.NE.0) GOTO 7130 X IF(OUTPUT.NE.' ') WRITE(60,7118) (BB(K),K,K=1,NVAR),YY X IF(OUTPUT.EQ.' ') PRINT 7118,(BB(K),K,K=1,NVAR),YY X 7118 FORMAT(4X,'CENSORSHIP',12(A7,I2,1X)) X DO 7119 I=1,NTOT X IF(OUTPUT.NE.' ') WRITE(60,7120) IIND(I), X + (X(J,I),J=1,NVAR),Y(I) X IF(OUTPUT.EQ.' ') PRINT 7120,IIND(I), X + (X(J,I),J=1,NVAR),Y(I) X 7119 CONTINUE X 7120 FORMAT(7X,I4,3X,10F10.3) X GOTO 7219 X 7130 IF(OUTPUT.NE.' ') WRITE(60,7133) X IF(OUTPUT.EQ.' ') PRINT 7133 X 7133 FORMAT(5X,' CENSORSHIP X Y') X DO 7134 I=1,NTOT X IF(OUTPUT.NE.' ') WRITE(60,7135) IIND(I),X(ICOL,I),Y(I) X IF(OUTPUT.EQ.' ') PRINT 7135,IIND(I),X(ICOL,I),Y(I) X 7134 CONTINUE X 7135 FORMAT(8X,I4,5X,2F10.3) X 7140 FORMAT(' ') XC X 7219 IF(OUTPUT .NE. ' ') CLOSE(UNIT=60) X PRINT * X PRINT * X PRINT *,' COMPUTATIONS FOR CORRELATION/REGRESSION' X PRINT *,' PROBLEMS ARE FINISHED' X 7220 PRINT * X WRITE(6,7230) X 7230 FORMAT('DO YOU WANT TO DO ANY OTHER ANALYSIS (Y/N) ? ') X READ(5,50) CHECK X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') IBACK=1 X IF(CHECK.EQ.'Y'.OR.CHECK.EQ.'y') RETURN X IF(CHECK.EQ.'N'.OR.CHECK.EQ.'n') STOP X GOTO 7220 X END X XC XC ********************************************************************** XC ********************* SUBROUTINE CHOL ******************************* XC ********************************************************************** XC X SUBROUTINE CHOL(A,N,U,NULLTY,NA,NU,IFAULT) XC XC * ALGORITHM AS 6 J.R.STATIST.SOC.C.(1968) VOL.17, NO.2 * XC * * XC * GIVEN A SYMMETRIC MATRIX OF ORDER N AS A LOWER TRIANGLE * XC * IN A( ), CALCULATE AN UPPER TRIANGLE, U( ), SUCH THAT * XC * UPRIME*U=A. U( ) MAY COINCIDE WITH A( ). A( ) MUST BE * XC * POSITIVE SEMIDEFINITE. * XC * ETA IS SET TO MULTIPLYING FACTOR DETERMINING THE * XC * EFFECTIVE ZERO FOR PIVOT. * XC * NULLTY IS RETURNED AS NO. OF EFFECTIVE ZERO PIVOTS. * XC * IFAULT IS RETURNED AS 1,IF N.LE.0, 2,IF A( ) IS NOT * XC * POSITIVE SEMI-DEFINITE WITHIN THE TOLERANCE BY ETA. * XC * OTHERWISE ZERO. * XC XC * NOTE : VARIABLES NA,NU, HAVE BEEN ADDED TO THE * XC * ARGUMENT LIST AND USED TO DIMENSION TO ARRAYS * XC * A AND U, RESPECTIVELY. (BY WOLYNETZ (1979)) * XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X X DIMENSION A(NA),U(NU) XC X DATA ETA /1.0E-9/ XC XC * THE VALUE OF ETA WILL DEPEND ON THE WORD LENGTH OF * XC * THE COMPUTER BEING USED. * XC X IFAULT=1 X IF(N.GT.0) THEN X IFAULT=2 X NULLTY=0 X J=1 X K=0 X X DO 10 ICOL=1,N X L=0 X X DO 11 IROW=1,ICOL X K=K+1 X W=A(K) X M=J X X DO 12 I=1,IROW X L=L+1 X IF(I.EQ.IROW) GOTO 13 X W=W-U(L)*U(M) X M=M+1 X 12 CONTINUE X X 13 IF(IROW.EQ.ICOL) GOTO 14 X IF(U(L).EQ.0.0) THEN X U(K) = 0.0 X ELSE X U(K)=W/U(L) X ENDIF X 11 CONTINUE X X 14 IF(DABS(W).GE.DABS(ETA*A(K))) THEN X IF(W.LT.0.0) GOTO 100 X U(K)=DSQRT(W) X ELSE X U(K)=0.0 X NULLTY=NULLTY+1 X ENDIF X J=J+ICOL X 10 CONTINUE X X IFAULT=0 X X ENDIF X 100 RETURN X END X XC XC ********************************************************************** XC ************************ SUBROUTINE COEFF *************************** XC ********************************************************************** XC X SUBROUTINE COEFF(I,X,IP,NTOT,ICOEFF,IA,IB,IC,ID,IE,IG,IH,IJ) XC XC * SUBROUTINE WHICH FINDS CONCORDANCE INFORMATION OF * XC * THE QUANTITY X(I). * XC * * XC * INPUT : X(I) : THE QUANTITIY TO BE EXAMINED * XC * IP(I) : CENSORED STATUS OF X(I) * XC * NTOT : NUMBER OF DATA * XC * OUTPUT : ICOEFF(I) : CONCORDANCE INFORMATION: * XC * FOR X(I) AND X(J) WITH IX(J), ICOEFF=-1 * XC * OTHERWISE, ICOEFF= 0 * XC XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X DIMENSION X(NTOT),IP(NTOT),ICOEFF(NTOT) XC X DO 100 J=1,NTOT X ICOEFF(J)=0 X IF(X(I).LT.X(J)) THEN X IF(IP(I).EQ.IA) GOTO 100 X IF(IP(J).EQ.ID) GOTO 100 X IF(IP(I).EQ.IB) GOTO 100 X IF(IP(J).EQ.IE) GOTO 100 X IF(IP(I).EQ.IC) GOTO 100 X IF(IP(J).EQ.IG) GOTO 100 X X ICOEFF(J)=1 X X ELSEIF(X(I).GT.X(J)) THEN X 50 IF(IP(I).EQ.ID) GOTO 100 X IF(IP(J).EQ.IA) GOTO 100 X IF(IP(I).EQ.IE) GOTO 100 X IF(IP(J).EQ.IB) GOTO 100 X IF(IP(I).EQ.IG) GOTO 100 X IF(IP(J).EQ.IC) GOTO 100 X X ICOEFF(J)=-1 X ENDIF X 100 CONTINUE X RETURN X END X XC XC ********************************************************************** XC ********************* SUBROUTINE COXREG ***************************** XC ********************************************************************** XC X SUBROUTINE COXREG(IND,XX,YY,NTOT,NVAR,OUTPUT,ICENS, X + RINFO,SCORE,FINFO,IL,IM,IP,Y,X, X + SWRK1,IWRK1,IWRK2,MVAR) XC XC * THIS PROGRAM COMPUTES A CORRELATION PROBABILITY ACCORDING * XC * TO COX'S (1972) PROPORTIONAL HAZARDS MODEL. * XC * ONLY ONE TYPE OF CENSORING (I.E. LOWER OR UPPER) * XC * IS ALLOWED IN Y, BUT UP TO NVAR INDEPENDENT VARIABLES CAN * XC * BE USED. THE HYPOTHESIS TESTED IS THE ABSENCE OF CORRELATION * XC * BETWEEN THE DEPENDENT VARIABLE AND INDEPENDENT VARIABLES. * XC * THEREFORE, THE REGRESSION COEFFICIENT IN COX MODEL BETA * XC * IS SET TO ZERO. * XC XC * NOTE NOTE NOTE: THE PROBABILITY CALCULATED MAY NOT BE ACCURATE * XC * WHEN THERE ARE A LARGE NUMBER OF TIED OBSERVATIONS (CF. * XC * R. G. MILLER, SURVIVAL ANALYSIS, 1981, PP. 136-7). * XC XC * * XC * INPUT IND(I) : INDICATOR OF CENSORING * XC * 0 : UNCENSORED DATA POINT * XC * 1 : Y(I) IS LOWER LIMIT * XC * -1 : Y(I) IS UPPER LIMIT * XC * XX(J,I) : INDEPENDENT VARIABLES (J=1,..NVAR) * XC * YY(I) : DEPENDENT VARIABLE * XC * NTOT : TOTAL NO. OF OBSERVATIONS * XC * NVAR : NO. OF INDEPENDENT VARIABLES * XC * * XC * WORK DF : DEGREE OF FREEDOM * XC * X(J,I) : =XX(J,I) * XC * Y(I) : =YY(I) * XC * IP(I) : =IND(I) * XC * IL(I) : INDICATOR OF TIES (# OF TIES) * XC * IM(I) : INDICATOR OF TIES (POSITION) * XC * RINFO(I): INFORMATION MATRIX AND ITS INVERSE * XC * MATRIX AFTER CALLING SUBROUTINE * XC * MATINV. * XC * SCORE(I): SCORE VECTOR * XC * * XC * OUTPUT CHI : GLOBAL CHI-SQUARE * XC * PROB : PROBABILITY OF CORRELATION * XC * * XC * SUBROUTINES * XC * SORT1, TIE, MATINV, PCHISQ * XC * * XC * REFERENCE: RUPERT G. MILLER JR., "SURVIVAL ANALYSIS", 1981, * XC * JOHN WILEY & SONS (NY:NY) * XC * * XC XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X X DIMENSION RINFO(MVAR,MVAR),SCORE(MVAR),FINFO(MVAR) X DIMENSION IND(NTOT),IL(NTOT),IM(NTOT),IP(NTOT),Y(NTOT),YY(NTOT) X DIMENSION X(MVAR,NTOT),XX(MVAR,NTOT) X DIMENSION SWRK1(MVAR), IWRK1(MVAR),IWRK2(MVAR) X CHARACTER*9 OUTPUT XC X DF=NVAR XC XC * SUBSTITUTE XX,YY,AND IND TO X, Y, IP TO AVOID ALTERATION OF * XC * THE ORIGINAL DATA * XC X DO 20 I=1,NTOT X IP(I)=IND(I) X Y(I)=YY(I) X IF(ICENS.EQ.-1) Y(I)=-YY(I) XC XC * IF THE OBSERVATION IS CENSORED, ADD A SMALL NUMBER TO AVOID TIES * XC * WITH DETECTED VALUE. * XC X IF(IP(I).NE.0) Y(I)=Y(I)*(1.0+FLOAT(ICENS)*0.0000001) XC X DO 10 J=1,NVAR X X(J,I)=XX(J,I) X IF(ICENS.EQ.-1) X(J,I)=-X(J,I) X 10 CONTINUE X 20 CONTINUE XC XC * SORT Y IN ASCENDING ORDER * XC X CALL SORT1(IP,X,Y,NTOT,NVAR,IL,SWRK1,MVAR) XC XC * CHECK TIED DATA POINTS AND GIVE THEM A SPECIAL FLAG. * XC X CALL TIE(IP,X,Y,NTOT,NVAR,IL,IM,SWRK1,MVAR) XC XC * COMPUTE SCORE VECTOR. DIMENSION IS NVAR * XC X DO 600 I=1,NVAR X SCORE(I)=0.0 XC X DO 500 J=1,NTOT X IF(IP(J).EQ.0) THEN X IF(IL(J).EQ.1) THEN X SUM=0.0 XC X DO 400 K=J,NTOT X SUM=SUM+X(I,K) X 400 CONTINUE XC X JJ=J+IM(J)-1 X XSUM=0.0 X DO 450 KL=J,JJ X XSUM=XSUM+X(I,KL) X 450 CONTINUE XC X DEN=REAL(NTOT+1-J) X SCORE(I)=SCORE(I)+XSUM-IM(J)*SUM/DEN X ENDIF X ENDIF X 500 CONTINUE X 600 CONTINUE XC XC * COMPUTE THE INFORMATION MATRIX. DIMENSION IS NVAR BY NVAR * XC X DO 1000 I=1,NVAR X DO 900 J=I,NVAR X RINFO(I,J)=0.0 XC X DO 800 K=1,NTOT X IF(IP(K).EQ.0) THEN X IF(IL(K).EQ.1) THEN X SUM1=0.0 X SUM2=0.0 X SUM3=0.0 XC X DO 700 L=K,NTOT X SUM1=SUM1+X(I,L) X SUM2=SUM2+X(J,L) X SUM3=SUM3+X(I,L)*X(J,L) X 700 CONTINUE X DEN=NTOT+1-K X RINFO(I,J)=RINFO(I,J)-REAL(IM(K)) X + *(SUM1*SUM2/DEN**2-SUM3/DEN) X ENDIF X ENDIF X 800 CONTINUE X RINFO(J,I)=RINFO(I,J) X 900 CONTINUE X 1000 CONTINUE XC XC * INVERT INFORMATION MATRX RINFO(I,J). THE INVERTED MATRIX * XC * IS STORED IN RINFO(I,J). * XC X CALL MATINV(RINFO,NVAR,DET,IWRK1,IWRK2,MVAR) XC XC * COMPUTE GLOBAL CHI-SQUARE: * XC * * XC * CHI = [SCORE(I)**T] X [RINFO(I,J)**-1] X [SCORE(J)] * XC * WHERE T IS TRANSVERSE. * XC X DO 1200 I=1,NVAR X FINFO(I)=0.0 X DO 1100 K=1,NVAR X FINFO(I)=FINFO(I)+RINFO(I,K)*SCORE(K) X 1100 CONTINUE X 1200 CONTINUE X CHI=0.0 XC X DO 1300 L=1,NVAR X CHI=CHI+FINFO(L)*SCORE(L) X 1300 CONTINUE X XC XC * GET THE REDUCED CHI-SQUARE * XC X RCHI=CHI/DF XC XC * COMPUTE CHI-SQUARE PROBABILITY * XC X PROB=PCHISQ(RCHI,NVAR) XC X IF(OUTPUT.EQ.' ') THEN X PRINT * X PRINT 1450 X PRINT 1400 X PRINT 1600,CHI X PRINT 1651,NVAR X PRINT 1650,PROB X PRINT * X ELSE X WRITE(60,1400) X WRITE(60,1450) X WRITE(60,1400) X WRITE(60,1600) CHI X WRITE(60,1651) NVAR X WRITE(60,1650) PROB X WRITE(60,1400) X ENDIF X 1400 FORMAT(' ') X 1450 FORMAT(5X,'CORRELATION TEST BY COX PROPORTIONAL HAZARD MODEL') X 1600 FORMAT(6X,' GLOBAL CHI SQUARE =',F9.3,' WITH ') X 1651 FORMAT(11X,I5,' DEGREES OF FREEDOM') X 1650 FORMAT(6X,' PROBABILITY =',F10.4) X X RETURN X END X XC XC ********************************************************************** XC ********************** SUBROUTINE DATA1 ****************************** XC ********************************************************************** XC XC * THIS SUBROUTINE'S PURPOSE IS TO ENABLE EASY, FOOL-PROOF * XC * KEYBOARD ENTRY OF INTEGER INPUT DATA. * XC X SUBROUTINE DATA1(INTEG) XC XC * VARIABLE DECLARATIONS * XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X CHARACTER*1 B(4) X INTEGER*4 CUTOFF,I1,I2,INTEG,TOTAL X INTEGER*4 N(4) XC XC XC * READ IN NUMBER AFTER WRITING PROMPT. FIELD SIZE = 4 * XC X X 3 READ(5,1) (B(I1),I1=1,4) X 1 FORMAT(4(A1)) X XC XC * ANALYZE DIGITS OF NUMBER * XC X DO 2 I1=1,4 X IF(B(I1).EQ.'0') N(I1)=0 X IF(B(I1).EQ.'1') N(I1)=1 X IF(B(I1).EQ.'2') N(I1)=2 X IF(B(I1).EQ.'3') N(I1)=3 X IF(B(I1).EQ.'4') N(I1)=4 X IF(B(I1).EQ.'5') N(I1)=5 X IF(B(I1).EQ.'6') N(I1)=6 X IF(B(I1).EQ.'7') N(I1)=7 X IF(B(I1).EQ.'8') N(I1)=8 X IF(B(I1).EQ.'9') N(I1)=9 XC IF((B(I1).EQ.' ').AND.(I1.EQ.1)) PRINT *,'ENTER NUMBER.' X IF((B(I1).EQ.' ').AND.(I1.EQ.1)) GOTO 3 XC X IF((B(I1).NE.'0').AND.(B(I1).NE.'1').AND.(B(I1).NE.'2').AND. X + (B(I1).NE.'3').AND.(B(I1).NE.'4').AND.(B(I1).NE.'5').AND. X + (B(I1).NE.'6').AND.(B(I1).NE.'7').AND.(B(I1).NE.'8').AND. X + (B(I1).NE.'9').AND.(B(I1).EQ.' ')) CUTOFF=I1 XC X IF((B(I1).NE.'0').AND.(B(I1).NE.'1').AND.(B(I1).NE.'2').AND. X + (B(I1).NE.'3').AND.(B(I1).NE.'4').AND.(B(I1).NE.'5').AND. X + (B(I1).NE.'6').AND.(B(I1).NE.'7').AND.(B(I1).NE.'8').AND. X + (B(I1).NE.'9').AND.(B(I1).EQ.' ')) GOTO 4 XC X IF((B(I1).NE.'0').AND.(B(I1).NE.'1').AND.(B(I1).NE.'2').AND. X + (B(I1).NE.'3').AND.(B(I1).NE.'4').AND.(B(I1).NE.'5').AND. X + (B(I1).NE.'6').AND.(B(I1).NE.'7').AND.(B(I1).NE.'8').AND. X + (B(I1).NE.'9').AND.(B(I1).NE.' ')) X + PRINT *,'POSITIVE INTEGER ONLY' XC X IF((B(I1).NE.'0').AND.(B(I1).NE.'1').AND.(B(I1).NE.'2').AND. X + (B(I1).NE.'3').AND.(B(I1).NE.'4').AND.(B(I1).NE.'5').AND. X + (B(I1).NE.'6').AND.(B(I1).NE.'7').AND.(B(I1).NE.'8').AND. X + (B(I1).NE.'9').AND.(B(I1).NE.' ')) GOTO 3 X 2 CONTINUE XC XC * TOTALS UP THE NUMBER * XC X CUTOFF=5 X 4 TOTAL=0 XC X LAST=CUTOFF-1 X DO 5 I2=1,LAST X TOTAL=TOTAL+N(I2)*(10**(CUTOFF-I2-1)) X 5 CONTINUE XC X INTEG=TOTAL X RETURN X END X XC XC ********************************************************************** XC ********************** SUBROUTINE DATA2 ***************************** XC ********************************************************************** XC X SUBROUTINE DATA2(B,J,IT,INTEG,LIND) XC XC * PURPOSE IS TO ENABLE EASY, FOOL-PROOF KEYBOARD ENTRY OF INTEGER * XC * INPUT DATA. * XC * INPUT * XC * B(4,IT) : INPUT IN CHARACTER FORM. THIS WILL BE TESTED * XC * AND CHANGED TO INTEGER * XC * J : J-TH INPUT, <=IT. * XC * IT : DIMENSION OF B * XC * OUTPUT * XC * INTEG : INTEGER OUTPUT * XC * LIND : INDICATOR OF READABILITY OF INPUT * XC * IF LIND=0, B IS SUCCESSFULLY CHANGED TO INTEG * XC * IF LIND=1, B CANNOT BE CONVERTED TO INTEGER * XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X CHARACTER*1 B(4,IT) X INTEGER CUTOFF,TOTAL X INTEGER N(4) XC XC * ANALYZE DIGITS OF NUMBER * XC X LIND=0 X DO 2 I1=1,4 X IF(B(I1,J).EQ.'0') GOTO 3 X IF(B(I1,J).EQ.'1') GOTO 3 X IF(B(I1,J).EQ.'2') GOTO 3 X IF(B(I1,J).EQ.'3') GOTO 3 X IF(B(I1,J).EQ.'4') GOTO 3 X IF(B(I1,J).EQ.'5') GOTO 3 X IF(B(I1,J).EQ.'6') GOTO 3 X IF(B(I1,J).EQ.'7') GOTO 3 X IF(B(I1,J).EQ.'8') GOTO 3 X IF(B(I1,J).EQ.'9') GOTO 3 X IF(B(I1,J).EQ.' ') GOTO 3 X LIND=1 X GOTO 6 XC X 3 IF(B(I1,J).EQ.'0') N(I1)=0 X IF(B(I1,J).EQ.'1') N(I1)=1 X IF(B(I1,J).EQ.'2') N(I1)=2 X IF(B(I1,J).EQ.'3') N(I1)=3 X IF(B(I1,J).EQ.'4') N(I1)=4 X IF(B(I1,J).EQ.'5') N(I1)=5 X IF(B(I1,J).EQ.'6') N(I1)=6 X IF(B(I1,J).EQ.'7') N(I1)=7 X IF(B(I1,J).EQ.'8') N(I1)=8 X IF(B(I1,J).EQ.'9') N(I1)=9 XC X IF((B(I1,J).EQ.' ').AND.(I1.EQ.1)) LIND=1 X IF((B(I1,J).EQ.' ').AND.(I1.EQ.1)) GOTO 6 XC X IF((B(I1,J).NE.'0').AND.(B(I1,J).NE.'1').AND.(B(I1,J).NE.'2') X + .AND.(B(I1,J).NE.'3').AND.(B(I1,J).NE.'4').AND. X + (B(I1,J).NE.'5').AND. (B(I1,J).NE.'6').AND.(B(I1,J).NE.'7') X + .AND.(B(I1,J).NE.'8').AND.(B(I1,J).NE.'9').AND.(B(I1,J).EQ.' ')) X + CUTOFF=I1 XC X IF((B(I1,J).NE.'0').AND.(B(I1,J).NE.'1').AND.(B(I1,J).NE.'2') X + .AND.(B(I1,J).NE.'3').AND.(B(I1,J).NE.'4').AND. X + (B(I1,J).NE.'5').AND.(B(I1,J).NE.'6').AND.(B(I1,J).NE.'7') X + .AND.(B(I1,J).NE.'8').AND.(B(I1,J).NE.'9') X + .AND.(B(I1,J).EQ.' ')) GOTO 4 XC X IF((B(I1,J).NE.'0').AND.(B(I1,J).NE.'1').AND.(B(I1,J).NE.'2') X + .AND.(B(I1,J).NE.'3').AND.(B(I1,J).NE.'4').AND.(B(I1,J).NE.'5') X + .AND.(B(I1,J).NE.'6').AND.(B(I1,J).NE.'7').AND.(B(I1,J).NE.'8') X + .AND.(B(I1,J).NE.'9').AND.(B(I1,J).NE.' ')) LIND=1 X 2 CONTINUE XC XC * TOTAL UP NUMBER * XC X CUTOFF=5 X 4 TOTAL=0 XC X LAST=CUTOFF-1 X DO 5 I2=1,LAST X TOTAL=TOTAL+N(I2)*(10**(CUTOFF-I2-1)) X 5 CONTINUE X INTEG=TOTAL X6 RETURN X END X XC XC ********************************************************************** XC ********************** SUBROUTINE DATREG **************************** XC ********************************************************************** XC X SUBROUTINE DATREG(NVAR,IND,X,Y,FILE,NTOT,ND,NC1,NC2,NC3,NC4, X + NC5,NC6,NC7,NC8,ICENS,NYC,NXC,NBC,NDIM,NDATA) XC XC * THIS SUBROUTINE READS DATA FROM THE FILE FOR CORRELATION/REGRESSION* XC * PROBLEMS * XC * * XC * INPUT NVAR : NUMBER OF THE INDEPENDENT VARIABLE * XC * FILE : NAME OF THE DATA FILE * XC * * XC * OUTPUT IND(1,I) : INDICATOR OF CENSORSHIP * XC * X(J,I) : INDEPENDENT VARIABLES * XC * Y(I) : DEPENDENT VARIABLES * XC * NTOT : NUMBER OF DATA POINTS * XC * ND : NUMBER OF DETECTED POINTS * XC * NC1 : NUMBER OF Y LOWER LIMITS * XC * NC2 : NUMBER OF X LOWER LIMITS * XC * NC3 : NUMBER OF DOUBLE LOWER LIMITS * XC * NC4 : NUMBER OF Y LOWER, X UPPER LIMITS * XC * NC5 : NUMBER OF Y UPPER LIMITS * XC * NC6 : NUMBER OF X UPPER LIMITS * XC * NC7 : NUMBER OF DOUBLE UPPER LIMITS * XC * NC8 : NUMBER OF Y UPPER, X LOWER LIMITS * XC * ICENS : INDICATOR OF CENSORING * XC * NYC : NC1+NC5 * XC * NXC : NC2+NC6 * XC * NBC : NC3+NC4+NC7+NC8 * XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) X X DIMENSION IND(NDIM, NDATA),X(NDIM, NDATA),Y(NDATA) X CHARACTER*9 FILE XC X OPEN(UNIT=40,FILE=FILE, STATUS='OLD', FORM='FORMATTED') XC X NTOT=1 X ND =0 X NC1 =0 X NC2 =0 X NC3 =0 X NC4 =0 X NC5 =0 X NC6 =0 X NC7 =0 X NC8 =0 XC X 10 READ(40,20,END=30) IND(1,NTOT),(X(J,NTOT),J=1,NVAR),Y(NTOT) X 20 FORMAT(I4,11F10.3) XC X IF(IND(1,NTOT).EQ.0) ND =ND+1 X IF(IND(1,NTOT).EQ.1) NC1=NC1+1 X IF(IND(1,NTOT).EQ.2) NC2=NC2+1 X IF(IND(1,NTOT).EQ.3) NC3=NC3+1 X IF(IND(1,NTOT).EQ.4) NC4=NC4+1 X IF(IND(1,NTOT).EQ.-1) NC5=NC5+1 X IF(IND(1,NTOT).EQ.-2) NC6=NC6+1 X IF(IND(1,NTOT).EQ.-3) NC7=NC7+1 X IF(IND(1,NTOT).EQ.-4) NC8=NC8+1 X NTOT=NTOT+1 X X IF(NTOT.GT.NDATA) THEN X WRITE(6,12) NDATA X WRITE(6,14) X 12 FORMAT(' ****WARNING!! THERE ARE MORE THAN',I5,' POINTS. ') X 14 FORMAT(' ARRAYS WILL OVERFLOW UNLESS DIMENSIONS ARE CHANGED') X STOP X ENDIF XC XC * THE NEXT BLOCK OF CODE CHECKS FOR LINES WHICH HAVE ALL INPUTS * XC * EQUAL TO ZERO. THIS WOULD BE THE EFFECT OF A BLANK LINE IN THE * XC * DATA FILE, AND COULD RESULT IN INCORRECT VALUES COMPUTED. * XC X K = 0 X IF(IND(1,NTOT-1) .NE. 0) K =1 X IF(Y(NTOT-1) .NE. 0.0) K = 1 X DO 21 J = 1, NVAR X IF(X(J,NTOT-1) .NE. 0.0) K = 1 X 21 CONTINUE X IF(K .EQ. 0) THEN X WRITE(6,22) X WRITE(6,24) NTOT X ENDIF X X 22 FORMAT(' ') X 24 FORMAT('WARNING: LINE ',I4,' IN THE DATA FILE MAY BE BLANK') X X GOTO 10 XC X 30 NTOT=NTOT-1 X NYC=NC1+NC5 X NXC=NC2+NC6 X NBC=NC3+NC4+NC7+NC8 X NLC=NC1+NC2+NC3+NC4 X NUC=NC5+NC6+NC7+NC8 X ICENS=0 X IF(NLC.EQ.0) ICENS=-1 X IF(NUC.EQ.0) ICENS=1 XC X CLOSE(UNIT=40) X RETURN X END X XC XC ********************************************************************** XC ********************* SUBROUTINE EM ********************************* XC ********************************************************************** XC X SUBROUTINE EM(IND,X,Y,NTOT,TOL,MAXITS,NVAR,IBET,ND,NC,OUTPUT, X + FILE,ALPHA,XX,Y2,W,WCEN,VCOV,WORK,SIGMAA,TOLA, X + LENG,LEGWRK,MVAR) XC XC * MAXIMUM LIKELIHOOD ESTIMATION IN A LINEAR MODEL * XC * FROM CONFINED AND CENSORED NORMAL DATA * XC * * XC * FORMAL PARAMETERS : * XC * * XC * NTOT INPUT : THE NUMBER OF OBSERVATIONS * XC * MPLONE INPUT : THE TOTAL NUMBER OF PARAMETERS TO BE * XC * ESTIMATED (I.E. NB+1) * XC * MAXITS INPUT : THE MAXIMUM NUMBER OF ITERATIONS ALLOWED * XC * IBET INPUT : IF THERE IS DATA WHICH IS CONFINED * XC * BETWEEN TWO VALUES, IBET=1 OTHERWISE 0 * XC * ALPHA INPUT : IF ALPHA(NVAR+2).LE.0.0, THE SUBROUTINE * XC * CALCULATES INITIAL PARAMETER ESTIMATES. * XC * IF >0.0, IT CONTAINS THE INITIAL * XC * ESTIMATE OF THE J-TH LOCATION PARAMETER * XC * FOR J=1,2,.....,NB. * XC * OUTPUT: THE MOST RECENT PARAMETER ESTIMATES * XC * BEFORE EXIT FROM THE SUBROUTINE. * XC * TOL INPUT : CONVERGENCE TO THE MAXIMUM LIKELIHOOD * XC * PARAMETER ESTIMATE IS REACHED IF THE * XC * DIFFERENCE BETWEEN CONSECUTIVE ESTIMATES * XC * OF THE J-TH PARAMETER IS LESS THAN TOL(J) * XC * FOR J=1,2,........,M. * XC * IND INPUT : INDICATOR OF CENSORED DATA * XC * XX INPUT : THE DESIGN MATRIX XX(I,J) CONTAINS THE * XC * COEFFICIENT OF THE J-TH LOCATION PARA- * XC * METER FOR I-TH OBSERVATION. * XC * Y INPUT : IF PP(I)=0, THE I-TH OBSERVATION IS * XC * COMPLETELY SPECIFIED IN Y(I); IF * XC * PP(I)=-1, Y(I) IS LEFT-CENSORED: IF * XC * PP(I)=1, Y(I) IS RIGHT-CENSORED: IF * XC * PP(I)=5, Y(I) IS CONFINED BETWEEN TWO * XC * VALUES. * XC * ROWX WORK : THE NUMBER OF ROWS OF X (=NTOT) * XC * COLX WORK : THE NUMBER OF COLUMNS OF X (=NB) * XC * W WORK : * XC * LENW WORK : = NB+NTOT * XC * WORK WORK : * XC * LENWRK WORK : = NB*NTOT * XC * LEG : DIMENSION SIZE >= LENW * XC * LEGWRK : DIMENSION SIZE >= LWNWRK * XC * MVAR : DIMENSION SIZE >= MPLONE * XC * VCOV OUTPUT: IF THE PROCEDURE CONVERGED TO THE MAX- * XC * IMUM LIKELIHOOD ESTIMATES, THE FIRST * XC * (NB+1)*(NB+1) POSITIONS CONTAIN AN ESTI- * XC * MATE OF THE VARIANCE-COVARIANCE MATRIX * XC * OF THESE ESTIMATES. * XC * SIGMAA OUTPUT: (J,J) COMPONENT CONTAINS THE STANDARD * XC * DEVIATION OF THE J-TH PARAMETER. * XC * IFAULT OUTPUT: FAILURE INDICATOR * XC * ICHECK OUTPUT: IF ERROR ANALYSIS IS NOT COMPLETED, * XC * ICHECK HAS A VALUE : 1 ,OTHERWISE 0 . * XC * * XC * SUBROUTINES * XC * EMALGO * XC * * XC * REF : M.S.WOLYNETZ AS 139 APL.STATIST.VOL.28 195 (1979) * XC * PLUS CORRECTIONS IN LATER ISSUES OF APPLIED STATISTICS.* XC * * XC * SUBROUTINE EMALGO IS COPIED DIRECTLY FROM WOLYNETZ, EXCEPT * XC * FOR A FEW CHANGES TO PRINT AND COMMENT STATEMENTS. * XC XC X IMPLICIT REAL*8 (A-H,O-Z), INTEGER(I-N) X INTEGER ROWX,COLX X X DIMENSION X(MVAR,NTOT),Y(NTOT),IND(NTOT),XX(NTOT,MVAR),Y2(NTOT) X DIMENSION W(LENG),WCEN(LENG),VCOV(LEGWRK),WORK(LEGWRK) X DIMENSION ALPHA(MVAR),SIGMAA(MVAR,MVAR),TOLA(MVAR) X CHARACTER*9 FILE,OUTPUT XC XC XC * INITIALIZATION * XC X MPLONE=NVAR+2 X NB=NVAR+1 X ROWX=NTOT X COLX=MPLONE X LENW=NTOT+MPLONE X LENWRK=NTOT*MPLONE XC X DO 10 I=1,MPLONE X TOLA(I)=TOL X 10 CONTINUE XC XC * * XC * IF THERE ARE DATA POINTS CONFINED BETWEEN TWO VALUES, READ THE * XC * INPUT DATA AGAIN. * XC * * X