R - 通过缺少值的 Dataframe 逐列循环 lm()答案

【问题标题】：R - Loop lm() columnwise through Dataframe with missing ValuesR - 通过缺少值的 Dataframe 逐列循环 lm()
【发布时间】：2020-08-28 10:06:15
【问题描述】：

我正在处理一个由至少三个变量（波长、辐照度、x）组成的数据框，我已经旋转了这些变量，以便每个波长都是一个新行，从而允许我在每个波长上运行 lm() 并提取系数，所以我可以看到 x 如何随波长和辐照度变化。

但是，我能够让它工作的唯一方法是在每个波长上显式运行 lm()。这对于较大的数据帧是不可行的，因为数据帧将具有数百个参数，这些参数会随着波长和辐照度的变化而变化。

我觉得这可以使用“应用”或编写循环来解决，但我没有任何运气让它们工作。

请参阅下面的问题示例。

我还是很新，所以任何指针表示赞赏

irr = rnorm(33, 10, 3)
wave = c(290, 290, 290, 300, 300, 300, 310, 310, 310, 320, 320, 320, 330, 330, 330, 340, 340, 340, 350, 350, 350, 360, 360, 360, 370, 370, 370, 380, 380, 380, 400, 400, 400)
x = rnorm(33, 50, 2)
df <- as.data.frame(cbind(wave, irr, x))
df_wide <- df %>%
  pivot_wider(names_from = "wave",
              values_from = "x")
"290_lm" <- lm(df_wide$`290` ~ df_wide$irr) 
"300_lm" <- lm(df_wide$`300` ~ df_wide$irr) #etc through each wavelength

## Attempt at loop

for (i in 2:(ncol(df_wide))){
  irr <- df_wide[2][i]
  lm_function <- paste(irr,
                       sep = "~")
  df_lm = lm(lm_function, 
             data = df_wide[2:12])
}

【问题讨论】：

我建议不要旋转，而是使用 lm 的 subset 选项。我也想知道线性模型中的因变量和自变量是什么。

标签： r loops apply lm

【解决方案1】：

当您使用长格式时，这可能会容易得多。只需使用lapply 子集您的数据。使用setNames，结果列表会得到好听的名字。

res <- setNames(lapply(unique(df$wave), function(w) 
  lm(x ~ irr, data=df[df$wave %in% w, ])),
  paste0("wave.", unique(df$wave)))
res
# $wave.290
# 
# Call:
#   lm(formula = x ~ irr, data = df[df$wave %in% w, ])
# 
# Coefficients:
#   (Intercept)          irr  
# 36.837        1.503  
# 
# 
# $wave.300
# 
# Call:
#   lm(formula = x ~ irr, data = df[df$wave %in% w, ])
# 
# Coefficients:
#   (Intercept)          irr  
# 54.3785      -0.5586 
# [...]

【讨论】：

【解决方案2】：

据我从您的描述中可以看出，您的问题与purrr::map 的示例相同，从而避免了扩展的需要。

library(dplyr)
library(purrr)

results_list <- 
  df %>% 
  split(.$wave) %>% 
  map(~ lm(x ~ irr, data = .x)) %>% 
  map(summary)

results_list$`350`
#> 
#> Call:
#> lm(formula = x ~ irr, data = .x)
#> 
#> Residuals:
#>      19      20      21 
#>  0.2924 -2.2947  2.0023 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)  
#> (Intercept)  52.7276     6.2200   8.477   0.0748 .
#> irr          -0.4977     0.6229  -0.799   0.5708  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.059 on 1 degrees of freedom
#> Multiple R-squared:  0.3897, Adjusted R-squared:  -0.2206 
#> F-statistic: 0.6385 on 1 and 1 DF,  p-value: 0.5708

根据您的数据

irr = rnorm(33, 10, 3)
wave = c(290, 290, 290, 300, 300, 300, 310, 310, 310, 320, 320, 320, 330, 330, 330, 340, 340, 340, 350, 350, 350, 360, 360, 360, 370, 370, 370, 380, 380, 380, 400, 400, 400)
x = rnorm(33, 50, 2)
df <- as.data.frame(cbind(wave, irr, x))

^{由reprex package (v0.3.0) 于 2020 年 5 月 12 日创建}

【讨论】：

【解决方案3】：

其他解决方案

library(tidyverse)
library(generics)
df %>% 
  group_by(wave) %>% 
  nest() %>% 
  mutate(model = map(data, ~ lm(x ~ irr, data = .x) %>% tidy)) %>% 
  select(-data) %>% 
  unnest(model)

【讨论】：

不错的方法，但我还必须加载包 generics。（使用 R 4.0.0 测试，上周的软件包）。
函数tidy或更准确地说是tidy.lm在包broom中找到：包generics只是另一个包装器。所以也可以使用library(broom)而不是library(generics)

【解决方案4】：

或者如下：

df <- data.frame(
  irr = rnorm(33, 10, 3),
  wave = c(290, 290, 290, 300, 300, 300, 310, 310, 310, 
           320, 320, 320, 330, 330, 330, 340, 340, 340, 350, 350, 350, 
           360, 360, 360, 370, 370, 370, 380, 380, 380, 400, 400, 400),
  x = rnorm(33, 50, 2)

)

mylm <- function(w) {
  m <- lm(x ~ irr, data = df, subset = (wave == w))
  ## outcomment the following if you just need the parameters
  # coef(m)
}

lapply(df$wave, mylm)

【讨论】：