使用mgcv的gaulss()拟合HGAM时出现逻辑值缺失错误

Hey there, let's break down why your gaulss() model is throwing that confusing error! The core issue is how location-scale families like gaulss() work in mgcv—they require you to explicitly define both the mean (location) and variance (scale) components of the model, and your original code only specifies the location part. This missing scale structure causes numerical instability during fitting, triggering the missing value where TRUE/FALSE needed error in the REML optimization loop.

What's Wrong with the Original Code? #

When you use family = gaulss() without a scale formula, mgcv might default to a constant scale, but combining this implicit setup with random effects (s(Plant, bs='re')) leads to unstable optimization. The error you see comes from the fitting loop hitting a missing value when evaluating convergence conditions.

How to Fix It #

You need to pass a list of two formulas to gam(): one for the location (mean of log(uptake)), and one for the scale (variance of the response). Here are working examples tailored to your CO2 dataset:

1. Constant Variance (Simplest Fix)

If you want variance to stay constant across observations (similar to a regular Gaussian model but using gaulss()), explicitly set the scale formula to ~1:

library(mgcv); library(datasets)
data <- datasets::CO2
data$Plant <- factor(data$Plant, ordered = FALSE)

# Working gaulss model with constant scale
CO2_modG_2 <- gam(
  list(
    log(uptake) ~ s(log(conc), k = 5, bs = 'tp') + s(Plant, k = 12, bs = 're'),
    ~ 1  # Explicit constant scale component
  ),
  data = data, 
  method = 'REML', 
  family = gaulss()
)

# Inspect results
summary(CO2_modG_2)

2. Variance Dependent on Random Effects

If you suspect variance of log(uptake) varies by plant (a common HGAM extension), add a random effect for Plant to the scale formula:

# gaulss model with variance varying by Plant
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')  # Scale tied to Plant random effect
  ),
  data = data, 
  method = 'REML', 
  family = gaulss()
)

3. Variance Dependent on a Smooth Term

You can also let variance change with conc by adding a smooth term to the scale formula:

# gaulss model with variance varying by conc
CO2_modG_4 <- gam(
  list(
    log(uptake) ~ s(log(conc), k = 5, bs = 'tp') + s(Plant, k = 12, bs = 're'),
    ~ s(log(conc), k = 5, bs = 'tp')  # Scale varies with conc
  ),
  data = data, 
  method = 'REML', 
  family = gaulss()
)

Key Takeaway #

Location-scale families in mgcv depend on explicit specification of both components to avoid numerical issues. Always use the list(location_formula, scale_formula) syntax when fitting models with gaulss(), scatreg(), or similar families—it clarifies your model intent and prevents the kind of error you encountered.

内容的提问来源于stack exchange，提问作者momcanally