nievergeltlab
7/13/2017 - 8:49 PM
How competing Cox models of risk of relapse can be compared using ROC analysis. This analysis uses an implementation of the permutation meth
How competing Cox models of risk of relapse can be compared using ROC analysis. This analysis uses an implementation of the permutation methods suggested for comparing competing models
###################################
### Load R libraries
#Contains survival analysis function
library(survival)
#Functions to build ROC curves from Cox models
library(risksetROC)
###################################
### Data loading
#Function to quickly read .csv files
rcsv <- function(...)
{
read.csv(..., header=T, stringsAsFactors=FALSE,na.strings=c("NA", "#N/A", "#NULL!"))
}
#Read the .csv files, of subject data and survival times
covs <- rcsv('miestimation/pgbd_mi_data_imputed-stata12.csv') #All predictor variables. Variables with imputations have 35 entries each.
surv <- rcsv('pgbdsurvivalroc.csv')
#Merge data based on ID.
#Note: all.y makes sure that all subjects in the surv dataframe are included here, whether or not they have data in the covs datasheet
#Note: I set it so that if two variables have the same name, the variables from the covs sheet will have the suffix _covs
dat <- merge(covs,surv,by="subject_id",all.y=TRUE,suffixes=c("_covs",""))
###################################
### Data handling (Make it usable for the analysis)
##Per talking with Peter, we just want to take an average of the multiply imputed variables to use as our data
#This is a list of all imputed variables.
#The scan function takes user input and saves it as a vector
#This is the lazy way of doing it because you can just copy/paste lists from .csv files
imputed_vars_cat <- scan(what="character",sep=",")
childhoodphysabuse_binary,chronicaffdis_bin,comorbidpersonality2,functmorbid_bin,migraines2,mixedepisodes_collapsed,suicidethobehav_collapsed
imputed_vars_cont <- scan(what="character",sep=",")
anxietybaseline_continuous,negative_leq,sheehanimpairment,sheehantotal,suicidethobehav_collapsed
#This loop will go through the list of imputed variables
#and for each subject, take the median value
#iv is an incrementor variable
for (iv in imputed_vars_cont)
{
#For an imputed variable, find the column number of all imputations
#This is done by grepping for the variable in the column names of the data
#note: I prepend a _ before iv, so the original, unimputed data does not come up in the grep result
var_location <- grep(paste("_",iv,sep=""),names(dat))
ivdat <- dat[,var_location]
#apply: take the input data (ivdat), and performs a function (in this case median) on every row (as indicated by the 1 value)
#Will return a vector of length
imputed_collapsed <- apply(ivdat,1,median)
#add the vector of
dat <- cbind(dat,imputed_collapsed)
names(dat)[dim(dat)[2]] <- paste(iv,"i",sep="_")
}
Mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
for (iv in imputed_vars_cat)
{
#For an imputed variable, find the column number of all imputations
#This is done by grepping for the variable in the column names of the data
#note: I prepend a _ before iv, so the original, unimputed data does not come up in the grep result
var_location <- grep(paste("_",iv,sep=""),names(dat))
ivdat <- dat[,var_location]
#apply: take the input data (ivdat), and performs a function (in this case median) on every row (as indicated by the 1 value)
#Will return a vector of length
imputed_collapsed <- apply(ivdat,1,Mode)
#add the vector of
dat <- cbind(dat,imputed_collapsed)
names(dat)[dim(dat)[2]] <- paste(iv,"i",sep="_")
}
#subset data to relevant covariates only
dat_sub <- subset(dat, select=c(time,event,age,sex,racecat,li_status_collapsed,anxietybaseline_continuous_i,
chronicaffdis_bin_i,migraines2_i, suicidethobehav_collapsed_i,
mixedepisodes_collapsed_i,comorbidpersonality2_i , sheehantotal_i,
negative_leq_i))
###################################
### Permutation and analysis code
#Note: As of Feb 6, 2017, the model reported in the paper is the 730 day, full model
### turn off resampling the phenotype to get the baseline analysis, which should be 0.75240930 0.67311803 0.07929127 if the seed is set to 17
do_not_permute = 0 #Set to 1 if you want to turn off permutations and run the default analysis
full_model = 1 # Set to 1 to include the sheehan and others in the full model (limited set of clinical predictors or a really expanded set of clinical predictors)
null_baseline = 1 # Set to 1 if you want to compare to a model with aboslutely no covariates
#Set number of permutations
nperm=10000
if (do_not_permute == 1)
{
nperm= 1
}
#Makes a matrix that contains number of repetitions dimension, ROC for full, reduced, and difference models, for 365 and 730 days
roc_outputs <- data.frame(matrix(data=NA,nrow=nperm,ncol=6))
names(roc_outputs) <- c("AUC_f_365","AUC_r_365","AUC_d_365","AUC_f_730","AUC_r_730","AUC_d_730")
#Print a status update for the first or every 1000th permutation
for (permnum in 1:nperm)
{
if(permnum %% 1000 == 0 | permnum == 1)
{
print(paste("permutation run", permnum))
}
##Permute the phenotype
#Set the seed for random sampling
set.seed(permnum)
#There are N = dim(dat_sub)[1] rows in the data. Sampling numbers 1:N without
#replacement N times will give us a randomly ordered list of row numbers
shuffled <- sample(1:dim(dat_sub)[1],dim(dat_sub)[1],replace=F)
#Make a permuted dataframe by loading the data in the shuffled order for the
#phenotype column (columns 1 and 2) and in the typical order for all other covariates
dat_sub_shuffle <- data.frame(dat_sub[shuffled,c(1:2)],dat_sub[,-c(1:2)]) #For use in permutation analysis
##Open issue: does it matter or not if I don't reset the row names?
#If you said to not permute, this uses the original data frame instead of the shuffled one
if (do_not_permute == 1)
{
dat_sub_shuffle <- dat_sub
}
#Randomly reorder the data, so that training/test sets contain different observations for each permutation
set.seed(permnum) #Note: set to 17 for the un-permuted data to recreate original results set.seed=17
if (do_not_permute == 1)
{
set.seed(17)
}
new_order <- sample(x=1:dim(dat_sub)[1], size=dim(dat_sub)[1], replace=F) #Defines the new ordering
dat_r <- dat_sub_shuffle[new_order,] #Accomplishes the reordering
#Determine which subjects have a survival event, and which do not
#the which function returns the row number of these subjects
cases <- which(dat_r$event == 1)
controls <- which(dat_r$event == 0)
#Make separate dataframes for cases and controls
#I do this because it's an easy way to insure that
#the training/test sets contain the correct proportion of cases
case_set <- dat_r[cases,]
control_set <- dat_r[controls,]
#Repeat the sequence 1...10
numberseq <- rep(1:10,dim(dat_r)[1])
#Cut the sequence down to the dimension of the datasets
#Now every subject will have been assigned a number from 1 to 10
#We will filter on this number to build training/test sets
control_numbers <- numberseq[1:length(controls)] #Give each control subject a set number
case_numbers <- numberseq[1:length(cases)] #Give each case subject a set number
###Perform cox regression on each subset of data
#Make a matrix to hold the predicted AUCs
#Make it into a data frame so variables can be named easily
pred_base_lps <- data.frame(matrix(data=NA,nrow=dim(dat_r)[1],ncol=2))
names(pred_base_lps) <- c("red","full")
row.names(pred_base_lps) <- paste("REM",1:dim(pred_base_lps)[1],sep="") # Have to remove the existing row names because we're going to merge on row names, which have been shuffled.
#Formula for basic model
if (null_baseline == 0)
{
red_model_formula <- formula(Surv(time,event) ~ age + sex + as.factor(racecat) + as.factor(li_status_collapsed))
}
if (null_baseline == 1)
{
red_model_formula <- formula(Surv(time,event) ~ 1)
}
#Formula for model with clinical predictors
#If choosing to include sheehan, negative leq, and chronicity
if (full_model ==1)
{
model_formula <- formula(Surv(time,event) ~ age + sex + as.factor(racecat) + as.factor(li_status_collapsed) + anxietybaseline_continuous_i +
migraines2_i + as.factor(suicidethobehav_collapsed_i) + mixedepisodes_collapsed_i +
comorbidpersonality2_i + sheehantotal_i + negative_leq_i + chronicaffdis_bin_i)
}
#If choosing to NOT include sheehan, negative leq, and chronicity
if (full_model ==0)
{
model_formula <- formula(Surv(time,event) ~ age + sex + as.factor(racecat) + as.factor(li_status_collapsed) + anxietybaseline_continuous_i +
migraines2_i + as.factor(suicidethobehav_collapsed_i) + mixedepisodes_collapsed_i +
comorbidpersonality2_i )
}
incstart <- 1 # Need to increment this matrix of AUCs based on the dimension of the test sets, which are not uniformly sized
for (i in 1:10)
{
train_set <- rbind(case_set[case_numbers != i,], control_set[control_numbers != i,]) #Get all cases and controls who aren't in the ith set. Combine the set of cases and controls. Regression will be done in these subjects
test_set <- rbind(case_set[case_numbers == i,], control_set[control_numbers == i,]) #Get the ith set to make the testing set
coxmod_base <- coxph(red_model_formula,data=train_set) #Fit the reduced cox model to the training set
coxmod_full <- coxph(model_formula, data= train_set) #Fit the full cox model to the training set
pred_base <- predict(coxmod_base,test_set,"lp") #obtain linear predictors for reduced model on the test set
pred_full <- predict(coxmod_full,test_set,"lp") #obtain linear predictors for full model on the test set
inc <- length(pred_base)
pred_base_lps[incstart:(incstart+inc-1),]$red <- pred_base
pred_base_lps[incstart:(incstart+inc-1),]$full <- pred_full
row.names(pred_base_lps)[incstart:(incstart+inc-1)] <- names(pred_full) #Set teh row names to match the existing data row names
incstart <- incstart + inc #Increment the counter
}
#Based on the row names, merge back with the original data
dmerge <- merge(pred_base_lps,dat_r,by="row.names")
#Now take the ROC from both models
#Do for each cutoff date
t1r365 <- CoxWeights(marker=dmerge$red,Stime = dmerge$time, status=dmerge$event, predict.time=365)
t1f365 <- CoxWeights(marker=dmerge$full,Stime = dmerge$time, status=dmerge$event, predict.time=365)
rocd365 <- t1f365$AUC-t1r365$AUC
t1r730 <- CoxWeights(marker=dmerge$red,Stime = dmerge$time, status=dmerge$event, predict.time=730)
t1f730 <- CoxWeights(marker=dmerge$full,Stime = dmerge$time, status=dmerge$event, predict.time=730)
rocd730 <- t1f730$AUC-t1r730$AUC
roc_outputs[permnum,] <- c(t1f365$AUC, t1r365$AUC, rocd365,t1f730$AUC, t1r730$AUC, rocd730)
}
##Handle the output depending on parameters
if (do_not_permute == 1 & full_model == 1 & null_baseline == 1)
{
t1r730_unpermuted_allpreds <- t1r730
t1f730_unpermuted_allpreds <- t1f730
t1r365_unpermuted_allpreds <- t1r365
t1f365_unpermuted_allpreds <- t1f365
roc_outputs_unpermuted_allpreds <- roc_outputs
save(t1r730_unpermuted_allpreds,t1f730_unpermuted_allpreds,t1r365_unpermuted_allpreds,t1f365_unpermuted_allpreds,roc_outputs_unpermuted_allpreds, file="roc_outputs_unpermuted_allpreds_nullbaseline_v1_feb6_2017.R")
}
table(0.7489 > roc_outputs_unpermuted_allpreds[,6])
if (do_not_permute == 0 & full_model == 1 & null_baseline == 1)
{
t1r730_unpermuted_allpreds <- t1r730
t1f730_unpermuted_allpreds <- t1f730
t1r365_unpermuted_allpreds <- t1r365
t1f365_unpermuted_allpreds <- t1f365
roc_outputs_unpermuted_allpreds <- roc_outputs
save(t1r730_unpermuted_allpreds,t1f730_unpermuted_allpreds,t1r365_unpermuted_allpreds,t1f365_unpermuted_allpreds,roc_outputs_unpermuted_allpreds, file="roc_outputs_permuted_allpreds_nullbaseline_v1_feb6_2017.R")
}
if (do_not_permute == 1 & full_model == 1 & null_baseline == 0)
{
t1r730_unpermuted_allpreds <- t1r730
t1f730_unpermuted_allpreds <- t1f730
t1r365_unpermuted_allpreds <- t1r365
t1f365_unpermuted_allpreds <- t1f365
roc_outputs_unpermuted_allpreds <- roc_outputs
save(t1r730_unpermuted_allpreds,t1f730_unpermuted_allpreds,t1r365_unpermuted_allpreds,t1f365_unpermuted_allpreds,roc_outputs_unpermuted_allpreds, file="roc_outputs_unpermuted_allpreds_v1_aug10_2016.R")
}
if (do_not_permute == 1 & full_model == 1 & null_baseline == 1)
{
t1r730_unpermuted_allpreds <- t1r730
t1f730_unpermuted_allpreds <- t1f730
t1r365_unpermuted_allpreds <- t1r365
t1f365_unpermuted_allpreds <- t1f365
roc_outputs_unpermuted_allpreds <- roc_outputs
save(t1r730_unpermuted_allpreds,t1f730_unpermuted_allpreds,t1r365_unpermuted_allpreds,t1f365_unpermuted_allpreds,roc_outputs_unpermuted_allpreds, file="roc_outputs_unpermuted_allpreds_nullbaseline_v1_feb6_2017.R")
}
if (do_not_permute == 1 & full_model == 0 & null_baseline == 0)
{
t1r730_unpermuted_lesspreds <- t1r730
t1f730_unpermuted_lesspreds <- t1f730
t1r365_unpermuted_lesspreds <- t1r365
t1f365_unpermuted_lesspreds <- t1f365
roc_outputs_unpermuted_lesspreds <- roc_outputs
save(t1r730_unpermuted_lesspreds,t1f730_unpermuted_lesspreds,t1r365_unpermuted_lesspreds,t1f365_unpermuted_lesspreds,roc_outputs_unpermuted_lesspreds, file="roc_outputs_unpermuted_lesspreds_v1_aug10_2016.R")
}
if (do_not_permute == 0 & full_model == 1 & null_baseline == 0)
{
t1r730_permed_allpreds <- t1r730
t1f730_permed_allpreds <- t1f730
t1r365_permed_allpreds <- t1r365
t1f365_permed_allpreds <- t1f365
roc_outputs_permed_allpreds <- roc_outputs
save(t1r730_permed_allpreds,t1f730_permed_allpreds,t1r365_permed_allpreds,t1f365_permed_allpreds,roc_outputs_permed_allpreds, file="roc_outputs_permed_allpreds_v1_aug10_2016.R")
}
if (do_not_permute == 0 & full_model == 0 & null_baseline == 0)
{
t1r730_permed_lesspreds <- t1r730
t1f730_permed_lesspreds <- t1f730
t1r365_permed_lesspreds <- t1r365
t1f365_permed_lesspreds <- t1f365
roc_outputs_permed_lesspreds <- roc_outputs
save(t1r730_permed_lesspreds,t1f730_permed_lesspreds,t1r365_permed_lesspreds,t1f365_permed_lesspreds,roc_outputs_permed_lesspreds, file="roc_outputs_permed_lesspreds_v1_aug10_2016.R")
}
load('roc_outputs_permed_allpreds_v1_aug10_2016.R')
load('roc_outputs_unpermuted_allpreds_v1_aug10_2016.R')
load('roc_outputs_permed_lesspreds_v1_aug10_2016.R')
load('roc_outputs_unpermuted_lesspreds_v1_aug10_2016.R')
pdf("PGBD_ROC_allpredictors_730days_v1_aug10_2016.pdf")
plot(t1f730_unpermuted_allpreds$FP,t1f730_unpermuted_allpreds$TP,type="l",lwd=2,cex.axis=1.45,main="",col="blue",xlab="False Positive Rate (1 - Specificity)",ylab="True Positive Rate (Sensitivity)")
lines(t1r730_unpermuted_allpreds$FP,t1r730_unpermuted_allpreds$TP,type="l",lwd=2,cex.axis=1.45,main="",col="red")
abline(0,1,col="black",lty=2)
pv <- table(roc_outputs_unpermuted_allpreds[,6] > roc_outputs_permed_allpreds[,6])[1] / length( roc_outputs_permed_allpreds[,6])
legend("topleft", legend=c(paste("Full model AUC =", round(t1f730_unpermuted_allpreds$AUC,3)),paste("Reduced model AUC =", round(t1r730_unpermuted_allpreds$AUC,3)),paste("P-value =", pv)),col=c("blue","red"),pch=c(18,18,NA))
dev.off()
pdf("PGBD_ROC_allpredictors_365days_v1_aug10_2016.pdf")
plot(t1f365_unpermuted_allpreds$FP,t1f365_unpermuted_allpreds$TP,type="l",lwd=2,cex.axis=1.45,main="",col="blue",xlab="False Positive Rate (1 - Specificity)",ylab="True Positive Rate (Sensitivity)")
lines(t1r365_unpermuted_allpreds$FP,t1r365_unpermuted_allpreds$TP,type="l",lwd=2,cex.axis=1.45,main="",col="red")
abline(0,1,col="black",lty=2)
pv <- table(roc_outputs_unpermuted_allpreds[,3] > roc_outputs_permed_allpreds[,3])[1] / length( roc_outputs_permed_allpreds[,3])
legend("topleft", legend=c(paste("Full model AUC =", round(t1f365_unpermuted_allpreds$AUC,3)),paste("Reduced model AUC =", round(t1r365_unpermuted_allpreds$AUC,3)),paste("P-value =", pv)),col=c("blue","red"),pch=c(18,18,NA))
dev.off()
pdf("PGBD_ROC_lesspredictors_730days_v1_aug10_2016.pdf")
plot(t1f730_unpermuted_lesspreds$FP,t1f730_unpermuted_lesspreds$TP,type="l",lwd=2,cex.axis=1.45,main="",col="blue",xlab="False Positive Rate (1 - Specificity)",ylab="True Positive Rate (Sensitivity)")
lines(t1r730_unpermuted_lesspreds$FP,t1r730_unpermuted_lesspreds$TP,type="l",lwd=2,cex.axis=1.45,main="",col="red")
abline(0,1,col="black",lty=2)
pv <- table(roc_outputs_unpermuted_lesspreds[,6] > roc_outputs_permed_lesspreds[,6])[1] / length( roc_outputs_permed_lesspreds[,6])
legend("topleft", legend=c(paste("Full model AUC =", round(t1f730_unpermuted_lesspreds$AUC,3)),paste("Reduced model AUC =", round(t1r730_unpermuted_lesspreds$AUC,3)),paste("P-value =", pv)),col=c("blue","red"),pch=c(18,18,NA))
dev.off()
pdf("PGBD_ROC_lesspredictors_365days_v1_aug10_2016.pdf")
plot(t1f365_unpermuted_lesspreds$FP,t1f365_unpermuted_lesspreds$TP,type="l",lwd=2,cex.axis=1.45,main="",col="blue",xlab="False Positive Rate (1 - Specificity)",ylab="True Positive Rate (Sensitivity)")
lines(t1r365_unpermuted_lesspreds$FP,t1r365_unpermuted_lesspreds$TP,type="l",lwd=2,cex.axis=1.45,main="",col="red")
abline(0,1,col="black",lty=2)
pv <- table(roc_outputs_unpermuted_lesspreds[,3] > roc_outputs_permed_lesspreds[,3])[1] / length( roc_outputs_permed_lesspreds[,3])
legend("topleft", legend=c(paste("Full model AUC =", round(t1f365_unpermuted_lesspreds$AUC,3)),paste("Reduced model AUC =", round(t1r365_unpermuted_lesspreds$AUC,3)),paste("P-value =", pv)),col=c("blue","red"),pch=c(18,18,NA))
dev.off()
#No Sheehan
table(roc_outputs_unpermuted_lesspreds[,1] > roc_outputs_lesspreds[,1])
table(roc_outputs_unpermuted_lesspreds[,2] > roc_outputs_lesspreds[,2])
table(roc_outputs_unpermuted_lesspreds[,3] > roc_outputs_lesspreds[,3])
pdf("roc_unpermuted_lesspreds_730days_v1_aug1_2016.pdf")
plot(t1f730$FP,t1f730$TP,type="l",lwd=2,cex.axis=1.45,main="",col="blue")
lines(t1r730$FP,t1r730$TP,type="l",lwd=2,cex.axis=1.45,main="",col="red")
abline(0,1,col="black",lty=2)
dev.off()
plot(t1f365$FP,t1f365$TP,type="l",lwd=2,cex.axis=1.45,main="",col="blue")
lines(t1r365$FP,t1r365$TP,type="l",lwd=2,cex.axis=1.45,main="",col="red")
abline(0,1,col="black",lty=2)
#Default AUC values without the
par(mfrow=c(2,2))
hist(roc_outputs[,1],freq=FALSE,xlim=c(0.5,.8),main="Full model AUCs")
abline(v= 0.75240930,col="red",lwd=2)
hist(roc_outputs[,2],freq=FALSE,xlim=c(0.5,.8),main="Reduced model AUCs")
abline(v= 0.67311803,col="red",lwd=2)
hist(roc_outputs[,3],freq=FALSE,xlim=c(-.2,.2),main="Full - Reduced AUCs")
abline(v= 0.07929127,col="red",lwd=2)
dev.off()
#730 days
#No Sheehan
table(roc_outputs_unpermuted_lesspreds[,4] > roc_outputs_permed_lesspreds[,4])
table(roc_outputs_unpermuted_lesspreds[,5] > roc_outputs_permed_lesspreds[,5])
table(roc_outputs_unpermuted_lesspreds[,6] > roc_outputs_permed_lesspreds[,6])
#With Sheehan
table(roc_outputs_unpermuted_allpreds[,4] > roc_outputs_permed_allpreds[,4])
table(roc_outputs_unpermuted_allpreds[,5] > roc_outputs_permed_allpreds[,5])
table(roc_outputs_unpermuted_allpreds[,6] > roc_outputs_permed_allpreds[,6])
#365 days
#No Sheehan
table(roc_outputs_unpermuted_lesspreds[,1] > roc_outputs_permed_lesspreds[,1])
table(roc_outputs_unpermuted_lesspreds[,2] > roc_outputs_permed_lesspreds[,2])
table(roc_outputs_unpermuted_lesspreds[,3] > roc_outputs_permed_lesspreds[,3])
#With Sheehan
table(roc_outputs_unpermuted_allpreds[,1] > roc_outputs_permed_allpreds[,1])
table(roc_outputs_unpermuted_allpreds[,2] > roc_outputs_permed_allpreds[,2])
table(roc_outputs_unpermuted_allpreds[,3] > roc_outputs_permed_allpreds[,3])
#Get correlation of measures
library(psych)
corrs <- corr.test(dat_sub[,c("anxietybaseline_continuous_i",
"chronicaffdis_bin_i","migraines2_i", "suicidethobehav_collapsed_i",
"mixedepisodes_collapsed_i","comorbidpersonality2_i" , "sheehantotal_i",
"negative_leq_i")],method="spearman")