mgcv：新しいデータが与えられたときの応答の予測分布を得る（負の2項の例）

Question

GAM（とGLM）では、条件付き尤度モデルに適合しています。したがって、モデルをフィッティングした後、新しい入力xとレスポンスyの場合、という特定の値の予測確率または密度を計算する必要があります。これはyです。たとえば、さまざまなモデルの検証データに対する適合性を比較するために、これを行うことができます。適合するGAMをmgcvに入れると便利でしょうか？それ以外の場合は、使用する密度の正確な形をどのように把握すれば、パラメータを適切に差し込むことができますか？具体的な例として

、負の二項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) 

をそして今、私はx=(.1,.2,.3,.4)与えy=7、たとえば、予測確率を計算します。

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) 
    } 

は今、あなたは一般的に

p.y(7) 
# 0.07810065 

---

使用し、y = 7の確率が、指定された新しいデータを条件にしたいため、数値積分することにより、このアプローチは簡単ではありません。例えば、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

注として経験的ベイズすなわち、すべての上記の方法は、推定された平滑化パラメータの条件です。 「完全ベイズ」のようなものをお望みなら、unconditional = TRUEをpredict()に設定してください。

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が一定となり、：おそらく一部の人々は、このような単純な解決策を想定している

備考2

ランダムではありません。これはではなく、予測分布です。予測分布はモデル推定の不確定性を統合するだろう。