stat_smooth gam не то же самое, что gam {mgcv}

Я использовал функцию stat_smooth в ggplot2, решил, что мне нужна «степень соответствия», и применил для этого gammgvc. Мне пришло в голову, что я должен проверить, что это одна и та же модель (stat_smooth vs mgvc gam), поэтому я использовал приведенный ниже код для проверки. Казалось бы, у них разные результаты, о чем свидетельствует график (График: stat_smoother gam (красный), mgcv gam (черный)). Однако я не знаю, почему у них разные результаты. Дело в том, что какой-то параметр по умолчанию у них отличается? Это что игра запускается с числовым x, а stat_smooth запускается с POSIXct x (если да - я не знаю, что с этим делать)? Похоже, что stat_smooth более плавный, но значения k такие же ...

Я думаю, что есть несколько сообщений о том, как строить гамма-выходы в ggplot2, но я действительно хотел бы знать, почему stat_smooth и mgcv вообще дают разные результаты. Я новичок в GAM (и R), поэтому вполне возможно, что мне не хватает чего-то простого. Однако перед тем, как спросить, я сделал поиск по этому форуму в Google.

Мои данные немного велики, чтобы их можно было легко передать, поэтому я использовал образец набора данных - я поместил исходный код в код, а также dput() под всем и мой sessionInfo() после этого.

Я пытался задать качественный вопрос, но это только второй. Всегда. Так что конструктивная критика приветствуется.

Спасибо!

library(readxl)
library(data.table)
library(ggplot2)
library(scales)
library(mgcv)

stackOF_data <- read_excel("mean-daily-flow-cumecs-vatnsdals.xlsx", sheet = "Data")
stackOF_data <- data.table(stackOF_data)
stackOF_data <- stackOF_data[,.(timeseries=as.POSIXct(Date,format("%Y-%m-%d")),mdf)]

a <- stackOF_data[,.(x=as.numeric(timeseries),y=mdf)]
a1 <- gam(y~s(x, k=100, bs="cs"),data=a) 
a2=data.table(gam_mdf= predict(a1,a))
a2=cbind(timeseries=stackOF_data$timeseries,a2)

# see if predict and actual are the same
p <- ggplot() + 
geom_line(data = a2, aes(x = timeseries, y = gam_mdf), size=1)+
scale_color_manual(values=c("black","magenta"))+
scale_y_continuous()+
scale_x_datetime(labels = date_format("%Y-%m-%d"), breaks = "1 month", minor_breaks = "1 week")+ 
theme(axis.text.x=element_text(angle=50, size=10,hjust=1))+
stat_smooth(data = stackOF_data, aes(x = (timeseries), y = mdf),method="gam", formula=y~s(x,k=100, bs="cs"), col="red", se=FALSE, size=1)
p

# data from: https://datamarket.com/data/set/235m/mean-daily-flow-cumecs-vatnsdalsa-river-1-jan-1972-31-dec-1974#!ds=235m&display=line&s=14l

> dput(a)
structure(list(x = c(126230400, 126316800, 126403200, 126489600, 
126576000, 126662400, 126748800, 126835200, 126921600, 127008000, 
127094400, 127180800, 127267200, 127353600, 127440000, 127526400, 
127612800, 127699200, 127785600, 127872000, 127958400, 128044800, 
128131200, 128217600, 128304000, 128390400, 128476800, 128563200, 
128649600, 128736000, 128822400, 128908800, 128995200, 129081600, 
129168000, 129254400, 129340800, 129427200, 129513600, 129600000, 
129686400, 129772800, 129859200, 129945600, 130032000, 130118400, 
130204800, 130291200, 130377600, 130464000, 130550400, 130636800, 
130723200, 130809600, 130896000, 130982400, 131068800, 131155200, 
131241600, 131328000, 131414400, 131500800, 131587200, 131673600, 
131760000, 131846400, 131932800, 132019200, 132105600, 132192000, 
132278400, 132364800, 132451200, 132537600, 132624000, 132710400, 
132796800, 132883200, 132969600, 133056000, 133142400, 133228800, 
133315200, 133401600, 133488000, 133574400, 133660800, 133747200, 
133833600, 133920000, 134006400, 134092800, 134179200, 134265600, 
134352000, 134438400, 134524800, 134611200, 134697600, 134784000, 
134870400, 134956800, 135043200, 135129600, 135216000, 135302400, 
135388800, 135475200, 135561600, 135648000, 135734400, 135820800, 
135907200, 135993600, 136080000, 136166400, 136252800, 136339200, 
136425600, 136512000, 136598400, 136684800, 136771200, 136857600, 
136944000, 137030400, 137116800, 137203200, 137289600, 137376000, 
137462400, 137548800, 137635200, 137721600, 137808000, 137894400, 
137980800, 138067200, 138153600, 138240000, 138326400, 138412800, 
138499200, 138585600, 138672000, 138758400, 138844800, 138931200, 
139017600, 139104000, 139190400, 139276800, 139363200, 139449600, 
139536000, 139622400, 139708800, 139795200, 139881600, 139968000, 
140054400, 140140800, 140227200, 140313600, 140400000, 140486400, 
140572800, 140659200, 140745600, 140832000, 140918400, 141004800, 
141091200, 141177600, 141264000, 141350400, 141436800, 141523200, 
141609600, 141696000, 141782400, 141868800, 141955200, 142041600, 
142128000, 142214400, 142300800, 142387200, 142473600, 142560000, 
142646400, 142732800, 142819200, 142905600, 142992000, 143078400, 
143164800, 143251200, 143337600, 143424000, 143510400, 143596800, 
143683200, 143769600, 143856000, 143942400, 144028800, 144115200, 
144201600, 144288000, 144374400, 144460800, 144547200, 144633600, 
144720000, 144806400, 144892800, 144979200, 145065600, 145152000, 
145238400, 145324800, 145411200, 145497600, 145584000, 145670400, 
145756800, 145843200, 145929600, 146016000, 146102400, 146188800, 
146275200, 146361600, 146448000, 146534400, 146620800, 146707200, 
146793600, 146880000, 146966400, 147052800, 147139200, 147225600, 
147312000, 147398400, 147484800, 147571200, 147657600, 147744000, 
147830400, 147916800, 148003200, 148089600, 148176000, 148262400, 
148348800, 148435200, 148521600, 148608000, 148694400, 148780800, 
148867200, 148953600, 149040000, 149126400, 149212800, 149299200, 
149385600, 149472000, 149558400, 149644800, 149731200, 149817600, 
149904000, 149990400, 150076800, 150163200, 150249600, 150336000, 
150422400, 150508800, 150595200, 150681600, 150768000, 150854400, 
150940800, 151027200, 151113600, 151200000, 151286400, 151372800, 
151459200, 151545600, 151632000, 151718400, 151804800, 151891200, 
151977600, 152064000, 152150400, 152236800, 152323200, 152409600, 
152496000, 152582400, 152668800, 152755200, 152841600, 152928000, 
153014400, 153100800, 153187200, 153273600, 153360000, 153446400, 
153532800, 153619200, 153705600, 153792000, 153878400, 153964800, 
154051200, 154137600, 154224000, 154310400, 154396800, 154483200, 
154569600, 154656000, 154742400, 154828800, 154915200, 155001600, 
155088000, 155174400, 155260800, 155347200, 155433600, 155520000, 
155606400, 155692800, 155779200, 155865600, 155952000, 156038400, 
156124800, 156211200, 156297600, 156384000, 156470400, 156556800, 
156643200, 156729600, 156816000, 156902400, 156988800, 157075200, 
157161600, 157248000, 157334400, 157420800, 157507200, 157593600, 
157680000), y = c(4.65, 4.65, 4.65, 4.48, 5.16, 5.52, 5.34, 5.34, 
4.82, 4.65, 4.48, 4.31, 4.31, 4.31, 4.14, 3.82, 3.98, 3.98, 4.31, 
5.71, 6.5, 6.3, 5.71, 5.71, 5.16, 4.65, 4.14, 3.98, 4.48, 4.48, 
4.31, 4.65, 4.31, 3.98, 3.98, 3.98, 3.98, 3.98, 3.98, 3.82, 3.67, 
3.67, 3.98, 3.98, 3.82, 3.82, 3.82, 4.14, 5.9, 4.48, 3.98, 3.98, 
3.82, 3.67, 3.67, 3.67, 4.65, 3.98, 4.31, 4.31, 3.67, 4.31, 6.1, 
7.3, 7.5, 7.5, 8.14, 10.8, 16.1, 14.8, 12.5, 9.9, 8.14, 6.9, 
6.1, 5.34, 5.16, 4.99, 4.99, 4.99, 4.99, 5.52, 6.3, 7.3, 6.9, 
5.9, 5.71, 5.71, 8.58, 31.5, 33.7, 18.4, 11.3, 16.1, 32.9, 45.3, 
54, 25.7, 18, 15.9, 15.6, 14.5, 15.9, 35.9, 37.5, 29.4, 27.5, 
30.1, 27.5, 30.8, 29.4, 22, 20.1, 35.9, 36.7, 32.9, 22, 18, 15.9, 
15.2, 14.8, 13, 12.7, 12.5, 11, 9.68, 8.8, 7.92, 7.3, 6.9, 7.3, 
10.3, 11, 11.3, 11.9, 12.5, 13.6, 12.2, 10.8, 9.9, 9.46, 8.8, 
7.5, 7.1, 7.71, 7.1, 6.1, 5.34, 5.34, 5.34, 5.52, 5.52, 6.3, 
6.7, 6.5, 5.9, 5.71, 5.9, 5.71, 5.52, 7.3, 7.5, 7.1, 7.3, 6.7, 
6.9, 7.3, 7.5, 10.8, 11.6, 8.58, 7.92, 7.1, 6.7, 6.5, 6.1, 5.9, 
5.9, 5.71, 5.52, 5.52, 5.52, 5.9, 5.9, 5.71, 5.52, 5.52, 5.34, 
5.34, 5.52, 6.5, 6.5, 5.71, 5.34, 5.16, 4.99, 4.82, 4.82, 4.99, 
4.82, 4.82, 4.82, 4.82, 4.82, 4.65, 4.48, 4.48, 4.31, 4.31, 4.14, 
4.14, 4.31, 4.48, 4.31, 4.31, 4.31, 4.99, 5.71, 6.3, 6.1, 6.1, 
5.9, 5.71, 5.52, 5.52, 5.52, 5.52, 5.52, 5.34, 5.34, 5.52, 5.52, 
5.52, 5.34, 5.34, 5.52, 5.34, 5.52, 5.52, 5.34, 5.34, 5.34, 5.34, 
5.71, 5.9, 6.3, 6.9, 7.5, 6.5, 6.1, 6.1, 5.9, 6.1, 6.1, 5.9, 
6.5, 6.5, 6.1, 5.9, 5.9, 5.71, 5.9, 5.9, 5.71, 4.99, 4.65, 5.16, 
5.34, 5.34, 4.65, 4.99, 5.71, 5.34, 5.34, 5.34, 5.34, 4.99, 5.34, 
5.34, 5.34, 5.34, 5.52, 5.34, 5.52, 5.71, 6.5, 7.71, 6.9, 6.5, 
6.7, 6.1, 5.9, 6.1, 5.9, 5.71, 7.92, 7.71, 7.1, 7.92, 5.34, 5.16, 
8.14, 10.1, 7.92, 7.3, 6.9, 6.9, 6.9, 8.58, 7.3, 6.9, 7.3, 6.3, 
5.16, 6.1, 5.52, 4.99, 5.34, 5.34, 5.34, 5.16, 5.71, 5.52, 5.52, 
5.16, 4.82, 5.52, 6.1, 5.9, 5.71, 5.52, 5.16, 4.99, 4.48, 4.82, 
5.16, 5.16, 5.16, 5.16, 5.16, 4.82, 4.65, 3.82, 4.14, 4.65, 4.65, 
4.31, 4.31, 5.16, 5.16, 5.16, 5.16, 5.16, 4.99, 4.65, 5.16, 5.16, 
5.16, 5.16, 5.16, 5.16, 5.16, 5.16, 5.34, 5.34)), .Names = c("x", 
"y"), row.names = c(NA, -365L), class = c("data.table", "data.frame"
), .internal.selfref = <pointer: 0x0000000005860788>)

> sessionInfo()
R version 3.2.2 (2015-08-14)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United         States.1252    LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                           LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] data.table_1.9.6 readxl_0.1.0     mgcv_1.8-7       nlme_3.1-121         scales_0.3.0     sos_1.3-8        brew_1.0-6       ggplot2_1.0.1   
[9] MASS_7.3-43     

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.1      lattice_0.20-33  digest_0.6.8     chron_2.3-47     grid_3.2.2       plyr_1.8.3       gtable_0.1.2     magrittr_1.5    
 [9] stringi_0.5-5    reshape2_1.4.1   Matrix_1.2-2     labeling_0.3     proto_0.3-10     tools_3.2.2      stringr_1.0.0    munsell_0.4.2   
[17] colorspace_1.2-6

Частичное решение

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

library(readxl)
library(data.table)
library(ggplot2)
library(scales)
library(mgcv)

stackOF_data <- read_excel("C:/Users/jel4049/Desktop/mean-daily-flow-cumecs-    vatnsdals.xlsx", sheet = "Data")
stackOF_data <- data.table(stackOF_data)

stackOF_data <- stackOF_data[,.(timeseries=as.POSIXct(Date,format("%Y-%m-%d")),mdf)]
a <- stackOF_data[,.(x=as.numeric(timeseries),y=mdf)]
a1 <- gam(y~s(x, k=100, bs="cs"),data=a) 
a2=data.table(gam_mdf = predict(a1,a))

preds <- predict(a1,se.fit=TRUE)
my_data <- data.frame(mu=preds$fit, low =(preds$fit - 1.96 * preds$se.fit), high = (preds$fit + 1.96 * preds$se.fit))


m <- ggplot()+
  geom_line(data = a2, aes(x=stackOF_data$timeseries, y=gam_mdf), size=1, col="blue")+
  geom_smooth(data=my_data,aes(ymin = low, ymax = high, x=stackOF_data$timeseries, y = mu), stat = "identity", col="green")
m

Теперь, по крайней мере, я знаю, что сводная информация и информация о качестве соответствия данных, которую я могу получить с помощью некоторых функций mgcv, соответствуют моим графикам.


person CJ9    schedule 09.10.2015    source источник
comment
@aosmith Спасибо за ответ. Выдает ошибку Error: Invalid input: time_trans works with objects of class POSIXct only Кажется, разница в требованиях gam и stat_smooth с точки зрения того, какую информацию о дате / времени они могут принять. Я не уверен, как сделать их такими же, и при этом у меня все еще есть мой график с читаемой осью x .... Рад помощи в этом, если у вас есть какие-то идеи! Спасибо.   -  person CJ9    schedule 12.10.2015
comment
@aosmith Я над этим работал. Я пытался:   -  person CJ9    schedule 14.10.2015
comment
@aosmith Я над этим работал. Я пробовал: w <- ggplot() + geom_line(data = a2, aes(x=a$x, y=gam_mdf), size=1, col="blue")+ geom_smooth(data = a, aes(x=x, y=y),method="gam", formula=y~s(x,k=100, bs="cs"), col="green",se=FALSE, size=1) w Это должно быть больше похоже на то, о чем вы думали, верно? Сюжет по-прежнему показывает разницу. Любые идеи? Спасибо.   -  person CJ9    schedule 14.10.2015
comment
@aosmith Спасибо, что заглянули еще раз. Я все еще хотел бы знать, почему они не совпадают. Однако я нашел работу, которую включил в свой ответ выше.   -  person CJ9    schedule 15.10.2015


Ответы (1)


Оказывается, различия, которые вы наблюдали, связаны с тем, что вы использовали значение по умолчанию для аргумента n в stat_smooth.

Со страницы справки:

n количество точек для оценки более плавного

Конечно, мне сразу не пришло в голову, что это означает, что n контролирует размер набора данных для newdata аргумента в predict и, следовательно, stat_smooth не использует исходный набор данных при прогнозировании. Но я читал этот хороший ответ на другой stat_smooth вопрос и понял, что чтобы понять, что происходит, мне нужно внимательнее посмотрите на stat_smooth прогнозы и прогнозы, сделанные вручную на основе подобранной gam модели.

Итак, используя набор данных из OP, который я назвал dat, мы можем проверить, что происходит.

График, когда k = 100, после подгонки модели через gam и добавления прогнозов в набор данных. Как вы заметили, синий (stat_smooth) и черный (прогнозы вручную) не совпадают.

dat$predgam = predict(gam(y ~ s(x, k = 100), data = dat))

(p1 = ggplot(dat, aes(x, y)) +
    geom_point() +
    geom_smooth(method = "gam", formula = y ~ s(x, k = 100)) +
    geom_line(aes(y = predgam)))

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

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

ggplot_build(p1)

Набор данных прогноза для stat_smooth является вторым в списке наборов данных.

ggplot_build(p1)$data[[2]]

Но посмотрите, сколько строк в этом наборе данных:

nrow(ggplot_build(p1)$data[[2]])
[1] 80

Значение по умолчанию для аргумента n - 80, но в вашем наборе данных 365 строк. Так что же произойдет, если вы измените n на 365? Я сделаю плавную линию более толстой, чтобы вы действительно могли ее видеть (синим цветом).

(p2 = ggplot(dat, aes(x, y)) +
    geom_point() +
    geom_smooth(method = "gam", formula = y ~ s(x, k = 100), n = 365, size = 2) +
    geom_line(aes(y = predgam)))

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

nrow(ggplot_build(p2)$data[[2]])
[1] 365

Если вы посмотрите на код функции predictdf, упомянутой в разделе «Подробности» на stat_smooth странице справки, вы увидите, что исходный набор данных не используется при построении прогнозов. Вместо этого последовательность составляется из исходной объясняющей переменной. Это может быть действительно важно при работе с небольшим набором данных, и вам нужно больше точек прогноза, чтобы линия выглядела гладкой. Однако в вашем случае исходный набор данных уже представляет собой красивую гладкую последовательность x, поэтому использование n = 365 дает те же прогнозы от stat_smooth, что и исходный набор данных.

Вы можете увидеть код для predictdf здесь.

person aosmith    schedule 16.10.2015
comment
Что здесь означает «k»: y ~ s (x, k = 100)? - person Kshitij15571; 11.11.2020
comment
@ Kshitij15571 См. Документацию для s(), чтобы увидеть описания аргументов, ?mgcv::s. Описание аргумента начинается с измерения основы, используемой для представления гладкого члена. - person aosmith; 11.11.2020
comment
Я считаю, что проблема может заключаться в построении графика в plot.gam() или plot(model) . Частичное решение @ CJ9 исправляет это, иначе кривая не будет хорошо соответствовать данным. См. Здесь тот же результат с использованием решения OP stackoverflow.com/a/67096214/2291642. Изменение n не повлияло на результат, по крайней мере, в моем случае. - person Bob; 14.04.2021