使用 gamm4 的 beta 和 quassi 系列模型选择

Question

我有两个响应符合 beta（也称为 betar）和 Poisson 系列，我正在研究使用 [=11= 拟合加法混合模型] 和准家族（计数数据过于分散）。

我知道我可以使用 mgcv 包中的 gamm 函数，它接受 beta 和准系列，但是我正在考虑它使用 PQL，并且 AIC 报告对比较模型没有用 - 这是我分析的主要 objective。

在计数响应的情况下，我知道 QAIC 已用于 ranking/comparing 过度分散的混合模型，但我找不到任何说明它适用于过度分散的 GAMM。

我知道这些可能是两个问题合而为一，但它们都有一个共同的主题，即大家庭的模型选择，并且可能有不同的解决方案。下面我为每个案例提供了可重现的示例。

##generate data
library(gamm4)
library(mgcv)
dat <- gamSim(1,n=400,scale=2)
dat<-subset(dat, select=c(x0,x1,x2,x3,f) )
dat$g <- as.factor(sample(1:20,400,replace=TRUE))#random factor
dat$yb<-runif(400)#yb ranges between 0-1 hence fitted with beta family
dat$f <- dat$f + model.matrix(~ g-1)%*%rnorm(20)*2
dat$yp <- rpois(400,exp(dat$f/7))#y2 is counts hence poisson family

#beta family example with gamm function (this works - however not sure if the subsequent model comparisons are valid!) 
m1b<-    gamm(yb~s(x0)+s(x1)+s(x2)+s(x3),family=betar(link='logit'),data=dat,random=list(g=~1))
m2b<-gamm(yb~s(x1)+s(x2)+s(x3),family=betar(link='logit'),data=dat,random=list(g=~1))
m3b<-gamm(yb~s(x0)+s(x2)+s(x3),family=betar(link='logit'),data=dat,random=list(g=~1))

#AIC to compare models  
AIC(m1b,m2b,m3b)

#try the same using gamm4 (ideally)- it obviously fails with beta family.
m<-gamm4(yb~s(x0)+s(x1)+s(x2)+s(x3),family=betar(link='logit'),data=dat,random = ~ (1|g)) 

##Example with quassi family - yp response is overdispersed count data (may not be overdispered in this example
#example using gamm function
m1p<-gamm(yp~s(x0)+s(x1)+s(x2)+s(x3),family = quasipoisson,data=dat,random=list(g=~1))
m2p<-gamm(yp~s(x1)+s(x2)+s(x3),family = quasipoisson,data=dat,random=list(g=~1))
m3p<-gamm(yp~s(x0)+s(x2)+s(x3),family = quasipoisson,data=dat,random=list(g=~1))

#AIC to compare models
AIC(m1p,m2p,m3p)

#again the example with using gamm4 function will not work as it doesnt accept quassi falimies 
m<-gamm4(yp~s(x0)+s(x1)+s(x2)+s(x3),family = quasipoisson,data=dat,random = ~ (1|g))

Answer 1

你在这里有很多问题，但我会尽力解决它们。基本上，您想用

拟合参数统计模型

随机效果（nlme、lme4）
来自指数族的分布... (MASS::glmmPQL, lme4::glmer)
...过度分散...
...或超出指数族的分布，例如 Beta 分布 (VGAM, betareg)
加法models/splines (splines) ...
...或惩罚回归样条，自动调整平滑项的复杂度
...使用真实似然模型而不是边际模型或 quasi-likelihood 模型（例如 GEE、PQL），因此您可以进行经典推理

上面的每个特定问题都会在 model-fitting 练习中增加 1 个或更多 "difficulty points" ... 通常一旦你的分数超过 +3 左右，你就必须找到一种方法来在一些你想要的事情上妥协或走捷径。您已经正确地将 gamm 和 gamm4 识别为做一些您想要的事情，但您无法得到所有的东西。一些建议：

过度分散

处理过度分散的一种方法是使用observation-level随机效应，例如

dat$obs <- factor(seq(nrow(dat)))
m <- gamm4(yp~s(x0)+s(x1)+s(x2)+s(x3),
           family = poisson,data=dat,random = ~ (1|g)+(1|obs))

另一种选择是自己调整过度分散，如果您认为这有道理，例如：

m0 <- gamm4(yp~s(x0)+s(x1)+s(x2)+s(x3),family = poisson,data=dat,random = ~ (1|g))

首先计算过度离散：

(phi <- sum(residuals(m0$gam,type="pearson")^2/df.residual(m0$gam)))
## [1] 1.003436

（如果我们用 m0$mer 重复这个练习，我们得到 0.9939696：结果几乎完全等于 1，因为我们首先从泊松分布生成数据......）

(qaic <- -2*logLik(m0$mer)/phi + 2*lme4:::npar.merMod(m0$mer))

N.B. 我 猜测从 gamm4 的各个组成部分构建可能性等是有意义的方法;使用风险自负

替代分布

glmmADMB 和 glmmTMB 包（都是 off-CRAN 但可通过 Google 找到 ...）都可以处理混合 Beta 模型。他们不能做惩罚回归样条，但你可以通过 splines::ns() 或 splines::bs() 使用常规样条（但你必须决定适当的复杂程度 - 也许你可以从初步猜测 gamm 或 mgcv 适合...)

library(glmmADMB)
library(splines)
m3b <- glmmadmb(yb~ns(x0,2)+ns(x1,2)+ns(x2,5)+ns(x3,2)+(1|g),
       family="beta",link="logit",data=dat)

glmmTMB 包可以原则上 这样做:

library(glmmTMB)
m2b <- glmmTMB(yb~ns(x0,2)+ns(x1,2)+ns(x2,5)+ns(x3,2)+(1|g),
       family=list(family="beta",link="logit"),data=dat)

但是这个包正在开发中，目前的结果集没有意义——所以我现在可能会犹豫是否要使用它。

使用 gamm4 的 beta 和 quassi 系列模型选择

Model selection with beta and quassi families using gamm4

r

mixed-models

gam

过度分散

替代分布