【问题标题】:R. Multivariate linear regression iterating column-wise over pairs of variablesR. 多元线性回归在变量对上逐列迭代
【发布时间】:2021-01-15 07:06:30
【问题描述】:

我有一个数据框 dfA(真正的有 1000 行和 400,000 列)。从第 6 列开始,变量名称是由 x 组成的“三元组”,带有 + 不同的前缀(GT_x、N_x、E_x),其中 x = rs1、rs7、rs300、rs502 等:

ID    SEX    PV    GAN    GAE    GT_rs1    N_rs1    E_rs1    GT_rs7    N_rs7    E_rs7    ...
2    0    7.8    0.3    0.4    0    1    1    1    0    2    ...
6    1    6.4    0.35    0.55    0    0    1    1    1    2    ...

这是我的数据的可重现示例:

dfA = data.frame(rbind(c("ID","SEX","PV","GAN","GAE","GT_rs1","N_rs1","E_rs1","GT_rs7","N_rs7","E_rs7"), 
                   c(2,0,7.8,0.3,0.4,0,1,1,1,0,2),
                   c(6,1,6.4,0.35,0.55,0,0,1,1,1,2)))
dfA = dfA %>% row_to_names(row_number = 1)

使用 R,我想运行形式的线性回归:

lm(PV ~ SEX + GAN + GT_x + N_x)

其中 x 是 rs1、rs7 等等。所以,我需要逐列迭代成对的变量。我想获得不同协变量(SEX、GAN、GT_x 和 N_x)的估计值、std.error、statistic 和 p.value。 SEX = 分类变量; PV, GAN = 定量变量; GT_x, N_x, E_x = 附加变量。

【问题讨论】:

  • 你能用dput(dfA)添加可重现的数据吗?拥有用于计算几个不同的lms(例如PV, SEX, GAN, GT_rs1, N_rs1, GT_rs7, N_rs7)的列以及足够的行来获得结果就足够了。
  • 模型中是否需要E_x
  • 嗨@GregorThomas。不,谢谢
  • 嗨@iago,我刚刚添加了可重现的数据。谢谢

标签: r linear-regression


【解决方案1】:

您可以通过将字符串粘贴在一起来构建公式 - 我们只需要知道您要粘贴在一起的字符串。

这应该有效 - 它未经测试,因为您共享的数据未与 dput 共享,因此不可复制/粘贴,并且它只有一组协变量,因此无法说明问题的复杂性。如果您有问题,请分享复制/粘贴数据来说明,我会尝试调试。

library(stringr)
library(dplyr)
library(broom)
# get all unique strings after underscores from your column names
suffix = str_extract(names(dfA), "_.*") %>% na.omit %>% unique
prefix = c("GT", "N")
base_formula = "PV ~ SEX + GAN +"
full_formula = paste(base_formula, paste0(prefix[1], suffix), "+", paste0(prefix[2], suffix))

mods = list()
for(i in seq_along(full_formula)) {
  mods[[suffix[i]]] = lm(as.formula(full_formula[i]), data = dfA)
}

stats = lapply(mods, tidy)
stats = bind_rows(stats, .id = "suffix")

【讨论】:

    【解决方案2】:

    这是一个在一个简单管道中使用purrr 的解决方案。

    您只需创建GT_xN_x 的列表即可使用。您可以使用一些正则表达式来做到这一点。

    library(purrr)
    
    nn <- names(df)
    pattern <- "^GT_|^N_"
    
    vars <- nn[grepl(pattern, nn)] # get the variables that start with GT_ and N_
    x <- sub(pattern, "", vars)    # get every x
    
    split(vars, x) %>%
     map(paste, collapse = " + ") %>% 
     sprintf("PV ~ SEX + GAN + %s", .) %>% 
     map(lm, data = df) %>% 
     map_dfr(broom::tidy, .id = "model")
    

    这会返回一个唯一的数据框。每个模型由列model 标识。 如果您更喜欢列表,只需将 map_dfr 替换为 map 并删除 .id


    在这里,我为您的数据创建了一个可重现的示例:

    set.seed(1)
    df <- data.frame(ID = 1:1000,
                     SEX = sample(0:1, 1000, replace = TRUE),
                     PV  = rnorm(1000),
                     GAN = rnorm(1000),
                     GAE = rnorm(1000))
    newcols <- unlist(lapply(c("GT_rs", "N_rs", "E_rs"), paste0, sample(100, 50)))
    df[newcols] <- replicate(50, rnorm(1000))
    
    df
    

    【讨论】:

    • 谢谢@Edo。正如 iago 指出的那样,您的解决方案的问题在于假设了一个算术序列,而事实并非如此。
    • 只需将 seq_len(n) 替换为您想要的任何数字向量。像 c(1,2,7,102,306)。无论如何它都会工作
    • 或者你不知道这些数字,你需要找到它们?
    • 谢谢@Edo。我不知道数字
    • 我编辑了答案。现在有了新的编辑,您可以看到您不再需要知道这些数字。即使缺少 GT 或 N 之一,它也能正常工作。
    【解决方案3】:

    由于 Edo 编辑了它的解决方案,我添加了它的变体:

    library(purrr)
    library(dplyr)
    library(broom)
    
    list("GT_rs", "N_rs") %>% 
        map(~dfA %>%  
                 select(matches(paste0(.x,"\\d+"))) %>% 
                 names %>% 
                 sub(pattern = .x, replacement = "")) %>% 
        reduce(intersect) %>% # until here we get the variables GT_rsx, N_rsx
        sprintf("PV ~ SEX + GAN + GT_rs%s + N_rs%s", ., .) %>%
        map(lm, data = dfA) %>%
        map_dfr(tidy, .id = "model") %>% 
        group_by(model) %>% 
        mutate(suffix = sub("N_rs", "", term[grepl("^N_rs\\d+$", term)]))
    

    【讨论】:

    • 谢谢@iago!。你的解决方案非常好。请问,在最后一步中,如何将“model”列的ID替换为N_rs ID?非常感谢
    • @Lucas 我更新了答案以获取带有 N_rs ID 的列 suffix
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2014-10-31
    • 2019-06-03
    • 1970-01-01
    • 1970-01-01
    • 2019-07-29
    • 2021-03-25
    相关资源
    最近更新 更多