\documentstyle[11pt] {report}
\marginparwidth 0pt
\oddsidemargin 0pt
\evensidemargin 0pt
\marginparsep 0pt
\topmargin=0.0in
\textwidth=6.25in
\textheight=8.25in
\pagestyle{plain}
\parskip=5pt
\parindent=30pt
\begin{document}
\thispagestyle{empty}
\centerline{\Huge\bf ASURV}
\bigskip
\centerline{\Huge\bf Astronomy SURVival Analysis }
\bigskip
\bigskip
\bigskip
\bigskip
\centerline{\LARGE\bf A Software Package for Statistical Analysis of }
\bigskip
\centerline{\LARGE\bf Astronomical Data Containing Nondetections}
\bigskip
\bigskip
\bigskip
\bigskip
\begin{tabbing}
xxxxxxxxxxx\= \kill
\>xxxxxxxxxxxxxxx\= \kill
{\Large\bf Developed by: } \\
\> \\
\> {\Large\bf Takashi Isobe } (Center for Space Research, MIT) \\
\> \\
\> {\Large\bf Michael LaValley } (Dept. of Statistics, Penn State)\\
\> \\
\>{\Large\bf Eric Feigelson }(Dept. of Astronomy \& Astrophysics, Penn State) \\
\> \\
\> \\
\> \\
{\Large\bf Available from: } \\
\> \\
\> code@stat.psu.edu \\
\> \\
\>or, \\
\> \\
\> Eric Feigelson \\
\> Dept. of Astronomy \& Astrophysics \\
\> Pennsylvania State University \\
\> University Park PA 16802 \\
\> (814) 865-0162 \\
\> Email: edf@astro.psu.edu (Internet)\\
\end{tabbing}
\bigskip
\bigskip
\bigskip
\bigskip
\bigskip
\centerline{\Large\bf Rev. 1.2, Summer 1992.}
\newpage
\centerline{\Large\bf TABLE OF CONTENTS}
\bigskip
\bigskip
\noindent{\large\bf
\begin{center}
\begin{verbatim}
1 Introduction ............................................ 3
2 Overview of ASURV ....................................... 4
2.1 Statistical Functions and Organization .......... 4
2.2 Cautions and caveats ............................ 6
3 How to run ASURV ........................................ 9
3.1 Data Input Formats .............................. 9
3.2 KMESTM instructions and example ................. 10
3.3 TWOST instructions and example .................. 11
3.4 BIVAR instructions and example .................. 13
4 Acknowledgements ........................................ 21
Appendices ................................................ 22
A1 Overview of survival analysis .................... 22
A2 Annotated Bibliography on Survival Analysis ...... 23
A3 How Rev 1.2 is Different From Rev 0.0 ............ 26
A4 Errors Removed in Rev 1.1 ........................ 28
A5 Errors Removed in Rev 1.2 ........................ 28
A6 Obtaining and Installing ASURV ................... 28
A7 User Adjustable Parameters ....................... 30
A8 List of subroutines used in ASURV Rev 1.2 ........ 32
\end{verbatim}
\end{center}
}
\bigskip
\bigskip
\bigskip
\centerline{\Large\bf NOTE}
\begin{small}
This program and its documentation are provided `AS IS' without
warranty of any kind. The entire risk as the results and performance
of the program is assumed by the user. Should the program prove
defective, the user assume the entire cost of all necessary correction.
Neither the Pennsylvania State University nor the authors of the program
warrant, guarantee or make any representations regarding use of, or the
results of the use of the program in terms of correctness, accuracy
reliability, or otherwise. Neither shall they be liable for any direct or
indirect damages arising out of the use, results of use or inability to
use this product, even if the University has been advised of the possibility
of such damages or claim. The use of any registered Trademark depicted
in this work is mere ancillary; the authors have no affiliation with
or approval from these Trademark owners.
\end{small}
\newpage
\centerline{\Large\bf 1 Introduction}
Observational astronomers frequently encounter the situation where they
observe a particular property (e.g. far infrared emission, calcium line
absorption, CO emission) of a previously defined sample of objects, but
fail to detect all of the objects. The data set then contains nondetections as
well as detections, preventing the use of simple and familiar statistical
techniques in the analysis.
A number of astronomers have recently recognized the existence
of statistical methods, or have derived similar methods, to deal with these
problems. The methods are collectively
called `survival analysis' and nondetections are called
`censored' data points. {\sl\bf ASURV} is a menu-driven stand-alone computer
package designed to assist astronomers in using methods from survival
analysis. Rev. 1.2 of {\sl\bf ASURV} provides the maximum-likelihood
estimator of the censored distribution,
several two-sample tests, correlation tests and linear regressions as
described in our papers in the {\it\bf Astrophysical Journal} (Feigelson
and Nelson, 1985; Isobe, Feigelson, and Nelson, 1986).
No statistical procedure can magically recover information that was
never measured at the telescope. However, there is frequently important
information implicit in the failure to detect some objects which can be
partially recovered under reasonable assumptions. We purposely provide
a variety of statistical tests - each
based on different models of where upper limits truly lie - so that the
astronomer can judge the importance of the different assumptions.
There are also reasons for {\bf not} applying these methods. We describe
some of their known limitations in {\bf \S 2.2}.
Because astronomers are frequently unfamiliar with the field of
statistics, we devote Appendix {\bf \S A1} to background
material. Both general issues concerning censored data and specifics
regarding the methods used in {\sl\bf ASURV} are discussed. More
mathematical presentations can be found in the references given in Appendix
{\bf \S A2}. Users of Rev 0.0, distributed between 1988 and 1991, are
strongly encouraged to
read Appendices {\bf \S A3-A5} to be familiar with the changes made in the
package. Appendices {\bf \S A6-A8} are needed only by code installers and
those needing to modify the I/O or array sizes.
Users wishing to get right down to work should read {\bf \S 2.1} to
find the appropriate program, and follow the instructions in {\bf \S 3}.
\bigskip
\bigskip
\centerline{\Large\bf Acknowledging ASURV}
We would appreciate that studies using this package include phrasing
similar to `... using {\sl\bf ASURV} Rev 1.2 ({\bf B.A.A.S.} reference),
which implements the methods presented in ({\bf Ap. J.} reference)', where
the {\bf B.A.A.S.} reference is the most recent published {\sl\bf ASURV}
Software Report (Isobe and Feigelson 1990; LaValley, Isobe and Feigelson
1992) and the {\bf Ap. J.} reference is Feigelson and Nelson (1985) for
univariate problems and Isobe,Feigelson and Nelson (1986) for bivariate
problems. Other general works appropriate for referral include the review
of survival methods for astronomers by Schmitt (1985), and the survival
analysis monographs by Miller (1981) or Lawless (1982).
\newpage
\centerline{\Large\bf 2 Overview of ASURV}
\medskip
\noindent{\large\bf 2.1 Statistical Functions and Organization}
The statistical methods for dealing with censored data might be
divided into a 2x2 grid: parametric $vs.$ nonparametric, and univariate $vs.$
bivariate. Parametric methods assume that the censored survival times
are drawn from a parent population with a known distribution function ($e.g.$
Gaussian, exponential), and this is the principal assumption of survival
models for industrial reliability. Nonparametric models make
no assumption about the
parent population, and frequently rely on maximum-likelihood techniques.
Univariate methods are devoted to determining the characteristics of the
population from which the censored sample is drawn, and comparing such
characteristics for two or more censored samples. Bivariate methods
measure whether the censored property of the sample is correlated with
another property of the objects and, if so, to perform a regression that
quantifies the relationship between the two variables.
We have chosen to concentrate on nonparametric models, since the
underlying distribution of astronomical populations is usually unknown.
The linear regression methods however, are either fully parametric
($e.g.$ EM algorithm regression) or semi-parametric
($e.g.$ Cox regression, Buckley-James regression).
Most of the methods are described in more detail in the astronomical
literature by Schmitt (1985), Feigelson and Nelson (1985) and Isobe et al.
(1986). The generalized Spearman's rho utilized in {\sl\bf ASURV}
Rev 1.2 is derived by Akritas (1990).
The program within {\sl\bf ASURV} giving univariate methods is
{\sl\bf UNIVAR}. Its first subprogram is {\sl\bf KMESTM}, giving the
Kaplan-Meier estimator for the distribution function of a randomly
censored sample. First derived in
1958, it lies at the root of non-parametric survival analysis. It is the
unique, self-consistent, generalized maximum-likelihood estimator for the
population from which the sample was drawn. When formulated in cumulative
form, it has analytic asymptotic (for large N) error bars. The median is
always well-defined, though the mean is not if the lowest point in the
sample is an upper limit. It is identical to the differential
`redistribute-to-the-right' algorithm, independently derived by Avni et al.
(1980) and others, but the differential form does not have simple analytic
error analysis.
The second univariate program is {\sl\bf TWOST}, giving a variety
of ways to test whether two censored samples are drawn from the same parent
population. They are mostly generalizations of standard tests for
uncensored data, such as the Wilcoxon and logrank nonparametric two-sample
tests. They differ in how the censored data are weighted or `scored' in
calculating the statistic. Thus each is more sensitive to differences at
one end or the other of the distribution. The Gehan and logrank tests are
widely used in biometrics, while some of the others are not. The tests
offered in Rev 1.2 differ significantly from those offered in Rev 0.0 and
details of the changes are in {\bf \S A3}.
{\sl\bf BIVAR} provides methods for bivariate data, giving three
correlation tests and three linear regression analyses. Cox hazard model
(correlation test), EM algorithm, and Buckley-James method (linear
regressions) can treat several independent variables if the dependent
variable contains only one kind of censoring ($i.e.$ upper or lower limits).
Generalized Kendall's tau (correlation test), Spearman's rho
(correlation test), and Schmitt's binned linear regression can treat
mixed censoring including censoring in the independent variable, but
can have only one independent variable.
{\sl\bf COXREG} computes the correlation probabilities by Cox's
proportional hazard model and {\sl\bf BHK} computes the generalized
Kendall's tau. {\sl\bf SPRMAN} computes correlation probabilities
based on a generalized version of Spearman's rank order correlation
coefficient. {\sl\bf EM} calculates linear regression
coefficients by the EM algorithm assuming a normal distribution of
residuals, {\sl\bf BJ} calculates the Buckley-James regression with
Kaplan-Meier residuals, and {\sl\bf TWOKM} calculates the binned
two-dimensional Kaplan-Meier distribution and associated linear
regression coefficients derived by Schmitt (1985). A bootstrapping
procedure provides error analysis for Schmitt's method in Rev 1.2. The
code for EM algorithm follows that given by
Wolynetz (1979) except minor changes. The code for Buckley-James method
is adopted from Halpern (Stanford Dept. of Statistics; private
communication). Code for the Kaplan-Meier estimator and some of the
two-sample tests was adapted from that given in Lee (1980). {\sl\bf COXREG},
{\sl\bf BHK}, {\sl\bf SPRMAN}, and the {\sl\bf TWOKM} code were developed
entirely by us.
Detailed remarks on specific subroutines can be found in the comments at
the beginning of each subroutine. The reference for the source of the code
for the subroutine is given there, as well as an annotated list of the
variables used in the subroutine.
\newpage
\noindent{\Large\bf 2.2 Cautions and Caveats}
The Kaplan-Meier estimator works with any underlying distribution
($e.g.$ Gaussian, power law, bimodal), but only if the censoring is `random'.
That is, the probability that the measurement of an object is censored can
not depend on the value of the censored variable. At first glance, this
may seem to be inapplicable to most astronomical problems: we detect the
brighter objects in a sample, so the distribution of upper limits always
depends on brightness. However, two factors often serve to randomize
the censoring distribution. First, the censored variable may not be
correlated with the variable by which the sample was initially
identified. Thus, infrared observations of a sample of radio bright
objects will be randomly censored if the radio and infrared emission are
unrelated to each other. Second, astronomical objects in a sample usually
lie at different distances, so that brighter objects are not always the
most luminous. In cases where the variable of interest is censored at
a particular value ($e.g.$ flux, spectral line equivalent width, stellar
rotation rate) rather than randomly censored, then parametric methods
appropriate to `Type 1' censoring should be used. These are described by
Lawless (1982) and Schneider (1986), but are not included in this package.
Thus, the censoring mechanism of each study MUST be understood
individually to judge whether the censoring is or is not likely to be
random. A good example of this judgment process is provided by
Magri et al. (1988). The appearance of the data, even if the upper limits
are clustered at one end of the distribution, is NOT a reliable measure. A
frequent (if philosophically distasteful) escape from the difficulty of
determining the nature of the censoring in a given experiment is to define
the population of interest to be the observed sample. The Kaplan-Meier
estimator then always gives a valid redistribution of the upper limits,
though the result may not be applicable in wider contexts.
The two-sample tests are somewhat less restrictive than the
Kaplan-Meier estimator, since they seek only to compare two samples rather
than determine the true underlying distribution. Because of this, the
two-sample tests do not require that the censoring patterns of the two samples
are random. The two versions of the Gehan test in {\sl\bf ASURV} assume
that the censoring patterns of the two samples are the same, but
the version with hypergeometric variance is more reliable in case of
different censoring patterns. The logrank test results appear to be
correct as long as the censoring patterns are not very different.
Peto-Prentice seems to be the test least affected by differences in
the censoring patterns. There is little known about the limitations of
the Peto-Peto test. These issues are discussed in Prentice and Marek (1979),
Latta (1981) and Lawless (1982).
Two of the bivariate correlation coefficients, generalized Kendall's tau
and Cox regression, are both known to be inaccurate when many tied values
are present in the data. Ties are particularly common when the data is
binned. Caution should be used in these cases. It is not known how the
third correlation method, generalized Spearman's rho, responds to ties in the
data. However, there is no reason to believe that it is more accurate than
Kendall's tau in this case, and it should also used be with caution in the
presence of tied values.
Cox regression, though widely used in biostatistical applications,
formally applies only if the `hazard rate' is constant; that is, the
cumulative distribution function of the censored variable falls
exponentially with increasing values. Astronomical luminosity functions,
in contrast, are frequently modeled by power law distributions. It is not
clear whether or not the results of a Cox regression are significantly
affected by this difference.
There are a variety of limitations to the three linear regression
methods -- {\sl\bf EM}, {\sl\bf BJ}, and {\sl\bf TWOKM} -- presented here.
First, only Schmitt's binned method implemented in {\sl\bf TWOKM} works when
censoring is present in both variables. Second, {\sl\bf\ EM} requires that
the residuals about the fitted line follow a Gaussian distribution.
{\sl\bf BJ} and {\sl\bf TWOKM} are less restrictive, requiring only that the
censoring distribution about the fitted line is random. Both we and
Schneider (1986) find little difference in the regression
coefficients derived from these two methods. Third, the Buckley-James
method has a defect in that the final solution occasionally oscillates
rather than converges smoothly. Fourth, there is considerable uncertainty
regarding the error analysis of the regression coefficients of all three
models. {\sl\bf\ EM} gives analytic formulae based on asymptotic normality,
while Avni and Tananbaum (1986) numerically calculate and examine the
likelihood surface. BJ gives an analytic formula for the slope only, but it
lies on a weak theoretical foundation. We now provide Schmitt's bootstrap
error analysis for {\sl\bf TWOKM}, although this may be subject to high
computational expense for large data sets. Users may thus wish to run
all methods and evaluate the uncertainty with caution. Fifth, Schmitt's
binned regression implemented in {\sl\bf TWOKM} has a number of drawbacks
discussed by Sadler et al. (1989), including slow or failed convergence
to the two-dimensional Kaplan-Meier distribution, arbitrary choice of
bin size and origin, and problems if either too many or too few bins are
selected. In our judgment, the most reliable linear regression method
may be the Buckley-James regression, and we suggest that Schmitt's regression
be reserved for problems with censoring present in both variables.
If we consider the field of survival analysis from a broader
perspective, we note a number of deficiencies with respect to censored
statistical problems in astronomy (see Feigelson, 1990). Most importantly,
survival analysis assumes the upper limits in a given experiment are
precisely known, while in astronomy they frequently represent n times
the RMS noise level in the experimental detector, where n = 2, 3, 5
or whatever. {\bf It is possible that all existing survival methods will
be inaccurate for astronomical data sets containing many points very close
to the detection limit.} Methods that combine the virtues of survival
analysis (which treat censoring) and measurement error models (which
treat known measurement errors in both uncensored and censored points)
are needed. See the discussion by Bickel and Ritov (1992) on this important
matter. A related deficiency is the absence of
weighted means or regressions associated with censored samples.
Survival analysis also has not yet produced any true multivariate
techniques, such as a Principal Components Analysis that permits
censoring. There also seems to be a dearth of nonparametric
goodness-of-fit tests like the Kolmogorov-Smirnov test.
Finally, we note that this package - {\sl\bf ASURV} - is not unique.
Nearly a dozen software packages for analyzing censored data have been
developed (Wagner and Meeker 1985). Four are part of large multi-purpose
commercial statistical software systems: SAS, SPSS, BMDP, and GLIM.
These packages are available on many university central computers. We have
found BMDP to be the most useful for astronomers (see Appendices to
Feigelson and Nelson 1985, Isobe et al. 1986 for instructions on its use).
It provides most of the functions in {\sl\bf KMESTM} and {\sl\bf TWOST}
as well as a full Cox regression, but no linear regressions. Other packages
(CENSOR, DASH, LIMDEP, STATPAC, STAR, SURVAN, SURVCALC, SURVREG) were written
at various universities, medical institutes, publishing houses and industrial
labs; they have not been evaluated by us.
\newpage
\centerline{\Large\bf 3 How to Run ASURV}
\medskip
\noindent{\large\bf 3.1 Data Input Formats}
{\sl\bf ASURV} is designed to be menu-driven. There are two basic input
structures: a data file, and a command file. The data file is assumed to
reside on disk. For each observed object, it contains the measured value
of the variable of interest and an indicator as to whether it is detected
or not. Listed below are the possible values that the {\bf censoring indicator
} can take on. Upper limits are most common in astronomical applications and
lower limits are most common in lifetime applications.
\begin{verbatim}
For univariate data: 1 : Lower Limit
0 : Detection
-1 : Upper Limit
For bivariate data: 3 : Both Independent and Dependent Variables are
Lower Limits
2 : Only independent Variable is a Lower Limit
1 : Only Dependent Variable is a Lower Limit
0 : Detected Point
-1 : Only Dependent Variable is an Upper Limit
-2 : Only Independent Variable is an Upper Limit
-3 : Both Independent and Dependent Variables are
Upper Limits
\end{verbatim}
The command file may either reside on disk, but is more frequently
given interactively at the terminal.
For the univariate program {\sl\bf UNIVAR}, the format of the data file
is determined in the subroutine {\sl\bf DATAIN}. It is currently set
at 10*(I4,F10.3) for {\sl\bf KMESTM}, where each column represents a
distinct subsample. For {\sl\bf TWOST}, it is set at I4,10*(I4,F10.3),
where the first column gives the sample identification number. Table 1
gives an example of a {\sl\bf TWOST} data file called 'gal2.dat'. It
shows a hypothetical study where six normal galaxies, six starburst galaxies
and six Seyfert galaxies were observed with the IRAS satellite. The variable
is the log of the 60-micron luminosity. The problem we address are the
estimation of the luminosity functions of each sample, and a determination
whether they are different from each other at a statistically significant
level. Command file input formats are given in subroutine {\sl\bf UNIVAR},
and inputs are parsed in subroutines {\sl\bf DATA1} and {\sl\bf DATA2}. All
data input and output formats can be changed by the user as described in
appendix {\bf \S A7}.
For the bivariate program {\sl\bf BIVAR}, the data file format is
determined by the subroutine {\sl\bf DATREG}. It is currently set
at I4,10F10.3. In most cases, only two variables are used with input
format I4,2F10.3. Table 2 shows an example of a bivariate censored problem.
Here one wishes to investigate possible relations between infrared and
optical luminosities in a sample of galaxies. {\sl\bf BIVAR} command file
input formats are sometimes a bit complex. The examples below illustrate
various common cases. The formats for the basic command inputs are given in
subroutine {\sl\bf BIVAR}. Additional inputs for the Spearman's rho
correlation, EM algorithm, Buckley-James method, and Schmitt's method are
given in subroutines {\sl\bf R3}, {\sl\bf R4}, {\sl\bf R5}, and {\sl\bf R6}
respectively.
The current version of {\sl\bf ASURV} is set up for data sets with
fewer than 500 data points and occupies about 0.46 MBy of core memory.
For problems with more than 500 points, the parameter values in the
subroutines {\sl\bf UNIVAR} and {\sl\bf BIVAR} must be changed as described
in appendix {\bf\S A7}. Note that for the generalized Kendall's tau
correlation coefficient, and possibly other programs, the computation time
goes up as a high power of the number of data points.
{\sl\bf ASURV} has been designed to work with data that can
contain either upper limits (right censoring) or lower limits (left
censoring). Most of these procedures assume restrictions on the
type of censoring present in the data. Kaplan-Meier estimation and
the two-sample tests presented here can work with data that has either
upper limits or lower limits, but not with data containing both. Cox
regression, the EM algorithm, and Buckley-James regression can use
several independent variables if the dependent variable
contains only one type of censored point (it can be either upper or lower
limits). Kendall's tau, Spearman's rho, and Schmitt's binned regression can
have mixed censoring, including censoring in the independent variable, but
they can only have one independent variable.
\bigskip
\bigskip
\noindent{\large\bf 3.2 KMESTM Instructions and Example}
Suppose one wishes to see the Kaplan-Meier estimator for the infrared
luminosities of the normal galaxies in Table 1. If one runs {\sl\bf ASURV}
interactively from the terminal, the conversation looks as follows:
\begin{verbatim}
Data type : 1 [Univariate data]
Problem : 1 [Kaplan-Meier]
Command file? : N [Instructions from terminal]
Data file : gal1.dat [up to 9 characters]
Title : IR Luminosities of Galaxies
[up to 80 characters]
Number of variables : 1 [if several variables in data file]
Name of variable : IR [ up to 9 characters]
Print K-M? : Y
Print out differential
form of K-M? : Y
25.0 [Starting point is set to 25]
5 [5 bins set]
2.0 [Bin size is set equal to 2]
Print original data? : Y
Need output file? : Y
Output file : gal1.out [up to 9 characters]
Other analysis? : N
\end{verbatim}
If an output file is not specified, the results will be sent to the
terminal screen.
If a command file residing on disk is preferred, run {\sl\bf ASURV}
interactively as above, specifying 'Y' or 'y' when asked if a command file is
available. The command file would then look as follows:
\begin{verbatim}
gal1.dat ... data file
IR Luminosities of Galaxies ... problem title
1 ... number of variables
1 ... which variable?
IR ... name of the variable
1 ... printing K-M (yes=1, no=0)
1 ... printing differential K-M (yes=1, no=0)
25.0 ... starting point of differential K-M
5 ... number of bins
2.0 ... bin size
1 ... printing data (yes=1, no=0)
gal1.out ... output file
\end{verbatim}
All inputs are read starting from the first column.
The resulting output is shown in Table 3. It shows, for example,
that the estimated normal galaxy cumulative IR luminosity function is 0.0
below log(L) = 26.9, 0.63 $\pm$ 0.21 for 26.9 $<$ log(L) $<$ 28.5,
0.83 $\pm$ 0.15 for 28.5 $<$ log(L) $<$ 30.1, and 1.00 above log(L) = 30.1.
The estimated mean for the sample is 27.8 $\pm$ 0.5. The 'C' beside two
values indicates these are nondetections; the Kaplan-Meier function does
not change at these values.
\bigskip
\bigskip
\noindent{\large\bf 3.3 TWOST Instructions and Example}
Suppose one wishes to see two sample tests comparing the subsamples
in Table 1. If one runs {\sl\bf ASURV} interactively from the terminal, the
conversation goes as follows:
\begin{verbatim}
Data type : 1 [Univariate data]
Problem : 2 [Two sample test]
Command file? : N [Instructions from terminal]
Data file : gal2.dat [up to 9 characters]
Problem Title : IR Luminosities of Galaxies
[up to 80 characters]
Number of variables : 1
[if the preceeding answer is more
than 1, the number of the variable
to be tested is now requested]
Name of variable : IR [up to 9 characters]
Number of groups : 3
Which combination ? : 0 [by the indicators in column one
1 of the data set]
Name of subsamples : Normal [up to 9 characters]
Starburst
Need detail print out ? : N
Print full K-M? : Y
Print out differential
form of K-M? : N
Print original data? : N
Output file? : Y
Output file name? : gal2.out [up to 9 characters]
Other analysis? : N
\end{verbatim}
A command file that performs the same operations goes as follows, after
answering 'Y' or 'y' where it asks for a command file:
\begin{verbatim}
gal2.dat ... data file
IR Luminosities of Galaxies ... title
1 ... number of variables
1 ... which variable?
IR ... name of the variable
3 ... number of groups
0 1 2 ... indicators of the groups
0 1 0 1 ... first group for testing
second group for testing
need detail print out ? (if Y:1, N:0)
need full K-M print out? (if Y:1, N:0)
0 ... printing differential K-M (yes=1, no=0)
0 ... print original data? (if Y:1, N:0)
Normal ... name of first group
Starburst ... name of second group
gal2.out ... output file
\end{verbatim}
The resulting output is shown in Table 4. It shows that the
probability that the distribution of IR luminosities of normal and
starburst galaxies are similar at the 5\% level in both the Gehan and Logrank
tests.
\bigskip
\bigskip
\noindent{\large\bf 3.4 BIVAR Instructions and Example}
Suppose one wishes to test for correlation and perform regression
between the optical and infrared luminosities for the galaxy sample in
Table 2. If one runs {\sl\bf ASURV} interactively from the terminal, the
conversation looks as follow:
\begin{verbatim}
Data type : 2 [Bivariate data]
Command file? : N [Instructions from terminal]
Title : Optical-Infrared Relation
[up to 80 characters]
Data file : gal3.dat [up to 9 characters]
Number of Indep. variable: 1
Which methods? : 1 [Cox hazard model]
another one ? : Y
: 4 [EM algorithm regression]
: N
Name of Indep. variable : Optical
Name of Dep. variable : Infrared
Print original data? : N
Save output ? : Y
Name of Output file : gal3.out
Tolerance level ? : N [Default : 1.0E-5]
Initial estimate ? : N
Iteration limits ? : Y
Max. iterations : 50
Other analysis? : N
\end{verbatim}
The user may notice that the above test run includes several input
parameters specific to the EM algorithm (tolerance level through maximum
iterations). All of the regression procedures, particularly Schmitt's
two-dimensional Kaplan-Meier estimator method that requires specification
of the bins, require such extra inputs.
A command file that performs the same operations goes as follows,
following the request for a command file name:
\begin{verbatim}
Optical-Infrared Relation .... title
gal3.dat .... data file
1 1 2 .... 1. number of independent variables
2. which independent variable
3. number of methods
1 4 .... method numbers (Cox, and EM)
Optical Infrared .... name of indep.& dep
variables
0 .... no print of original data
gal3.out .... output file name
1.0E-5 .... tolerance
0.0 0.0 0.0 0.0 .... estimates of coefficients
50 .... iteration
\end{verbatim}
The resulting output is shown in Table 5. It shows that the
correlation between optical and IR luminosities is significant at a
confidence level P $<$ 0.01\%, and the linear fit is
$log(L_{IR})\alpha(1.05 \pm 0.08)*log(L_{OPT})$. It is important to run all
of the methods in order to get a more complete understanding of the
uncertainties in these estimates.
If you want to use Buckley-James method, Spearman's rho, or Schmitt's
method from a command file, you need the next extra inputs:
\begin{verbatim}
(for B-J method)
1.0e-5 tolerance
30 iteration
(for Spearman's Rho)
1 Print out the ranks used in computation;
if you do not want, set to 0
(for Schmitt's)
10 10 bin # for X & Y
10 if you want to set the binning
information, set it to the positive
integer; if not, set to 0.
1.e-5 tolerance
30 iteration
0.5 0.5 bin size for X & Y
26.0 27.0 origins for X & Y
1 print out two dim KM estm;
if you do not need, set to 0.
100 # of iterations for bootstrap error
analysis; if you do not want error
analysis, set to 0
-37 negative integer seed for random number
generator used in error analysis.
\end{verbatim}
\newpage
\centerline{\Large\bf Table 1}
\bigskip
\bigskip
\centerline{\large\bf Example of UNIVAR Data Input for TWOST}
\bigskip
\centerline{\large\bf IR,Optical and Radio Luminosities of Normal,}
\centerline{\large\bf Starburst and Seyfert Galaxies}
\begin{verbatim}
____
0 0 28.5 |
0 0 26.9 |
0 -1 29.7 Normal galaxies
0 -1 28.1 |
0 0 30.1 |
0 -1 27.6 ____
1 -1 29.0 |
1 0 29.0 |
1 0 30.2 Starburst galaxies
1 -1 32.4 |
1 0 28.5 |
1 0 31.1 ____
2 0 31.9 |
2 -1 32.3 Seyfert galaxies
2 0 30.4 |
2 0 31.8 ____
| | |
#1 #2 #3
---I---I---------I--
Column # 4 8 17
Note : #1 : indicator of subsamples (or groups)
If you need only K-M estimator, leave out this indicator.
#2 : indicator of censoring
#3 : variable (in this case, the values are in log form)
\end{verbatim}
\newpage
\centerline{\Large\bf Table 2}
\bigskip
\bigskip
\centerline{\large\bf Example of Bivariate Data Input for MULVAR}
\bigskip
\centerline{\large\bf Optical and IR luminosity relation of IRAS galaxies}
\begin{verbatim}
0 27.2 28.5
0 25.4 26.9
-1 27.2 29.7
-1 25.9 28.1
0 28.8 30.1
-1 25.3 27.6
-1 26.5 29.0
0 27.1 29.0
0 28.9 30.2
-1 29.9 32.4
0 27.0 28.5
0 29.8 31.1
0 30.1 31.9
-1 29.7 32.3
0 28.4 30.4
0 29.3 31.8
| | |
#1 #2 #3
_________ ______
Optical IR
---I---------I---------I--
Column # 4 13 22
Note : #1 : indicator of censoring
#2 : independent variable (may be more Than one)
#3 : dependent variable
\end{verbatim}
\newpage
\centerline{\Large\bf Table 3}
\bigskip
\centerline{\large\bf Example of KMESTM Output}
\bigskip
\begin{small}
\begin{verbatim}
KAPLAN-MEIER ESTIMATOR
TITLE : IR Luminosities of Galaxies
DATA FILE : gal1.dat
VARIABLE : IR
# OF DATA POINTS : 6 # OF UPPER LIMITS : 3
VARIABLE RANGE KM ESTIMATOR ERROR
FROM 0.000 TO 26.900 1.000
FROM 26.900 TO 28.500 0.375 0.213
27.600 C
28.100 C
FROM 28.500 TO 30.100 0.167 0.152
29.700 C
FROM 30.100 ONWARDS 0.000 0.000
SINCE THERE ARE LESS THAN 4 UNCENSORED POINTS,
NO PERCENTILES WERE COMPUTED.
MEAN= 27.767 +/- 0.515
DIFFERENTIAL KM ESTIMATOR
BIN CENTER D
26.000 3.750
28.000 1.250
30.000 1.000
32.000 0.000
34.000 0.000
-------
SUM = 6.000
(D GIVES THE ESTIMATED DATA POINTS IN EACH BIN)
INPUT DATA FILE: gal1.dat
CENSORSHIP X
0 28.500
0 26.900
-1 29.700
-1 28.100
0 30.100
-1 27.600
\end{verbatim}
\end{small}
\newpage
\centerline{\Large\bf Table 4}
\bigskip
\centerline{\large\bf Example of TWOST Output}
\bigskip
\begin{small}
\begin{verbatim}
******* TWO SAMPLE TEST ******
TITLE : IR Luminosities of Galaxies
DATA SET : gal2.dat
VARIABLE : IR
TESTED GROUPS ARE Normal AND Starburst
# OF DATA POINTS : 12, # OF UPPER LIMITS : 5
# OF DATA POINTS IN GROUP I : 6
# OF DATA POINTS IN GROUP II : 6
GEHAN`S GENERALIZED WILCOXON TEST -- PERMUTATION VARIANCE
TEST STATISTIC = 1.652
PROBABILITY = 0.0986
GEHAN`S GENERALIZED WILCOXON TEST -- HYPERGEOMETRIC VARIANCE
TEST STATISTIC = 1.687
PROBABILITY = 0.0917
LOGRANK TEST
TEST STATISTIC = 1.814
PROBABILITY = 0.0696
PETO & PETO GENERALIZED WILCOXON TEST
TEST STATISTIC = 1.730
PROBABILITY = 0.0837
PETO & PRENTICE GENERALIZED WILCOXON TEST
TEST STATISTIC = 1.706
PROBABILITY = 0.0881
KAPLAN-MEIER ESTIMATOR
DATA FILE : gal2.dat
VARIABLE : Normal
# OF DATA POINTS : 6 # OF UPPER LIMITS : 3
VARIABLE RANGE KM ESTIMATOR ERROR
FROM 0.000 TO 26.900 1.000
FROM 26.900 TO 28.500 0.375 0.213
27.600 C
28.100 C
FROM 28.500 TO 30.100 0.167 0.152
29.700 C
FROM 30.100 ONWARDS 0.000 0.000
SINCE THERE ARE LESS THAN 4 UNCENSORED POINTS,
NO PERCENTILES WERE COMPUTED.
MEAN= 27.767 +/- 0.515
KAPLAN-MEIER ESTIMATOR
DATA FILE : gal2.dat
VARIABLE : Starburst
# OF DATA POINTS : 6 # OF UPPER LIMITS : 2
VARIABLE RANGE KM ESTIMATOR ERROR
FROM 0.000 TO 28.500 1.000
FROM 28.500 TO 29.000 0.600 0.219
29.000 C
FROM 29.000 TO 30.200 0.400 0.219
FROM 30.200 TO 31.100 0.200 0.179
FROM 31.100 ONWARDS 0.000 0.000
32.400 C
PERCENTILES
75 TH 50 TH 25 TH
17.812 28.750 29.900
MEAN= 29.460 +/- 0.460
\end{verbatim}
\end{small}
\newpage
\centerline{\Large\bf Table 5}
\bigskip
\bigskip
\centerline{\large\bf Example of BIVAR Output}
\bigskip
\begin{center}
\begin{verbatim}
CORRELATION AND REGRESSION PROBLEM
TITLE IS Optical-Infrared Relation
DATA FILE IS gal3.dat
INDEPENDENT DEPENDENT
Optical AND Infrared
NUMBER OF DATA POINTS : 16
UPPER LIMITS IN Y X BOTH MIX
6 0 0 0
CORRELATION TEST BY COX PROPORTIONAL HAZARD MODEL
GLOBAL CHI SQUARE = 18.458 WITH
1 DEGREES OF FREEDOM
PROBABILITY = 0.0000
LINEAR REGRESSION BY PARAMETRIC EM ALGORITHM
INTERCEPT COEFF : 0.1703 +/- 2.2356
SLOPE COEFF 1 : 1.0519 +/- 0.0793
STANDARD DEVIATION : 0.3687
ITERATIONS : 27
\end{verbatim}
\end{center}
\newpage
\centerline{\Large\bf 4 Acknowledgements}
The production and distribution of {\sl\bf ASURV Rev 1.2} was
principally funded at Penn State by
a grant from the Center for Excellence in Space Data and Information
Sciences (operated by the Universities Space Research Association in
cooperation with NASA), NASA grants NAGW-1917 and NAGW-2120, and
NSF grant DMS-9007717. T.I. was supported by NASA grant NAS8-37716.
We are grateful to Prof. Michael Akritas (Dept. of Statistics, Penn
State) for advice and guidance on mathematical issues, and
to Drs. Mark Wardle (Dept. of Physics and Astronomy, Northwestern
University), Paul Eskridge (Harvard-Smithsonian Center for Astrophysics),
Eric Smith (Laboratory for Astronomy and Solar Physics, NASA-Goddard
Space Flight Center) and Eric Jensen (Wisconsin)
for generous consultation and assistance on the coding.
We would also like to thank Drs. Peter Allan (Rutherford Appleton Laboratory),
Simon Morris (Carnegie Observatories), Karen Strom (UMASS), and Marco
Lolli (Bologna) for their help in debugging {\sl\bf ASURV Rev 1.0}.
\newpage
\bigskip
\centerline{\Large\bf APPENDICES}
\bigskip
\bigskip
\noindent{\large\bf A1 Overview of Survival Analysis }
Statistical methods dealing with censored data have a long and
confusing history. It began in the 17th century with the development of
actuarial mathematics for the life insurance industry and demographic
studies. Astronomer Edmund Halley was one of the founders. In the
mid-20th century, this application grew along with biomedical and clinical
research into a major field of biometrics. Similar (and sometimes
identical) statistical methods were being developed in two other fields:
the theory of reliability, with industrial and engineering applications;
and econometrics, with attempts to understand the relations between
economic forces in groups of people. Finally, we find the same mathematical
problems and some of the same solutions arising in modern astronomy as
outlined in {\bf \S 1} above.
Let us give an illustration on the use of survival analysis in these
disparate fields. Consider a linear regression problem. First, an
epidemiologist wishes to determine how the human longevity or `survival time'
(dependent variable) depends on the number of cigarettes smoked per day
(independent variable). The experiment lasts 10 years, during which some
individuals die and others do not. The survival time of the living
individuals is only known to be greater than their age when the experiment
ends. These values are therefore `right censored data points'. Second,
an engine manufacturing company engines wishes to know the average time
between breakdowns as a function of engine speed to determine the optimal
operating range. Test engines are set running until 20\% of them fail;
the average `lifetime' and dependence on speed is then calculated with
80\% of the data points right-censored. Third, an astronomer observes a sample
of nearby galaxies in the far-infrared to determine the relation between
dust and molecular gas. Half have detected infrared luminosities but
half are upper limits (left-censored data points). The astronomer then seeks
the relationship between infrared luminosities and the CO traces of molecular
material to investigate star formation efficiency. The CO observations may
also contain censored values.
Astronomers have dealt with their upper limits in a number of
fashions. One is to consider only detections in the analysis; while
possibly acceptable for some purposes (e.g. correlation), this will
clearly bias the results in others (e.g. estimating luminosity functions).
A second procedure considers the ratio of detected objects to observed
objects in a given sample. When plotted in a range of luminosity bins,
this has been called the `fractional luminosity function' and has been
frequently used in extragalactic radio astronomy. This is mathematically
the same as actuarial life tables. But it is sometimes used incorrectly
in comparing different samples: the detected fraction clearly depends on
the (usually uninteresting) distance to the objects as well as their
(interesting) luminosity. Also, simple $sqrt$N error bars do not
apply in these fractional luminosity functions, as is frequently assumed.
A third procedure is to take all of the data, including the exact
values of the upper limits, and model the properties of the parent
population under certain mathematical constraints, such as maximizing
the likelihood that the censored points follow the same distribution as the
uncensored points. This technique, which is at the root of much of modern
survival analysis, was independently developed by astrophysicists (Avni et al.
1980; Avni and Tananbaum 1986) and is now commonly used in observational
astronomy.
The advantage accrued in learning and using survival analysis methods
developed in biometrics, econometrics, actuarial and reliability
mathematics is the wide range of useful techniques developed in these
fields. In general, astronomers can have faith in the mathematical validity
of survival analysis. Issues of great social importance (e.g.
establishing Social Security benefits, strategies for dealing with the
AIDS epidemic, MIL-STD reliability standards) are based on it. In detail,
however, astronomers must study the assumptions underlying each method and
be aware that some methods at the forefront of survival analysis that may
not be fully understood.
\S {\bf A2} below gives a bibliography of selected works concerning
survival analysis statistical methods. We have listed some recent
monographs from the biometric and reliability field that we have found to
be useful (Kalbfleisch and Prentice 1980; Lee 1980; Lawless 1982; Miller
1981; Schneider 1986), as well as one from econometrics (Amemiya 1985).
Papers from the astronomical literature dealing with these methods include
Avni et al. (1980), Schmitt (1985), Feigelson and Nelson (1985), Avni and
Tananbaum (1986), Isobe et al. (1986), and Wardle and Knapp (1986). It is
important to recognize that the methods presented in these papers and in
this software package represent only a small part of the entire body of
statistical methods applicable to censored data.
\bigskip
\bigskip
\noindent{\large\bf A2 Annotated Bibliography}
\begin{description}
\item [] Akritas, M. ``Aligned Rank Tests for Regression With Censored Data'',
Penn State Dept. of Statistics Technical Report, 1989. \\
{\it Source for ASURV's generalized Spearman's rho computation.}
\item [] Amemiya, T. {\bf Advanced Econometrics} (Harvard U. Press:Cambridge
MA) 1985. \\
{\it Reviews econometric `Tobit' regression models including censoring.}
\item [] Avni, Y., Soltan, A., Tananbaum, H. and Zamorani, G. ``A Method for
Determining Luminosity Functions Incorporating Both Flux Measurements
and Flux Upper Limits, with Applications to the Average X-ray to Optical
Luminosity Ration for Quasars", {\bf Astrophys. J.} 235:694 1980. \\
{\it The discovery paper in the astronomical literature for the
differential Kaplan-Meier estimator.}
\item [] Avni, Y. and Tananbaum, H. ``X-ray Properties of Optically Selected
QSOs", {\bf Astrophys. J.} 305:57 1986. \\
{\it The discovery paper in the astronomical literature of the linear
regression with censored data for the Gaussian model.}
\item [] Bickel, P.J and Ritov, Y. ``Discussion: `Censoring in
Astronomical Data Due
to Nondetections' by Eric D. Feigelson'', in {\bf Statistical Challenges
in Modern Astronomy}, eds. E.D. Feigelson and G.J. Babu, (Springer-Verlag:
New York) 1992. \\
{\it A discussion of the possible inadequacies of survival analysis for
treating low signal-to-noise astronomical data.}
\item [] Brown, B .J. Jr., Hollander, M., and Korwar, R. M. ``Nonparametric
Tests of Independence for Censored Data, with Applications to Heart
Transplant Studies" from {\bf Reliability and Biometry}, eds. F. Proschan
and R. J. Serfling (SIAM: Philadelphia) p.327 1974.\\
{\it Reference for the generalized Kendall's tau correlation coefficient.}
\item [] Buckley, J. and James, I. ``Linear Regression with Censored Data",
{\bf Biometrika} 66:429 1979.\\
{\it Reference for the linear regression with Kaplan-Meier residuals.}
\item [] Feigelson, E. D. ``Censored Data in Astronomy'', {\bf Errors,
Bias and Uncertainties in Astronomy}, eds. C. Jaschek and F. Murtagh,
(Cambridge U. Press: Cambridge) p. 213 1990.\\
{\it A recent overview of the field.}
\item [] Feigelson, E. D. and Nelson, P. I. ``Statistical Methods for
Astronomical Data with Upper Limits: I. Univariate Distributions",
{\bf Astrophys. J.} 293:192 1985.\\
{\it Our first paper, presenting the procedures in UNIVAR here.}
\item [] Isobe, T., Feigelson, E. D., and Nelson, P. I. ``Statistical Methods
for Astronomical Data with Upper Limits: II. Correlation and Regression",
{\bf Astrophys. J.} 306:490 1986.\\
{\it Our second paper, presenting the procedures in BIVAR here.}
\item [] Isobe, T. and Feigelson, E. D. ``Survival Analysis, or What To Do with
Upper Limits in Astronomical Surveys", {\bf Bull. Inform. Centre Donnees
Stellaires}, 31:209 1986.\\
{\it A precis of the above two papers in the Newsletter of Working Group
for Modern Astronomical Methodology.}
\item [] Isobe, T. and Feigelson, E. D. ``ASURV'', {\bf Bull. Amer. Astro.
Society}, 22:917 1990.\\
{\it The initial software report for ASURV Rev 1.0.}
\item [] Kalbfleisch, J. D. and Prentice, R. L. {\bf The Statistical Analysis
of Failure Time Data} (Wiley:New York) 1980.\\
{\it A well-known monograph with particular emphasis on Cox regression.}
\item [] Latta, R. ``A Monte Carlo Study of Some Two-Sample Rank Tests With
Censored Data'', {\bf Jour. Amer. Stat. Assn.}, 76:713 1981. \\
{\it A simulation study comparing several two-sample tests available in
ASURV.}
\item [] LaValley, M., Isobe, T. and Feigelson, E.D. ``ASURV'', {\bf Bull.
Amer. Astro. Society} 1992
{\it The new software report for ASURV Rev. 1.1.}
\item [] Lawless, J. F. {\bf Statistical Models and Methods for Lifetime Data}
(Wiley:New York) 1982.\\
{\it A very thorough monograph, emphasizing parametric models and
engineering applications.}
\item [] Lee, E. T. {\bf Statistical Methods for Survival Data Analysis}
(Lifetime Learning Pub.:Belmont CA) 1980.\\
{\it Hands-on approach, with many useful examples and Fortran programs.}
\item [] Magri, C., Haynes, M., Forman, W., Jones, C. and Giovanelli, R.
``The Pattern of Gas Deficiency in Cluster Spirals: The Correlation of
H I and X-ray Properties'', {\bf Astrophys. J.} 333:136 1988. \\
{\it A use of bivariate survival analysis in astronomy, with a
good discussion of the applicability of the methods.}
\item [] Millard, S. and Deverel, S. ``Nonparametric Statistical Methods for
Comparing Two Sites Based on Data With Multiple Nondetect Limits'',
{\bf Water Resources Research}, 24:2087 1988. \\
{\it A good introduction to the two-sample tests used in ASURV, written
in nontechnical language.}
\item [] Miller, R. G. Jr. {\bf Survival Analysis} (Wiley, New York) 1981.\\
{\it A good introduction to the field; brief and informative lecture notes
from a graduate course at Stanford.}
\item [] Prentice, R. and Marek, P. ``A Qualitative Discrepancy Between
Censored Data Rank Tests'', {\bf Biometrics} 35:861 1979. \\
{\it A look at some of the problems with the Gehan two-sample test, using
data from a cancer experiment.}
\item[] Sadler, E. M., Jenkins, C. R. and Kotanyi, C. G..
``Low-Luminosity Radio Sources in Early-Type Galaxies'',
{\bf Mon. Not. Royal Astro. Soc.} 240:591 1989. \\
{\it An astronomical application of survival analysis, with
useful discussion of the methods in the Appendices.}
\item [] Schmitt, J. H. M. M. ``Statistical Analysis of Astronomical Data
Containing Upper Bounds: General Methods and Examples Drawn from X-ray
Astronomy", {\bf Astrophys. J.} 293:178 1985.\\
{\it Similar to our papers, a presentation of survival analysis for
astronomers. Derives {\sl\bf TWOKM} regression for censoring in both
variables.}
\item [] Schneider, H. {\bf Truncated and Censored Samples from Normal
Populations} (Dekker: New York) 1986.\\
{\it Monograph specializing on Gaussian models, with a good discussion of
linear regression.}
\item [] Wagner, A. E. and Meeker, W. Q. Jr. ``A Survey of Statistical Software
for Life Data Analysis", in {\bf 1985 Proceedings of the Statistical
Computing Section}, (Amer. Stat. Assn.:Washington DC), p. 441.\\
{\it Summarizes capabilities and gives addresses for other software
packages.}
\item [] Wardle, M. and Knapp, G. R. ``The Statistical Distribution of the
Neutral-Hydrogen Content of S0 Galaxies", {\bf Astron. J.} 91:23 1986.\\
{\it Discusses the differential Kaplan-Meier distribution and its error
analysis.}
\item [] Wolynetz, M. S. ``Maximum Likelihood Estimation in a Linear Model
from Confined and Censored Normal Data", {\bf Appl. Stat.} 28:195
1979.\\
{\it Published Fortran code for the EM algorithm linear regression.}
\end{description}
\bigskip
\bigskip
\noindent{\Large\bf A3 Rev 1.1 Modifications and Additions}
Each of the three major areas of {\sl\bf ASURV}; the KM
estimator, the two-sample tests, and the bivariate methods have been
updated in going from Rev 0.0 to Rev 1.1.
\noindent{\large\bf A3.1 KMESTM}
In the Survival Analysis literature, the value of the survival function
at a x-value is the probability that a given observation will be at least as
large as that x-value. As a result, the survival curve starts with a
value of one and declines to zero as the x-values increase. The
KM estimate of the survival curve should mirror this behavior by
starting at one and declining to zero as more and more of the observations
are passed. In Rev 0.0, the KM estimate for right-censored (lower limits)
data was given in that form, but the KM estimate for left-censored
(upper limits) data started at zero and increased to one. As a result,
the x-values where jumps in the KM estimate occurred were correct, and the
jumps were of the correct height, but the reported survival value at most
points were (in a strict sense) incorrect. In Rev 1.1, this has been
corrected so that the KM estimate will move in the proper direction for both
upper limits and lower limits data.
A differential, or binned, Kaplan-Meier estimate has also been added to
the package. This allows the user to find the number of points falling into
specified bins along the X-axis according to the Kaplan-Meier estimated
survival curve. However, astronomers are strongly encouraged to use the
integral KM for which analytic error analysis is available. There is
{\bf no known error analysis} for the differential KM.
\noindent{\large\bf A3.2 TWOST}
In Rev 0.0 the code for two-sample tests relied heavily on code published
in {\bf Statistical Methods for Survival Data Analysis} by E.T. Lee. Since the
publication of this book in 1980, there have been several articles and
simulation studies done on the various two-sample tests. Lee's book
uses a permutation variance for its tests, which assumes that both
groups being considered are subject to the same censoring pattern. Tests
with hypergeometric variance form seem to be more `robust' against differences
in the censoring patterns, and some statistical software packages
($e.g.$ SAS) have replaced permutation variances with hypergeometric variances.
We have also realized that Rev 0.0 presented Lee's Peto-Peto test, rather
than the Peto-Prentice test described in Feigelson and Nelson (1985) and the
Rev. 0.0 {\sl\bf ASURV} manual.
In light of these developments we have modified the two-sample tests
calculated in {\sl\bf ASURV}. Rev 1.1 offers two versions of
Gehan's test: one with permutation variance (which will match the Gehan's
test value from Rev 0.0) and one with hypergeometric variance. The logrank
test now uses a hypergeometric variance, but the Peto-Peto test still uses a
permutation variance. The Cox-Mantel test has been removed as it was very
similar to the logrank test, and the Peto-Prentice test has been added. The
Peto-Prentice test uses an asymptotic variance form that has been shown to
do very well in simulation studies (Latta, 1981).
{\sl\bf ASURV} now automatically does all the two-sample tests,
instead of asking the user to specify which tests to run. These tests
are not very time consuming and the user will do well to consider
the results of all the tests. If the p-values differ significantly,
then the Peto-Prentice test is probably the most reliable (Prentice
and Marek, 1979). The two-sample tests all use different, but
reasonable, weightings of the data points, so large discrepancies
between the results of the tests indicates that caution should be
used in drawing conclusions based on this data.
\noindent{\large\bf A3.3 BIVAR}
The bivariate methods have been extended in two ways. First,
standard error estimates for the slope and intercept are offered for
Schmitt's method of linear regression. These error estimates are
based upon the bootstrap, a statistical technique which randomly
resamples the set of data points, with replacement, many times
and then runs Schmitt's procedure on the artificial data sets created by the
resampling. Two hundred iterations is considered sufficient to get reasonable
estimates of the standard error of the estimate in most situations.
However, this might be computationally intensive for a large data set.
The bivariate methods have also been extended by an
additional correlation routine, a generalized Spearman's rho procedure
(Akritas 1990). The usual Spearman's rho correlation estimate for uncensored
data is simply the correlation between the ranks of the independent and
dependent variables. In order to extend the procedure to censored data, the
Kaplan-Meier estimate of the survival curve is used to assign ranks to the
observations. Consequently, the ranks assigned to the observations may not
be whole numbers. Censored points are assigned half (for left-censored)
or twice (for right-censored) the rank that they would have had were
they uncensored. This method is based on preliminary findings and
has not been carefully scrutinized either theoretically or empirically. It is
offered in the code to serve as a less computer intensive approximation to the
Kendall's tau correlation, which becomes computationally infeasible for large
data sets (n $>$ 1000). {\it The generalized Spearman's Rho routine is not
dependable for small data sets (n $<$ 30)}. In that situation the generalized
Kendall's tau routine should be relied upon. It should also be noted that
the test statistics reported by {\sl\bf ASURV} for Kendall's tau and Spearman's
rho are not directly comparable. The test statistic reported for Spearman's
rho is an estimated correlation, and the test statistic given for Kendall's
tau is an estimated function of the correlation. It is the reported
probabilities that should be compared.
\noindent{\large\bf A3.4 Interface, Outputs and Manual}
The screen-keyboard interface has been streamlined somewhat. For
example. the user is now provided all two-sample tests without requesting
them individually. New inputs have been introduced for the new programs
(differential Kaplan-Meier function, the generalized Spearman's rho, and
error analysis for Schmitt's regression). The printed outputs of the
programs have been clarified in several places where ambiguities were
reported. For example, the Kaplan-Meier estimator now specifies whether the
change occurs before or after a given value, the meaning of the correlation
probabilities is now given explicitly, and warnings are printed when there
is good reason to suspect the results of a test are unreliable.
The Users Manual has been reorganized so that material that is not actually
needed to operate the program has been located in several Appendices.
\bigskip
\bigskip
\noindent{\Large\bf A4 Errors Removed in Rev 1.1}
Since {\sl\bf ASURV} was released in November 1991, several errors
have been discovered in the package by users and have been reported to us.
In revision 1.1, all the bugs that have been reported are corrected. We
have also taken the opportunity to correct some subtle programming errors that
we came across in the code. The major errors were:
\begin{itemize}
\item The command files gal1.com and gal2.com, provided in asurv.etc did not
match the new input formats. This caused the output files created by
the command files not to match the examples in the manual.
\item The character variable containing the name of the data file was
not defined properly in the subroutine which prints out the results of
Kaplan-Meier estimation.
\item In subroutine TWOST, the variable WRK6 was wrongly specified as WWRK6.
\item Various problems with truncation and integer arithmetic were found
and removed from the code.
\item To help VAX/VMS users, a carriage control line was added to {\sl\bf
ASURV}. For information on this addition, look in {\bf \S A7}
\end{itemize}
\bigskip
\noindent{\Large\bf A5 Errors removed in Rev 1.2}
After releasing {\sl\bf ASURV} Rev 1.1 in March 1992, we determined that
Rev. 1.1, and all previous versions, contained an inconsistency in the
way the Kaplan-Meier routine treated certain data sets. The problem occurred
when multiple upper limits were tied at the smallest data point, or when
multiple lower limits were tied at the largest data point. Since such an
event would be very unlikely in the biomedical setting, and {\sl\bf ASURV}
was initially modeled on biomedical software, no contingency for such a
situation was contained in the package. However, this type of censoring
occurs commonly in astronomical applications, so the package has been
modified to reflect this. The Kaplan-Meier routine in Rev 1.2 temporarily
changes all such tied censored points to detections so they can be
treated consistently.
{\it Other reported bugs:}
Two bugs are present in the univariate
two-sample tests. First, the generalized Wilcoxon and Peto-Prentice
univariate two-sample tests may give different results of one switches
the sample identifiers when ties are present. Second, in certain unlikely
circumstances when the user runs multiple tests without exiting ASURV, slight
changes in test results may be seen. When compiling in Sun Solaris,
three warnings due to unreferenced statement labels are reported (subroutines
GAMMA, PCHISQ and SORT2). These were reported and corrected in February 1997.
A crash is reported in Schmitt's linear regression method with many bin
subdivisions when using Sun's Fortran compiler under Solaris 2.6, which does not
show up using GNU Fortran on a PC. Crashes are also seen in the bootstrap
error analysis of Schmitt's regression for many bins.
A compiler error is encountered with g77 Fortran under the linux
operating system. The problem arises because an ASURV subroutine is
named STAT, which is also an internal function name. The problem is fixed
by renaming STAT to ASSTAT in the source code (subroutines GEHAN, PWLCXN
and STAT). This problem was reported and fixed in October 2000.
\bigskip
\noindent{\Large\bf A6 Obtaining and Installing ASURV}
This program is available as a stand-alone package to any member of
the astronomical community without charge. We provide the FORTRAN source code,
not the executable files. We prefer to distribute it by electronic mail.
Scientists with network connectivity should send their request to:
\begin{center}
\begin{verbatim}
code@stat.psu.edu
\end{verbatim}
\end{center}
\noindent specifying {\sl\bf `ASURV Rev. 1.1'} and providing both email and
regular mail addresses. ASCII versions of the code can also be mailed, if
necessary, on a 3 1/2-inch double-sided high-density MS-DOS diskette, or a
1/2 inch 9-track tape (1600 BPI, written on a UNIX machine). For any
questions regarding the package, contact:
\begin{center}
\begin{verbatim}
Prof. Eric Feigelson
Department of Astronomy & Astrophysics
Pennsylvania State University
University Park PA 16802
Telephone: (814) 865-0162
Telex: 842-510
Email: edf@astro.psu.edu
\end{verbatim}
\end{center}
The package consists of 59 subroutines of Fortran totally about 1/2 MBy. It
is completely self-contained, requiring no external libraries or programs. The
code is distributed in ten files: a brief READ.ME file; six files
called asurv1.f-asurv6.f containing the source code; two documentation files,
with the users manual in ASCII and LaTeX; and one file called asurv.etc
containing test datasets, test outputs and a subroutine list. We request
that recipients not
distribute the package themselves beyond their own institution. Rev. 0.0
of {\sl\bf ASURV} has been incorporated into the larger astronomical software
system SDAS/IRAF, distributed by the Space Telescope Science Institute.
Installation requires removing the email headers, Fortran
compilation, linking, and executing. We have written the code
consistent with Fortran 77 conventions. It has been successfully ported
to a Sun SPARCstation under UNIX, a DEC VAX under VMS, a personal computer
under MS-DOS using Microsoft FORTRAN, and an IBM mainframe under VM/CMS.
It can also be ported to a Macintosh with small format changes (see \S
{\bf A7} below).
\footnote{SPARCstation is a trademark of Sun Microsystems, Inc; UNIX is a
trademark of AT\& T Corporation; VAX and VMS are trademarks of Digital
Equipment Corporation; MS-DOS and Microsoft FORTRAN are trademarks of
Microsoft Corporation; VM/CMS is a trademark of International
Business Machines Corporation; and Macintosh is a trademark of Apple
Corporation.}
When {\sl\bf ASURV} is compiled using Microsoft FORTRAN the
user will notice several warnings from the compiler about labels used across
blocks and formal arguments not being used. The user should not be alarmed,
these warnings are caused by the compilation of the subroutines separately
and do not reflect program errors.
All of the functions have been tested against textbook formulae,
published examples, and/or commercial software packages. However, some
methods are not widely used by researchers in other fields, and their behavior
is not well-documented. We would appreciate hearing about any difficulties or
unusual behavior encountered when running the code. A bug report form for
{\sl\bf ASURV} can be found in the asurv.etc file.
\bigskip
\bigskip
\noindent{\Large\bf A7 User Adjustable Parameters}
{\sl\bf ASURV} is initially set to handle data sets of up to
500 points, with up to four variables, however a user may wish to
consider even larger data sets. With this in mind, we have modified
Rev 1.1 to be easy for a user to adjust to a given sample size.
The sizes of all the arrays in {\sl\bf ASURV} are controlled by
two parameter statements; one is in {\sl\bf UNIVAR} and one is in
{\sl\bf BIVAR}. Both statements are currently of the form:
\begin{verbatim}
PARAMETER(MVAR=4,NDAT=500,IBIN=50)
\end{verbatim}
MVAR is the number of variables allowed in a data set
and NDAT is the number of observations allowed in a data set. To work
with larger data sets, it is only necessary to change the values MVAR and/or
NDAT in the two parameter statements. If MVAR is set to be greater
than ten, then the data input formats should also be changed to read more
than ten variables from the data file. Clearly, adjusting either MVAR or NDAT
upwards will increase the memory space required to run {\sl\bf ASURV}.
The following table provides listings of format statements that a user
might wish to change if data values do not match the default formats:
\begin{verbatim}
Input Formats:
Problem Subroutine Default Format
------- ---------- --------------
Kaplan-Meier DATAIN Statement 30 FORMAT(10(I4,F10.3))
Two Sample tests
DATAIN 40 FORMAT(I4,10(I4,F10.3))
Correlation/Regression
DATREG 20 FORMAT(I4,11F10.3)
Output Formats:
Problem Subroutine Default Format
------- ---------- --------------
Kaplan-Meier KMPRNT Statements 550 [For Uncensored Points]
555 [For Censored Points]
850 [For Percentiles]
1000 [For Mean]
KMDIF 680 [For Differential KM]
Two Sample Tests TWOST Statements 612, 660, 665,780
Problem Subroutine Default Format
------- ---------- --------------
Correlation/Regression
Cox Regression COXREG Statements 1600, 1651, 1650
Kendall's Tau BHK 2005, 2007
Spearman's Rho SPRMAN 230, 240, 250, 997
EM Algorithm EM 1150, 1200, 1250, 1300, 1350
Buckley-James BJ 1200, 1250, 1300, 1350
Schmitt TWOKM 1710, 1715, 1720, 1790, 1795, 1800,
1805, 1820
\end{verbatim}
Macintosh users who have Microsoft Fortran may have difficulty with the
statements that read from the keyboard and write to the screen. Throughout
{\sl\bf ASURV} we have used statements of the form,
\begin{verbatim}
WRITE(6,format) READ(5,format).
\end{verbatim}
Macintosh users may wish to replace the 6 and 5 in statements of the above
type by `*'.
VAX/VMS users may notice that the first character of some statements
sent to the screen is not printed. This can be corrected by removing the 'C'
at the start of the following statement, in the {\sl\bf ASURV} subroutine,
\begin{verbatim}
C OPEN(6,CARRIAGECONTROL='LIST',STATUS='OLD')
\end{verbatim}
Finally, if the data have values which extend beyond three decimal places,
then the user should reduce the value of `CONST' in the subroutine
{\sl\bf CENS} to have at least two more decimal places than the data.
\newpage
\noindent{\Large\bf A8 List of subroutines}
The following is an alphabetic listing of all the subroutines
used in {\sl\bf ASURV} and their respective byte lengths.
\bigskip
\bigskip
\centerline{\large\bf Subroutine List}
\bigskip
\begin{center}
\begin{tabular}{|ll|ll|ll|} \hline
Subroutine & Bytes & Subroutine & Bytes & Subroutine & Bytes \\ \hline
&&&&& \\
AARRAY &4863 & FACTOR &1207 & REARRN & 3057 \\
AGAUSS &1345 & GAMMA &1520 & REGRES & 5505 \\
AKRANK &5887 & GEHAN &3477 & RMILLS & 1448 \\
ARISK &2675 & GRDPRB &11995 & SCHMIT & 8258\\
ASURV &7448 & KMADJ &3758 & SORT1 & 2512\\
BHK &6936 & KMDIF &5058 & SORT2 & 2214\\
BIN &15950 & KMESTM &7570 & SPEARHO & 2021\\
BIVAR &28072 & KMPRNT &8363 & SPRMAN & 8743\\
BJ &5750 & LRANK &6237 & STAT & 1859\\
BUCKLY &9830 & MATINV &3495 & SUMRY & 2324\\
CENS &2225 & MULVAR &15444 & SYMINV & 3007\\
CHOL &2672 & PCHISQ &2257 & TIE & 2784\\
COEFF &1948 & PETO &8054 & TWOKM & 16447\\
COXREG &8032 & PLESTM &6122 & TWOST & 15877\\
DATA1 &2650 & PWLCXN &2158 & UNIVAR & 27677\\
DATA2 &3470 & R3 &1629 & UNPACK & 1702\\
DATAIN &2962 & R4 &6142 & WLCXN & 6358\\
DATREG &4102 & R5 &3058 & XDATA & 1826\\
EM &12581 & R6 &9654 & XVAR & 5319\\
EMALGO &13183 & RAN1 &1432 & &\\ \hline
\end{tabular}
\end{center}
\end{document}