r - 使用 bam 的零膨胀模型 (ziP) 中的错误

Question

我正在尝试使用 bam 来运行以下广义加性模型:

m <- bam(result ~ factor(city) + factor(year) + lnpopulation + s(lnincome), data=full_df, na.action=na.omit, family=ziP(theta = NULL, link = "identity",b=0))

但出现以下错误:

Error in bam(result : extended families not supported by bam

bam 的文档提到了以下内容:

This is a family object specifying the distribution and link to use in fitting etc. See glm and family for more details. The extended families listed in family.mgcv can also be used.

family.mgcv 确实包含 ziP。我究竟做错了什么？任何指导将不胜感激。谢谢!

从 r-help 重新发布。

此致

米卢

Answer 1

错误消息表明您正在使用的 mgcv 版本不支持这些由 mgcv 的 gam()< 支持的扩展系列 函数。

谢天谢地，您的问题的解决方案很简单，因为 Simon Wood 现在已经从版本 1.8-19 实现了此功能(bam() 中的扩展系列)，如 ChangeLog 中所述:

1.8-19

** bam() now accepts extended families (i.e. nb, tw, ocat etc)

当前版本是 1.8-22，它修复了一些与您正在寻找的功能相关的错误，因此请确保您更新到最新版本。

这是一个修改自 ?ziP 的例子

 ## function to simulated zip data

 rzip <- function(gamma,theta= c(-2,.3)) {
 ## generate zero inflated Poisson random variables, where 
 ## lambda = exp(gamma), eta = theta[1] + exp(theta[2])*gamma
 ## and 1-p = exp(-exp(eta)).
   y <- gamma; n <- length(y)
   lambda <- exp(gamma)
   eta <- theta[1] + exp(theta[2])*gamma
   p <- 1- exp(-exp(eta))
   ind <- p > runif(n)
   y[!ind] <- 0
   np <- sum(ind)
   ## generate from zero truncated Poisson, given presence...
   y[ind] <- qpois(runif(np,dpois(0,lambda[ind]),1),lambda[ind])
   y
 } 

 library('mgcv')

 ## Simulate some ziP data...
 set.seed(1);n<-400
 dat <- gamSim(1,n=n)
 dat$y <- rzip(dat$f/4-1)

 b <- bam(y ~ s(x0) + s(x1) + s(x2) + s(x3),
          family = ziP(), data = dat)

这为我提供了以下拟合模型:

> summary(b)

Family: Zero inflated Poisson(-1.855,1.244) 
Link function: identity 

Formula:
y ~ s(x0) + s(x1) + s(x2) + s(x3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.97426    0.04988   19.53   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
        edf Ref.df      F p-value    
s(x0) 2.396  2.989  2.336  0.0759 .  
s(x1) 2.784  3.464 77.217  <2e-16 ***
s(x2) 7.397  8.317 59.364  <2e-16 ***
s(x3) 1.235  1.428  0.269  0.5888    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Deviance explained =   69%
fREML = 593.94  Scale est. = 1         n = 400