Как я уже упоминал в комментариях выше, вам действительно лучше делать это за пределами ggplot и вместо этого делать это с полной моделью двух сглаживаний, из которых вы можете вычислить неопределенности разницы и т. Д..
Это, в основном, короткая версия сообщения в блоге , которое я написал год или около того.
Примерные данные ОП
set.seed(1) # make example reproducible
n <- 5000 # set sample size
df <- data.frame(x= rnorm(n), g= factor(rep(c(0,1), n/2))) # generate data
df$y <- NA # include y in df
df$y[df$g== 0] <- df$x[df$g== 0]**2 + rnorm(sum(df$g== 0))*5 # y for group g= 0
df$y[df$g== 1] <-2 + df$x[df$g== 1]**2 + rnorm(sum(df$g== 1))*5 # y for g= 1 (with intercept 2)
Начните с подгонки модели для данных примера:
library("mgcv")
m <- gam(y ~ g + s(x, by = g), data = df, method = "REML")
Здесь я подгоняю GAM с гладко-множительным взаимодействием (бит by), и для этой модели нам нужно также включить g как параметрический эффект, поскольку групповые сглаживания центрированы около 0, поэтому нам нужно включить групповые средние значения в параметрическую часть модели.
Далее нам нужна сетка данных вдоль xпеременная, по которой мы будем оценивать разницу между двумя оцененными сглаживаниями:
pdat <- with(df, expand.grid(x = seq(min(x), max(x), length = 200),
g = c(0,1)))
pdat <- transform(pdat, g = factor(g))
, затем мы используем эти данные прогнозирования для генерации матрицы Xp, которая представляет собой матрицу, отображающую значения ковариатаes к значениям базисного расширения для сглаживаний;мы можем манипулировать этой матрицей, чтобы получить сглаживание различий, которое мы хотим:
xp <- predict(m, newdata = pdat, type = "lpmatrix")
Следующий код, чтобы определить, какие строки и столбцы в xp принадлежат сглаживаниям для соответствующих уровней g;поскольку в модели есть только два уровня и только один сглаженный член, это совершенно тривиально, но для более сложных моделей это необходимо, и важно правильно настроить имена компонентов, чтобы биты grep() работали.
## which cols of xp relate to splines of interest?
c1 <- grepl('g0', colnames(xp))
c2 <- grepl('g1', colnames(xp))
## which rows of xp relate to sites of interest?
r1 <- with(pdat, g == 0)
r2 <- with(pdat, g == 1)
Теперь мы можем различать строки xp для пары уровней, которые мы сравниваем
## difference rows of xp for data from comparison
X <- xp[r1, ] - xp[r2, ]
Поскольку мы фокусируемся на разнице, нам нужно обнулить все столбцы, а несвязан с выбранной парой сглаживаний, которая включает в себя любые параметрические члены.
## zero out cols of X related to splines for other lochs
X[, ! (c1 | c2)] <- 0
## zero out the parametric cols
X[, !grepl('^s\\(', colnames(xp))] <- 0
(В этом примере эти две строки делают одно и то же, но в более сложных примерах необходимы оба.)
Теперь у нас есть матрица X, которая содержит разницу между двумя базисными разложениями для пары сглаживаний, которые нас интересуют, но чтобы получить это в терминах подгоночных значений ответа y, нам нужно умножитьэта матрица по вектору коэффициентов:
## difference between smooths
dif <- X %*% coef(m)
Теперь dif содержит разницу между двумя сглаживаниями.
We cснова использовать X и ковариационную матрицу модельных коэффициентов для вычисления стандартной ошибки этой разности, а затем 95% (в данном случае) доверительного интервала для разности оценок.
## se of difference
se <- sqrt(rowSums((X %*% vcov(m)) * X))
## confidence interval on difference
crit <- qt(.975, df.residual(m))
upr <- dif + (crit * se)
lwr <- dif - (crit * se)
Обратите внимание, что здесьпри вызове vcov() мы используем эмпирическую байесовскую ковариационную матрицу, но не ту, которая была исправлена для выбора параметров гладкости.Функция, которую я покажу в ближайшее время, позволяет вам учитывать эту дополнительную неопределенность с помощью аргумента unconditional = TRUE.
Наконец, мы собираем результаты и график:
res <- data.frame(x = with(df, seq(min(x), max(x), length = 200)),
dif = dif, upr = upr, lwr = lwr)
ggplot(res, aes(x = x, y = dif)) +
geom_ribbon(aes(ymin = lwr, ymax = upr, x = x), alpha = 0.2) +
geom_line()
Это дает
Что согласуется с оценкой, показывающей, что модель со сглаживанием на уровне группы не обеспечивает существенно лучшего соответствия, чем модель с другими групповыми средствами, но только один общий сглаживательв x:
r$> m0 <- gam(y ~ g + s(x), data = df, method = "REML")
r$> AIC(m0, m)
df AIC
m0 9.68355 30277.93
m 14.70675 30285.02
r$> anova(m0, m, test = 'F')
Analysis of Deviance Table
Model 1: y ~ g + s(x)
Model 2: y ~ g + s(x, by = g)
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 4990.1 124372
2 4983.9 124298 6.1762 73.591 0.4781 0.8301
Заключение
В сообщении, которое я упомянул, есть функция, которая превращает описанные выше шаги в простую функцию, smooth_diff():
smooth_diff <- function(model, newdata, f1, f2, var, alpha = 0.05,
unconditional = FALSE) {
xp <- predict(model, newdata = newdata, type = 'lpmatrix')
c1 <- grepl(f1, colnames(xp))
c2 <- grepl(f2, colnames(xp))
r1 <- newdata[[var]] == f1
r2 <- newdata[[var]] == f2
## difference rows of xp for data from comparison
X <- xp[r1, ] - xp[r2, ]
## zero out cols of X related to splines for other lochs
X[, ! (c1 | c2)] <- 0
## zero out the parametric cols
X[, !grepl('^s\\(', colnames(xp))] <- 0
dif <- X %*% coef(model)
se <- sqrt(rowSums((X %*% vcov(model, unconditional = unconditional)) * X))
crit <- qt(alpha/2, df.residual(model), lower.tail = FALSE)
upr <- dif + (crit * se)
lwr <- dif - (crit * se)
data.frame(pair = paste(f1, f2, sep = '-'),
diff = dif,
se = se,
upper = upr,
lower = lwr)
}
Используя эту функцию, мы можем повторить весь анализ и отобразить разницу следующим образом:
out <- smooth_diff(m, pdat, '0', '1', 'g')
out <- cbind(x = with(df, seq(min(x), max(x), length = 200)),
out)
ggplot(out, aes(x = x, y = diff)) +
geom_ribbon(aes(ymin = lower, ymax = upper, x = x), alpha = 0.2) +
geom_line()
Я не буду показывать график здесь, поскольку он идентичен показанному выше, за исключением меток оси.