使用 family = gaulss() 时 mgcv GAM 中的错误消息答案

【问题标题】：Error message in mgcv GAM when using family = gaulss()使用 family = gaulss() 时 mgcv GAM 中的错误消息
【发布时间】：2019-08-19 14:31:10
【问题描述】：

我正在尝试使用高斯位置尺度模型系列在 mgcv 中使用分层通用加法模型。但是，这一系列函数在尝试拟合时会引发一个神秘的错误。

遵循最近一篇论文 (https://peerj.com/articles/6876/) 中关于 HGAM 的指导以及 lmvar 包文章 (https://cran.r-project.org/web/packages/lmvar/vignettes/Intro.html) 中注释的指导。


library(mgcv); library(datasets)

# Using the CO2 dataset for illustration

data <- datasets::CO2

# Changing the 'Plant' factor from ordered to unordered

data$Plant <- factor(data$Plant, ordered = FALSE)

# This model fits with no errors

CO2_modG_1 <- gam(log(uptake) ~ s(log(conc), k = 5, bs = 'tp') + s(Plant, k = 12, bs = 're'), data = data, method = 'REML', family = 'gaussian')

# This model fails with error

CO2_modG_2 <- gam(log(uptake) ~ s(log(conc), k = 5, bs = 'tp') + s(Plant, k = 12, bs = 're'), data = data, method = 'REML', family = gaulss())

这会返回错误信息：


Error in while (sqrt(mean(dlb/(dlb + lami * dS * rm)) * mean(dlb)/mean(dlb +  : 
  missing value where TRUE/FALSE needed

【问题讨论】：

标签： r gam mgcv

【解决方案1】：

当使用gaulss() 系列拟合高斯位置尺度模型时，您必须将两个公式对象的列表传递给gam()。

您没有说要如何对模型的比例组件进行建模，所以这里有一个示例，它应该等效于 gaussian() 系列，其中我们有一个比例常数项和一个您提供的线性预测器为平均值。

CO2_modG_2 <- gam(list(log(uptake) ~ s(log(conc), k = 5, bs = 'tp') + 
                         s(Plant, k = 12, bs = 're'), 
                       ~ 1), 
                  data = data, method = 'REML', family = gaulss())

如果你想让每个植物都有自己的差异，那么我们可以在第二个公式中添加项：

CO2_modG_3 <- gam(list(log(uptake) ~ s(log(conc), k = 5, bs = 'tp') + 
                         s(Plant, k = 12, bs = 're'), 
                       ~ s(Plant, k = 12, bs = 're')), 
                  data = data, method = 'REML', family = gaulss())

重要的是，您必须为该族提供一个包含两个公式对象的列表，并且第二个公式应该只有一个波浪号加上公式的右侧，以指定比例参数的线性预测器。

list(response ~ x1 + s(x2), # location linear predictor, left & right hand sided
              ~ x1 + s(x3)  # scale linear predictor, right-hand side only
    )

因此，根据我上面展示的第一个示例，如果您希望在线性预测器中使用常数项 only 来表示比例，则需要使用 R 符号表示截距或常数项； ~ 1

list(response ~ x1 + s(x2), # location linear predictor, left & right hand sided
              ~ 1           # scale linear predictor, constant
    )

【讨论】：