正如Gustavo Thomas 等人在使用 GAMLSSs 的纵向多级实验分析一文中所解释的那样，整体和特定于预测器的蠕虫图具有“不同的形状表示模型中的不同不足”的特征： https ://arxiv.org/pdf /1810.03085.pdf。
灵活回归和平滑一书的第 12.4 节：在 R 中使用 GAMLSS。里格比等人。值得一读，因为它全面介绍了如何解释蠕虫图。该部分以这些陈述结束：“一般来说，建立一个没有不合适区域的模型可能并不总是可能的。” “无论如何，当使用具有许多不合适区域的模型来支持结论时，需要格外小心。”。然而，校准被提到作为一种解决方案，以最大限度地减少失配。
如何纠正模型不匹配取决于蠕虫图中检测到的问题的性质。如果这些问题表明需要考虑您的连续预测变量之一的非线性效应以改善模型拟合，那么您将需要对该预测变量的效应进行非线性建模而不是线性建模。（其他类型的校正可能涉及给定模型中的预测变量和随机效应，为响应变量指定不同类型的分布，省略或包括模型各个部分的预测变量，转换预测变量等）
请注意，根据 cs() 函数的帮助文件：
函数 scs() 与函数 cs() 的不同之处在于它允许对平滑参数进行交叉验证，这与固定有效自由度 df 的 cs() 不同。请注意，推荐的平滑函数现在是函数 pb()，它允许使用局部最大似然估计平滑参数。函数 pb() 基于 Eilers 和 Marx (1996) 的惩罚 beta 样条 (P-splines)。
因此，您可能需要考虑在模型中使用 pb() 而不是 cs()。
附录：
这是一些 R 代码，用于为模型生成数据，其中二次拟合比线性甚至平滑拟合更好。它将帮助您对蠕虫图的外观建立一些直觉。数据是根据https://www.theanalysisfactor.com/r-tutorial-4/生成的。
14, 15, 16, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 29, 30),
Outcome = c(126.6, 101.8, 71.6, 101.6, 68.1, 62.9, 45.5, 41.9,
46.3, 34.1, 38.2, 41.7, 24.7, 41.5, 36.6, 19.6,
22.8, 29.6, 23.5, 15.3, 13.4, 26.8, 9.8, 18.8, 25.9, 19.3)),
.Names = c("Time", "Outcome"),
row.names = c(1L, 2L, 3L, 5L, 7L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 19L, 20L, 21L, 22L, 23L, 25L, 26L, 27L, 28L, 29L, 30L, 31L),
class = "data.frame")
Data
数据的标题如下所示：
Time Outcome
1 0 126.6
2 1 101.8
3 2 71.6
5 4 101.6
7 6 68.1
9 8 62.9```
The plot of the Outcome variable versus the predictor variable Time can be obtained with:
```library(ggplot2)
theme_set(theme_bw())
ggplot(Data, aes(x = Time, y = Outcome)) +
geom_point(size=3, colour="dodgerblue")
现在，在 gamlss 框架内为这些数据拟合 3 个可能的模型：
linear.model <- gamlss(Outcome ~ Time, data = Data, family=NO)
quadratic.model <- gamlss(Outcome ~ Time + I(Time^2), data = Data, family=NO)
smooth.model <- gamlss(Outcome ~ pb(Time), data = Data, family=NO)
summary(linear.model)
summary(quadratic.model)
summary(smooth.model)
比较 3 个拟合模型的（广义）AIC 值：
GAIC(linear.model, quadratic.model, smooth.model)
二次模型是“赢家”，因为它具有最小的 AIC 值：
df AIC
quadratic.model 4.000000 197.0357
smooth.model 5.251898 197.8349
linear.model 3.000000 219.0893
现在为时间预测器构建蠕虫图：
wp(linear.model, xvar=Time)
wp(quadratic.model, xvar=Time)
wp(smooth.model, xvar=Time)
线性模型拟合的蠕虫图显示了一些失配问题：
二次和平滑模型拟合的蠕虫图看起来比线性模型拟合的蠕虫图好一点。
我们还可以直接根据时间预测器绘制模型残差：
Data$linear.model.residuals <- residuals(linear.model)
Data$quadratic.model.residuals <- residuals(quadratic.model)
Data$smooth.model.residuals <- residuals(smooth.model)
plot1 <- ggplot(Data, aes(x = Time, y = linear.model.residuals)) +
geom_point(size=3, colour="darkgrey") +
geom_hline(yintercept = 0, linetype=2, colour="red") +
ggtitle("Linear Model Residuals vs. Time") +
coord_cartesian(ylim=c(-3,3))
plot2 <- ggplot(Data, aes(x = Time, y = quadratic.model.residuals)) +
geom_point(size=3, colour="darkgrey") +
geom_hline(yintercept = 0, linetype=2, colour="red") +
ggtitle("Quadratic Model Residuals vs. Time") +
coord_cartesian(ylim=c(-3,3))
plot3 <- ggplot(Data, aes(x = Time, y = smooth.model.residuals)) +
geom_point(size=3, colour="darkgrey") +
geom_hline(yintercept = 0, linetype=2, colour="red") +
ggtitle("Smooth Model Residuals vs. Time") +
coord_cartesian(ylim=c(-3,3))
library(cowplot)
plot_grid(plot1, plot2, plot3, ncol=3)
最后这些图使我们更容易辨别线性模型的残差中存在二次模式，这需要在模型中加以考虑。
如果您愿意，您可以将线性模型的残差与时间图分开，并检查对应于相应蠕虫图中使用的时间间隔划分的图部分：
w.linear <- wp(linear.model, xvar=Time, main="Given: Time")
w.linear
在 w.linear 的 R 输出的 $classes 部分中报告了 Time 观察值范围划分的切点：
> w.linear
$classes
[,1] [,2]
[1,] -0.5 8.5
[2,] 8.5 15.5
[3,] 15.5 24.5
[4,] 24.5 30.5
$coef
[,1] [,2] [,3] [,4]
[1,] 0.6061177 0.79644473 0.26190049 -0.29589027
[2,] -1.0467772 -0.54040972 0.08504976 -0.05550396
[3,] -0.1400464 -0.64524770 -0.15331613 0.02095304
[4,] 0.7161490 -0.03070935 -0.08930395 -0.19956330
这些切点是 -0.5、8.5、15.5、24.5 和 30.5。我们可以绘制残差与时间的关系，并仅为“中间”切点绘制垂直线：
plot11 <- ggplot(Data, aes(x = Time, y = linear.model.residuals)) +
geom_point(size=3, colour="darkgrey") +
geom_hline(yintercept = 0, linetype=2, colour="red") +
ggtitle("Linear Model Residuals vs. Time") +
coord_cartesian(ylim=c(-3,3)) +
geom_vline(xintercept = w.linear$classes[1,2],
colour="blue", linetype=3, size=1.5) +
geom_vline(xintercept = w.linear$classes[2,2],
colour="blue", linetype=3, size=1.5) +
geom_vline(xintercept = w.linear$classes[3,2],
colour="blue", linetype=3, size=1.5)
plot11
这使我们可以放大特定的时间间隔并确定模型拟合在这些间隔中如何分解：
