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

#r #ggplot2 #plot #regression #non-linear-regression

#r #ggplot2 #график #регрессия #нелинейная регрессия

Вопрос:

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

 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 не сработает, я переключу свое форматирование.

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

Ответ №1:

Первое, что я хотел бы сделать здесь, это преобразовать числа, которые вы рассматриваете как даты, в фактические даты. Если вы этого не сделаете, это 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"))