Как вывести линию предсказания модели перед geom_rect - PullRequest
Новая
       25

Как вывести линию предсказания модели перед geom_rect

0 голосов
/ 01 апреля 2020

Я исследую влияние времени с момента пожара на разнообразие видов. Я пытаюсь сделать график, который имеет разные цвета в разное время с момента пожара. Однако размещение цветов на графике привело к исчезновению линии прогнозирования модели. Мне интересно, есть ли какой-нибудь способ вывести строку перед geom_rect?

Загруженные пакеты: 

library(voxel)
library(gamm4)
library(ggplot2)

Мои данные: 

data <- read.csv('StacksOverflow.csv')

 structure(list(Lscape = c(158L, 158L, 158L, 158L, 158L, 158L), 
    TSF = c(5, 5, 5, 18.5, 5, 18.5), VegtypeNew = structure(c(1L, 
    1L, 1L, 2L, 1L, 1L), .Label = c("spinsandplain", "woodlndsandplain"
    ), class = "factor"), FF = c(2L, 2L, 2L, 1L, 2L, 1L), ThreeYearRain = c(913.799997, 
    913.799997, 913.799997, 913.799997, 913.799997, 913.799997
    ), Div = c(2.2629743, 1.9630117, 1.7336569, 1.2816843, 2.4155056, 
    1.4240443), triodia_low = c(19L, 6L, 21L, 32L, 11L, 32L)), row.names = c(NA, 
6L), class = "data.frame")

Расширенные данные: 

structure(list(Lscape = c(158L, 158L, 158L, 158L, 158L, 158L, 
158L, 158L, 201L, 201L, 201L, 201L, 201L, 201L, 201L, 201L, 235L, 
235L, 235L, 235L, 235L, 235L, 235L, 235L, 237L, 237L, 237L, 237L, 
237L, 237L, 237L, 237L, 254L, 254L, 254L, 254L, 254L, 254L, 254L, 
254L, 287L, 287L, 287L, 287L, 287L, 287L, 287L, 287L, 304L, 304L, 
304L, 304L, 304L, 304L, 304L, 304L, 311L, 311L, 311L, 311L, 311L, 
311L, 311L, 311L, 312L, 312L, 312L, 312L, 312L, 312L, 312L, 312L, 
323L, 323L, 323L, 323L, 323L, 323L, 323L, 323L, 326L, 326L, 326L, 
326L, 326L, 326L, 326L, 326L, 327L, 327L, 327L, 327L, 327L, 327L, 
327L, 327L, 337L, 337L, 337L, 337L, 337L, 337L, 337L, 337L, 355L, 
355L, 355L, 355L, 355L, 355L, 355L, 355L, 370L, 370L, 370L, 370L, 
370L, 370L, 370L, 370L, 379L, 379L, 379L, 379L, 379L, 379L, 379L, 
379L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 414L, 414L, 
414L, 414L, 414L, 414L, 414L, 414L, 435L, 435L, 435L, 435L, 435L, 
435L, 435L, 435L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 
438L, 438L, 438L, 438L, 438L, 438L, 438L, 438L, 447L, 447L, 447L, 
447L, 447L, 447L, 447L, 447L, 452L, 452L, 452L, 452L, 452L, 452L, 
452L, 452L), TSF = c(5, 5, 5, 18.5, 5, 18.5, 18.5, 18.5, 11.5, 
4.5, 0.5, 20, 11.5, 0.5, 1, 4.5, 1, 1, 4.5, 5, 4.5, 2, 5, 1, 
6, 6, 4.5, 6, 14.5, 17, 4.5, 6, 1, 1, 7, 4.5, 2, 2, 7, 7, 20, 
4, 3.5, 4, 3.5, 3.5, 11.5, 20, 6, 0.5, 5, 6, 6, 0.5, 7, 6, 3.5, 
3.5, 3.5, 11.5, 11.5, 1, 1, 11.5, 1, 1, 4, 1, 1, 4, 1, 10.5, 
7, 17.5, 0.5, 0.5, 0.5, 17.5, 7, 0.5, 18, 1.5, 3.5, 18, 18, 5, 
3.5, 18.5, 14.5, 1.5, 7, 1.5, 7, 7, 7, 7, 10.5, 1.5, 0, 1.5, 
7, 3, 7, 10.5, 0.5, 20, 0.5, 2, 2, 1.5, 2, 3, 20, 1, 1.5, 10.5, 
17, 1.5, 1.5, 10.5, 3, 1, 1, 1, 4.5, 1, 6.5, 1, 10, 1.5, 12.5, 
1.5, 1.5, 1.5, 1.5, 1.5, 2, 7, 12.5, 2, 7, 2, 2, 2, 1.5, 18.5, 
18.5, 1.5, 7, 1.5, 1.5, 5, 12.5, 6.5, 1.5, 1.5, 1.5, 1.5, 1.5, 
1.5, 6.5, 6.5, 1.5, 6.5, 18.5, 6.5, 7, 1.5, 1, 7, 7, 1, 7, 1, 
7.5, 7.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5), VegtypeNew = structure(c(1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 
2L, 2L, 1L, 2L, 2L, 2L, 2L), .Label = c("spinsandplain", "woodlndsandplain"
), class = "factor"), FF = c(2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 
3L, 3L, 4L, 2L, 2L, 5L, 3L, 4L, 5L, 5L, 2L, 2L, 4L, 3L, 5L, 4L, 
5L, 4L, 4L, 5L, 3L, 3L, 4L, 5L, 5L, 3L, 4L, 6L, 5L, 5L, 3L, 4L, 
1L, 5L, 3L, 4L, 4L, 4L, 2L, 1L, 2L, 2L, 3L, 4L, 4L, 3L, 2L, 3L, 
3L, 3L, 4L, 3L, 3L, 3L, 4L, 2L, 6L, 6L, 6L, 6L, 5L, 2L, 7L, 3L, 
2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 3L, 3L, 4L, 4L, 3L, 4L, 3L, 3L, 
4L, 5L, 2L, 3L, 3L, 2L, 2L, 2L, 2L, 4L, 3L, 6L, 4L, 4L, 3L, 3L, 
4L, 0L, 2L, 4L, 3L, 2L, 3L, 3L, 0L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 
3L, 2L, 5L, 6L, 3L, 3L, 3L, 3L, 2L, 4L, 3L, 4L, 5L, 4L, 3L, 3L, 
6L, 4L, 3L, 5L, 5L, 5L, 5L, 4L, 3L, 2L, 1L, 2L, 2L, 3L, 2L, 2L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 3L, 2L, 2L, 2L, 3L, 
2L, 2L, 2L, 3L, 2L, 4L, 2L, 3L, 4L, 5L, 5L, 3L, 4L, 3L, 3L, 4L
), ThreeYearRain = c(913.799997, 913.799997, 913.799997, 913.799997, 
913.799997, 913.799997, 913.799997, 913.799997, 938.899988, 938.899988, 
938.899988, 938.600004, 938.899988, 938.899988, 938.899988, 938.499989, 
930.700005, 932.800001, 930.700005, 930.700005, 932.800001, 930.700005, 
930.700005, 932.800001, 932.699991, 934.799987, 934.799987, 934.799987, 
932.699991, 934.799987, 932.699991, 934.799987, 896.99999, 896.99999, 
908.999991, 908.999991, 908.999991, 908.999991, 910.199991, 898.399994, 
928.000009, 939.800006, 935.500004, 928.000009, 928.000009, 923.700007, 
931.499996, 931.499996, 866.200004, 866.200004, 867.000003, 867.000003, 
867.000003, 867.300002, 867.000003, 869.3, 926.800003, 926.800003, 
926.800003, 926.800003, 933.800003, 934.600006, 934.600006, 934.600006, 
924.2, 925.100002, 924.2, 924.2, 924.2, 922.2, 924.2, 924.7, 
974.799995, 983.500006, 983.500006, 983.500004, 983.500006, 974.799995, 
974.799994, 983.500006, 839.1, 839.1, 839.1, 839.100001, 839.300001, 
839.1, 838.699999, 839.100001, 839.100001, 838.699999, 842.300004, 
842.300004, 842.900006, 842.300004, 842.900006, 842.300004, 936.900014, 
936.900014, 936.900014, 932.999984, 933.099983, 932.999984, 936.900014, 
936.900014, 870.499995, 870.499995, 877.399998, 877.399998, 876.099997, 
876.099997, 876.099997, 859.199997, 957.199982, 966.299982, 955.699998, 
955.699998, 957.199982, 955.699998, 955.699998, 956.299985, 852.2, 
852.2, 852.600006, 852.500001, 852.500001, 852.500001, 852.600006, 
852.500001, 906.700011, 904.700001, 912.600007, 912.600007, 914.600007, 
906.700001, 906.399998, 914.600007, 925.599982, 933.299992, 933.299992, 
933.299992, 933.299992, 926.500012, 935.899994, 935.199992, 916.800001, 
916.100001, 916.800001, 916.400003, 918.700003, 904.100001, 916.800001, 
918.700003, 899.1, 904.100001, 906.000003, 903.400002, 904.100001, 
903.400002, 906.000003, 906.000003, 905.7, 903.099999, 903.099999, 
905.7, 912.199994, 893.200002, 905.399999, 904.999998, 933.700012, 
933.700012, 933.700012, 933.700012, 933.700012, 932.30001, 932.300008, 
932.300008, 878.500006, 878.500006, 878.500006, 879.300004, 879.300004, 
879.300004, 879.300004, 873.200008), Div = c(2.2629743, 1.9630117, 
1.7336569, 1.2816843, 2.4155056, 1.4240443, 1.5178948, 0.8993031, 
1.2022801, 1.9287665, 2.0237769, 2.004871, 1.5020684, 2.1776591, 
2.093787, 2.3139276, 2.7244402, 2.7026829, 1.6644725, 2.0696347, 
1.9561853, 2.6018987, 2.5800017, 2.1867866, 2.4144821, 1.7389892, 
2.1427451, 1.6544538, 1.8651966, 1.7569776, 1.8257533, 1.4048204, 
2.7384914, 2.9344488, 2.2306909, 2.5085619, 1.8874836, 2.3431509, 
1.8401602, 1.8620274, 1.8038997, 2.5909049, 2.2265328, 2.0882065, 
2.4737837, 2.2995223, 1.4231311, 2.0577752, 1.6463134, 2.1464331, 
2.2636437, 2.0992589, 1.7666974, 1.835061, 1.7732171, 2.0813243, 
1.865505, 2.0200607, 1.2510612, 1.021761, 0.8111482, 0.2617645, 
2.0282081, 1.1145976, 2.2596683, 2.3517629, 1.9424972, 1.9191269, 
1.4222035, 2.6007698, 2.0071984, 1.9049132, 1.073374, 0.9576897, 
1.6273043, 1.7701581, 0.6890092, 1.5764456, 0.384906, 1.5099996, 
1.6713486, 2.5483064, 2.2033185, 2.0798843, 1.9082985, 2.1580972, 
1.6952798, 1.6303402, 1.9461221, 1.4116405, 1.5347693, 2.6924921, 
1.727278, 1.9384415, 1.6659585, 1.612819, 1.6592884, 2.7129796, 
0, 2.7098898, 1.3785924, 2.7635218, 1.1481271, 1.8597007, 2.2191531, 
1.088549, 2.431015, 1.3702099, 2.1018035, 2.3442348, 2.3599146, 
2.789816, 1.8340235, 1.0606126, 2.5852679, 1.7791063, 1.2273106, 
2.2432636, 2.5642458, 1.3306642, 2.6771856, 1.5062567, 2.0903266, 
2.0398412, 2.4821503, 0.5979376, 1.479214, 1.9188301, 1.2267089, 
2.4491421, 1.5366949, 2.516592, 2.4084849, 2.4385928, 2.549348, 
2.7090074, 2.3337573, 1.8982968, 1.7956341, 2.3752386, 1.6587394, 
2.6663039, 2.4853204, 1.9325793, 2.4431141, 1.6976331, 0.8791745, 
2.6625573, 1.9596877, 1.9287565, 2.4590816, 2.4963942, 1.8767916, 
1.3954333, 2.5155936, 2.2327274, 2.6613726, 2.580748, 2.3142567, 
2.2280879, 1.7925025, 1.663008, 2.3488945, 2.0746398, 1.7050203, 
2.0108246, 1.7317251, 2.4936515, 0.9556999, 1.3716151, 2.0694067, 
1.4944032, 1.0984774, 1.2868726, 1.6429103, 1.3720737, 1.8037795, 
1.8745583, 1.8921264, 1.8320377, 1.201682, 1.8489571, 1.798546, 
0.8486856), triodia_low = c(19L, 6L, 21L, 32L, 11L, 32L, 16L, 
29L, 17L, 20L, 0L, 24L, 37L, 0L, 3L, 29L, 4L, 2L, 31L, 28L, 20L, 
12L, 6L, 6L, 26L, 28L, 27L, 32L, 37L, 26L, 15L, 27L, 2L, 1L, 
19L, 5L, 13L, 10L, 33L, 14L, 25L, 22L, 15L, 34L, 15L, 7L, 36L, 
25L, 25L, 0L, 25L, 4L, 21L, 0L, 33L, 16L, 16L, 15L, 22L, 25L, 
25L, 0L, 0L, 18L, 2L, 0L, 26L, 0L, 0L, 7L, 0L, 13L, 28L, 35L, 
0L, 0L, 0L, 31L, 29L, 0L, 14L, 5L, 14L, 11L, 12L, 16L, 21L, 26L, 
22L, 7L, 23L, 10L, 23L, 17L, 19L, 7L, 27L, 3L, 0L, 2L, 29L, 14L, 
30L, 12L, 0L, 35L, 0L, 29L, 4L, 5L, 14L, 15L, 33L, 0L, 3L, 21L, 
34L, 0L, 2L, 28L, 16L, 0L, 1L, 0L, 11L, 0L, 32L, 0L, 27L, 2L, 
28L, 3L, 0L, 4L, 1L, 6L, 14L, 27L, 25L, 12L, 7L, 10L, 16L, 9L, 
4L, 15L, 40L, 2L, 18L, 5L, 3L, 6L, 1L, 33L, 2L, 5L, 12L, 4L, 
7L, 3L, 17L, 30L, 5L, 7L, 17L, 15L, 16L, 9L, 0L, 26L, 16L, 0L, 
24L, 1L, 27L, 32L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), class = "data.frame", row.names = c(NA, 
-184L))

Модель: 

m1b <-gamm4(Div~TSF+FF+s(triodia_low, k=6)+VegtypeNew+ThreeYearRain,random=~(1|Lscape),data=data)

Печать: 

p <-plotGAMM (m1b,smooth.cov="triodia_low",groupCovs = NULL,orderedAsFactor=T,rawOrFitted="raw",plotCI=T,grouping = NULL)
p + labs(x= "Years since fire") +
  labs(y="Species diversity (H')") +
  theme_bw() +
  theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor = element_blank()) +
  theme(panel.border = element_blank()) +
  theme(axis.line = element_line(colour="black")) +
  theme(axis.title = element_text(size=22)) +    # face="bold"
  theme(axis.ticks = element_line()) +
  scale_x_continuous(breaks = seq(0,40,by=5)) +
  geom_rect(aes(xmin=0,xmax=0.6,ymin=-Inf,ymax=Inf),alpha=0.002, fill="coral4") +
  geom_rect(aes(xmin=.6,xmax=1,ymin=-Inf,ymax=Inf),alpha=0.002, fill="gold") +
  geom_rect(aes(xmin=1,xmax=5,ymin=-Inf,ymax=Inf),alpha=0.002, fill="darkred") +
  geom_rect(aes(xmin=5,xmax=10,ymin=-Inf,ymax=Inf),alpha=0.002,fill="chocolate") +
  geom_rect(aes(xmin=10,xmax=40,ymin=-Inf,ymax=Inf),alpha=0.002,fill="orangered") +
  theme(axis.text = element_text(size=18, colour="black")) +
  theme(text = element_text(family = "Arial")) +
  theme(legend.position= "none")

Сюжет:

enter image description here

Любая помощь будет принята с благодарностью:)

1 Ответ

0 голосов
/ 07 апреля 2020

Используя версию для разработки gratia , вы можете воспроизвести показанный вами график с помощью нескольких простых вызовов для создания данных, предсказать et c.

Чтобы установить версию для разработки gratia do 

# install.packages("remotes")
remotes::install_github("gavinsimpson/gratia")

После установки вы можете создать объект, подходящий для черчения, используя: 

library('mgcv')
library('gamm4')
library('gratia')
library('ggplot2')
library('dplyr')

## model fit
m1b <- gamm4(Div ~ TSF + FF + s(triodia_low, k=6) + VegtypeNew + ThreeYearRain,
             random = ~ (1|Lscape), data = df)

## data to predict at
new_df <- data_slice(m1b, var1 = 'triodia_low', n = 100)
## predict and cast to a tibble
pred_df <- as_tibble(predict(m1b[["gam"]], new_df, se.fit = TRUE))
## add to the data we're predicting at
pred_df <- bind_cols(new_df, pred_df)
## grab the inverse link of the model (not needed here, but is for non-Normal fits)
ilink <- inv_link(m1b)
## create the upper and lower credible interval
pred_df <- mutate(pred_df,
                  lwr = ilink(fit - (2 * se.fit)),
                  upr = ilink(fit + (2 * se.fit)),
                  fit = ilink(fit))

Сам сюжет может быть создан с помощью: 

ggplot(pred_df, aes(x = triodia_low, y = fit)) +
    labs(x = "Years since fire", y = "Species diversity (H')") +
    theme_bw() +
    theme(panel.grid.major = element_blank()) +
    theme(panel.grid.minor = element_blank()) +
    theme(panel.border = element_blank()) +
    theme(axis.line = element_line(colour="black")) +
    theme(axis.title = element_text(size=22)) +
    theme(axis.ticks = element_line()) +
    scale_x_continuous(breaks = seq(0,40,by=5)) +
    geom_rect(aes(xmin=0,   xmax=0.6, ymin=-Inf, ymax=Inf), alpha=0.01, fill="coral4") +
    geom_rect(aes(xmin=0.6, xmax=1,   ymin=-Inf, ymax=Inf), alpha=0.01, fill="gold") +
    geom_rect(aes(xmin=1,   xmax=5,   ymin=-Inf, ymax=Inf), alpha=0.01, fill="darkred") +
    geom_rect(aes(xmin=5,   xmax=10,  ymin=-Inf, ymax=Inf), alpha=0.01, fill="chocolate") +
    geom_rect(aes(xmin=10,  xmax=40,  ymin=-Inf, ymax=Inf), alpha=0.01, fill="orangered") +
    geom_point(data = df, mapping = aes(x = triodia_low, y = Div)) +
    geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.4) +
    geom_line()

Большая часть кода ggplot принадлежит вам, но теперь, когда у нас есть полный контроль, мы можем поместить слои данных на передний план, оставив эти слои до конца.

Обратите внимание, что я не уверен, что это отличный сюжет. Отображаемая функция, которую вы показываете, зависит от других ковариат в наборе данных. Здесь, как и в случае voxel::plotGAMM(), мы прогнозируем по модели и, следовательно, мы должны предоставить что-то для других ковариат. Следуя voxel::plotGAMM и mgcv::vis.gam, мы фиксируем другие ковариаты, не показанные на

  1. значении данных наблюдения, ближайшем к медиане (для непрерывных переменных), или
  2. модальная категория (для множителя параметри c слагаемых)

Итак, результирующая цифра соответствует условию этих значений. В частности это для уровня spinsandplain VegTypeNew. Таким образом, это немного вводит в заблуждение.

...