DATA MINING Desktop Survival Guide by Graham Williams
Desktop Survival Project Home Introduction Getting Started The Business Problem Data Loading Data Exploring Data Interactive Graphics Statistical Tests Models Network Analysis Text Mining Decision Trees Random Forests Boosting Bagging Support Vector Machine Linear Regression Neural Network Naive Bayes Survival Analysis Evaluation and Deployment Transforming Data Deployment Troubleshooting Issues Moving into R R Getting Help Data Graphics in R Understanding Data Preparing Data Issues Evaluating Models Reporting Fraud Analysis Archetype Analysis Algorithms Bayes Classifier K-Nearest Neighbours Linear Models Open Products AlphaMiner Borgelt Data Mining Suite KNime R Rattle Weka Closed Products Clementine Equbits Foresight GhostMiner InductionEngine ODM Enterprise Miner Statistica Data Miner TreeNet Virtual Predict Installing Rattle Projects Bibliography Index Preface Goals Organisation Features Audience Typographical Conventions A Note on Languages Currency Acknowledgements Introduction Data Mining The Business Problem Types of Analysis Data Mining Applications A Framework for Modelling Agile Data Mining R Rattle Why R and Rattle? Data Preparation Number of Algorithms Repeatability Performance Open Source Data Mining Business Case Sample Business Case Pros and Cons Books on R Getting Started Initial Interaction with R Quitting Rattle and R First Contact Loading a Dataset Building a Model Understanding Our Data Evaluating the Model Evaluating the Model Interacting with R Interacting with Rattle Projects Toolbar Menus Interacting with Plots Keyboard Navigation Summary Command Summary The Business Problem Solar Panel Efficiency Water Collection Others Other Business Problems Fraud Detection Loan Approval Documenting the Business Problem Summary Resources Exercises Data Data Nomenclature Loading Data into Rattle CSV Data Datasets Reading Direct from URL Play Golf Weather Data Other Data Sources ARFF Data ODBC Sourced Data Setting Up a Data Source Name Netezza Setup Teradata Setup R Data R Dataset Data Entry Data Tab Options in Rattle Sampling Data Variable Roles Automatic Role Identification Weights Calculator Manipulating Data Loading Data CSV Data Locating and Loading Data Loading the File CSV Options Basic Data Summary ARFF Data ODBC Sourced Data R Dataset R Data Library Data Options Sampling Data Variable Roles Automatic Role Identification Weights Calculator Command Summary Exploring Data Summarising Data Summary Describe Basics Kurtosis Skewness Missing Exploring Distributions Box Plot Histogram Cumulative Distribution Plot Benford's Law Other Digits Stratified Benford Plots Bar Plot Dot Plot Mosaic Plot GGobi Scatterplot Data Viewer Brushing Identify Multivariate Outliers Other Options Quality Plots Using R Further GGobi Documentation Correlation Analysis Hierarchical Correlation Principal Components Single Variable Overviews Interactive Graphics Interactive Visualisations Latticist GGobi Scatterplot Multiple Plots Brushing Other Plots Data Viewer Brushing Identify Multivariate Outliers Other Options Quality Plots Using R Further GGobi Documentation Documenting Interactive Explorations Code Review Chapter Exercises Command Summary Statistical Tests Documenting Interactive Explorations Code Review Further Resources Chapter Exercises Command Summary Models A Framework for Modelling Descriptive Analytics Predictive Analytics Documenting Models Summary Code Review Exercises Further Resources Command Summary Cluster Analysis Summary Clusters Basic Clustering Hot Spots Alternative Clustering Other Cluster Examples KMeans Export KMeans Clusters Discriminant Coordinates Plot Number of Clusters Hierarchical Clusters Other Cluster Algorithms Association Analysis Summary Overview Algorithm Usage Read Transactions file format sep cols rm.duplicates Summary Apriori data parameter appearance control Inspect Examples Video Marketing Survey Data Other Examples Resources and Further Reading Basket Analysis General Rules Network Analysis Documenting Interactive Explorations Code Review Chapter Exercises Command Summary Text Mining Application to Text Text Mining with R Decision Trees Knowledge Representation Search Heuristic Measures Tutorial Example Rattle R Tuning Parameters Min Split (Rarg[]minsplit) Min Bucket (minbucket) Priors (prior) Loss Matrix Complexity (cp) Other Options Simple Example Convert Tree to Rules Predicting Salary Group Issues Summary Code Review Iris Wine Exercises Resources Command Summary Random Forests Formalities Tutorial Example Tuning Parameters Number of Trees Sample Size Number of Variables Summary Overview Algorithm Usage Random Forest importance classwt Examples Resources and Further Reading Summary Overview Example Algorithm Resources and Further Reading Boosting Formalities Tutorial Example Tuning Parameters Summary Overview AdaBoost Algorithm Examples Step by Step Using gbm Extensions and Variations Alternating Decision Tree Resources and Further Reading Documenting Code Review Further Resources Chapter Exercises Command Summary Bootstrapping Summary Usage Further Information Summary Overview Example Algorithm Resources and Further Reading Bagging Support Vector Machine Formalities Tutorial Example Tuning parameters Examples Resources and Further Reading Overview Examples Resources and Further Reading Linear Regression Linear Regression Formalities Tutorial Example Tuning parameters Generalized Regression Formalities Tutorial Example Tuning parameters Logistic Regression Formalities Tutorial Example Tuning parameters Discussion Probit Regression Formalities Tutorial Example Tuning Parameters Multinomial Regression Formalities Tutorial Example Tuning Parameters Neural Network Formalities Tutorial Example Tuning parameters Documenting Code Review Further Resources Chapter Exercises Command Summary Naive Bayes Summary Code Review Resources Exercises Command Summary Survival Analysis Sample Data Simple Lung Descriptive Analysis Regression survreg Simple Lung coxph Simple Lung Apply to New Data More Input Variables Decision Tree Example from Singer and Willett Other Approaches Design Package Random Survival Forests Prediction on Test Data Evaluation The Evaluate Tab Confusion Matrix Measures Graphical Measures Risk Charts Cost Curves Lift ROC Curves Area Under Curve Precision versus Recall Sensitivity versus Specificity Predicted versus Observed Scoring Documenting Interactive Explorations Code Review Chapter Exercises Command Summary Transforming Data Rescale Data Recenter Scale [0,1] Rank Median/MAD Peer Relativity Profiling Index Impute Zero/Missing Mean/Median/Mode Constant Remap Binning Indicator Variables Join Categorics Math Transforms Outliers Cleanup Delete Ignored Delete Selected Delete Missing Delete Obs with Missing Other Transformations Removing Duplicates Command Summary Deployment Documenting Deployment Code Review Chapter Exercises Command Summary Troubleshooting Cairo A factor has new levels Issues Model Selection Overfitting Imbalanced Classification Sampling Cost Based Learning Model Deployment and Interoperability SQL PMML XML for Data Bibliographic Notes Documenting Code Review Chapter Exercises Command Summary Moving into R Interacting with R Basic Command Line Windows, Icons, Mouse, Pointer--WIMP The Current Rattle State Samples Projects The Rattle Log Further Tuning Models Emacs and ESS Documenting Code Review Chapter Exercises Command Summary R Evaluation Exercises Assignment Libraries and Packages Searching for Objects Package Management Information About a Package Testing Package Availability Packages and Namespaces Basic Programming in R Principles Folders and Files Flow Control If Statement For Loop Functions Apply Methods Objects System Running System Commands System Parameters Misc Internet Memory Management Memory Usage Garbage Collection Errors Frivolous Sudoku Further Resources Using R Specific Purposes Survey Analysis Getting Help R Documentation Data Data Types Numbers Strings Building Strings Splitting Strings Substitution Trim Whitespace Evaluating Strings Logical Dates and Times Space Data Structures Vectors Arrays Lists Sets Matricies Exercises Data Frames Accessing Columns Removing Columns Exercises General Manipulation Factors Elements Rows and Columns Finding Index of Elements Partitions Head and Tail Reverse a List Sorting Unique Values Loading Data Interactive Responses Interactive Data Entry Available Datasets The Iris Dataset CSV Data Used In The Book The Wine Dataset The Cardiac Arrhythmia Dataset The Adult Survey Dataset Foreign Formats Stata Data Conversions Reading Variable Width Data Saving Data Formatted Output Automatically Generate Filenames Reading a Large File Manipulating Data Manipulating Data As SQL Using SQLite ODBC Data Database Connection Excel Access Clipboard Data Spatial Data Simple Map A Density Map Overlays and Point in Polygon Other Data Formats Fixed Width Data Global Positioning System Documenting a Dataset Common Data Problems Graphics in R Basic Plot Controlling Axes Arrow Axes Legends and Points Tables Within Plots Colour Labels in Plots Axis Labels Legend Labels Within Plots Maths in Labels Multiple Plots MatPlot Multiple Plots Using ggplot2 Using GGPlot Networks Symbols Other Graphic Elements Making an Animation Animated Mandelbrot Adding a Logo to a Graphic Graphics Devices Setup Screen Devices Multiple Devices File Devices Multiple Plots Copy and Print Devices Graphics Parameters Plotting Region Locating Points on a Plot Scientific Notation and Plots Understanding Data Single Variable Overviews Textual Summaries Multiple Line Plots Separate Line Plots Pie Chart Fan Plot Stem and Leaf Plots Histogram Barplot Trellis Histogram Histogram Uneven Distribution Bump Chart Density Plot Basic Histogram Basic Histogram with Density Curve Practical Histogram Multiple Variable Overviews Scatterplot Scatterplot with Marginal Histograms Multi-Dimension Scatterplot Correlation Plot Colourful Correlations Fluctuation Plot Heat Map Projection Pursuit RADVIZ Parallel Coordinates Categoric and Numeric Measuring Data Distributions Textual Summaries Boxplot Multiple Boxplots Boxplot by Class Tuning a Boxplot Boxplot Using Lattice Boxplot Using ggplot Violin Plot What Distribution Miscellaneous Plots Line and Point Plots Matrix Data Multiple Plots Aligned Plots Probability Scale Network Plot Sunflower Plot Stairs Plot Graphing Means and Error Bars Bar Charts With Segments Bar Plot With Means 3d Bar Plot Stacks Versus Lines Multi-Line Title Mathematics Plots for Normality Basic Bar Chart Bar Chart Displays Multiple Dot Plots Alternative Multiple Dot Plots 3D Plot Clustered Box Plot Perspective Plots Star Plot Residuals Plot Waterfall Plots Dates and Times Simple Time Series Multiple Time Series Plot Time Series Plot Time Series with Axis Labels Grouping Time Series for Box Plot Time Series Heatmap Textual Summaries Stem and Leaf Plots Histogram Barplot Density Plot Basic Histogram Basic Histogram with Density Curve Practical Histogram Correlation Plot Colourful Correlations Measuring Data Distributions Textual Summaries Boxplot Multiple Boxplots Boxplot by Class Box and Whisker Plot Box and Whisker Plot Clustered Box Plot Further Resources Map Displays Further Resources Preparing Data Data Selection and Extraction Training and Test Datasets Data Cleaning Review Data Selectively Changing Vector Values Replace Indices By Names Missing Values Remove Levels from a Factor Variable Manipulations Remove Columns Reorder Columns Remove Non-Numeric Columns Remove Variables with no Variance Cleaning the Wine Dataset Cleaning the Cardiac Dataset Cleaning the Survey Dataset Imputation Nearest Neighbours Multiple Imputation Data Linking Simple Linking Record Linkage Data Transformation Aggregation Sum of Columns Pivot Tables Normalising Data Binning Interpolation Variable Selection Classification Classification Classification Issues Incremental or Online Modelling Model Tuning Tuning rpart Unbalanced Classification Building Models Temporal Analysis Evaluation Basics Basic Measures Cross Validation Graphical Performance Measures Lift The ROC Curve Other Examples 10 Fold Cross Validation Area Under Curve Calibration Curves Reporting Generating Open Document Format Getting Started with odfWeave OpenOffice.org Macro Support Generating HTML Generating PDF with LATEX Configuration Figure Sizes Fraud Analysis Archetype Analysis

Example from Singer and Willett

From http://www.ats.ucla.edu/stat/R/examples/alda/ch14.htm and Copyrighted (c) by UCLA Academic Technology Services.

Obtain the data:

> library(survival)
> rearrest<-read.table("http://www.ats.ucla.edu/stat/examples/alda/rearrest.csv", 
                       sep=",", header=T)

Separate data into two datasets, grouped on personal. Next, run an empty Cox model and use the survfit function to generate the event times and estimated survival probabilities at the event times for both datasets. These probabilities are plotted, using a step function. The plot represents the estimated survivor function for the recidivism data stratified by the predictor personal.

> rearrest0 <- subset(rearrest, personal == 0)
> rearrest1 <- subset(rearrest, personal == 1)
> f14.1.0 <- summary(survfit(Surv(rearrest0$months, abs(rearrest0$censor - 1)) ~ 1))
> f14.1.1 <- summary(survfit(Surv(rearrest1$months, abs(rearrest1$censor - 1)) ~ 1))
> s.hat.0 <- f14.1.0[1][[1]]
> time.0 <- f14.1.0[2][[1]]
> s.hat.1 <- f14.1.1[1][[1]]
> time.1 <- f14.1.1[2][[1]]
> s.hat.steps.0 <- stepfun(time.0, c(1, s.hat.0))
> s.hat.steps.1 <- stepfun(time.1, c(1, s.hat.1))
> plot(s.hat.steps.0, do.points = FALSE, xlim = c(0, 36), ylim = c(0,1),
    ylab = "Estimated Survival", xlab = "Months after release", main = "")
> lines(s.hat.steps.1, do.points = FALSE, xlim = c(0,36), lty = 2)
> legend("bottomleft", c("Personal = 0", "Personal = 1"), lty = c(1, 2))

Now graph the corresponding cumulative hazard stratified by personal, calculated by taking the negative log of the survival probabilities.

> h.hat.0 <- -log(s.hat.0)
> h.hat.1 <- -log(s.hat.1)
> h.hat.steps.0 <- stepfun(time.0, c(0, h.hat.0))
> h.hat.steps.1 <- stepfun(time.1, c(0, h.hat.1))
> plot(h.hat.steps.0, do.points = FALSE, xlim = c(0, 36), ylim = c(0,1.4),
    ylab = "Estimated Cumulative Hazard", xlab = "Months after release", main = "")
> lines(h.hat.steps.1, do.points = FALSE, xlim = c(0,36), lty = 2)
> legend("bottomright", c("Personal = 0", "Personal = 1"), lty = c(1, 2))

Introduce a function, smooth, to take the time and estimated survival probabilities from survfit and aggregates the hazard over a given width of time. This function uses the same smoothing algorithm the book authors used. Using one of R's kernel smoothing functions (like ksmooth) would not exactly reproduce the graphs as they appear in the text. In this instance, the bandwidth of the graph in the text appears to be 8.

> smooth<- function(width, time, survive){
    n=length(time)
    lo=time[1] + width
    hi=time[n] - width
    npt=50
    inc=(hi-lo)/npt
    s=lo+t(c(1:npt))*inc
    slag = c(1, survive[1:n-1])
    h=1-survive/slag
    x1 = as.vector(rep(1, npt))%*%(t(time))
    x2 = t(s)%*%as.vector(rep(1,n)) 
    x = (x1 - x2) / width
    k=.75*(1-x*x)*(abs(x)<=1)
    lambda=(k%*%h)/width
    smoothed= list(x = s, y = lambda)
    return(smoothed)
  }
> bw1.0 <- smooth(8, time.0, s.hat.0)
> bw1.1 <- smooth(8, time.1, s.hat.1)
> plot(bw1.0$x, bw1.0$y, type = "l", xlim = c(0, 36), ylim = c(0, .08),
    xlab = "Months after release", ylab = "Smoothed hazard")
> points(bw1.1$x, bw1.1$y, type = "l", lty = 2)
> legend("bottomright", c("Personal = 0", "Personal = 1"), lty = c(1, 2))

Figure 14.2, page 508. To generate the first in this set of graphs, we created a step function from the log of the cumulative hazard calculated in the previous example.

> l.hat.0 <- log(h.hat.0)
> l.hat.1 <- log(h.hat.1)
> log.h.steps.0<- stepfun(time.0, c(l.hat.0[1], l.hat.0))
> log.h.steps.1<- stepfun(time.1, c(l.hat.1[1], l.hat.1))
> par(mfrow=c(1,1))
> plot(log.h.steps.0, do.points = FALSE, xlim = c(0, 36), ylim = c(-6.0,1),
    ylab = "Log(Cumulative Hazard)", xlab = "Months after release", main = "")
> lines(log.h.steps.1, do.points = FALSE, xlim = c(0,36), lty = 2)
> points(c(-5, 36), c(0,0), type = "l")
> legend("bottomright", c("Personal = 0", "Personal = 1"), lty = c(1, 2))

To generate the second and third graphs in this set, we first ran a Cox model with personal as a predictor. We then created one data frame where personal is 0 and one data frame where personal is 1 and used the model to obtain estimates of the hazard rate for each value of personal (using the proportional hazards assumption). We then plotted the fitted curves from the model and overlaid the non-parametric estimates for the log cumulative hazard at the times of events on the second graph. The third graph is based on the cumulative hazard rather than the log cumulative hazard.

> attach(rearrest)
> f14.2 <- coxph(Surv(months, abs(censor - 1))~personal, method = "efron")
> personal.0 <- data.frame(personal=0)
> personal.1 <- data.frame(personal=1)
> s.base <- survfit(f14.2, newdata = personal.0, type = "kaplan-meier")
> s.personal <- survfit(f14.2, newdata = personal.1, type = "kaplan-meier")
> h.base <- -log(s.base$surv)
> h.personal <- -log(s.personal$surv)
> l.h.base <- log(h.base)
> l.h.personal <- log(h.personal)
> plot(s.base$time, l.h.base, type = "l", lty = 1, ylim = c(-6, 1),
    xlab = "Months after release", ylab = "Cumulative Hazard")
> points(s.personal$time, l.h.personal, type = "l", lty = 2)
> points(time.0, l.hat.0, pch = 3)
> points(time.1, l.hat.1, pch = 4)
> legend("bottomright", c("Personal = 0", "Personal = 1"), lty = c(1, 2))

> plot(s.base$time, h.base, type = "l", lty = 1, ylim = c(0,1.5),
    xlab = "Months after release", ylab = "Cumulative Hazard")
> points(s.personal$time, h.personal, type = "l", lty = 2)
> points(time.0, h.hat.0, pch = 3)
> points(time.1, h.hat.1, pch = 4)
> legend("bottomright", c("Personal = 0", "Personal = 1"), lty = c(1, 2))

Table 14.1, page 525. To generate the parameter estimates and model statistics seen in this table, we ran Cox models with the predictors specified in the table.

> tab14.1A <- coxph(Surv(months, abs(censor - 1))~personal)
> summary(tab14.1A)

Call:
coxph(formula = Surv(months, abs(censor - 1)) ~ personal)

  n= 194 

           coef exp(coef) se(coef)     z Pr(>|z|)  
personal 0.4787    1.6140   0.2025 2.364   0.0181 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

         exp(coef) exp(-coef) lower .95 upper .95
personal     1.614     0.6196     1.085       2.4

Rsquare= 0.027   (max possible= 0.994 )
Likelihood ratio test= 5.32  on 1 df,   p=0.02104
Wald test            = 5.59  on 1 df,   p=0.01807
Score (logrank) test = 5.69  on 1 df,   p=0.01702

> tab14.1B <- coxph(Surv(months, abs(censor - 1))~property)
> summary(tab14.1B)

Call:
coxph(formula = Surv(months, abs(censor - 1)) ~ property)

  n= 194 

           coef exp(coef) se(coef)    z Pr(>|z|)    
property 1.1946    3.3022   0.3493 3.42 0.000626 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

         exp(coef) exp(-coef) lower .95 upper .95
property     3.302     0.3028     1.665     6.548

Rsquare= 0.08   (max possible= 0.994 )
Likelihood ratio test= 16.2  on 1 df,   p=0.00005684
Wald test            = 11.7  on 1 df,   p=0.000626
Score (logrank) test = 13.14  on 1 df,   p=0.0002896

> tab14.1C <- coxph(Surv(months, abs(censor - 1))~cage)
> summary(tab14.1C)

Call:
coxph(formula = Surv(months, abs(censor - 1)) ~ cage)

  n= 194 

         coef exp(coef) se(coef)      z  Pr(>|z|)    
cage -0.06813   0.93414  0.01563 -4.359 0.0000131 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

     exp(coef) exp(-coef) lower .95 upper .95
cage    0.9341      1.071     0.906    0.9632

Rsquare= 0.112   (max possible= 0.994 )
Likelihood ratio test= 22.97  on 1 df,   p=0.000001646
Wald test            = 19  on 1 df,   p=0.00001306
Score (logrank) test = 19.18  on 1 df,   p=0.00001191

> tab14.1D <- coxph(Surv(months, abs(censor - 1))~personal + property + cage)
> summary(tab14.1D)

Call:
coxph(formula = Surv(months, abs(censor - 1)) ~ personal + property + 
    cage)

  n= 194 

             coef exp(coef) se(coef)      z Pr(>|z|)    
personal  0.56914   1.76674  0.20521  2.773  0.00555 ** 
property  0.93579   2.54922  0.35088  2.667  0.00765 ** 
cage     -0.06671   0.93546  0.01678 -3.976  0.00007 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

         exp(coef) exp(-coef) lower .95 upper .95
personal    1.7667     0.5660    1.1817    2.6415
property    2.5492     0.3923    1.2816    5.0708
cage        0.9355     1.0690    0.9052    0.9667

Rsquare= 0.182   (max possible= 0.994 )
Likelihood ratio test= 38.96  on 3 df,   p=0.00000001770
Wald test            = 29.02  on 3 df,   p=0.000002214
Score (logrank) test = 30.3  on 3 df,   p=0.000001195

Table 14.2, page 533. To generate this table, we wrote a function risk.score that takes an ID as an argument and returns a vector containing the information for the given ID (including the risk score). We then called this function for the set of IDs that appear in the table and combined the vectors of information from each function call into a table.

> risk.score <- function(ID.num){
    temp <- subset(rearrest, id==ID.num)
    log.score <- tab14.1D$coef[1]*temp$personal + tab14.1D$coef[2]*temp$property + tab14.1D$coef[3]*temp$cage
    score <- exp(log.score)
    day <- temp$months * (365/12)
    table.row <- cbind(temp$id, temp$personal, temp$property, temp$cage, score, day, temp$months, temp$censor)
    colnames(table.row) <- c("ID", "personal", "property", "centered.age", "risk.score", "day", "months", "censor")
    rownames(table.row) <- c()
    return(table.row)}
> tab.14.2 <- rbind(
    risk.score(22),
    risk.score(8),
    risk.score(187),
    risk.score(26),
    risk.score(5),
    risk.score(130),
    risk.score(106),
    risk.score(33))
> tab.14.2

      ID personal property centered.age risk.score
[1,]  22        0        0    0.2577249  0.9829532
[2,]   8        1        1   22.4507434  1.0071815
[3,] 187        1        0   -7.2001806  2.8561843
[4,]  26        0        1   -7.3014811  4.1491205
[5,]   5        1        1   -7.1645886  7.2637815
[6,] 130        0        1   22.3905107  0.5723740
[7,] 106        0        0   16.2029679  0.3392704
[8,]  33        1        0   27.0612841  0.2904806
            day     months censor
[1,]   51.96441  1.7084189      1
[2,]   18.98700  0.6242300      1
[3,] 1095.00000 36.0000000      1
[4,]   71.95072  2.3655031      0
[5,]    8.99384  0.2956879      0
[6,]  485.66735 15.9671458      1
[7,]  355.75633 11.6960986      0
[8,]   84.94182  2.7926078      1

Figure 14.4, page 538. To generate these plots, we first reran the fourth model from table 14.1 and created two data frames: one representing the baseline (all predictors equal to zero) and one representing the average (all predictors equal to their sample means). We then used the model to predict survival and cumulative hazard rates for both the baseline and average case and plotted the results in the first four graphs. For the fifth and sixth graphs, we again used the smooth function described earlier in this page to generate kernel smoothed values of the hazard function.

> tab14.1D <- coxph(Surv(months, abs(censor - 1))~personal + property + cage)
> base <- data.frame(personal=0, property=0, cage=0)
> ss <- survfit(tab14.1D)
> ss.base <- survfit(tab14.1D, newdata = base)
> avg <- data.frame(personal=mean(personal), property=mean(property), cage=0)
> ss.avg <- survfit(tab14.1D, newdata = avg)
> par(mfrow=c(1,2))
> plot(ss.base$time, ss.base$surv, type = "l", ylim = c(0,1),
    xlab = "Months after release", ylab = "Estimated Survival")
> points(ss.base$time, ss.base$surv, pch = 20)
> plot(ss.avg$time, ss.avg$surv, type = "l", ylim = c(0,1),
    xlab = "Months after release", ylab = "Estimated Survival")
> points(ss.avg$time, ss.avg$surv, pch = 20)

> plot(ss.base$time, -log(ss.base$surv), type = "l", ylim = c(0,1.5),
    xlab = "Months after release", ylab = "Estimated Cumulative Hazard")
> points(ss.base$time, -log(ss.base$surv), pch = 20)
> plot(ss.avg$time, -log(ss.avg$surv), type = "l", ylim = c(0,1.5),
    xlab = "Months after release", ylab = "Estimated Cumulative Hazard")
> points(ss.avg$time, -log(ss.avg$surv), pch = 20)

> smooth.base <- smooth(8, ss.base$time, ss.base$surv)
> smooth.avg<- smooth(8, ss.avg$time, ss.avg$surv)
> plot(smooth.base$x, smooth.base$y, type = "l", xlim = c(0, 36), ylim = c(0, .08),
    xlab = "Months after release", ylab = "Smoothed hazard")
> plot(smooth.avg$x, smooth.avg$y, type = "l", xlim = c(0, 36), ylim = c(0, .08),
    xlab = "Months after release", ylab = "Smoothed hazard")

Figure 14.5, page 541.To generate these plots, we first reran the fourth model from table 14.1 and created four data frames, all with age variable set to its mean (zero, since age is centered). One data frame had personal and property values both set to 0, one data frame had personal set to 1 and property to 0, one data frame had property set to 1 and personal set to 0, and the fourth data frame had both personal and property set to zero. We then used the model to predict survival (first graph), cumulative hazard (second graph), and log cumulative hazard (third graph) for each of the four scenarios represented in a data frame and plotted each set of estimates.

> tab14.1D <- coxph(Surv(months, abs(censor - 1))~personal + property + cage)
> base <- data.frame(personal=0, property=0, cage=0)
> personal.only <- data.frame(personal=1, property=0, cage=0)
> property.only <- data.frame(personal=0, property=1, cage=0)
> both <- data.frame(personal=1, property=1, cage=0)
> s.base <- survfit(tab14.1D, newdata = base)
> s.personal <- survfit(tab14.1D, newdata = personal.only)
> s.property <- survfit(tab14.1D, newdata = property.only)
> s.both <- survfit(tab14.1D, newdata = both)
> par(mfrow=c(1,1))
> plot(s.base$time, s.base$surv, type = "l", lty = 1, ylim = c(0,1),
    xlab = "Months after release", ylab = "Estimated Survival")
> points(s.personal$time, s.personal$surv, type = "l", lty = 2)
> points(s.property$time, s.property$surv, type = "l", lty = 3)
> points(s.both$time, s.both$surv, type = "l", lty = 4)
> legend("bottomleft", c("Personal = 0, Property = 0", "Personal = 1, Property = 0", 
  "Personal = 0, Property = 1", "Personal = 1, Property = 1"), lty = c(1, 2, 3, 4))

> h.base <- -log(s.base$surv)
> h.personal <- -log(s.personal$surv)
> h.property <- -log(s.property$surv)
> h.both <- -log(s.both$surv)
> plot(s.base$time, h.base, type = "l", lty = 1, ylim = c(0,2),
    xlab = "Months after release", ylab = "Estimated Cumulative Hazard")
> points(s.personal$time, h.personal, type = "l", lty = 2)
> points(s.property$time, h.property, type = "l", lty = 3)
> points(s.both$time, h.both, type = "l", lty = 4)
> legend("topleft", c("Personal = 0, Property = 0", "Personal = 1, Property = 0", 
  "Personal = 0, Property = 1", "Personal = 1, Property = 1"), lty = c(1, 2, 3, 4))

> l.h.base <- log(h.base)
> l.h.personal <- log(h.personal)
> l.h.property <- log(h.property)
> l.h.both <- log(h.both)
> plot(s.base$time, l.h.base, type = "l", lty = 1, ylim = c(-7,1),
    xlab = "Months after release", ylab = "Log(Cumulative Hazard)")
> points(s.personal$time, l.h.personal, type = "l", lty = 2)
> points(s.property$time, l.h.property, type = "l", lty = 3)
> points(s.both$time, l.h.both, type = "l", lty = 4)
> legend("bottomright", c("Personal = 0, Property = 0", "Personal = 1, Property = 0", 
  "Personal = 0, Property = 1", "Personal = 1, Property = 1"), lty = c(1, 2, 3, 4))

Support further development through the purchase of the PDF version of the book.
The PDF version is a formatted comprehensive draft book (with over 800 pages).
Brought to you by Togaware. This page generated: Sunday, 22 August 2010