Построение кривых полиномиальной регрессии в R

Я провел полиномиальную регрессию для данных, которые я включаю, от квадратичного до септического, но я застрял, пытаясь построить эти регрессионные кривые на моем графике рассеяния. Я прошу помощи в создании кода, который будет работать для каждого полиномиального порядка.

Time <- Mangan_2008_total$YearMonthDay
Abundance <- Mangan_2008_total$`Total Nymph Aa`

Выше приведен код, который я использовал, чтобы изолировать данные, которые я включил в этот пост, от более крупного набора данных. Ниже приведены данные для справки.

(В настоящее время даты указаны в формате григорианского календаря. Я планирую в какой-то момент преобразовать их в календарь юлианских дней.)

Time
1   20080301
2   20080316
3   20080402
4   20080416
5   20080428
6   20080514
7   20080527
8   20080608
9   20080627
10  20080709
11  20080722
12  20080806
13  20080818
14  20080901
15  20080915
16  20080930
17  20081013
18  20081029
19  20081110
20  20081124

Abundance
1   0
2   0
3   26
4   337
5   122
6   232
7   190
8   381
9   148
10  201
11  69
12  55
13  35
14  15
15  6
16  1
17  0
18  1
19  0
20  0

Я скомпилировал эти данные во фрейм данных для упрощения манипуляций:

Mangan_2008_nymph <- data.frame(Time, Abundance)

Вот код диаграммы рассеяния, которую я сделал в ggplot:

Nymph_2008_Plot <- ggplot(Mangan_2008_nymph, aes(Time, Abundance)) + 
  geom_point(size=4, col='red') + ggtitle("2008 Amblyomma americanum Abundance") + 
  xlab("Time") + ylab("Nymph Abundance")
Nymph_2008_Plot

Вот код, который я использовал для регрессионного анализа (чтобы запустить полиномиальную регрессию более высокого порядка, я просто поменял местами 2 (значение степени) на соответствующий полиномиальный порядок):

Quadratic_2008_Nymph <- lm(Abundance ~ poly(Time, 2))
summary(Quadratic_2008_Nymph)

Вот где я застреваю. Как мне нанести на график кривые полиномиальной регрессии? Если есть способ сделать это с помощью формата ggplot, который будет предпочтительнее. Если построение этих полиномиальных кривых на графике ggplot не сработает, я переключу форматирование.

Заранее спасибо и прокомментируйте, если мне нужно уточнить / предоставить дополнительную информацию.


r regression plot ggplot2 non-linear-regression
person Ben G.    schedule 04.12.2020    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Первое, что я хотел бы сделать здесь, это преобразовать числа, которые вы рассматриваете как даты, в настоящие даты. Если вы этого не сделаете, lm даст неверный результат; Например, строки 1 и 2 вашего фрейма данных представляют данные с интервалом в 15 дней (20080316 - 20080301 = 15), но тогда строки 2 и 3 разнесены на 17 дней, но при регрессии они будут различаться на 86 дней (20080402 - 20080316 = 86). Это приведет к тому, что lm выдаст бессмысленные результаты.

Хорошей практикой является выработка привычки как можно раньше менять числа или символьные строки, представляющие данные даты и времени, на фактические даты и время. Это значительно упрощает остальную часть вашего кода. В вашем случае это было бы просто:

Mangan_2008_nymph$Time <- as.Date(as.character(Mangan_2008_nymph$Time), "%Y%m%d")

Для самого графика вы можете добавить линии полиномиальной регрессии, используя geom_smooth. Вы указываете метод lm и формулу (в терминах x и y, не в терминах имен переменных). Также неплохо сопоставить линию с эстетикой цвета, чтобы она отображалась в легенде. Сделайте это один раз для каждого полиномиального порядка, который вы хотите добавить.

library(ggplot2)

ggplot(Mangan_2008_nymph, aes(Time, Abundance)) + 
  geom_point(size = 4, col = 'red') + 
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), aes(color = "2"), se = F) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 3), aes(color = "3"), se = F) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 4), aes(color = "4"), se = F) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 5), aes(color = "5"), se = F) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 6), aes(color = "6"), se = F) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 7), aes(color = "7"), se = F) +
  labs(x = "Time", y = "Nymph abundance", color = "Degree") +
  ggtitle("2008 Amblyomma americanum Abundance")

Лично я считаю, что это приводит к тому, что сюжет становится немного загроможденным, и я бы подумал об отказе от пары полиномиальных уровней, чтобы его было легче читать. Вы также можете легко добавить несколько стилистических поправок:

ggplot(Mangan_2008_nymph, aes(Time, Abundance)) + 
  geom_point(size = 2, col = 'gray50') + 
  geom_smooth(method = "lm", formula = y ~ poly(x, 3), aes(color = "3"), se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 5), aes(color = "5"), se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 7), aes(color = "7"), se = FALSE) +
  scale_color_brewer(palette = "Set1") +
  labs(x = "Time", y = "Nymph abundance", color = "Degree") +
  ggtitle("2008 Amblyomma americanum Abundance") +
  theme_classic() +
  theme(panel.grid.major = element_line())

введите описание изображения здесь

Используемые данные

Mangan_2008_nymph <- structure(list(Time = c(20080301L, 20080316L, 20080402L, 
20080416L, 20080428L, 20080514L, 20080527L, 20080608L, 20080627L, 20080709L, 
20080722L, 20080806L, 20080818L, 20080901L, 20080915L, 20080930L, 
20081013L, 20081029L, 20081110L, 20081124L), Abundance = c(0L, 
0L, 26L, 337L, 122L, 232L, 190L, 381L, 148L, 201L, 69L, 55L, 
35L, 15L, 6L, 1L, 0L, 1L, 0L, 0L)), class = "data.frame", row.names = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", 
"14", "15", "16", "17", "18", "19", "20"))
person Allan Cameron    schedule 04.12.2020