【问题标题】:R: Efficient formula creation within `lm`R:在“lm”中高效创建公式
【发布时间】:2023-03-06 08:20:01
【问题描述】:

我希望为我正在对遗传数据进行的一些迭代建模编写通用代码。

这是我的数据框的一个子集:

> head(exprTarget)
      patient CDR Diagnosis diag_key  UNC93B1     CTSC     PLEK   LGALS9      GRN    CYTH4     C1QA     C1QC     C1QB   LAPTM5     CTSS    FCER1G   ALOX5AP
16955   16955   2       MCI        1 2.468387 3.306170 1.669025 2.197085 4.817537 2.303606 3.126281 3.537686 4.077572 4.660030 2.960342 1.0880424 2.0820685
16365   16365   5        AD        2 2.312767 3.205852 1.276787 1.942052 4.924718 2.461212 2.641784 3.592875 3.758567 4.215387 2.536174 0.9872809 0.7559553
17155   17155   5        AD        2 3.276758 4.039103 2.482880 3.347225 5.465345 2.990894 6.004585 6.108294 6.762214 5.708623 4.358901 2.5924355 3.6172763
17135   17135   5        AD        2 2.245509 3.056953 1.877469 2.083920 4.492934 1.827284 2.584534 3.012729 3.369049 3.892801 2.990098 0.7350252 1.1568519
16625   16625   4        AD        2 2.575806 3.978674 2.060418 2.327522 4.981906 2.685569 4.694788 4.725954 5.460863 5.260811 4.021172 2.5871655 3.3241311
16295   16295   4        AD        2 3.107424 3.701104 2.880653 2.880653 5.115831 2.723281 4.224342 4.717155 5.110232 5.031450 3.980189 2.0809520 1.9699207

我正在尝试使用diag_key 作为我的响应变量,其右侧的所有列都充当预测变量,即所需的公式是:

lm(diag_key ~ . - patient - CDR - Diagnosis, data= exprTarget)

我想让它更通用。具体来说,我希望能够只传递将我的临床注释与基因表达数据分开的列号,在上面的示例中,这将是第 4 列 diag_key,尽管对于不同的实现它可能会有所不同。

我目前的目标是仅使用此信息重新创建上述公式。这是我目前的尝试,注意response对应于分隔列号,即上例中的4:

 clinical<- colnames(exprTarget)[1:(response-1)]

 lm(as.formula(paste("exprTarget[,response] ~ . ",clinical, sep = "-")), data= exprTarget)

此命令生成公式exprTarget[, 4] ~ . - patient,所以很明显问题是在使用paste 将公式的开头与我要删除的列连接起来时,只粘贴了列表中的第一项。

非常感谢您在导航方面的任何帮助。

【问题讨论】:

    标签: r automation paste lm


    【解决方案1】:

    可以这样避免构造公式:

    nc <- ncol(exprTarget)
    fm <- lm(exprTarget[response:nc])
    

    但如果你真的想要这个公式:

    formula(fm)
    ## diag_key ~ UNC93B1 + CTSC + PLEK + LGALS9 + GRN + CYTH4 + C1QA + 
    ##  C1QC + C1QB + LAPTM5 + CTSS + FCER1G + ALOX5AP
    

    注意:上面使用的可重现形式的输入是:

    Lines <- "
          patient CDR Diagnosis diag_key  UNC93B1     CTSC     PLEK   LGALS9      GRN    CYTH4     C1QA     C1QC     C1QB   LAPTM5     CTSS    FCER1G   ALOX5AP
    16955   16955   2       MCI        1 2.468387 3.306170 1.669025 2.197085 4.817537 2.303606 3.126281 3.537686 4.077572 4.660030 2.960342 1.0880424 2.0820685
    16365   16365   5        AD        2 2.312767 3.205852 1.276787 1.942052 4.924718 2.461212 2.641784 3.592875 3.758567 4.215387 2.536174 0.9872809 0.7559553
    17155   17155   5        AD        2 3.276758 4.039103 2.482880 3.347225 5.465345 2.990894 6.004585 6.108294 6.762214 5.708623 4.358901 2.5924355 3.6172763
    17135   17135   5        AD        2 2.245509 3.056953 1.877469 2.083920 4.492934 1.827284 2.584534 3.012729 3.369049 3.892801 2.990098 0.7350252 1.1568519
    16625   16625   4        AD        2 2.575806 3.978674 2.060418 2.327522 4.981906 2.685569 4.694788 4.725954 5.460863 5.260811 4.021172 2.5871655 3.3241311
    16295   16295   4        AD        2 3.107424 3.701104 2.880653 2.880653 5.115831 2.723281 4.224342 4.717155 5.110232 5.031450 3.980189 2.0809520 1.9699207"
    exprTarget <- read.table(text = Lines)
    response <- 4
    

    【讨论】:

      【解决方案2】:

      您可以为此使用reformulate。我基于响应的列名,因为这似乎比直接使用列索引更安全,但您当然可以只使用列索引。

      resp_col = "diag_key"
      idx = match(resp_col, names(exprTarget))
      
      my_formula = reformulate(names(exprTarget)[(idx+1):ncol(exprTarget)], response=names(exprTarget)[idx])
      
      diag_key ~ UNC93B1 + CTSC + PLEK + LGALS9 + GRN + CYTH4 + C1QA + 
          C1QC + C1QB + LAPTM5 + CTSS + FCER1G + ALOX5AP
      

      你可以把它打包成一个函数:

      lm_form = function(data, resp_col) {
      
        idx = match(resp_col, names(data))
      
        form = reformulate(names(data)[(idx+1):ncol(data)], response=names(data)[idx])
      
        lm(form, data=data)
      }
      
      my_model = lm_form(exprTarget, "diag_key")
      

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 2020-10-03
        • 2013-06-27
        • 2012-12-10
        • 2020-03-04
        • 2015-06-25
        • 1970-01-01
        相关资源
        最近更新 更多