set.seed(1298) # ------------------------------------------------------------------------------ # Index # 0. Preliminaries # 1. Reading in data and preprocessing # 2. Model development # 2.1 Selecting a random side to use for HC # 2.2 Correcting miRNA expression levels for different variables using linear regression # 2.3 Box-Cox # 2.4 Principal component dimensionality reduction # 3. Model fitting and evaluation (via cross validation): individual-level analysis # 4. Breast-level analysis # ------------------------------------------------------------------------------ # Preliminaries library(glmnet) #library(survival) #library(penalized) getSensSpec <- function(obs,pred,weights=rep(1L,length(obs)),threshold=NULL,dir="<"){ if(is.null(threshold)) threshold <- c(-Inf,sort(unique(pred)),Inf) n <- length(threshold) x <- matrix(rep(threshold,each=length(pred)),ncol=n) obsn <- matrix(weights*!obs,nrow=length(obs),ncol=n) obs <- matrix(weights*obs,nrow=length(obs),ncol=n) pred <- matrix(pred,nrow=length(pred),ncol=n) # see ?pROC::roc under "Thresholds" for details on semantics # sens = Pr(pred>threshold|obs==1) # spec = Pr(pred<=threshold|obs==0) sens <- if(dir==">") colSums(obs*(pred<=x))/colSums(obs) else colSums(obs*(pred>x))/colSums(obs) spec <- if(dir==">") 1-colMeans((obsn)*(pred<=x))/colMeans(obsn) else 1-colMeans((obsn)*(pred>x))/colMeans(obsn) data.frame(threshold,sens,spec) } table_wtd <- function(...,weights=NULL,relative=F){ x <- list(...) n <- length(x[[1]]) if(is.null(weights)) weights <- rep(1L,n) if(any(lengths(x)!=length(weights))) stop("all arguments must have the same length.") p <- do.call(paste,x) o <- order(p) q <- p[o] weights <- weights[o] r <- rev(q) s <- rev(c(TRUE,r[-1L]!=r[-n])) u <- q[s] w <- cumsum(weights) y <- w[s] f <- c(y[1L],diff(y)) out <- cbind(do.call(data.frame,x)[match(u,p),,drop=F],frequency=f) nms <- as.list(match.call())[-1L][seq_along(x)] nms[names(nms)!=""] <- names(nms)[names(nms)!=""] names(out) <- c(as.character(nms),"Frequency") if(relative) out$Frequency <- out$Frequency/sum(out$Frequency) return(out) } rank_wtd <- function(x,weights=NULL,ties.method="average"){ tab <- table_wtd(x,weights=weights) tab <- tab[order(tab$x),,drop=F] switch(ties.method,stop("'ties.method' not supported"), average={ r <- with(tab,cumsum(Frequency)-.5*(Frequency-1)) r[match(x,tab$x)] }, min={ r <- with(tab,cumsum(Frequency)-(Frequency-1L)) r[match(x,tab$x)] }, max={ r <- with(tab,cumsum(Frequency)) r[match(x,tab$x)] } ) } auc_wtd <- function(x,y,weights=rep(1L,length(x))){ # cf. Hmisc::somers2(x,y)['C'] if(any(!{y %in% 0:1})) stop('"y" must be binary') n <- sum(weights*y) m <- sum(weights*!y) r <- rank_wtd(x,weights=weights,ties.method="average") (sum(weights*y*r)-n*(n+1)/2)/(n*m) } largestNoLarger <- function(x,y,smallestNoSmaller=FALSE){ # for every element y[i] of y, # identify j such that x[j] is the (last) largest element in x that is no larger than y[i] if(smallestNoSmaller){ return(largestNoLarger(-x,-y)) # identify j such that x[j] is the (last) smallest element in x that is no larger than y[i] } n <- length(x) m <- length(y) group <- c(rep(0L,n),rep(1L,m)) id <- seq_along(group) z <- c(x,y) o <- order(z) i <- match(id,o) group <- group[o] id <- id[o] z <- z[o] a <- cumsum(!group) b <- match(a,a) w <- id[b][i][group[i]==1L] w[w>n] <- NA return(w) } interpolate <- function(x,y,xout){ # assume y~x monotonic; x,y submitted s.t. x in increasing order w <- largestNoLarger(x,xout,F) yd <- ifelse(w==length(x),rev(y)[1],y[w+1])-y[w] xd <- ifelse(w==length(x),max(x)+1,x[w+1])-x[w] yout <- y[w]+(xout-x[w])*yd/xd return(yout) } boxcox <- function(x,nonzeroLambda=T,lambda=NULL){ #https://www.r-bloggers.com/2022/10/box-cox-transformation-in-r/ if(is.null(lambda)){ b <- MASS::boxcox(x~1,plot=F,lambda=seq(-2,2,length=1000)) if(!nonzeroLambda){ ci <- range(b$x[b$y>=max(b$y)-1/2*qchisq(.95,1)]) if(min(ci)<0&035] <- 35 } # Recalculate the DCT values for(cl in colsCT_T) data[[paste0("D",cl)]] <- data[[cl]]-data[["CT_miR_99b_5p_T"]] for(cl in colsCT_CL) data[[paste0("D",cl)]] <- data[[cl]]-data[["CT_miR_99b_5p_CL"]] colsDCT_T <- paste0("D",colsCT_T) colsDCT_CL <- paste0("D",colsCT_CL) colsDCT <- c(colsDCT_T,colsDCT_CL) # ------------------------------------------------------------------------------ # 2. Model development develop <- function( data, response="HC_BC_HR", formula=model1, noiseReduction=T,#cloudy & RNA input correctie boxCox=T, PCAoption=4, # PCAoption==1-3: not used for final analysis # PCAoption==4: include all miRNAs, after noise red. and box-cox but add least number of PC's that explain >= 80% of variance, rather than add fixed number of PC's LASSO=T ){ # 3.1 Selecting a random side to use for HC # -------------------------------------------- newvar <- gsub("_T$","",colsDCT_T) newvar <- newvar[!grepl("miR_(126|99b|151a|18a)",newvar)] healthy_control <- !data[[outcome]] w <- sample(1:2,nrow(data),replace=T) chosen_side <- as.data.frame(matrix(nrow=nrow(data),ncol=length(newvar))) colnames(chosen_side) <- newvar for(var in newvar){ x <- data[paste0(var,c("_T","_CL"))] chosen_side[[var]] <- ifelse(healthy_control,ifelse(is.na(x[cbind(1:nrow(data),w)]),3-w,w),1) data[[var]] <- x[cbind(1:nrow(data),chosen_side[[var]])] } colsDCT <- newvar # -------------------------------------------- # 2.2 Correcting miRNA expression levels for different variables using linear regression # -------------------------------------------- depvar <- colsDCT df <- data[depvar] colnames(df) <- depvar df$Subject_nr <- data$Subject_nr Cloudiness <- with(data,cbind(TandHC_cloudiness,CL_cloudiness)) RNA_input <- with(data,cbind(TandHC_RNA_input,CL_RNA_input)) noise_reduction_model_fits <- vector("list",length(depvar)) names(noise_reduction_model_fits) <- depvar if(noiseReduction){ for(var in depvar){ w <- cbind(1:nrow(data),chosen_side[[var]]) df$Cloudiness <- Cloudiness[w] df$RNA_input <- RNA_input[w] noise_reduction_model_fits[[var]] <- if(!noiseReduction) lm(df[[var]]~1,data=df) else lm(df[[var]]~as.factor(Cloudiness)+as.factor(RNA_input),data=df) data[[var]] <- if(!noiseReduction) data[[var]] else resid(noise_reduction_model_fits[[var]]) } } colsRES_DCT <- colsDCT # -------------------------------------------- # 2.3 Box-Cox # -------------------------------------------- boxcox_trans <- vector("list",length(colsRES_DCT)) names(boxcox_trans) <- colsRES_DCT for(cl in colsRES_DCT){ if(!boxCox){ boxcox_trans[[cl]]$m <- -Inf boxcox_trans[[cl]]$k <- 1 boxcox_trans[[cl]]$lambda <- 1 next } x <- data[[cl]] if(all(is.na(x))){ boxcox_trans[[cl]]$m <- 0 boxcox_trans[[cl]]$k <- 0 } else if(any(x<=0,na.rm=T)){ m <- min(abs(x),na.rm=T)/2 k <- abs(min(x,na.rm=T))+m if(k==0) k <- 1 boxcox_trans[[cl]]$k <- k boxcox_trans[[cl]]$m <- m } else{ boxcox_trans[[cl]]$k <- 0 boxcox_trans[[cl]]$m <- 0 } lambda <- if(all(is.na(x))) NA else attr(boxcox(x+boxcox_trans[[cl]]$k),"lambda") boxcox_trans[[cl]]$lambda <- lambda } if(boxCox){ for(cl in colsRES_DCT){ data[[cl]] <- as.numeric(boxcox(data[[cl]]+boxcox_trans[[cl]]$k,lambda=boxcox_trans[[cl]]$lambda)) } } colsRES_BC_DCT <- colsRES_DCT # -------------------------------------------- # 2.4 Principal component dimensionality reduction # -------------------------------------------- selectPCA <- function(pca,thresh=0.8,returnWhichOnly=F){ vv <- diag(cov(pca)) tt <- sum(vv) n <- which(cumsum(vv/tt)>=thresh)[1L] if(returnWhichOnly) return(seq_len(n)) out <- pca[,seq_len(n),drop=FALSE] attr(out,"components") <- seq_len(n) attr(out,"cum_var_explained") <- cumsum(vv/tt) return(out) } if(PCAoption==1){ x <- data[colsRES_BC_DCT[!grepl("miR_(125|145|21|29c)",colsRES_BC_DCT)]] pca <- prcomp(x,scale.=T) pca_orig <- pca$x data <- cbind(data,pca_orig) rot <- pca$rotation contribution <- sweep(abs(rot),2,colSums(abs(rot)),"/") rownames(contribution) <- colnames(x) contribution } else if(PCAoption==2){ x <- data[colsRES_BC_DCT] pca <- prcomp(x,scale.=T) pca_orig <- pca$x data <- cbind(data,pca_orig) rot <- pca$rotation contribution <- sweep(abs(rot),2,colSums(abs(rot)),"/") rownames(contribution) <- colnames(x) contribution PComponents <- selectPCA(pca$x,thresh=0.999999999999999,returnWhichOnly=T) form <- as.character(formula) form[3] <- gsub("(~|\\+)( |)PC[:0-9:]+","",form[3]) form[3] <- paste0(form[3]," + ",paste(paste0("PC",PComponents),collapse=" + ")) formula <- as.formula(paste(form[2],form[1],form[3])) } else if(PCAoption==3){ x <- data[colsRES_BC_DCT[!grepl("miR_(125|145|21|29c)",colsRES_BC_DCT)]] pca <- prcomp(x,scale.=T) pca_orig <- pca$x data <- cbind(data,pca_orig) rot <- pca$rotation contribution <- sweep(abs(rot),2,colSums(abs(rot)),"/") rownames(contribution) <- colnames(x) contribution PComponents <- selectPCA(pca$x,thresh=.8,returnWhichOnly=T) form <- as.character(formula) form[3] <- gsub("(~|\\+)( |)PC[:0-9:]+","",form[3]) form[3] <- paste0(form[3]," + ",paste(paste0("PC",PComponents),collapse=" + ")) formula <- as.formula(paste(form[2],form[1],form[3])) } else if(PCAoption==4){ x <- data[colsRES_BC_DCT] pca <- prcomp(x,scale.=T) pca_orig <- pca$x data <- cbind(data,pca_orig) rot <- pca$rotation contribution <- sweep(abs(rot),2,colSums(abs(rot)),"/") rownames(contribution) <- colnames(x) contribution PComponents <- selectPCA(pca$x,thresh=0.8,returnWhichOnly=T) form <- as.character(formula) form[3] <- gsub("(~|\\+)( |)PC[:0-9:]+","",form[3]) form[3] <- paste0(form[3]," + ",paste(paste0("PC",PComponents),collapse=" + ")) formula <- as.formula(paste(form[2],form[1],form[3])) } else contribution <- NULL # -------------------------------------------- # MLE with/without LASSO # -------------------------------------------- if(LASSO){ lambdas <- exp(seq(2,-10,length=100)) mm <- model.matrix(formula,data=data)[,-1,drop=F] f <- cv.glmnet(x=mm,y=data[[response]],family="binomial",type.measure="auc", alpha=1,nfolds=10,lambda=lambdas,keep=T) # b <- coef(f,s="lambda.min") } else{ mm <- model.matrix(formula,data=data)[,-1,drop=F] fit <- glm(formula,family=binomial,data=data) b <- coef(fit) } predict <- function(newdata){ newvar <- gsub("_T$","",colsDCT_T) newvar <- newvar[!grepl("miR_(126|99b|151a|18a)",newvar)] healthy_control <- !newdata[[outcome]] w <- sample(1:2,nrow(newdata),replace=T) chosen_side <- as.data.frame(matrix(nrow=nrow(newdata),ncol=length(newvar))) colnames(chosen_side) <- newvar for(var in newvar){ x <- newdata[paste0(var,c("_T","_CL"))] chosen_side[[var]] <- ifelse(healthy_control,ifelse(is.na(x[cbind(1:nrow(newdata),w)]),3-w,w),1) newdata[[var]] <- x[cbind(1:nrow(newdata),chosen_side[[var]])] } # Noise reduction depvar <- colsDCT Cloudiness <- with(newdata,cbind(TandHC_cloudiness,CL_cloudiness)) RNA_input <- with(newdata,cbind(TandHC_RNA_input,CL_RNA_input)) for(var in depvar){ w <- cbind(1:nrow(newdata),chosen_side[[var]]) df <- data.frame(side=w) df$Cloudiness <- Cloudiness[w] df$RNA_input <- RNA_input[w] p <- if(!noiseReduction) 0 else stats::predict(noise_reduction_model_fits[[var]],newdata=df,type="response") o <- newdata[[var]] newdata[[var]] <- o-p } # Box-Cox for(cl in colsRES_DCT){ newdata[[cl]] <- boxcox(pmax(boxcox_trans[[cl]]$m,newdata[[cl]]+boxcox_trans[[cl]]$k),lambda=boxcox_trans[[cl]]$lambda) } # PCA pca_orig <- stats:::predict.prcomp(pca,newdata=newdata) newdata <- cbind(newdata,pca_orig) current.na.action <- options('na.action') options(na.action='na.pass') x <- model.matrix(formula,data=newdata) options(na.action=current.na.action$na.action) out <- plogis(as.vector(x%*%b)) attr(out,"model.matrix") <- x return(out) } attr(predict,"scalingParamPCA") <- list(mean=pca$center,sd=pca$scale) attr(predict,"coefPCA") <- pca$rotation attr(predict,"contributions") <- contribution if(PCAoption%in%c(3,4)) attr(predict,"numberOfPComponents") <- length(PComponents) attr(predict,"coef") <- b attr(predict,"scaled_coef") <- b[-1]*apply(mm,2,sd) return(predict) } plotContributionsPC <- function(x=attr(fit1,"contributions")*100,which=1:4){ y <- x[,which,drop=F] convert <- function(x,source_range,target_range){ source_range <- range(source_range) target_range <- range(target_range) min(target_range)+(x-min(source_range))*diff(target_range)/diff(source_range) } barplot <- function(i){ z <- y[,i] xlower=1:m+(1:m-1)*ysep; xupper=xlower-1 xrange=c(0,100); yrange=c(0,m+(m-1)*ysep) frame <- data.frame(x=c(0,1)+(i-1)*(1+xsep),y=c(0,m+(m-1)*ysep)) bars <- lapply(seq_along(z), function(j)data.frame(x=c(xlower[j],xupper[j])[c(1,1,2,2,1)],y=c(0,z[j])[c(1,2,2,1,1)],w=j)) bars <- do.call(rbind,bars) bars <- data.frame( x=convert(bars$y,xrange,frame$x), y=convert(bars$x,yrange,frame$y), w=bars$w ) bars <- split(bars,bars$w) for(j in seq_along(bars)) with(bars[[j]],polygon(x,y,col="white")) lbl <- paste0(formatC(z,format="f",digits=1),"%") a <- z+10>=60 txt <- data.frame(x=convert(1:m+(1:m-1)*ysep-.5,yrange,frame$y), y=convert(ifelse(a,z-5,z+5),xrange,frame$x)) txt$adj <- ifelse(a,1,0) txt$labels <- lbl for(j in 1:nrow(txt)) with(as.data.frame(txt[j,]),text(y=x,x=y,labels=labels,adj=c(adj,0.5),cex=.7)) } par0 <- par(no.readonly=T) par(mar=c(5,4,4,1)+.1-c(2,-1,3.5,0.5)) # default c(5,4,4,1)+.1 par(mgp=c(3,1,0)-c(1.3,.6,0)) # default c(3,1,0) par(tcl=-.3) # default -.5 n <- length(which) m <- nrow(y) xsep <- .01*n ysep <- .01*m plot(1,1,xlim=c(0,n+(n-1)*xsep),ylim=c(0,m+(m-1)*ysep),type="n",xlab="",ylab="",main="",font.lab=2,axes=F) axis(2,at=1:m+(1:m-1)*ysep-.5,labels=F,lwd=0,lwd.ticks=1) mtext(rownames(y),side=2,at=1:m+(1:m-1)*ysep-.5,srt=90,las=1,adj=1,line=.5) for(i in 1:n){ axis(1,at=c(0,1)+(i-1)*(1+xsep),lwd=0,lwd.ticks=1,labels=F) mtext(0,at=0+(i-1)*(1+xsep),adj=0,line=.5,side=1) mtext(100,at=1+(i-1)*(1+xsep),adj=1,line=.5,side=1) mtext(paste0("PC",i),at=.5+(i-1)*(1+xsep),adj=0.5,line=1.5,side=1,font=2) #abline(v=c(0,1)+(i-1)*(1+xsep),lty=3) polygon(x=(c(0,1)+(i-1)*(1+xsep))[c(1,2,2,1,1)],y=par()$usr[3:4][c(1,1,2,2,1)],col="lightgrey",border="transparent") } box() for(i in 1:n) barplot(i) } # ------------------------------------------------------------------------------ # 3. Model fitting and evaluation (via cross validation): individual-level analysis # model evaluation evaluate <- function(testdata,predict,response="HC_BC_HR",plot=F){ pred <- predict(testdata) obs <- testdata[[response]]*1 roc <- with(testdata, getSensSpec(obs,pred,dir="<")) auc <- with(testdata,auc_wtd(x=pred,y=obs)) if(plot){ par0 <- par(no.readonly=T) par(mar=c(5,4,4,1)+.1-c(2,1,3.5,0.5)) par(mgp=c(3,1,0)-c(1.3,.6,0)) par(tcl=-.3) plot(1,1,xlim=0:1,ylim=0:1,type="n",xlab="1 - Specificity",ylab="Sensitivity",main="",font.lab=2) with(roc, lines(1-spec, sens)) abline(0,1,lty=3) return(list(AUC=auc,ROC=roc,par=par(),par0=par0,fit=predict)) } else return(list(AUC=auc,ROC=roc,fit=predict)) } # cross validation crossValidate <- function(data,nFolds=10,nRepeats=20,develop, evaluate,response="HC_BC_HR",trace=T,group=1:nrow(data),...){ CV_output <- vector("list",nRepeats) names(CV_output) <- paste0("Iteration",seq_len(nRepeats)) id <- match(group,unique(group)) n <- max(id) for(iteration in seq_len(nRepeats)){ folds <- rep(seq_len(nFolds),length=n) folds <- sample(folds,n) CV_output[[iteration]] <- vector("list",nFolds) names(CV_output[[iteration]]) <- paste0("Fold",seq_len(nFolds)) for(fold in seq_len(nFolds)){ test_id <- which(folds==fold) train_id <- which(folds!=fold) test <- data[id%in%test_id,,drop=F] train <- data[id%in%train_id,,drop=F] f <- develop(train,response=response,...) CV_output[[iteration]][[fold]] <- evaluate(test,predict=f,response=response) if(trace){ cat("\rIteration ",iteration,"; fold ",fold," ",sep="") flush.console() } } } invisible(CV_output) } plotContributionsPC <- function(x=attr(fit1,"contributions")*100,which=1:4){ y <- x[,which,drop=F] convert <- function(x,source_range,target_range){ source_range <- range(source_range) target_range <- range(target_range) min(target_range)+(x-min(source_range))*diff(target_range)/diff(source_range) } barplot <- function(i){ z <- y[,i] xlower=1:m+(1:m-1)*ysep; xupper=xlower-1 xrange=c(0,100); yrange=c(0,m+(m-1)*ysep) frame <- data.frame(x=c(0,1)+(i-1)*(1+xsep),y=c(0,m+(m-1)*ysep)) bars <- lapply(seq_along(z), function(j)data.frame(x=c(xlower[j],xupper[j])[c(1,1,2,2,1)],y=c(0,z[j])[c(1,2,2,1,1)],w=j)) bars <- do.call(rbind,bars) bars <- data.frame( x=convert(bars$y,xrange,frame$x), y=convert(bars$x,yrange,frame$y), w=bars$w ) bars <- split(bars,bars$w) for(j in seq_along(bars)) with(bars[[j]],polygon(x,y,col="white")) lbl <- paste0(formatC(z,format="f",digits=1),"%") a <- z+10>=60 txt <- data.frame(x=convert(1:m+(1:m-1)*ysep-.5,yrange,frame$y), y=convert(ifelse(a,z-5,z+5),xrange,frame$x)) txt$adj <- ifelse(a,1,0) txt$labels <- lbl for(j in 1:nrow(txt)) with(as.data.frame(txt[j,]),text(y=x,x=y,labels=labels,adj=c(adj,0.5),cex=.7)) } par0 <- par(no.readonly=T) par(mar=c(5,4,4,1)+.1-c(2,-1,3.5,0.5)) par(mgp=c(3,1,0)-c(1.3,.6,0)) par(tcl=-.3) n <- length(which) m <- nrow(y) xsep <- .01*n ysep <- .01*m plot(1,1,xlim=c(0,n+(n-1)*xsep),ylim=c(0,m+(m-1)*ysep),type="n",xlab="",ylab="",main="",font.lab=2,axes=F) axis(2,at=1:m+(1:m-1)*ysep-.5,labels=F,lwd=0,lwd.ticks=1) mtext(rownames(y),side=2,at=1:m+(1:m-1)*ysep-.5,srt=90,las=1,adj=1,line=.5) for(i in 1:n){ axis(1,at=c(0,1)+(i-1)*(1+xsep),lwd=0,lwd.ticks=1,labels=F) mtext(0,at=0+(i-1)*(1+xsep),adj=0,line=.5,side=1) mtext(100,at=1+(i-1)*(1+xsep),adj=1,line=.5,side=1) mtext(paste0("PC",i),at=.5+(i-1)*(1+xsep),adj=0.5,line=1.5,side=1,font=2) polygon(x=(c(0,1)+(i-1)*(1+xsep))[c(1,2,2,1,1)],y=par()$usr[3:4][c(1,1,2,2,1)],col="lightgrey",border="transparent") } box() for(i in 1:n) barplot(i) return(NULL) } ### Models ### --------------------------------------------------------- # Fit 1 miRNA <- gsub("_T$","",grep("DCT_miR_(.*?)T",colnames(data),value=T)) miRNA <- miRNA[!grepl("miR_(126|99b|151a|18a)",miRNA)] model1 <- as.formula(paste0("HC_BC_HR~",paste(miRNA,collapse="+"),"+Age")) fit1 <- develop(data,formula=model1,boxCox=T,noiseReduction=T) attr(fit1,"scaled_coef") #coefficients: # (Intercept) 0.69 # DCT_miR_125a_5p 0.39 # DCT_miR_145_5p 0.48 # DCT_miR_148a_3p -0.36 # DCT_miR_153_3p -0.32 # DCT_miR_155_5p -0.07 # DCT_miR_16_5p 0.19 # DCT_miR_181a_5p 0.60 # DCT_miR_19a_3p 0.71 # DCT_miR_205_5p -0.56 # DCT_miR_21_5p -0.11 # DCT_miR_221_3p 3.00 # DCT_miR_222_3p -2.25 # DCT_miR_29c_5p -0.78 # DCT_miR_30b_5p 0.57 # DCT_miR_320a_3p -0.62 # DCT_miR_339_5p -0.64 # DCT_miR_374b_5p -0.61 # DCT_miR_425_5p 0.48 # DCT_miR_92a_3p -0.42 # Age 0.03 contr <- attr(fit1,"contributions")*100 rownames(contr) <- gsub("_","-",gsub("^(.*?)miR_","",rownames(contr))) apparent1 <- evaluate(testdata=data,predict=fit1,plot=T) apparent1$AUC cv1 <- crossValidate(data,evaluate=evaluate, develop=develop,formula=model1,boxCox=T,noiseReduction=T) mean(sapply(unlist(cv1,recursive=F),function(x)x$AUC)) # mean AUC across folds and CV iterations (cv_auc <- mean(sapply(unlist(cv1,recursive=F),function(x)x$AUC))) # mean AUC across folds and CV iterations (cv_auc2 <- median(sapply(unlist(cv1,recursive=F),function(x)x$AUC))) # median AUC across folds and CV iterations (cv_auc3 <- quantile(sapply(unlist(cv1,recursive=F),function(x)x$AUC),probs=.025)) # 2.5th percentile (cv_auc4 <- quantile(sapply(unlist(cv1,recursive=F),function(x)x$AUC),probs=.975)) # 97.5th percentile ROCs <- lapply(unlist(cv1,recursive=F),function(x)x$ROC) poolROCs <- function(ROCs,FUN=mean,...){ x <- seq(0,1,length=1001) y <- do.call(rbind,lapply(ROCs,function(z)interpolate(z$spec,z$sens,xout=x))) data.frame(spec=c(0,x,1),sens=c(1,unname(apply(y,2,FUN,...,simplify=T)),0)) } with(poolROCs(ROCs,median),lines(1-spec,sens,col="black",lty=2)) #Calculate the 95% CI for the apparent AUC calculate_auc_ci <- function(auc, data, response, confidence_level = 0.95) { n1 <- sum(data[[response]] == 1) n2 <- sum(data[[response]] == 0) Q1 <- auc / (2 - auc) Q2 <- 2 * auc^2 / (1 + auc) SE <- sqrt((auc * (1 - auc) + (n1 - 1) * (Q1 - auc^2) + (n2 - 1) * (Q2 - auc^2)) / (n1 * n2)) Z <- qnorm(1 - (1 - confidence_level) / 2) lower_ci <- auc - Z * SE upper_ci <- auc + Z * SE return(c(lower_ci, upper_ci)) } apparent_auc <- apparent1$AUC ci <- calculate_auc_ci(apparent_auc, data, response = "HC_BC_HR") ci_lower <- formatC(ci[1], format = "f", digits = 3) ci_upper <- formatC(ci[2], format = "f", digits = 3) print(ci) #legend lgd <- c( paste0("Apparent performance (AUC = ", formatC(apparent_auc, format = "f", digits = 3), ", 95% CI [", ci_lower, ", ", ci_upper, "])"), paste0("Median across CV folds and iterations (AUC = ", formatC(cv_auc2, format = "f", digits = 3), ", 2.5% =", formatC(cv_auc3, format = "f", digits = 3), "; 97.5% = ", formatC(cv_auc4, format = "f", digits = 3), ")") ) legend("bottomright",lty=c(1,2),col=c("black","black"), lgd,cex=1,bty="n") # Explore model across folds/repeats tmp <- lapply(unlist(cv1,recursive=F),function(x)attr(x$fit,"coef")[-1,"s1"]) vars <- unique(unlist(lapply(tmp,names))) tmp <- lapply(tmp,function(x){ out <- x[match(vars,names(x))] names(out) <- vars out }) tmp <- do.call(rbind,tmp) tmp[is.na(tmp)] <- 0 inclusion_fraction <- colMeans(tmp != 0) a <- inclusion_fraction * 100 sort(a) # Fit 2 # Without age miRNA <- gsub("_T$","",grep("DCT_miR_(.*?)T",colnames(data),value=T)) miRNA <- miRNA[!grepl("miR_(126|99b|151a|18a)",miRNA)] model1 <- as.formula(paste0("HC_BC_HR~",paste(miRNA,collapse="+"))) fit1 <- develop(data,formula=model1,boxCox=T,noiseReduction=T) attr(fit1,"scaled_coef") contr <- attr(fit1,"contributions")*100 rownames(contr) <- gsub("_","-",gsub("^(.*?)miR_","",rownames(contr))) apparent1 <- evaluate(testdata=data,predict=fit1,plot=T) apparent1$AUC cv1 <- crossValidate(data,evaluate=evaluate, develop=develop,formula=model1,boxCox=T,noiseReduction=T) mean(sapply(unlist(cv1,recursive=F),function(x)x$AUC)) # mean AUC across folds and CV iterations (cv_auc <- mean(sapply(unlist(cv1,recursive=F),function(x)x$AUC))) # mean AUC across folds and CV iterations (cv_auc2 <- median(sapply(unlist(cv1,recursive=F),function(x)x$AUC))) # median AUC across folds and CV iterations (cv_auc3 <- quantile(sapply(unlist(cv1,recursive=F),function(x)x$AUC),probs=.025)) # 2.5th percentile (cv_auc4 <- quantile(sapply(unlist(cv1,recursive=F),function(x)x$AUC),probs=.975)) # 97.5th percentile ROCs <- lapply(unlist(cv1,recursive=F),function(x)x$ROC) poolROCs <- function(ROCs,FUN=mean,...){ x <- seq(0,1,length=1001) y <- do.call(rbind,lapply(ROCs,function(z)interpolate(z$spec,z$sens,xout=x))) data.frame(spec=c(0,x,1),sens=c(1,unname(apply(y,2,FUN,...,simplify=T)),0)) } with(poolROCs(ROCs,median),lines(1-spec,sens,col="black",lty=2)) #Calculate the 95% CI for the apparent AUC calculate_auc_ci <- function(auc, data, response, confidence_level = 0.95) { n1 <- sum(data[[response]] == 1) n2 <- sum(data[[response]] == 0) Q1 <- auc / (2 - auc) Q2 <- 2 * auc^2 / (1 + auc) SE <- sqrt((auc * (1 - auc) + (n1 - 1) * (Q1 - auc^2) + (n2 - 1) * (Q2 - auc^2)) / (n1 * n2)) Z <- qnorm(1 - (1 - confidence_level) / 2) lower_ci <- auc - Z * SE upper_ci <- auc + Z * SE return(c(lower_ci, upper_ci)) } apparent_auc <- apparent1$AUC ci <- calculate_auc_ci(apparent_auc, data, response = "HC_BC_HR") ci_lower <- formatC(ci[1], format = "f", digits = 3) ci_upper <- formatC(ci[2], format = "f", digits = 3) print(ci) #legend lgd <- c( paste0("Apparent performance (AUC = ", formatC(apparent_auc, format = "f", digits = 3), ", 95% CI [", ci_lower, ", ", ci_upper, "])"), paste0("Median across CV folds and iterations (AUC = ", formatC(cv_auc2, format = "f", digits = 3), ", 2.5% =", formatC(cv_auc3, format = "f", digits = 3), "; 97.5% = ", formatC(cv_auc4, format = "f", digits = 3), ")") ) legend("bottomright",lty=c(1,2),col=c("black","black"), lgd,cex=1,bty="n") # Explore model across folds/repeats tmp <- lapply(unlist(cv1,recursive=F),function(x)attr(x$fit,"coef")[-1,"s1"]) vars <- unique(unlist(lapply(tmp,names))) tmp <- lapply(tmp,function(x){ out <- x[match(vars,names(x))] names(out) <- vars out }) tmp <- do.call(rbind,tmp) tmp[is.na(tmp)] <- 0 inclusion_fraction <- colMeans(tmp != 0) a <- inclusion_fraction * 100 sort(a) # Fit 3 #PCA on 80% of variance PCs model3 <- as.formula(paste0("HC_BC_HR~Age")) fit5 <- develop(data,formula=model3,boxCox=T,noiseReduction=T,PCAoption=4) contr <- attr(fit5,"contributions")*100 rownames(contr) <- gsub("_","-",gsub("^(.*?)miR_","",rownames(contr))) plotContributionsPC(contr,which=1:5) apparent5 <- evaluate(data,predict=fit5,plot=T) apparent5$AUC cv5 <- crossValidate(data,develop=develop,formula=model3,boxCox=T,noiseReduction=T,PCAoption=4,evaluate=evaluate) (cv_auc <- mean(sapply(unlist(cv5,recursive=F),function(x)x$AUC))) # mean AUC across folds and CV iterations (cv_auc2 <- median(sapply(unlist(cv5,recursive=F),function(x)x$AUC))) # median AUC across folds and CV iterations (cv_auc3 <- quantile(sapply(unlist(cv5,recursive=F),function(x)x$AUC),probs=.025)) # 25th percentile (cv_auc4 <- quantile(sapply(unlist(cv5,recursive=F),function(x)x$AUC),probs=.975)) # 75th percentile ROCs <- lapply(unlist(cv5,recursive=F),function(x)x$ROC) poolROCs <- function(ROCs,FUN=mean,...){ x <- seq(0,1,length=1001) y <- do.call(rbind,lapply(ROCs,function(z)interpolate(z$spec,z$sens,xout=x))) data.frame(spec=c(0,x,1),sens=c(1,unname(apply(y,2,FUN,...,simplify=T)),0)) } with(poolROCs(ROCs,median),lines(1-spec,sens,col="black",lty=2)) #Calculate the 95% CI for the apparent AUC calculate_auc_ci <- function(auc, data, response, confidence_level = 0.95) { n1 <- sum(data[[response]] == 1) n2 <- sum(data[[response]] == 0) Q1 <- auc / (2 - auc) Q2 <- 2 * auc^2 / (1 + auc) SE <- sqrt((auc * (1 - auc) + (n1 - 1) * (Q1 - auc^2) + (n2 - 1) * (Q2 - auc^2)) / (n1 * n2)) Z <- qnorm(1 - (1 - confidence_level) / 2) lower_ci <- auc - Z * SE upper_ci <- auc + Z * SE return(c(lower_ci, upper_ci)) } apparent_auc <- apparent5$AUC ci <- calculate_auc_ci(apparent_auc, data, response = "HC_BC_HR") ci_lower <- formatC(ci[1], format = "f", digits = 3) ci_upper <- formatC(ci[2], format = "f", digits = 3) print(ci) #legend lgd <- c( paste0("Apparent performance (AUC = ", formatC(apparent_auc, format = "f", digits = 3), ", 95% CI [", ci_lower, ", ", ci_upper, "])"), paste0("Median across CV folds and iterations (AUC = ", formatC(cv_auc2, format = "f", digits = 3), ", 2.5% =", formatC(cv_auc3, format = "f", digits = 3), "; 97.5% = ", formatC(cv_auc4, format = "f", digits = 3), ")") ) legend("bottomright",lty=c(1,2),col=c("black","black"), lgd,cex=1,bty="n") contr <- lapply(unlist(cv5,recursive=F),function(x)attr(x$fit,"contributions")) for(i in seq_len(length(contr)-1)) contr[[i+1]] <- contr[[i]]+contr[[i+1]] contr <- contr[[length(contr)]]/length(contr) rownames(contr) <- gsub("_","-",gsub("^(.*?)miR_","",rownames(contr))) plotContributionsPC(contr*100,which=1:5) # average across CV folds/repeats # Explore model across folds/repeats tmp <- lapply(unlist(cv5,recursive=F),function(x)attr(x$fit,"coef")[-1,"s1"]) vars <- unique(unlist(lapply(tmp,names))) tmp <- lapply(tmp,function(x){ out <- x[match(vars,names(x))] names(out) <- vars out }) tmp <- do.call(rbind,tmp) tmp[is.na(tmp)] <- 0 inclusion_fraction <- colMeans(tmp != 0) a <- inclusion_fraction * 100 sort(a) # ------------------------------------------------------------------------------ # 4. BC vs CL within BC patients dataBC <- data[data[[outcome]],,drop=F] dataBC <- dataBC[dataBC$BC_left_or_right%in%c("Left","Right"),,drop=F] dataBC$cluster <- dataBC$Subject_nr n <- nrow(dataBC) dataBC <- rbind(dataBC,dataBC) dataBC$BC <- c(rep(1,n),rep(0,n)) tmp <- gsub("_T$","",colsDCT_T) colsDCT <- tmp[!grepl("miR_(126|99b|151a|18a)",tmp)] colsDCT_T <- paste0(colsDCT,"_T") colsDCT_CL <- paste0(colsDCT,"_CL") for(i in seq_along(colsDCT)){ dataBC[[colsDCT[i]]] <- ifelse(dataBC$BC,dataBC[[colsDCT_T[i]]],dataBC[[colsDCT_CL[i]]]) } isNA <- is.na(dataBC[colsDCT]) colSums(isNA) (excl <- unique(dataBC$Subject_nr[rowSums(isNA)>1])) dataBC <- dataBC[!dataBC$Subject_nr%in%excl,,drop=F] dataBC$Cloudiness <- ifelse(dataBC$BC,dataBC$TandHC_cloudiness,dataBC$CL_cloudiness) dataBC$RNA_input <- ifelse(dataBC$BC,dataBC$TandHC_RNA_input,dataBC$CL_RNA_input) develop2 <- function( data, response="BC", formula=model1, boxCox=T, noiseReduction=T, PCAoption=2, # PCAoption==1: not used for analysis # PCAoption==2: include all, after noise red. and box-cox LASSO=T ){ # -------------------------------------------- # 3.2 Correcting miRNA expression levels for different variables using linear regression # -------------------------------------------- depvar <- colsDCT df <- data[depvar] colnames(df) <- depvar df$Subject_nr <- data$Subject_nr Cloudiness <- data$Cloudiness RNA_input <- data$RNA_input noise_reduction_model_fits <- vector("list",length(depvar)) names(noise_reduction_model_fits) <- depvar if(noiseReduction){ for(var in depvar){ df$Cloudiness <- Cloudiness df$RNA_input <- RNA_input noise_reduction_model_fits[[var]] <- if(!noiseReduction) lm(df[[var]]~1,data=df) else lm(df[[var]]~as.factor(Cloudiness)+as.factor(RNA_input),data=df) data[[var]] <- if(!noiseReduction) data[[var]] else resid(noise_reduction_model_fits[[var]]) } } colsRES_DCT <- colsDCT # -------------------------------------------- # 3.3 Box-Cox # -------------------------------------------- boxcox_trans <- vector("list",length(colsRES_DCT)) names(boxcox_trans) <- colsRES_DCT for(cl in colsRES_DCT){ if(!boxCox){ boxcox_trans[[cl]]$m <- -Inf boxcox_trans[[cl]]$k <- 1 boxcox_trans[[cl]]$lambda <- 1 next } x <- data[[cl]] if(all(is.na(x))){ boxcox_trans[[cl]]$m <- 0 boxcox_trans[[cl]]$k <- 0 } else if(any(x<=0,na.rm=T)){ m <- min(abs(x),na.rm=T)/2 k <- abs(min(x,na.rm=T))+m if(k==0) k <- 1 boxcox_trans[[cl]]$k <- k boxcox_trans[[cl]]$m <- m } else{ boxcox_trans[[cl]]$k <- 0 boxcox_trans[[cl]]$m <- 0 } lambda <- if(all(is.na(x))) NA else attr(boxcox(x+boxcox_trans[[cl]]$k),"lambda") boxcox_trans[[cl]]$lambda <- lambda } if(boxCox){ for(cl in colsRES_DCT){ data[[cl]] <- as.numeric(boxcox(data[[cl]]+boxcox_trans[[cl]]$k,lambda=boxcox_trans[[cl]]$lambda)) } } colsRES_BC_DCT <- colsRES_DCT # -------------------------------------------- # 3.4 Principal component dimensionality reduction # -------------------------------------------- selectPCA <- function(pca,thresh=0.8,returnWhichOnly=F){ vv <- diag(cov(pca)) tt <- sum(vv) n <- which(cumsum(vv/tt)>=thresh)[1L] if(returnWhichOnly) return(seq_len(n)) out <- pca[,seq_len(n),drop=FALSE] attr(out,"components") <- seq_len(n) attr(out,"cum_var_explained") <- cumsum(vv/tt) return(out) } if(PCAoption==1){ x <- data[colsRES_BC_DCT[!grepl("miR_(125|145|21|29c)",colsRES_BC_DCT)]] pca <- prcomp(x,scale.=T) pca_orig <- pca$x data <- cbind(data,pca_orig) rot <- pca$rotation contribution <- sweep(abs(rot),2,colSums(abs(rot)),"/") rownames(contribution) <- colnames(x) contribution } else if(PCAoption==2){ x <- data[colsRES_BC_DCT] pca <- prcomp(x,scale.=T) pca_orig <- pca$x data <- cbind(data,pca_orig) rot <- pca$rotation contribution <- sweep(abs(rot),2,colSums(abs(rot)),"/") rownames(contribution) <- colnames(x) contribution } else contribution <- NULL # -------------------------------------------- # MLE with/without LASSO # -------------------------------------------- pclogit <- function(lambda,dat=data){ mm <- model.matrix(formula,data=dat)[,-1,drop=F] out <- penalized::penalized(Surv(rep(1,nrow(dat)),dat[[response]])~strata(cluster), penalized=mm,lambda1=lambda,model="cox",standardize=TRUE,data=dat,trace=F) return(out) } if(LASSO){ lambdas <- exp(seq(2,-10,length=100)) mm <- model.matrix(formula,data=data)[,-1,drop=F] findLambda <- function(dev=pclogit,nFolds=10){ n <- nrow(data) CV_output <- vector("list",nFolds) names(CV_output) <- paste0("Fold",seq_len(nFolds)) id <- match(data$cluster,unique(data$cluster)) n <- max(id) folds <- rep(seq_len(nFolds),length=n) folds <- sample(folds,n) for(fold in seq_len(nFolds)){ test_id <- which(folds==fold) train_id <- which(folds!=fold) test <- data[id%in%test_id,,drop=F] train <- data[id%in%train_id,,drop=F] f <- lapply(lambdas,function(lambda) dev(lambda,train)) x <- lapply(f,function(g) (mm[test_id,,drop=F]%*%coef(g,which="all"))[,]) auc <- sapply(seq_along(x),function(i)with(test,auc_wtd(x=x[[i]],y=test[[response]]))) lambda <- lambdas[which.max(auc)] CV_output[[fold]] <- lambda } lambda <- mean(unlist(CV_output)) return(lambda) } lambda <- findLambda(pclogit) fit <- pclogit(lambda) b <- coef(fit,which="all") } else{ fit <- pclogit(0) mm <- model.matrix(formula,data=data)[,-1,drop=F] b <- coef(fit,which="all") } predict <- function(newdata){ # Noise reduction depvar <- colsDCT Cloudiness <- with(newdata,Cloudiness) RNA_input <- with(newdata,RNA_input) for(var in depvar){ df <- data.frame(id=newdata$Subject_nr) df$Cloudiness <- Cloudiness df$RNA_input <- RNA_input p <- if(!noiseReduction) 0 else stats::predict(noise_reduction_model_fits[[var]],newdata=df,type="response") o <- newdata[[var]] newdata[[var]] <- o-p } # Box-Cox for(cl in colsRES_DCT){ newdata[[cl]] <- boxcox(pmax(boxcox_trans[[cl]]$m,newdata[[cl]]+boxcox_trans[[cl]]$k),lambda=boxcox_trans[[cl]]$lambda) } # PCA pca_orig <- stats:::predict.prcomp(pca,newdata=newdata) newdata <- cbind(newdata,pca_orig) current.na.action <- options('na.action') options(na.action='na.pass') x <- model.matrix(formula,data=newdata)[,-1,drop=F] options(na.action=current.na.action$na.action) out <- plogis(as.vector(x%*%b)) attr(out,"model.matrix") <- x return(out) } attr(predict,"scalingParamPCA") <- list(mean=pca$center,sd=pca$scale) attr(predict,"coefPCA") <- pca$rotation attr(predict,"contributions") <- contribution attr(predict,"coef") <- b attr(predict,"scaled_coef") <- b*apply(mm,2,sd) return(predict) } # Model A miRNA <- paste(colsDCT,collapse="+") modelA <- as.formula(paste0("BC~",miRNA)) fitA <- develop2(dataBC,formula=modelA) apparentA <- evaluate(dataBC,fitA,response="BC",plot=T) cvA <- crossValidate(dataBC,develop=develop2,response="BC",formula=modelA,boxCox=T,noiseReduction=T,LASSO=T,evaluate=evaluate,group=dataBC$cluster) (cv_auc <- mean(sapply(unlist(cvA,recursive=F),function(x)x$AUC))) # mean AUC across folds and CV iterations (cv_auc2 <- median(sapply(unlist(cvA,recursive=F),function(x)x$AUC))) # median AUC across folds and CV iterations (cv_auc3 <- quantile(sapply(unlist(cvA,recursive=F),function(x)x$AUC),probs=.025)) # 2.5th percentile (cv_auc4 <- quantile(sapply(unlist(cvA,recursive=F),function(x)x$AUC),probs=.975)) # 97.5th percentile ROCs <- lapply(unlist(cvA,recursive=F),function(x)x$ROC) poolROCs <- function(ROCs,FUN=mean,...){ x <- seq(0,1,length=1001) y <- do.call(rbind,lapply(ROCs,function(z)interpolate(z$spec,z$sens,xout=x))) data.frame(spec=c(0,x,1),sens=c(1,unname(apply(y,2,FUN,...,simplify=T)),0)) } with(poolROCs(ROCs,median),lines(1-spec,sens,col="black",lty=2)) #Calculate the 95% CI for the apparent AUC calculate_auc_ci <- function(auc, data, response, confidence_level = 0.95) { n1 <- sum(data[[response]] == 1) n2 <- sum(data[[response]] == 0) Q1 <- auc / (2 - auc) Q2 <- 2 * auc^2 / (1 + auc) SE <- sqrt((auc * (1 - auc) + (n1 - 1) * (Q1 - auc^2) + (n2 - 1) * (Q2 - auc^2)) / (n1 * n2)) Z <- qnorm(1 - (1 - confidence_level) / 2) lower_ci <- auc - Z * SE upper_ci <- auc + Z * SE return(c(lower_ci, upper_ci)) } apparent_auc <- apparentA$AUC ci <- calculate_auc_ci(apparent_auc, data, response = "HC_BC_HR") ci_lower <- formatC(ci[1], format = "f", digits = 3) ci_upper <- formatC(ci[2], format = "f", digits = 3) print(ci) #legend lgd <- c( paste0("Apparent performance (AUC = ", formatC(apparent_auc, format = "f", digits = 3), ", 95% CI [", ci_lower, ", ", ci_upper, "])"), paste0("Median across CV folds and iterations (AUC = ", formatC(cv_auc2, format = "f", digits = 3), ", 2.5% =", formatC(cv_auc3, format = "f", digits = 3), "; 97.5% = ", formatC(cv_auc4, format = "f", digits = 3), ")") ) legend("bottomright",lty=c(1,2),col=c("black","black"), lgd,cex=1,bty="n") contr <- attr(fitA,"contributions")*100 # Explore model across folds/repeats --> does not work... tmp <- lapply(unlist(cvA,recursive=F),function(x)attr(x$fitA,"coef")[-1,"s1"]) vars <- unique(unlist(lapply(tmp,names))) tmp <- lapply(tmp,function(x){ out <- x[match(vars,names(x))] names(out) <- vars out }) tmp <- do.call(rbind,tmp) tmp (1-colMeans(is.na(tmp))) tmp[is.na(tmp)] <- 0 colMeans(tmp!=0)*100 contr <- lapply(unlist(cvA,recursive=F),function(x)attr(x$fitA,"contributions")) for(i in seq_len(length(contr)-1)) contr[[i+1]] <- contr[[i]]+contr[[i+1]] contr <- contr[[length(contr)]]/length(contr) rownames(contr) <- gsub("_","-",gsub("^(.*?)miR_","",rownames(contr)))