Прогнозы на графике с одновременным интервалом от игрового диапазона в диапазоне сглаженной переменной

Question

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

library(mgcv)

mod <- gam(dv_value ~ age_grps + period.f + s(born_adult), data = dat, contrasts = list(age_grps = contr.sum, period.f = contr.sum))

Сначала я вычисляю предсказанные значения по всему диапазону переменной born_adult с одновременныминтервал, который, кажется, работает довольно хорошо:

rmvn <- function(n, mu, sig) { 
  L <- mroot(sig)
  m <- ncol(L)
  t(mu + L %*% matrix(rnorm(m*n), m, n))
}

Vb <- vcov(mod)

pred <- predict(mod, se.fit = TRUE)

se.fit <- pred$se.fit

N <- 10000

BUdiff <- rmvn(N, mu = rep(0, nrow(Vb)), sig = Vb)

Cg <- predict(mod, type = "lpmatrix")
simDev <- Cg %*% t(BUdiff)

absDev <- abs(sweep(simDev, 1, se.fit, FUN = "/"))

masd <- apply(absDev, 2L, max)

crit <- quantile(masd, prob = 0.95, type = 8)

predData <- transform(cbind(data.frame(pred), dat),
                      uprP = fit + (crit * se.fit),
                      lwrP = fit - (crit * se.fit))

Однако, пытаясь построить результаты, я получаю действительно странный график:

ggplot() +
  geom_ribbon(aes(x = born_adult, ymin = lwrP, ymax = uprP), data = predData, alpha = 0.2, fill = "red")

https://www.dropbox.com/s/uskj9oyq8ud3zx2/plot1.png?dl=0

Но, при обработке моих управляющих переменных, я получаю правильные прогнозы для отдельных «срезов» моих данных:

ggplot() +
  geom_ribbon(aes(x = born_adult, ymin = lwrP, ymax = uprP), data = predData, alpha = 0.2, fill = "red") + 
  facet_wrap(vars(period.f, age_grps))

https://www.dropbox.com/s/yju68yl8kes8mp1/plot2.png?dl=0

Я также попытался предсказатьНовый смоделированный набор данных использует ту же структуру, что и мои данные, однако проблема осталась прежней. Есть ли возможность показать «средние» прогнозы по всему диапазону моей независимой сглаженной переменной, не обращая внимания на контрольные переменные? Я полагаю, что это может сработать, взяв средние прогнозы, сгруппированные по значениям переменной born_adult: predData <- group_by(born_adult) %>% summarize(fit = mean(fit)) Однако я не знаю, как взять среднее значение одновременных интервалов для единичных прогнозов.

И последнее, но не менее важное: вот небольшое подмножество данных, которые я использую:

dat <- structure(list(dv_value = c(0.8, 0.8, 0.4, 0.8, 1, 0.6, 0.6, 
1, 0.8, 1, 1, 1, 1, 0.4, 0.8, 0.8, 1, 0.4, 1, 0.6, 1, 0.8, 0.6, 
0, 0.6, 0.8, 0.8, 1, 0.8, 0.8, 0.8, 1, 1, 1, 0.8, 1, 0.6, 1, 
0.6, 0.8, 0.8, 0.8, 0.6, 1, 1, 1, 0.6, 1, 1, 1, 0.8, 1, 0.6, 
0.6, 1, 1, 0.8, 0.6, 0.8, 0.6, 1, 0.8, 0.8, 0.6, 0.8, 0.8, 1, 
1, 0.8, 0.8, 0.8, 1, 1, 0.6, 1, 1, 1, 1, 1, 1, 0.6, 0.8, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 0.6, 1, 0.6, 0.6, 0.6, 0.8, 0.8, 0.8, 0.8, 
1, 0.4, 0.8, 1, 1, 1, 1, 0.4, 1, 1, 0.6, 1, 1, 0.4, 0.6, 0.8, 
1, 1, 0.6, 1, 1, 0.6, 1, 0.8, 0.8, 1, 0.8, 1, 0.8, 1, 0.6, 0.8, 
1, 0.8, 0.6, 0.6, 1, 0.8, 0.6, 1, 0.6, 1, 0.6, 0.8, 1, 0.6, 1, 
0.8, 0.8, 0.8, 1, 1, 1, 1, 0.2, 1, 0.6, 1, 0.8, 0.8, 1, 0.6, 
1, 0.4, 1, 0.8, 0.8, 0.4, 1, 1, 0.8, 0.8, 0.8, 1, 0.8, 0.6, 0.6, 
0.4, 0.2, 1, 0.8, 0.4, 1, 1, 0.8, 1, 0.8, 0.6, 1, 1, 1, 0.8, 
1, 0.6, 0.8, 0.8, 1, 1, 0.8, 1), age_grps = structure(c(1L, 3L, 
3L, 3L, 1L, 2L, 3L, 3L, 2L, 2L, 2L, 3L, 2L, 2L, 3L, 3L, 3L, 2L, 
3L, 3L, 3L, 2L, 3L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 3L, 2L, 3L, 3L, 
2L, 1L, 2L, 2L, 2L, 3L, 3L, 2L, 1L, 2L, 3L, 1L, 2L, 3L, 2L, 3L, 
2L, 2L, 3L, 3L, 2L, 3L, 3L, 2L, 3L, 1L, 3L, 1L, 2L, 3L, 2L, 3L, 
2L, 3L, 2L, 3L, 3L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 3L, 2L, 2L, 1L, 
2L, 3L, 2L, 3L, 3L, 2L, 3L, 1L, 3L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 
3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 3L, 2L, 2L, 3L, 3L, 2L, 2L, 
3L, 3L, 1L, 2L, 1L, 3L, 2L, 2L, 3L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 
3L, 3L, 3L, 2L, 3L, 2L, 3L, 3L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 
2L, 2L, 3L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 3L, 1L, 1L, 
2L, 3L, 2L, 3L, 3L, 3L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 
2L, 1L, 3L, 3L, 3L, 3L, 2L, 1L, 2L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 
2L, 3L, 2L, 3L, 3L, 2L), .Label = c("1", "2", "3"), class = "factor"), 
    period.f = structure(c(9L, 9L, 6L, 5L, 10L, 2L, 3L, 6L, 13L, 
    5L, 2L, 2L, 13L, 6L, 7L, 13L, 3L, 7L, 5L, 9L, 5L, 7L, 9L, 
    10L, 7L, 13L, 3L, 13L, 6L, 2L, 10L, 6L, 9L, 9L, 9L, 13L, 
    6L, 7L, 5L, 13L, 3L, 13L, 6L, 10L, 13L, 3L, 7L, 2L, 3L, 9L, 
    10L, 2L, 6L, 6L, 2L, 7L, 6L, 5L, 13L, 2L, 13L, 2L, 3L, 9L, 
    13L, 9L, 7L, 10L, 2L, 13L, 2L, 13L, 10L, 7L, 7L, 9L, 3L, 
    6L, 5L, 5L, 9L, 7L, 13L, 2L, 3L, 6L, 6L, 2L, 13L, 10L, 13L, 
    13L, 10L, 13L, 6L, 5L, 2L, 5L, 6L, 6L, 13L, 7L, 13L, 7L, 
    13L, 13L, 13L, 9L, 13L, 3L, 13L, 13L, 10L, 3L, 10L, 7L, 13L, 
    7L, 5L, 3L, 13L, 9L, 5L, 10L, 2L, 6L, 6L, 2L, 13L, 13L, 13L, 
    9L, 6L, 10L, 5L, 13L, 13L, 7L, 6L, 6L, 7L, 7L, 6L, 3L, 2L, 
    9L, 2L, 5L, 9L, 9L, 2L, 13L, 10L, 13L, 9L, 10L, 2L, 6L, 7L, 
    6L, 2L, 5L, 13L, 5L, 3L, 9L, 7L, 13L, 7L, 3L, 9L, 7L, 9L, 
    3L, 2L, 7L, 2L, 3L, 7L, 7L, 6L, 3L, 5L, 9L, 9L, 10L, 6L, 
    6L, 10L, 2L, 10L, 6L, 6L, 5L, 13L, 3L, 13L, 3L, 3L, 2L), .Label = c("1991", 
    "1992", "1993", "1994", "1995", "1996", "1998", "2000", "2002", 
    "2005", "2008", "2014", "2018"), class = "factor"), born_adult = c(1994, 
    1953, 1937, 1944, 1996, 1977, 1944, 1953, 2001, 1976, 1963, 
    1950, 1978, 1984, 1938, 1969, 1928, 1977, 1943, 1945, 1951, 
    1968, 1959, 1971, 1978, 1998, 1951, 1976, 1951, 1987, 1950, 
    1969, 1955, 1946, 1981, 2008, 1968, 1975, 1957, 1942, 1950, 
    1978, 1993, 1986, 1974, 1982, 1960, 1948, 1953, 1943, 1980, 
    1963, 1943, 1944, 1958, 1953, 1937, 1971, 1971, 1983, 1954, 
    1984, 1979, 1952, 1984, 1946, 1959, 1949, 1979, 1953, 1947, 
    1980, 1979, 1996, 1973, 1964, 1952, 1955, 1948, 1980, 1961, 
    1994, 1991, 1949, 1979, 1947, 1941, 1955, 1962, 2004, 1974, 
    1993, 1976, 1994, 1994, 1974, 1976, 1990, 1946, 1947, 1961, 
    1941, 1991, 1986, 1983, 1983, 1988, 1953, 1990, 1965, 1961, 
    1971, 1979, 1977, 1956, 1948, 2015, 1973, 1988, 1935, 2004, 
    1983, 1948, 1993, 1976, 1960, 1959, 1980, 1968, 1968, 1970, 
    1940, 1949, 1964, 1941, 2005, 1959, 1954, 1969, 1988, 1959, 
    1989, 1971, 1975, 1989, 1980, 1953, 1955, 1959, 1972, 1986, 
    1988, 1974, 1981, 1998, 2001, 1959, 1970, 1960, 1944, 1986, 
    1984, 2000, 1946, 1978, 1930, 1952, 1956, 1979, 1982, 1969, 
    1980, 1961, 1973, 1951, 1979, 1982, 1970, 1974, 1998, 1944, 
    1941, 1950, 1948, 1978, 1999, 1955, 1930, 1961, 1942, 1962, 
    1980, 1983, 1974, 1992, 1949, 2003, 1949, 1949, 1976)), row.names = c(NA, 
-200L), class = c("tbl_df", "tbl", "data.frame"))

Любая помощь очень ценится!

Answer 1

Спасибо Гэвину за предоставление решения для многогранных прогнозов! Тем не менее, для того, чтобы сформировать график прогнозов для наблюдаемых значений, я полагаю, что может существовать обходной путь с применением методов, предложенных King et al. 2001, чтобы получить прогнозы по диапазону независимой переменной на одном графике.

Основная проблема заключается в том, что прогнозы для когорт варьируются в зависимости от их значений на других ковариатах, что делает график таким волнистым. Чтобы решить эту проблему, мы можем предположить, что общий прогноз для данной когорты ненаблюдаем, но может быть вменен из прогнозов и их стандартных ошибок, которые мы получаем из функции mgcv::predict.gam. Используя методы, описанные на странице 53 в King et al. В статье 2001 года мы можем рассчитать общие прогнозируемые значения с их соответствующей стандартной ошибкой.

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

и

Чтобы реализовать эти вычисления в R, мы можем просто использовать некоторые dplyr функции:

predData <- transform(cbind(data.frame(pred)))
predBornAdult <- predData %>% 
group_by(born_adult) %>% 
mutate(m = n(),
       mean_fit = mean(fit),
       S_sq = (fit - mean_fit)^2/(m - 1)) %>%
replace_na(list(S_sq = 0)) %>% # For those cohorts, where we only have one prediction
summarize(fit = mean(fit),
          S_sq = mean(S_sq),
          se.fit2 = mean(se.fit^2) + S_sq,
          se.fit = sqrt(se.fit2)) %>% 
ungroup() %>% 
mutate(uprP = fit + crit * se.fit,
       lwrP = fit - crit * se.fit) %>% 
select(born_adult, fit, uprP, lwrP)

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

ggplot(predBornAdult, aes(x=born_adult, 
                     y=fit)) + 
geom_errorbar(aes(ymin = lwrP, ymax = uprP)) +
geom_point(size = 1)

Поскольку мы всегда хотим легко определять тренды как для прогнозов, так и для их неопределенности, теперь мы можем добавить geom_smooth, чтобы визуализировать потенциальную основутренды:

ggplot(predBornAdult, aes(x=born_adult, y=fit)) + 
geom_errorbar(aes(ymin = lwrP, ymax = uprP), alpha = 0.2) +
geom_point(alpha = 0.2, size = 1) + 
geom_smooth(aes(y = fit), se = F, alpha = 0.5) + 
geom_smooth(aes(y = lwrP), se = F, alpha = 0.5, linetype = "solid", size = 0.5) + 
geom_smooth(aes(y = uprP), se = F, alpha = 0.5, linetype = "solid", size = 0.5)

Answer 2

Я думаю, что это просто проблема из-за того, что все данные перепутаны, потому что вы прогнозируете для наблюдений.

Этот график основан на ваших, но я строю подгонянные значения, от mgcv::predict.gam()и вы сразу увидите проблему:

Красная и синяя линии - это верхний и нижний интервалы, соответственно, а черная линия соответствуетзначение от mgcv::predict.gam(). Поскольку последним не манипулировали вообще, я склонен полагать, что интервалы здесь.

Это то, что вы получите, если вы используете точечные / поперечные функции вероятных интервалов:

, которые, помимо того, чтобы быть более узкими, демонстрируют то же поведение.

Если вы просто пытаетесь получить предсказания и одновременные интервалы в диапазонеborn_adult для каждой комбинации двух факторных переменных, затем вы должны создать новые данные для прогнозирования, повторяющие последовательность значений born_adult для всех комбинаций age_grps и period.f. Вот что для 50 значений born_adult - подобранное сглаживание в основном линейное, поэтому даже 50 является избыточным, но интервалы сглаживаются при увеличении n - используя expand.grid():

pdat <- with(dat, expand.grid(
  born_adult = seq(min(born_adult), max(born_adult), length = 50),
  age_grps = unique(age_grps),
  period.f = unique(period.f)))

Затем, повторяя ваш код, но добавляя newdata = pdat к одновременным вычислениям интервалов, мы получаем их для наших данных прогнозирования, а не исходных данных

Vb <- vcov(mod)
pred2 <- predict(mod, newdata = pdat, se.fit = TRUE)
N <- 10000
BUdiff <- rmvn(N, mu = rep(0, nrow(Vb)), sig = Vb)
Cg <- predict(mod, newdata = pdat, type = "lpmatrix")
simDev <- Cg %*% t(BUdiff)
absDev <- abs(sweep(simDev, 1, pred2$se.fit, FUN = "/"))
masd <- apply(absDev, 2L, max)
crit2 <- quantile(masd, prob = 0.95, type = 8)

Здесь я создаю то же самое, что и ваш predData, ноЯ cbind() на pdat вместо исходных данных и добавляю поточечные интервалы просто как проверку:

predData2 <- transform(cbind(data.frame(pred2), pdat),
                       uprP = fit + (crit2 * se.fit),
                       lwrP = fit - (crit2 * se.fit),
                       uprCI = fit + (2 * se.fit),
                       lwrCI = fit - (2 * se.fit))

, которая при построении с использованием

ggplot(predData2) +
  geom_ribbon(aes(x = born_adult, ymin = lwrP, ymax = uprP),
              alpha = 0.2, fill = "red") +
  geom_ribbon(aes(x = born_adult, ymin = lwrCI, ymax = uprCI),
              alpha = 0.2, fill = "red") +
  geom_line(aes(x = born_adult, y = fit)) +
  facet_wrap(vars(period.f, age_grps))

дает следующее:

Если вы хотите использовать его только для наблюдаемых комбинаций age_grps и period.f, которые необходимо создатьданные прогноза несколько иные, но общая идея все же применима. (Или вы можете просто сделать то, что я сделал, а затем удалить все строки, где комбинация age_grps и period.f не является одной из наблюдаемых комбинаций.)