mgcv：獲得給定新數據的響應預測分佈（負二項式示例）

Question

在GAM（和GLM，就此而言）中，我們擬合了條件似然模型。因此，在擬合模型之後，對於新輸入x和響應y，我應該能夠計算給定x的特定值y的預測概率或密度。例如，我可能想要這樣做來比較驗證數據上各種模型的適合性。有沒有一種方便的方法可以在mgcv中使用合適的GAM來做到這一點？否則，我如何計算出使用的密度的確切形式，以便適當地插入參數？

作爲一個具體例子，考慮一個負二項分佈GAM：

## From ?negbin 
library(mgcv) 
set.seed(3) 
n<-400 
dat <- gamSim(1,n=n) 
g <- exp(dat$f/5) 

## negative binomial data... 
dat$y <- rnbinom(g,size=3,mu=g) 

## fit with theta estimation... 
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=nb(),data=dat) 

，現在我想計算的，也就是說，y=7給出x=(.1,.2,.3,.4)預測概率。

Answer 1

是的。 mgcv正在做（經驗）貝葉斯估計，所以你可以獲得預測分佈。對於你的例子，這裏是如何。

# prediction on the link (with standard error) 
l <- predict(b, newdata = data.frame(x0 = 0.1, x1 = 0.2, x2 = 0.3, x3 = 0.4), se.fit = TRUE) 

# Under central limit theory in GLM theory, link value is normally distributed 
# for negative binomial with `log` link, the response is log-normal 
p.mu <- function (mu) dlnorm(mu, l[[1]], l[[2]]) 

# joint density of `y` and `mu` 
p.y.mu <- function (y, mu) dnbinom(y, size = 3, mu = mu) * p.mu(mu) 

# marginal probability (not density as negative binomial is discrete) of `y` (integrating out `mu`) 
# I have carefully written this function so it can take vector input 
p.y <- function (y) { 
    scalar.p.y <- function (scalar.y) integrate(p.y.mu, lower = 0, upper = Inf, y = scalar.y)[[1]] 
    sapply(y, scalar.p.y) 
    } 

現在既然你想的y = 7概率，條件上指定新的數據，使用

p.y(7) 
# 0.07810065 

---

一般情況下，這種方法用數值積分是不容易的。例如，如果其他鏈接功能如sqrt()用於負二項式，則響應的分佈不是那麼直截了當（雖然也不難推導出來）。

現在我提供一個基於抽樣的方法或蒙特卡洛方法。這與貝葉斯過程最爲相似。

N <- 1000 # samples size 

set.seed(0) 

## draw N samples from posterior of `mu` 
sample.mu <- b$family$linkinv(rnorm(N, l[[1]], l[[2]])) 

## draw N samples from likelihood `Pr(y|mu)` 
sample.y <- rnbinom(1000, size = 3, mu = sample.mu) 

## Monte Carlo estimation for `Pr(y = 7)` 
mean(sample.y == 7) 
# 0.076 

---

備註1

注意，作爲經驗貝葉斯，所有上述方法估計的平滑參數爲條件的。如果你想要像「完整貝葉斯」那樣的東西，請在predict()中設置unconditional = TRUE。

備註2

也許有些人假定溶液這樣簡單：

mu <- predict(b, newdata = data.frame(x0 = 0.1, x1 = 0.2, x2 = 0.3, x3 = 0.4), type = "response") 
dnbinom(7, size = 3, mu = mu) 

這樣的結果是在迴歸係數有條件的（假定固定不確定性），從而mu變得固定和不是隨機的。這是而不是預測分佈。預測分佈將整合模型估計的不確定性。