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

Question

我正在尝试使用高斯位置尺度模型系列在 mgcv 中使用分层通用加性模型。然而，这一系列函数在尝试拟合时会抛出一个神秘错误。

遵循最近一篇论文 (https://peerj.com/articles/6876/) along with the guidance from the notes in the lmvar package write-up (https://cran.r-project.org/web/packages/lmvar/vignettes/Intro.html) 中关于 HGAM 的指导。


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

此returns错误信息：


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

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

因此，根据我上面显示的第一个示例，如果您想要线性预测变量中的常数项仅，您需要使用 R 的符号表示截距或常数项； ~ 1

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

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

Error message in mgcv GAM when using family = gaulss()

r

gam

mgcv