十一、timeROC

作者: 白米饭睡不醒 | 来源:发表于2022-01-26 09:51 被阅读0次

    1.准备输入数据

    load("TCGA-KIRC_sur_model.Rdata")
    

    2.构建lasso回归模型

    输入数据是表达矩阵(仅含tumor样本)和对应的生死。

    x=t(exprSet)
    y=meta$event
    library(glmnet)
    cv_fit <- cv.glmnet(x=x, y=y, nlambda = 1000,alpha = 1)
    model_lasso_min <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.min)
    choose_gene_mi n=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]
    length(choose_gene_min)
    #> [1] 39
    

    3.模型预测和评估

    3.1自己预测自己

    lasso.prob <- predict(cv_fit, newx=x , s=cv_fit$lambda.min )
    re=cbind(y ,lasso.prob)
    head(re)
    #>                              y         1
    #> TCGA-A3-3307-01A-01T-0860-13 0 0.1107463
    #> TCGA-A3-3308-01A-02R-1324-13 0 0.3985909
    #> TCGA-A3-3311-01A-02R-1324-13 1 0.2875707
    #> TCGA-A3-3313-01A-02R-1324-13 1 0.3884061
    #> TCGA-A3-3316-01A-01T-0860-13 0 0.3288315
    #> TCGA-A3-3317-01A-01T-0860-13 0 0.3793098
    

    3.2 time-ROC

    new_dat=meta
    library(timeROC)
    library(survival)
    library(survminer)
    new_dat$fp=as.numeric(lasso.prob[,1])
    with(new_dat,
         ROC <<- timeROC(T=time,#结局时间 
                         delta=event,#生存结局 
                         marker=fp,#预测变量 
                         cause=1,#阳性结局赋值,比如死亡与否
                         weighting="marginal",#权重计算方法,marginal是默认值,采用km计算删失分布
                         times=c(60,100),#时间点,选取5年(60个月)和8年生存率,小的年份写前面
                         ROC = TRUE,
                         iid = TRUE)
    )
    auc_60 = ROC$AUC[[1]]
    auc_100 = ROC$AUC[[2]]
    dat = data.frame(tpr_60 = ROC$TP[,1],
                     fpr_60 = ROC$FP[,1],
                     tpr_100 = ROC$TP[,2],
                     fpr_100 = ROC$FP[,2])
    library(ggplot2)
    
    ggplot() + 
      geom_line(data = dat,aes(x = fpr_60, y = tpr_60),color = "blue") + 
      geom_line(data = dat,aes(x = fpr_100, y = tpr_100),color = "red")+
      geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
      theme_bw()+
      annotate("text",x = .75, y = .25,
               label = paste("AUC of 60 = ",round(auc_60,2)),color = "blue")+
      annotate("text",x = .75, y = .15,label = paste("AUC of 100 = ",round(auc_100,2)),color = "red")+
      scale_x_continuous(name  = "fpr")+
      scale_y_continuous(name = "tpr")
    

    4.切割数据构建模型并预测

    4.1 切割数据

    用R包caret切割数据,生成的结果是一组代表列数的数字,用这些数字来给表达矩阵和meta取子集即可。

    library(caret)
    set.seed(12345679)
    sam<- createDataPartition(meta$event, p = .5,list = FALSE)
    train <- exprSet[,sam]
    test <- exprSet[,-sam]
    train_meta <- meta[sam,]
    test_meta <- meta[-sam,]
    

    4.2 切割后的train数据集建模

    和上面的建模方法一样。

    #计算lambda
    x = t(train)
    y = train_meta$event
    cv_fit <- cv.glmnet(x=x, y=y, nlambda = 1000,alpha = 1)
    #构建模型
    model_lasso_min <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.min)
    #挑出基因
    choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]
    length(choose_gene_min)
    #> [1] 18
    

    4.3 模型预测

    用训练集构建模型,预测测试集的生死,注意newx参数变了。

    lasso.prob <- predict(cv_fit, newx=t(test), s=cv_fit$lambda.min)
    re=cbind(test_meta$event ,lasso.prob)
    

    4.4 time-ROC

    new_dat = test_meta
    library(timeROC)
    library(survival)
    library(survminer)
    new_dat$fp=as.numeric(lasso.prob[,1])
    with(new_dat,
         ROC <<- timeROC(T=time,#结局时间 
                         delta=event,#生存结局 
                         marker=fp,#预测变量 
                         cause=1,#阳性结局赋值,比如死亡与否
                         weighting="marginal",#权重计算方法,marginal是默认值,采用km计算删失分布
                         times=c(60,100),#时间点,选取5年(60个月)和8年生存率
                         ROC = TRUE,
                         iid = TRUE)
    )
    auc_60 = ROC$AUC[[1]]
    auc_100 = ROC$AUC[[2]]
    dat = data.frame(tpr_60 = ROC$TP[,1],
                     fpr_60 = ROC$FP[,1],
                     tpr_100 = ROC$TP[,2],
                     fpr_100 = ROC$FP[,2])
    library(ggplot2)
    
    ggplot() + 
      geom_line(data = dat,aes(x = fpr_60, y = tpr_60),color = "blue") + 
      geom_line(data = dat,aes(x = fpr_100, y = tpr_100),color = "red")+
      geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
      theme_bw()+
      annotate("text",x = .75, y = .25,
               label = paste("AUC of 60 = ",round(auc_60,2)),color = "blue")+
      annotate("text",x = .75, y = .15,label = paste("AUC of 100 = ",round(auc_100,2)),color = "red")+
      scale_x_continuous(name  = "fpr")+
      scale_y_continuous(name = "tpr")
    

    *生信技能树课程笔记

    相关文章

      网友评论

        本文标题:十一、timeROC

        本文链接:https://www.haomeiwen.com/subject/ucgwhrtx.html