R中的循环逻辑回归答案

【问题标题】：loop logistic regression in RR中的循环逻辑回归
【发布时间】：2022-01-10 18:36:34
【问题描述】：

我对 R 非常陌生。我有一个包含超过 50 万个基因位点的庞大数据集。我正在努力通过我的自变量循环逻辑回归。数据只有 10 行（10 名患者），但有 500K+ 列（遗传位点和表型）。

我需要循环运行每个基因位点一次，每个模型就像表型 = 每个基因位点 + 年龄 + 性别。（遗传位点是每个模型中唯一不同的变量）。

如何创建循环？我试过了：

logistic_regression <- lapply(mydata [,1:n], function(x) glm(outcome ~ geneticsite[,x] + age + gender , family=binomial,data=mydata))

我不断收到一条警告消息，提示找不到 geneticsite。有人可以帮我吗？谢谢！

更多详情：这是一项针对非常罕见表型的初步研究，因此我们只能测试 10 名患者。我们希望将 OR（遗传位点）和 P 值（遗传位点）作为输出。

【问题讨论】：

重组你的数据，让每个基因位点都是它自己的变量，然后你可以使用 GROUP_BY 和 tidymodels 在一次调用中进行回归难道不是有意义的吗？
我认为这是您需要的答案：coderedirect.com/questions/523195/…
首先一个包含 10 个观测值和 2-3 个预测变量的多变量模型可能不是很稳定，另外你打算用 500k 模型做什么？如果您有每个站点变量名称的向量，您可以执行以下操作：sites <- c('site1', 'site2'); lr <- lapply(sites, function(x) glm(reformulate(c('age', 'gender', x), 'outcome'), family = binomial, data = mydata))

标签： r loops logistic-regression

【解决方案1】：

这是我使用的一些示例数据

mydata <- 
  data.frame(
    GeneticSite = c(1,1,1,1,1,1,1,1),
    GeneticSite1 = c(2,3,2,3,2,3,2,3),
    Age = c(17, 52, 19, 18, 62, 53, 41, 24),
    Gender = c(1,2,1,1,2,2,2,1),
    Outcome = c(1,1,1,1,0,0,0,1)
  )

然后我使用library(dplyr) 对数据进行更长时间的旋转，以便我们可以按所有遗传位点进行分组。然后我创建线性回归函数 Age + Gender 并使用 summary 来获取我用来计算 OR 值的所有系数。

Linear_regression <- mydata %>%
  pivot_longer(starts_with("GeneticSite"),
               names_to = ".value",
               names_pattern = "(^GeneticSite)") %>% 
  group_by(GeneticSite) %>% 
  do(Reg = lm(Outcome ~ Age + Gender, data = .)) %>%
  summarise(GeneticSite,rsq = summary(Reg)$r.squared,
            Intercept_Estimate = summary(Reg)$coefficients[1],
            Age_Estimate = summary(Reg)$coefficients[2],
            Gender_Estimate = summary(Reg)$coefficients[3],
            Intercept_Std_Error = summary(Reg)$coefficients[4],
            Age_Std_Error = summary(Reg)$coefficients[5],
            Gender_Std_Error = summary(Reg)$coefficients[6],
            Intercept_t_value = summary(Reg)$coefficients[7],
            Age_t_value = summary(Reg)$coefficients[8],
            Gender_t_value = summary(Reg)$coefficients[9],
            Intercept_p_value = summary(Reg)$coefficients[10],
            Age_p_value = summary(Reg)$coefficients[11],
            Gender_p_value = summary(Reg)$coefficients[12]) %>%
  mutate(OR_Estimate = exp(Intercept_Estimate)/(1+exp(Intercept_Estimate))) %>% 
  mutate(OR_Age = exp(Age_Estimate)/(1+exp(Age_Estimate))) %>% 
  mutate(OR_Gender = exp(Age_Estimate)/(1+exp(Age_Estimate)))

这就是我知道如何标记每一列的方式

> summary(y)$coefficients

             Estimate     Std. Error   t value     Pr(>|t|)
(Intercept)  1.750000e+00  0.3331066  5.253574e+00 0.0001558944
x$Age        9.688771e-18  0.0151608  6.390673e-16 1.0000000000
x$Gender    -7.500000e-01  0.5211766 -1.439052e+00 0.1737764924

示例输出：

"GeneticSite","rsq","Intercept_Estimate","Age_Estimate","Gender_Estimate","Intercept_Std_Error","Age_Std_Error","Gender_Std_Error","Intercept_t_value","Age_t_value","Gender_t_value","Intercept_p_value","Age_p_value","Gender_p_value","OR_Estimate","OR_Age","OR_Gender"
"1",1,0.6,1.75,-9.13466064460698e-18,-0.749999999999999,0.537118251352767,0.0244460541141712,0.840372000724438,3.25812797385401,-3.73666056777305e-16,-0.892461908956351,0.0224922711346427,1,0.413031579096883,0.85195280196831,0.5,0.5
"2",2,1,2,9.08261087251563e-18,-1,4.20114922415409e-16,1.48859041349867e-17,5.45878506026351e-16,4760602143102169,0.610148418944103,-1831909461464908,1.33726733138161e-16,0.651229019050892,3.47517050246851e-16,0.880797077977882,0.5,0.5
"3",3,0.351351351351351,1.21621621621622,-0.0270270270270271,0.351351351351353,2.0286480005184,0.162162162162162,5.15550724280418,0.599520575233073,-0.166666666666667,0.0681506852389265,0.656182728286001,0.894863086577493,0.956680908817521,0.771396988381907,0.493243654508354,0.493243654508354

【讨论】：

由于geneticSite 是单独的列不起作用，需要将pivot_longer() 添加到IMO 解决方案中。
哦，真的，谢谢！
这给出了相关性和相应的 p 值，对吗？我正在寻找 OR 和 P 值，特别是在遗传学位点和结果之间（这是一个二进制变量 1 对 0）。谢谢！