Dezmințirea mitului GAMLSS | R-bloggeri

Modelele aditive generalizate pentru locație, scară și formă (GAMLSS) sunt populare pentru estimarea ecuațiilor de referință, de exemplu, pentru măsurătorile funcției pulmonare. O lucrare recentă a susținut că „dezamăgește” mitul GAMLSS în acest context, arătând că modelele liniare pe bucăți mult mai simple funcționează la fel de bine. O nouă scrisoare către editor respinge aceste afirmații.

Citare

Robert A. Rigby, Mikis D. Stasinopoulos, Achim Zeileis, Sanja Stanojevic, Gillian Heller, Fernanda de Bastiani, Thomas Kneib, Andreas Mayr, Reto Stauffer, Nikolaus Umlauf (2025). „Scrisoare către editor care respinge „Deconfirmarea mitului GAMLSS: simplitatea domnește în diagnosticarea funcției pulmonare”. Urmează în Medicina respiratorie. Preprint disponibil la arXiv.org E-Print Arhiva arXiv:2512.09179 (stat.AP). doi:10.48550/arXiv.2512.09179

Abstract

Citim cu interes articolul de mai sus al lui Zavorsky (2025, Respiratory Medicine) referitor la ecuațiile de referință pentru testarea funcției pulmonare. Autorul compară un model aditiv generalizat pentru locație, scară și formă (GAMLSS), care este standardul adoptat de Global Lung Function Initiative (GLI), cu un model de regresie liniară segmentată (SLR), pentru variabilele funcției pulmonare. Autorul prezintă o comparație interesantă; cu toate acestea, există unele probleme fundamentale cu abordarea. Salutăm această oportunitate de a discuta problemele pe care le ridică. Susținerea autorului este că (1) SLR oferă „precizii de predicție la egalitate cu GAMLSS”; și (2) ecuațiile modelului GAMLSS sunt „complicate și necesită tabele spline suplimentare”, în timp ce SLR este „mai simplu, mai parsimonios și mai accesibil unui public mai larg”. Cu respect, nu suntem de acord cu ambele puncte.

Citiți scrisoarea completă ›

Replicare

Zavorsky furnizează datele pentru analiza sa în NHANES_2007_2012_Only_Acceptable_Spirometry_Values.csv.

Pe baza codului lui Zavorsky, am putut replica majoritatea analizelor sale și am putut reevalua modelele. Codul de replicare pentru scrisoarea către editor este disponibil în undebunking.R.

Pentru ilustrare, o pereche de modele pentru una dintre variabilele de răspuns (în subgrupul masculin) este ajustată și evaluată mai jos.

Ilustrare

Pentru a ilustra fluxul de lucru pe care l-am folosit pentru evaluarea regresiei liniare segmentate (SLR) și a modelelor aditive generalizate pentru locație, scară și formă (GAMLSS) pentru datele de spirometrie, luăm în considerare variabila de răspuns FEV1 pentru subgrupul masculin. VEMS este volumul expirator forțat într-o secundă în litri, adică volumul de aer care poate fi suflat forțat în prima secundă, după inspirația completă.

După descărcarea datelor (vezi linkul de mai sus) putem citi CSV, calcula vârsta pătrată (pe care Zavorsky o folosește oarecum confuz în loc de vârstă în SLR) și selectăm subsetul masculin.

d <- "NHANES_2007_2012_Only_Acceptable_Spirometry_Values.csv" |>
  read.csv() |>
  transform(Age2 = Age^2) |>
  subset(Sex == "Male")

Apoi încărcăm pachetele de care avem nevoie, toate disponibile de la CRAN, cu excepția pachetului topmodels care este încă pe R-Forge.

library("segmented")
library("gamlss")
library("distributions3")
library("topmodels")

Pentru a activa fluxul de lucru pentru prognozele probabilistice și evaluările modelului de la distributions3 şi topmodels pachete pentru modelele SLR montate de segmented avem nevoie de o prodist() metoda de extragere/predire distribuția de probabilitate corespunzătoare. Acest lucru nu funcționează din cutie deoarece segmented folosește doar o singură variație constantă, în timp ce Zavorsky folosește o varianță constantă pe bucăți. Funcția de mai jos implementează acest lucru pentru cazul special al unui segment lm() cu exact un punct de întrerupere.

prodist.segmented <- function(object, newdata = NULL, ...) {
  ## in-sample and out-of-sample data
  fitdata <- model.frame(object)
  if (is.null(newdata)) newdata <- fitdata

  ## predicted mean (m) and standard deviation (s)
  m <- predict(object, newdata = newdata)
  psi <- object$psi(, "Est.")
  seg <- as.numeric(fitdata((object$nameUV$Z)) > psi) + 1
  s <- sqrt(tapply(residuals(object)^2, seg, mean))
  s <- s(as.numeric(newdata((object$nameUV$Z)) > psi) + 1)

  ## Normal distribution object
  distributions3::Normal(mu = m, sigma = s)
}

Apoi am încadrat modelele SLR și GAMLSS în specificația (oarecum complicată) pe care Zavorsky a luat-o în considerare pentru lucrarea sa.

lm_fev1_m <- lm(Baseline_FEV1_L ~ Age2 + Height + I(Height^2) + Weight, data = d)
seg_fev1_m <- segmented(lm_fev1_m, seg.Z = ~ Age2, psi = 20^2)

gam_fev1_m <- gamlss(Baseline_FEV1_L ~ log(Height) + log(Age) + cs(Age, df = 3),
  sigma.formula = ~ log(Age) + cs(Age, df = 3), family = BCCGo(mu.link = "log"), 
  data = d, method = mixed(50, 500),
  control = gamlss.control(trace = TRUE, maxit = 500))

Pentru a arăta că SLR cu variația constantă pe bucăți nu oferă o potrivire satisfăcătoare pentru mai multe grupuri de vârstă, oferim vizualizarea de mai jos. Acesta descrie proporțiile empirice ale observațiilor sub limita inferioară de 5% a normalului (LLN) la diferite grupuri de vârstă (intervale de 7,5 ani). Linia roșie cu triunghiuri corespunde modelului SLR, iar linia albastră cu cercuri corespunde GAMLSS. Intervalele corespunzătoare de 95% sunt adăugate cu gri în fiecare panou.

În timp ce proporția de observații sub LLN corespunde aproximativ cu nivelul nominal de 5% pentru GAMLSS în toate intervalele de vârstă, SLR se descurcă destul de slab. Pentru tinerii (cu vârsta cuprinsă între 5 și 12,5 ani) LLN este mult prea mic, astfel încât foarte puține observații scad efectiv sub LLN. În schimb, pentru persoanele cu vârsta cuprinsă între 12,5 și 20 de ani, LLN este prea mare și, prin urmare, prea multe observații sunt sub LLN. Cu alte cuvinte, SLR-ul nu ar reuși să indice valori extrem de scăzute ale FEV1, fie producând prea multe, fie prea puține cazuri.

Folosind procast() functia de la topmodels prezicem pentru fiecare persoană LLN corespunzătoare de 5% și apoi cumulăm proporțiile depășite în intervale de vârstă de 7,5 ani.

## proportion of observations below LLN
d_lln <- rbind(
  data.frame(
    proportion = d$Baseline_FEV1_L < procast(gam_fev1_m, type = "quantile", at = 0.05, drop = TRUE),
    fit = "gamlss",
    age = d$Age),
  data.frame(
    proportion = d$Baseline_FEV1_L < procast(seg_fev1_m, type = "quantile", at = 0.05, drop = TRUE),
    fit = "segmented",
    age = d$Age)
)

## set up age intervals
breaks <- seq(5, 80, by = 7.5)
mid <- (head(breaks, -1) + tail(breaks, -1))/2
d_lln <- transform(d_lln,
  age = cut(age, breaks, include.lowest = TRUE)
)

## aggregate in age groups
a_lln <- d_lln |>
  aggregate(x = proportion ~ fit + age, FUN = mean) |>
  transform(nobs = aggregate(proportion ~ fit + age, data = d_lln, FUN = length)$proportion) |>
  transform(alpha = 0.05) |>
  transform(se = sqrt(alpha * (1 - alpha)/nobs)) |>
  transform(lower = alpha + qnorm(0.025) * se) |>
  transform(upper = alpha + qnorm(0.975) * se) |>
  transform(age = mid(age))

Pentru a afișa aceste valori agregate într-un afișaj atrăgător, folosim pachetul recent tinyplot.

library("tinyplot")
pal <- palette.colors()(6:7)
tinytheme("clean2", grid.lwd = 1, grid.lty = 2)
tinyplot(
  proportion ~ age | fit,
  data = a_lln,
  type = "o",
  pch = c(19, 17),
  lwd = 1.5,
  cex = 1.5,
  col = pal,
  ylim = c(0.01, 0.09),
  yaxl = "%",
  xlab = "Age (mid points of 7.5 year intervals)",
  ylab = "Proportion of observations below LLN"
)
tinyplot_add(
  alpha ~ age | fit,
  ymin = a_lln$lower,
  ymax = a_lln$upper,
  type = "ribbon",
  col = "gray"
)
tinyplot_add()

Pentru a evidenția în mod clar că nu numai LLN-ul de 5% este specificat greșit, ci și că, în general, SLR oferă o potrivire probabilistică slabă pentru FEV1, folosim o diagramă QQ cu tendințe (cunoscută și sub numele de diagramă de vierme) bazată pe scorurile z (cunoscute și ca reziduuri cuantile). Folosind qqrplot() functia de la topmodels aceasta poate fi configurată după cum urmează.

qqrplot(seg_fev1_m, detrend = TRUE, confint = "polygon",
  col = pal(2), pch = 17, main = "")
qqrplot(gam_fev1_m, detrend = TRUE, plot = FALSE) |>
  points(col = pal(1), pch = 19)

În scrisoare, diagramele QQ sunt utilizate fără declinare și sunt considerate separat pentru diferite grupuri de vârstă.

Consultați codul de replicare complet (legat mai sus) pentru mai multe detalii.