使用map对每个列执行回归,并整理R中的结果

vx6bjr1n  于 2023-06-27  发布在  其他
关注(0)|答案(2)|浏览(67)

我有一个纵向研究的数据集,其中有10个结果变量和16个自变量。但为了简化我的问题,我编辑了数据集,现在它包含:

  • 4列作为结果变量
  • 1列作为自变量(TimePointNum2)
  • 6例受试者(受试者)
  • 每例受试者随访5次(时间点编号2)

数据集可以通过here下载。如果您不满意,我已经使用dput()为dataframe m2创建了对象(参见下面的数据)

library(tidyverse)
library(broom)
library(lme4)
library(nlme)
m2 <- read_csv('question.csv')

我的目标是:

  • 针对 Dataframe 的每个结果变量(MUFA、SFA、体重、体脂)运行(迭代或循环)回归模型,其中每个结果的自变量是TimePointNum2。

所以我运行了很多次lmer(这是违反DRY规则的)

mod_SFA_num <- lmer(SFA ~ TimePointNum2 + (1| Subject), data = m2)
mod_MUFA_num <- lmer(MUFA ~ TimePointNum2 + (1| Subject), data = m2)
mod_BodyFat_num <- lmer(BodyFat ~ TimePointNum2 + (1| Subject), data = m2)
mod_weight_num <- lmer(Weight ~ TimePointNum2 + (1| Subject), data = m2)

我看到了一些使用purrr::map的建议,最相似的方法是

mtcars %>% 
  names() %>% 
  paste('mpg ~', .) %>% 
  map(as.formula) %>% 
  map(lm, data = mtcars) %>%
  map_df(~broom::tidy(.))

受上面例子的启发,所以我运行这些代码:

m2 %>%
  names() %>%
  paste(., '~ TimePoint + (1| Subject)') %>%
  map(as.formula) %>%
  map(lmer, data = m2) %>%
  map(~broom::tidy(.))

但是我收到了一个错误(很可能是第5行的lmer

Error in `map()`:
ℹ In index: 1.
Caused by error in `mkRespMod()`:
! response must be numeric
Run `rlang::last_trace()` to see where the error occurred.

我的期望输出类似于下表:
| 模型||价值|Std.Error| DF| t值|p值|
| - -----|- -----|- -----|- -----|- -----|- -----|- -----|
| | | | | | | |
| 身体脂肪|(截距)|28.4583| 0.68122|三百七十五|41.7758| 0|
| | 时间点编号2| 0.06894| 0.02536|三百七十五|2.71861| 0.0069|
| 重量|(截距)|16.4157| 0.42985|三百七十五|38.1895| 0|
| | 时间点编号2| 0.08064| 0.01945|三百七十五|4.14603| 0|
| SFA|(截距)|67.1564| 1.79525|三百七十四|37.4079| 0|
| | 时间点编号2| 0.73767| 0.48301|三百七十四|1.52724| 0.1275|
| PUFA|(截距)|12.8161| 0.85069|三百六十五|十五点零六五六|0|
| | 时间点编号2|-0.5978| 0.23454|三百六十五|-2.5491| 0.0112|
我的问题是:
1.我们怎样才能纠正这个错误呢?
1.如何使用purrr::map对每个感兴趣的列进行循环回归?
1.如何生成所需的表(如上所示,可能使用tidy())?
注:tidyverse友好方法优先
数据

m2 <- structure(
  list(
    Subject = c(
      "B001",
      "B001",
      "B001",
      "B001",
      "B001",
      "B002",
      "B002",
      "B002",
      "B002",
      "B002",
      "B003",
      "B003",
      "B003",
      "B003",
      "B003",
      "B004",
      "B004",
      "B004",
      "B004",
      "B004",
      "B005",
      "B005",
      "B005",
      "B005",
      "B005",
      "B006",
      "B006",
      "B006",
      "B006",
      "B006"
    ),
    Weight = c(
      68.5,
      69.1,
      70.2,
      70,
      70.1,
      56.7,
      56.1,
      55.8,
      56.5,
      55.2,
      62.9,
      62.5,
      62,
      62.7,
      63.3,
      63.9,
      62.6,
      62.1,
      61.5,
      62.6,
      53.5,
      52,
      51.8,
      52.4,
      52.9,
      77,
      77.1,
      77.8,
      79.2,
      78.7
    ),
    BodyFat = c(
      21.8,
      21.1,
      21.7,
      22.2,
      22,
      34,
      33.3,
      33,
      32.4,
      32.8,
      36.9,
      36.6,
      36.1,
      36.8,
      36.9,
      33.6,
      32.5,
      32.7,
      31.9,
      32.2,
      12.2,
      10.6,
      11,
      11.5,
      11.6,
      23,
      22.5,
      22.6,
      23.1,
      23.9
    ),
    SFA = c(
      62.895,
      84.705,
      83.49,
      66.64,
      72.19,
      56.195,
      93.945,
      92.635,
      88.51,
      83.505,
      92.67,
      90.81,
      83.37,
      90.205,
      84.195,
      88.065,
      53.69,
      93.14,
      52.57,
      4.505,
      63.505,
      92.59,
      80.87,
      89.125,
      67.305,
      71.41,
      89.1,
      91.175,
      80.635,
      89.425
    ),
    MUFA = c(
      14.455,
      11.71,
      12.58,
      21.135,
      26.175,
      13.285,
      4.91,
      6.08,
      6.735,
      12.745,
      7.325,
      7.605,
      15.735,
      7.985,
      12.32,
      7.92,
      42.045,
      4.57,
      24.305,
      0,
      18.585,
      5.955,
      14.815,
      8.775,
      20.295,
      11.485,
      6.745,
      6.32,
      10.805,
      6.8
    ),
    TimePointNum2 = c(
      1,
      2,
      3,
      4,
      5,
      1,
      2,
      3,
      4,
      5,
      1,
      2,
      3,
      4,
      5,
      1,
      2,
      3,
      4,
      5,
      1,
      2,
      3,
      4,
      5,
      1,
      2,
      3,
      4,
      5
    )
  ),
  row.names = c(NA, -30L),
  class =  "data.frame"
)
nimxete2

nimxete21#

要从混合效应模型获得整洁的输出,您需要使用tidy() from {broom.mixed}。您可以按照nest-unnest工作流在每列上Map回归模型,并从中获得整洁的结果。

library(tidyverse)
library(broom.mixed)
library(lme4)
library(nlme)

m2 %>%
  pivot_longer(cols = -c("Subject", "TimePointNum2"), 
               names_to = "outcome", values_to = "value") %>% 
  mutate(outcome_fct = as.factor(outcome)) %>% 
  nest(data = -outcome_fct) %>% 
  mutate(
    data = map(data, ~ pivot_wider(.x, names_from = "outcome", values_from = "value")),
    formula = map(data, ~ reformulate("TimePointNum2 + (1| Subject)", 
                                      response = tail(names(.x), n = 1))),
    fit = map2(data, formula, ~ lmer(formula = .y, data = .x)),
    tidied = map(fit, broom.mixed::tidy)
  ) %>% 
  select(outcome_fct, tidied) %>% 
  unnest(tidied)

#> # A tibble: 16 × 7
#>    outcome_fct effect   group    term            estimate std.error statistic
#>    <fct>       <chr>    <chr>    <chr>              <dbl>     <dbl>     <dbl>
#>  1 Weight      fixed    <NA>     (Intercept)      63.4       3.77      16.8  
#>  2 Weight      fixed    <NA>     TimePointNum2     0.0583    0.0977     0.597
#>  3 Weight      ran_pars Subject  sd__(Intercept)   9.20     NA         NA    
#>  4 Weight      ran_pars Residual sd__Observation   0.757    NA         NA    
#>  5 BodyFat     fixed    <NA>     (Intercept)      26.6       3.87       6.87 
#>  6 BodyFat     fixed    <NA>     TimePointNum2    -0.0483    0.0703    -0.687
#>  7 BodyFat     ran_pars Subject  sd__(Intercept)   9.45     NA         NA    
#>  8 BodyFat     ran_pars Residual sd__Observation   0.545    NA         NA    
#>  9 SFA         fixed    <NA>     (Intercept)      83.0       8.13      10.2  
#> 10 SFA         fixed    <NA>     TimePointNum2    -1.74      2.29      -0.760
#> 11 SFA         ran_pars Subject  sd__(Intercept)   7.16     NA         NA    
#> 12 SFA         ran_pars Residual sd__Observation  17.7      NA         NA    
#> 13 MUFA        fixed    <NA>     (Intercept)      11.8       3.62       3.26 
#> 14 MUFA        fixed    <NA>     TimePointNum2     0.189     1.09       0.173
#> 15 MUFA        ran_pars Subject  sd__(Intercept)   0        NA         NA    
#> 16 MUFA        ran_pars Residual sd__Observation   8.45     NA         NA
cygmwpex

cygmwpex2#

# select the columns which are double type which aren't TimePointNum2
variables <- m2 %>%
    select_if(is.double) %>%
    colnames() %>%
    .[.!= "TimePointNum2"]
    
    

variables %>%
    paste("~ TimePointNum2 + (1| Subject)") %>%    # create a formula for each variable
    map(lmer, data = m2) %>%          # fit a linear mixed model for each variable
    map(broom::tidy) %>%            # tidy the results
    bind_rows() %>%                 # bind the results into a single data frame
    mutate(variable = rep(variables, each = 4)) %>% # add the variable name as a column
    select(variable, everything()) # move the variable column to the start

# A tibble: 16 × 7
   variable effect   group    term            estimate std.error statistic
   <chr>    <chr>    <chr>    <chr>              <dbl>     <dbl>     <dbl>
 1 Weight   fixed    NA       (Intercept)      63.4       3.77      16.8  
 2 Weight   fixed    NA       TimePointNum2     0.0583    0.0977     0.597
 3 Weight   ran_pars Subject  sd__(Intercept)   9.20     NA         NA    
 4 Weight   ran_pars Residual sd__Observation   0.757    NA         NA    
 5 BodyFat  fixed    NA       (Intercept)      26.6       3.87       6.87 
 6 BodyFat  fixed    NA       TimePointNum2    -0.0483    0.0703    -0.687
 7 BodyFat  ran_pars Subject  sd__(Intercept)   9.45     NA         NA    
 8 BodyFat  ran_pars Residual sd__Observation   0.545    NA         NA    
 9 SFA      fixed    NA       (Intercept)      83.0       8.13      10.2  
10 SFA      fixed    NA       TimePointNum2    -1.74      2.29      -0.760
11 SFA      ran_pars Subject  sd__(Intercept)   7.16     NA         NA    
12 SFA      ran_pars Residual sd__Observation  17.7      NA         NA    
13 MUFA     fixed    NA       (Intercept)      11.8       3.62       3.26 
14 MUFA     fixed    NA       TimePointNum2     0.189     1.09       0.173
15 MUFA     ran_pars Subject  sd__(Intercept)   0        NA         NA    
16 MUFA     ran_pars Residual sd__Observation   8.45     NA         NA

相关问题