R gfam 分组家庭

说明

与 gam 或 bam 一起使用的系列，允许单变量响应向量由来自多个不同分布的变量组成。响应变量以 2 列矩阵的形式提供，其中第一列包含响应观测值，第二列索引其来源的分布(族)。 gfam 将族列表作为其单个参数。

对于通过共享某些组件的模型链接的不同源的数据进行建模非常有用。不共享的平滑模型组件通常使用 by 变量进行处理(请参阅 gam.models )。

用法

gfam(fl)

参数

fl	家庭名单。这些可以是继承自 family 或 extended.family 的任何系列，可与 gam 一起使用，前提是它们通常不需要矩阵响应变量。

细节

gfam 的每个组件函数使用列表 fl 中提供的族来获取该族数据子集所需的数量，并适当地组合结果。例如，它通过计算应用于其数据子集的族特定偏差和每个族的导数并对它们求和，提供模型的总偏差(两倍负对数似然)及其导数。其他量的计算方法相同。

正则指数族计算的数量与扩展族不同，因此 gfam 在内部将这些族产生的结果转换为 extended.family 形式。

显然，尺度参数必须针对每个族单独处理，并被视为要估计的参数，就像其他extended.family非位置分布参数一样。这又是在内部处理的。此要求是始终生成 extended.family 的部分原因，即使 fl 的所有元素都是标准指数族也是如此。因此，平滑参数估计总是通过 REML 或 NCV 进行。

请注意，零偏差目前是通过假设每个系列的单个参数模型而不是仅一个参数来计算的，这可能会稍微降低解释的偏差。还要注意，残差检查可能应该通过按家庭分解残差来完成。因此，不提供函数来促进 qq.gam 的残差检查。

对响应量表的预测需要提供家族索引向量以及响应名称，作为新预测数据的一部分。但是，诸如 ocat 等通常为预测类型 "response" 生成矩阵预测的系列，在成为 gfam 的一部分时将无法执行此操作。

gfam 依赖于 Wood、Pya 和 Saefken (2016) 中的方法。

值

类 extended.family 的对象。

例子

library(mgcv)
## a mixed family simulator function to play with...
sim.gfam <- function(dist,n=400) {
## dist can be norm, pois, gamma, binom, nbinom, tw, ocat (R assumed 4)
## links used are identitiy, log or logit.
  dat <- gamSim(1,n=n,verbose=FALSE)
  nf <- length(dist) ## how many families
  fin <- c(1:nf,sample(1:nf,n-nf,replace=TRUE)) ## family index
  dat[,6:10] <- dat[,6:10]/5 ## a scale that works for all links used
  y <- dat$y;
  for (i in 1:nf) {
    ii <- which(fin==i) ## index of current family
    ni <- length(ii);fi <- dat$f[ii]
    if (dist[i]=="norm") {
      y[ii] <- fi + rnorm(ni)*.5
    } else if (dist[i]=="pois") {
      y[ii] <- rpois(ni,exp(fi))
    } else if (dist[i]=="gamma") {
      scale <- .5
      y[ii] <- rgamma(ni,shape=1/scale,scale=exp(fi)*scale)
    } else if (dist[i]=="binom") {
      y[ii] <- rbinom(ni,1,binomial()$linkinv(fi))
    } else if (dist[i]=="nbinom") {
      y[ii] <- rnbinom(ni,size=3,mu=exp(fi))
    } else if (dist[i]=="tw") {
      y[ii] <- rTweedie(exp(fi),p=1.5,phi=1.5)
    } else if (dist[i]=="ocat") {
      alpha <- c(-Inf,1,2,2.5,Inf)
      R <- length(alpha)-1
      yi <- fi
      u <- runif(ni)
      u <- yi + log(u/(1-u)) 
      for (j in 1:R) {
        yi[u > alpha[j]&u <= alpha[j+1]] <- j
      }
      y[ii] <- yi
    }
  }
  dat$y <- cbind(y,fin)
  dat
} ## sim.gfam

## some examples

dat <- sim.gfam(c("binom","tw","norm"))
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),
         family=gfam(list(binomial,tw,gaussian)),data=dat)
predict(b,data.frame(y=1:3,x0=c(.5,.5,.5),x1=c(.3,.2,.3),
        x2=c(.2,.5,.8),x3=c(.1,.5,.9)),type="response",se=TRUE)
summary(b)
plot(b,pages=1)

## set up model so that only the binomial observations depend
## on x0...

dat$id1 <- as.numeric(dat$y[,2]==1)
b1 <- gam(y~s(x0,by=id1)+s(x1)+s(x2)+s(x3),
         family=gfam(list(binomial,tw,gaussian)),data=dat)
plot(b1,pages=1) ## note the CI width increase	 

作者

Simon N. Wood simon.wood@r-project.org

参考

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and model selection for general smooth models. Journal of the American Statistical Association 111, 1548-1575 doi:10.1080/01621459.2016.1180986