在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)。
是的。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变得固定而不是随机。这不是预测分布。预测分布将整合模型估计的不确定性。