EXAMPLES OF NEURAL NETWORK SURVIVAL MODELS The neural network software can also be used to model survival data. In this application, the network models the "hazard function", which determines the instantaneous risk of "death". The hazard may depend on the values of inputs (covariates) for a particular case. It may be constant through time (for each case), or it may change over time. The data may be censored - ie, for some cases, it may be known only that the individual survived for at least a certain time. The details of these models are described in model-spec.doc. Command files for the examples below are found in the ex-surv directory. The sub-directory "pbc" contains files for the examples used in my talk on "Survival analysis using a Bayesian neural network" at the 2001 Joint Statistical Meetings in Atlanta. Survival models with constant hazard. I generated 700 cases of synthetic data in which the hazard function is constant through time, but depends on one input, which was randomly generated from a standard normal distribution. When the input was x, the hazard was h(x) = 1+(x-1)^2, which makes the distribution of the time of death be exponential with mean 1/h(x). For the first 200 cases, the time of death was censored at 0.5; if the actual time of death was later than this, the time recorded was minus the censoring time (ie, -0.5). The first 500 of these cases were used for training, and the remaining 200 for testing. The network used to model this data had one input (x), and one output, which was interpreted as log h(x). One layer of eight hidden units was used. The following commands specify the network, model, and data source: > net-spec clog.net 1 8 1 / ih=0.05:1 bh=0.05:1 ho=x0.05:1 bo=100 > model-spec clog.net survival const-hazard > data-spec clog.net 1 1 / cdata@1:500 . cdata@-1:500 . We can now train the network using the same general methods as for other models. Here is one set of commands that work reasonably well: > net-gen clog.net fix 0.5 > mc-spec clog.net repeat 10 heatbath hybrid 100:10 0.2 > net-mc clog.net 1 > mc-spec clog.net repeat 4 sample-sigmas heatbath hybrid 500:10 0.4 > net-mc clog.net 100 This takes 8.1 seconds on the system used (see Ex-test-system.doc). Predictions for test cases can now be made using net-pred. For example, the following command produces the 10% and 90% quantiles of the predictive distribution for time of death, based on iterations from 25 onward. These predictions can be compared to the actual times of death (the "targets"): > net-pred itq clog.net 25: Number of iterations used: 76 Case Inputs Targets 10% Qnt 90% Qnt 1 -0.01 0.29 0.05 1.08 2 -1.43 0.04 0.02 0.37 3 -0.68 0.03 0.03 0.57 4 0.47 0.11 0.07 1.59 5 0.63 0.97 0.08 1.73 ( middle lines omitted ) 196 1.39 0.21 0.08 1.83 197 0.04 0.07 0.05 1.13 198 1.66 0.56 0.08 1.67 199 1.33 1.48 0.08 1.80 200 -1.41 0.08 0.02 0.39 Since there is only one input for this model, we can easily look at the posterior distribution of the log of the hazard function with the following command: > net-eval clog.net 25:%5 / -3 3 100 | plot This plots sixteen functions that are drawn from the posterior distribution of the log of the hazard function. These functions mostly match the actual function, which is log(1+(x-1)^2), for values of x up to about one. There is considerable uncertainty for values of x above one, however. Some of the functions from the posterior are a good match for the true function in this region, but most are below the true function. Note that the data here are fairly sparse. Survival models with piecewise-constant hazard. One can also define models in which the hazard function depends on time, as well, perhaps, on various inputs. To demonstrate this, I generated synthetic data for which the hazard function was given by h(x,t) = (1+(x-1)^2) * t, which makes the cumulative distribution function for the time of death of an individual for whom the covariate is x be 1-exp(-(1+(x-1)^2)*t^2/2). I generated 1000 cases (none of which were censored), and used the first 700 for training and the remaining 300 for testing. The network for this problem has two input units. The first input is set to the time, t, the second to the value of the covariate, x. There is a single output unit, whose value will be interpreted as the log of the hazard function. If we decide to use one layer of eight hidden units, we might specify the network as follows: > net-spec vlog.net 2 8 1 / ih=0.05:1:1 bh=0.05:1 ho=x0.05:1 bo=100 The prior for the input-hidden weights uses separate hyperparameters for the two input units (representing time, t, and the covariate, x), allowing the smoothness of the hazard function to be different for the two. In particular, if the hazard is actually constant, the time input should end up being almost ignored. (This can be demonstrated by applying this model to the data used in the previous example.) The hazard function determines the likelihood by means of its integral from time zero up to the observed time of death, or the censoring time. Computing this exactly would be infeasible if the hazard function was defined to be the network output when the first input is set to a given time. To allow these integrals to be done in a reasonable amount of compute time, the hazard function is instead piecewise constant with respect to time (but continuous with respect to the covariates). The time points where the hazard function changes are given in the model-spec command. For example, we can use the following specification: > model-spec vlog.net survival pw-const-hazard 0.05 0.1 0.2 0.35 \ 0.5 0.7 1.0 1.5 The eight numbers above are the times at which the hazard function changes. The log of the value of the hazard before the first time (0.05) is found from the output of the network when the first input is set to this time (and other inputs are set to the covariates for the case); similarly for the log of the value of the hazard function after the last time point (1.5). The log of the hazard for other pieces is found by setting the first input of the network to the average of the times at the ends of the piece. For instance, the log of the value of the hazard function in the time range 0.7 to 1.0 is found by setting the first input to (0.7+1.0)/2 = 0.85. Alternatively, if the "log" option is used, everything is done in terms of the log of the time (see model-spec.doc for details). Note that the number of time points specified is limited only by the greater amount of computation required when there are many of them. The hyperparameters controlling the weight priors will be able to control "overfitting" even if the hazard function has many pieces. The data specification for this problem says there is just one input (the covariate) and one target (the time of death, or censoring time): > data-spec vlog.net 1 1 / vdata@1:700 . vdata@-1:700 . Note that data-spec should be told there is only one input even though net-spec is told there are two inputs - the extra input for the net is the time, not an input read from the data file. The following commands sample from the posterior distribution: > net-gen vlog.net fix 0.5 > mc-spec vlog.net repeat 10 heatbath hybrid 100:10 0.2 > net-mc vlog.net 1 > mc-spec vlog.net repeat 4 sample-sigmas heatbath hybrid 500:10 0.4 > net-mc vlog.net 100 This takes 84 seconds on the system used (see Ex-test-system.doc). The net-pred command can be used to make predictions for test cases. Here is a command to display the median of the predictive distribution for time of death, based on iterations from 25 on, along with the covariates (inputs) and the actual times of death (targets) for the test cases: > net-pred itd vlog.net 25: Number of iterations used: 76 Case Inputs Targets Medians |Error| 1 0.00 0.07 0.84 0.7661 2 -0.22 0.42 0.76 0.3317 3 -2.24 0.04 0.39 0.3485 4 -0.81 0.26 0.60 0.3335 5 1.02 3.12 1.05 2.0715 ( middle lines omitted ) 296 -0.81 0.28 0.59 0.3139 297 1.19 1.60 1.05 0.5540 298 -1.53 0.57 0.45 0.1161 299 0.99 2.67 1.06 1.6090 300 0.28 1.13 0.93 0.1987 Average abs. error guessing median: 0.32963+-0.01648 The predictive medians can be plotted against the inputs as follows: > net-pred idb vlog.net 25: / vdata@-1:700 | plot-points (Note that by specifying only the test inputs, without the targets, and using the 'b' option, the output piped to plot-points has just the input and median.) The fairly small amount of scatter seen in this plot is a result of computing the medians by Monte Carlo, which produces some random error. The exact median time of death as predicted by the model is a smooth function of the input. These predictions can be compared to the true median time of death for individuals with the given covariates, as determined by the way the data was generated, which are in the file "vmedians". Here is a command to compare them on the test data: > (tail -300 vmedians; echo " " > net-pred idb vlog.net 25: / vdata@-1:700) | plot-points The result is in vlog-medians.png. The log of the true hazard function for this problem is the sum of a function of x only and a function of t only - what is called a "proportional hazards" model. A network model specified as follows can discover this: > net-spec v2log.net 2 8 8 1 / ih0=0.05:1:1 bh0=0.05:1 \ ih1=0.05:1:1 bh1=0.05:1 \ ho1=x0.05:1 ho0=x0.05:1 bo=100 The input-hidden weights for the two hidden layers have hierarchical priors that allow one of them to (mostly) look at only the time input and the other to look at only the covariate. There are no connections between the hidden layers. Final output can therefore be an additive function of time and the covariate, if the hyperparameters are set in this way. The remaining commands can be the same as above, except that the stepsize needs to be reduced to get a reasonable acceptance rate (the number of leapfrog steps is increased to compensate): > model-spec v2log.net survival pw-const-hazard 0.05 0.1 0.2 0.35 \ 0.5 0.7 1.0 1.5 > data-spec v2log.net 1 1 / vdata@1:700 . vdata@-1:700 . > net-gen v2log.net fix 0.5 > mc-spec v2log.net repeat 10 heatbath hybrid 100:10 0.1 > net-mc v2log.net 1 > mc-spec v2log.net repeat 4 sample-sigmas heatbath hybrid 1000:10 0.25 > net-mc v2log.net 400 Running this takes 19 minutes on the system used (Ex-test-system.doc). The model does in fact discover that an additive form for the log hazard is appropriate, as can be seen by examining the hyperparameters with commands such as the following: > net-plt t h1@ v2log.net | plot > net-plt t h3@ v2log.net | plot One can see from this that the Markov chain spends most of its time in regions of the hyperparameter space in which one of the hidden layers looks more at the first input (time) than the second input, and the other hidden layer looks more at the second input than the first.