Fitting a Survival Model to Simulated Data

This section details how the language for statistical programming ’R’ (R Development Core Team 2011) can be used to fit a Cox proportional hazards model to simulated data.

The following Rscript code was used to generate the results presented in Section 2.2.4:

  #!/usr/bin/Rscript
  library(survival)

  # Set parameter values
  n_components <- 1000
  n_hours_observed <- 45
  failure_rate <- 1/60
  a1 <- 0.05
  a2 <- 0.1

  # Generate simulated data
  observations <- matrix(0,n_components*n_hours_observed,5)
  colnames(observations) <- c(’start’,’stop’,’z1’,’z2’,’failure’)
  row <- 0
  for (c in 1:n_components) {
      for (h in 1:n_hours_observed) {
       row <- row+1
       start <- h-1
   stop <- h
   z1 <- abs(rnorm(1,0,10))
   z2 <- runif(1,0,5)
   if (failure_rate*exp(a1*z1+a2*z2) >= runif(1)) {
      # failure occured
      observations[row,] = c(start,stop,z1,z2,1)
      break
   } else {
      observations[row,] = c(start,stop,z1,z2,0)
   }
      }
  }
  observations <- as.data.frame(observations)

  # Fit a Cox proportional hazards model to the data
  fitted_model <- coxph(Surv(start,stop,failure) ~ z1 + z2,
                          data=observations[1:row,])

  # Print summary of fitted model
  summary(fitted_model)

  # Create a data frame with covariate values to use when
  # plotting the estimated survival function
  attach(observations)
  observations.eliminated_covariates <- data.frame(z1=0,z2=0)
  detach()

  # Plot the survival function
  png(filename="data_fitting_example.png",width=800,height=600)
  plot(survfit(fitted_model, newdata=observations.eliminated_covariates),
       xlab="Hours of observation",
       ylab="Proportion not failed")
  dev.off()

3.1 Fitting a Survival Model to Simulated Data