Регрессия лёсса в каждой группе с помощью dplyr::group_by()
Хорошо, я машу своим белым флагом.
Я пытаюсь вычислить регрессию лесса на моем наборе данных.
Я хочу, чтобы Лесс вычислила различный набор точек, которые изображены как плавные линии для каждой группы.
Проблема в том, что при вычислении лёсса экранируется функция 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 как столбец в отдельной строке, я сталкиваюсь с проблемами, потому что я вызвал dplyr::range ранее.
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
3 ответа
Это аккуратный 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.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))))
Обратите внимание, что для выполнения лессов требуется минимальное количество наблюдений (~4? Я точно не помню). Кроме того, это займет некоторое время для запуска, поэтому протестируйте часть данных, чтобы убедиться, что они работают правильно.
К сожалению, описанные выше подходы в моем случае не сработали. Таким образом, я реализовал предсказание Лесса в виде обычной функции, которая работала очень хорошо. В приведенном ниже примере данные содержатся вdf
кадр данных, пока мы группируем поdf$profile
и хотите вписать предсказание Лесса вdf$daily_sum
ценности.