Доверительные интервалы прогнозируемых вероятностей на основе порядковой регрессии, построенной с помощью mgcv::gam(..., family = ocat(R = ...))

Question

Я хочу получить доверительные интервалы предсказанных вероятностей из модели порядковой регрессии, построенной с mgcv::gam(..., family = ocat(R = ...)), но я не могу понять, как. Мои попытки ниже.

Я создаю пример набора данных и строю порядковую регрессию с gam(), Обратите внимание, что x отсортировано

library("mgcv")
set.seed(1)
d <- data.frame(
  y = sample(1:5, size = 200, replace = TRUE),
  x = sort(runif(200))
)
d.gam <- gam(y ~ x, family = ocat(R = 5), data = d)

Если я использую predict.gam(..., type = "response")Я получаю матрицу прогнозируемой вероятности каждой категории для каждого наблюдения. поскольку x был включен как линейный термин, предсказанные вероятности также являются монотонными (например, вероятность категории 1 постоянно уменьшается, а вероятность категории 4 постоянно увеличивается). Все идет нормально.

> d.response <- predict(d.gam, d, se = TRUE, type = "response")
> d.response$fit[c(1:3, 100:103, 198:200), ]

         [,1]      [,2]      [,3]      [,4]      [,5]
1   0.1533162 0.2327623 0.2160818 0.2333859 0.1644538
2   0.1529262 0.2324397 0.2160737 0.2336929 0.1648675
3   0.1528949 0.2324138 0.2160730 0.2337176 0.1649007
100 0.1452906 0.2259199 0.2157159 0.2397503 0.1733233
101 0.1451200 0.2257698 0.2157034 0.2398865 0.1735203
102 0.1451026 0.2257544 0.2157021 0.2399005 0.1735405
103 0.1451008 0.2257528 0.2157020 0.2399019 0.1735425
198 0.1342714 0.2158035 0.2144608 0.2486086 0.1868556
199 0.1342414 0.2157748 0.2144561 0.2486328 0.1868948
200 0.1341483 0.2156856 0.2144414 0.2487081 0.1870167

Это, однако, не относится к стандартным ошибкам.

> d.response$se.fit[c(1:3, 100:103, 198:200), ] 

          [,1]       [,2]         [,3]       [,4]       [,5]
1   0.03015269 0.02490338 0.0005903245 0.02372872 0.03191767
2   0.02918116 0.02417559 0.0006422390 0.02298270 0.03101629
3   0.02910391 0.02411755 0.0006462272 0.02292332 0.03094436
100 0.01574518 0.01384974 0.0011449970 0.01257287 0.01816704
101 0.01566419 0.01379678 0.0011579659 0.01251146 0.01810748
102 0.01565671 0.01379206 0.0011593446 0.01250581 0.01810230
103 0.01565595 0.01379159 0.0011594855 0.01250525 0.01810178
198 0.03108344 0.02975641 0.0048968388 0.02510754 0.04062915
199 0.03115134 0.02982819 0.0049153610 0.02516275 0.04073214
200 0.03136282 0.03005193 0.0049732667 0.02533468 0.04105333

Прежде всего, я даже не уверен, действительно ли это в масштабе вероятности. Если это так, то почему шкала, по-видимому, различается по категориям (например, SE категории 3 намного меньше, чем SE категории 5)? Кроме того, SE понижается сначала и затем повышается как x увеличивается. Поскольку данные предположительно равномерно распределены по xЯ ожидаю, что аналогичные SE вместе x, Почему это не так?

Поскольку в описанной выше процедуре, скорее всего, что-то не так, я подумал, что следует вывести SE на основе линейного предиктора (predict(..., type = "link")). Но, в отличие от вышеизложенного, predict(..., type = "link") возвращает одно значение для каждого наблюдения, и я не смог выяснить взаимосвязь между линейным предиктором и вероятностью каждой категории.

> d.link <- predict(d.gam, d, se = TRUE, type = "link")
> head(d.link$fit)

        1         2         3         4         5         6 
0.7088250 0.7118324 0.7120737 0.7124732 0.7143695 0.7146253 

> head(d.link$se.fit)
        1         2         3         4         5         6 
0.2322826 0.2252680 0.2247092 0.2237855 0.2194252 0.2188400 

Итак, мои вопросы следующие:

1. Являются ли значения, возвращаемые predict.gam(..., type = "response")$se.fit в шкале вероятностей? Если нет, то что они?
2. Могу ли я получить вероятность каждой категории на основе вывода predict.gam(..., type = "link")? Если так, то как?
3. Что наиболее важно, как я мог вычислить доверительные интервалы предсказанной вероятности каждой категории в каждом наблюдении?

Заранее спасибо!