使用 dplyr::group_by() 对每个组进行黄土回归答案

【问题标题】：loess regression on each group with dplyr::group_by()使用 dplyr::group_by() 对每个组进行黄土回归
【发布时间】：2018-10-14 05:50:48
【问题描述】：

好的，我正在挥舞我的白旗。

我正在尝试在我的数据集上计算黄土回归。

我希望 loess 计算一组不同的点，为每个组绘制一条平滑线。

问题是黄土计算逃逸dplyr::group_by函数，所以黄土回归是在整个数据集上计算的。

互联网搜索让我相信这是因为 dplyr::group_by 不应该以这种方式工作。

我只是不知道如何在每个组的基础上进行这项工作。

以下是我尝试失败的一些示例。

test2 <- test %>% 
  group_by(CpG) %>% 
  dplyr::arrange(AVGMOrder) %>% 
  do(broom::tidy(predict(loess(Meth ~ AVGMOrder, span = .85, data=.))))

> test2
# A tibble: 136 x 2
# Groups:   CpG [4]
   CpG            x
   <chr>      <dbl>
 1 cg01003813 0.781
 2 cg01003813 0.793
 3 cg01003813 0.805
 4 cg01003813 0.816
 5 cg01003813 0.829
 6 cg01003813 0.841
 7 cg01003813 0.854
 8 cg01003813 0.866
 9 cg01003813 0.878
10 cg01003813 0.893

这个可行，但我不知道如何将结果应用于原始数据框中的列。我想要的结果是第 x 列。如果我将 x 作为列应用到单独的行中，我会遇到问题，因为我之前调用了 dplyr::arrange。

test2 <- test %>% 
  group_by(CpG) %>% 
  dplyr::arrange(AVGMOrder) %>% 
  dplyr::do({
    predict(loess(Meth ~ AVGMOrder, span = .85, data=.))
  })

这只是失败并出现以下错误。

“错误：结果 1、2、3、4 必须是数据帧，而不是数字”

而且它仍然没有作为dplyr::mutate 的新列应用

fems <- fems %>% 
  group_by(CpG) %>% 
  dplyr::arrange(AVGMOrder) %>% 
  dplyr::mutate(Loess = predict(loess(Meth ~ AVGMOrder, span = .5, data=.)))

这是我的第一次尝试，大部分类似于我想做的事情。问题是这个是对整个数据帧而不是每个 CpG 组执行黄土预测。

我真的被困在这里了。我在网上读到 purr 包可能会有所帮助，但我无法弄清楚。

数据如下所示：

> head(test)
    X geneID        CpG                                        CellLine       Meth AVGMOrder neworder Group SmoothMeth
1  40     XG cg25296477 iPS__HDF51IPS14_passage27_Female____165.592.1.2 0.81107210         1        1     5  0.7808767
2  94     XG cg01003813 iPS__HDF51IPS14_passage27_Female____165.592.1.2 0.97052120         1        1     5  0.7927130
3 148     XG cg13176022 iPS__HDF51IPS14_passage27_Female____165.592.1.2 0.06900448         1        1     5  0.8045080
4 202     XG cg26484667 iPS__HDF51IPS14_passage27_Female____165.592.1.2 0.84077890         1        1     5  0.8163997
5  27     XG cg25296477  iPS__HDF51IPS6_passage33_Female____157.647.1.2 0.81623880         2        2     3  0.8285259
6  81     XG cg01003813  iPS__HDF51IPS6_passage33_Female____157.647.1.2 0.95569240         2        2     3  0.8409501

独特（测试$CpG） [1] “cg25296477” “cg01003813” “cg13176022” “cg26484667”

所以，明确地说，我想对我的数据框中的每个唯一 CpG 进行黄土回归，将生成的“回归 y 轴值”应用于与原始 y 轴值 (Meth) 匹配的列。

我的实际数据集有几千个这样的 CpG，而不仅仅是四个。

https://docs.google.com/spreadsheets/d/1-Wluc9NDFSnOeTwgBw4n0pdPuSlMSTfUVM0GJTiEn_Y/edit?usp=sharing

【问题讨论】：

您看过 R for Data Science 的 Many Models 章节吗？它完成了一个非常相似的练习
我去看看。谢谢。
所以您希望将黄土预测值作为数据集中的附加列？我认为您可以将第一个示例中的do(broom::tidy...) 更改为do(x = broom::tidy...)，或使用broom::augment。将测试我何时可以制作一些数据或您是否提供一些数据
现在正在尝试。还添加了谷歌表格链接以测试数据框
这可能对您有所帮助：stackoverflow.com/a/49616753/8583393

标签： r dplyr purrr loess broom

【解决方案1】：

这是一种使它工作的整洁的 Tidyverse 方法：

library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)

models <- fems %>%
        tidyr::nest(-CpG) %>%
        dplyr::mutate(
                # Perform loess calculation on each CpG group
                m = purrr::map(data, loess,
                               formula = Meth ~ AVGMOrder, span = .5),
                # Retrieve the fitted values from each model
                fitted = purrr::map(m, `[[`, "fitted")
        )

# Apply fitted y's as a new column
results <- models %>%
        dplyr::select(-m) %>%
        tidyr::unnest()

# Plot with loess line for each group
ggplot(results, aes(x = AVGMOrder, y = Meth, group = CpG, colour = CpG)) +
        geom_point() +
        geom_line(aes(y = fitted))

【讨论】：

如果您能解释一下嵌套和映射发生了什么，那就太好了。我很惊讶第一张地图以data 作为参数，而不是一个点或数据框的实际名称。不过这对我有用，所以谢谢！
@MokeEire 默认情况下，nest 将分配 data 作为新列的名称。 nest 基本上从每个组中制作一个数据框，以提供一系列较小的数据框。因此，通过将 data 作为其第一个参数，map 将依次获取每个嵌套数据帧并在每个数据帧上计算 loess。
如何将其应用于不同的群体？但到不同的群体和不同的列？我有很多指标需要去趋势。所以必须适合许多非参数样条。
当所有Meth == 0.8 组的CpG 时，您将如何提取x 值？我知道cg13176022 会导致NA，但我有类似的情况，我需要为每个组的建模y 值找到x 值。好奇predict 在哪里发挥作用。谢谢。

【解决方案2】：

您可能已经想通了 - 但如果没有，这里有一些帮助。

基本上，您需要向 predict 函数提供您想要预测的值的 data.frame（向量也可以工作，但我没有尝试过）。

所以对于你的情况：

fems <- fems %>% 
  group_by(CpG) %>% 
  arrange(CpG, AVGMOrder) %>% 
  mutate(Loess = predict(loess(Meth ~ AVGMOrder, span = .5, data=.),
    data.frame(AVGMOrder = seq(min(AVGMOrder), max(AVGMOrder), 1))))

注意，loess 需要最少数量的观察来运行（~4？我记不清了）。此外，这需要一段时间才能运行，因此请使用您的数据片段进行测试，以确保其正常工作。

【讨论】：