Comparing Approaches

Here are three different approaches to producing an interactive plot of a Cox proportional hazards regression model. Interactivity can be very useful for exploring survival data because (given a model) we can examine how a curve changes when covariates are varied. In many cases, we are interested in comparing survival curves for different groups, such as treated individuals versus non-treated. We will show here how much effort it requires to create an interactive grouped survival curve figure using healthvis, shiny, and d3. The examples below use the veteran data set in the survival package in R.

Healthvis

Figure: http://healthviz.appspot.com/display/hs_20001

Code:

# Create a coxph object (survival package)
cobj # Plot the object
survivalVis(cobj, data=veteran, plot.title="Veteran Survival Data", group="trt",
group.names=c("Treatment", "No Treatment"), line.col=rainbow_hcl(2,start=50,end=270))

Using the survivalVis() function, we specify the covariate we want to group the data with (in this case, “trt” or treatment). We pass in the regression object (much as one would using a familiar plot() command) and set other options such as the plot’s title, group names, and colors to use for the lines. Once this function is run, a plot is generated on the Google App Engine server. This plot can then be perma-linked (as it is above) and/or embedded into a web page or other container. We do not embed the figure here because 1) wordpress does not allow embedding of iframes 2) the shiny example will require a link, and we want the comparison to be similar.

shiny

Figure: http://glimmer.rstudio.com/prpatil/survival_ex/

Code:

server.R

library(shiny)
library(survival)
library(colorspace)

veteran$trt veteran$prior cobj tmp1 tmp2$trt
shinyServer(function(input, output) {
	tmpdata 		tmp1$age 		tmp1$celltype 		tmp1$prior 		rbind(tmp1, tmp2)
	})

	output$survPlot 		plot(survfit(cobj, newdata=tmpdata()), xlab="Days", ylab="Survival Probability", col=rainbow_hcl(2))
		})
})

ui.R

library(shiny)
library(survival)

veteran$trt veteran$prior cobj vars
shinyUI(pageWithSidebar(

  headerPanel("Veteran Survival Data"),

 sidebarPanel(
	numericInput("age","Age",0,min=34,max=81),
	selectInput("celltype", "Cell Type", as.list(levels(veteran$celltype))),
	selectInput("prior", "Prior", list(0,10))),

  mainPanel(plotOutput("survPlot"))
))

Deploying the shiny version of this figure requires a bit more code, although there is more opportunity to customize. We should note that if we wanted to group by a different covariate, we’d have to either rewrite the panels and reactivity actions or add the ability to dynamically update the UI. It is also possible to incorporate javascript (such as the d3 code below), but this requires an even finer understanding of shiny (that we are still gaining…).

d3

Code:

    this.visualize = function() {

        // Initialize baseline hazard function
        var data = update_hazard(this.init_data, this.coef, this.covar);

        // create xAxis
        var xAxis = d3.svg.axis().scale(this.x).orient('bottom');

        // Add the x-axis.
        this.vis.append('svg:g')
                .attr('class', 'x axis')
                .attr('transform', 'translate(0,' + (this.h-30) + ')')
                .call(xAxis);

	this.vis.append('text')
		.attr('class', 'x label')
		.attr('x', this.w/2-12)
		.attr('y', this.h+10)
		.text('Time');

        // create left yAxis
        var yAxis = d3.svg.axis().scale(this.y).ticks(6).orient('left');

        // Add the y-axis to the left
        this.vis.append('svg:g')
                .attr('class', 'y axis')
                .attr('transform', 'translate(-5,0)')
                .call(yAxis);

	this.vis.append('text')
		.attr('class', 'y label')
		.attr('x', -200)
		.attr('y', -30)
		.attr('transform', 'rotate(-90)')
		.text('Survival');

	if(typeof this.colors == "string"){
	        var colors = [this.colors];
	} else {
		var colors = this.colors;
	}

        // Line drawer
        var x = this.x;
        var y = this.y;
        this.line = d3.svg.line()
                          .x(function(d){return x(d.time);})
                          .y(function(d){return y(d.haz);})
                          .interpolate('step-after');

        // Add path layer
        this.vis.selectAll('.line')
                .data(data)
                .enter().append('path')
                 .attr('class', 'line')
		 .style('stroke', function(d,i){return colors[i];})
                 .attr('d', this.line);

	// Add legend, if there are groups

	var group_names = this.group_names;

	if(!(group_names == "")){
		var legend = this.vis.append('g')
			  .attr('class', 'legend')
			  .attr('x', this.w - 165)
			  .attr('y', 125)
			  .attr('height', 200)
			  .attr('width', 200);

		legend.selectAll('rect')
		   .data(colors).enter().append('rect')
		  .attr('x', this.w - 165)
		  .attr('y', function(d,i){return i*20;})
		  .attr('width', 10)
		  .attr('height', 10)
		  .style('fill', function(d) { return d; });

		legend.selectAll('text')
		   .data(group_names).enter().append('text')
		  .attr('x', this.w - 145)
		  .attr('y', function(d,i){return i*20 + 10;})
		  .text(function(d) { return d; });
	}
    };

    this.update_covar = function(newcov){
        for (var i=0; i<this.coef_names.length; i++) {
            this.covar[this.coef_names[i]]=0;
        }

        for(var j=0; j < newcov.length; j++){
            if(this.vtype[newcov[j].name] == 'continuous'){
                this.covar[newcov[j].name] = parseFloat(newcov[j].value);
            } else {
                this.covar[(newcov[j].name+newcov[j].value)]=1;
            }
        }
    };

    this.update = function(newcov) {
        this.update_covar(newcov);

        var tmp = update_hazard(this.init_data, this.coef, this.covar);
        this.vis.selectAll('path.line')
            .data(tmp)
            .transition().duration(1800).delay(100).ease('elastic')
            .attr('width', 0)
            .attr('d', this.line);

        var x = this.x;
        var y = this.y;
    };
}

This is an excerpt from the d3 code required to generate the figure seen in the healthvis example (actually requires around 60 lines). This does not include the fitting of the proportional hazards model, which is done in R. We feel that this approach probably requires the most effort, which was part of the motivation behind developing the healthvis package.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s