Доверительные интервалы прогнозируемых вероятностей на основе порядковой регрессии, построенной с помощью 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. Самое главное, как я могу рассчитать доверительные интервалы прогнозируемой вероятности каждой категории в каждом наблюдении?

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