如何用hnp包为nls拟合的非线性模型生成半正态图

使用R语言stats包的nls()函数拟合非线性模型，相关代码与数据集如下：

模型拟合代码 #

nonlinear_func <- function(x, beta_0, beta_1) {
  return (beta_0 * (1 - exp(beta_1 * x)) + 1.30)
}

nonlinear_model <- nls(y ~ nonlinear_func(x, beta_0, beta_1), 
            start = list(beta_0 = 16, beta_1 = -0.2), 
            data = data.set)

数据集 #

data.set <- structure(list(x = c(13.05, 6.05, 13.21, 9.55, 18.14, 9.55, 14.48, 
15.28, 9.87, 15.92, 12.41, 12.41, 12.57, 15.12, 10.66, 16.87, 
12.57, 15.92, 9.71, 15.92, 17.35, 6.37, 11.94, 11.14, 8.91, 13.05, 
17.67, 10.66, 17.19, 7, 10.82, 11.62, 16.71, 18.3, 11.78, 12.89, 
10.82, 9.23, 14.32, 7.64, 5.09, 15.44, 10.35, 8.91, 14.32, 13.21, 
8.91, 15.6, 14.16, 15.28, 12.57), y = c(16.8, 11.3, 16.9, 11.8, 
17.5, 15.8, 20, 19.1, 15.8, 18.8, 18.6, 18.4, 18.4, 18.7, 15.2, 
18.3, 15.7, 17.5, 14.8, 16.7, 19.8, 10.3, 18.6, 14.4, 14.3, 17.8, 
21, 18.8, 18.9, 13.7, 17.6, 17.5, 19.6, 18.8, 15.3, 17, 15.9, 
13.3, 17.3, 13.6, 9.3, 17.7, 14.2, 14.9, 18.4, 18.2, 14.3, 19.7, 
18.6, 18.1, 15.5)), row.names = c(17L, 18L, 19L, 20L, 21L, 22L, 
23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 
36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 
49L, 50L, 51L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 
12L, 13L, 14L, 15L, 16L), class = "data.frame")

尝试使用hnp包绘制带模拟包络的半正态图，编写了自定义函数但生成的图无模拟包络：

d.fun <- function(obj) resid(obj)

s.fun <- function(n, obj) {
}

f.fun <- function(data) {
  nls(y ~ nonlinear_func(x, beta_0, beta_1), 
      start = list(beta_0 = 16, beta_1 = -0.2), 
      data = data.set)
}
library(hnp)
hnp(nonlinear_model, newclass = TRUE, 
diagfun = d.fun, simfun = s.fun, 
fitfun = f.fun, data = data.set)

要求仅使用hnp包，不得使用car包的qqPlot函数。

出现无模拟包络的问题有两个核心原因：

1. s.fun函数为空：该函数负责生成模拟响应值，是构建模拟包络的核心，空函数无法提供模拟残差；
2. f.fun函数硬编码数据集：没有使用传入的模拟数据集参数，导致每次拟合都用原始数据，无法生成不同的模拟残差。

修改后的完整代码如下：

# 定义残差提取函数：从模型对象中提取残差
d.fun <- function(obj) resid(obj)

# 定义模拟函数：生成n组符合模型假设的模拟响应值
s.fun <- function(n, obj) {
  # 获取模型拟合值和残差标准差
  mu <- fitted(obj)
  sigma <- sqrt(sum(resid(obj)^2) / df.residual(obj))
  # 生成n组正态分布模拟残差，计算模拟响应值
  replicate(n, mu + rnorm(length(mu), mean = 0, sd = sigma))
}

# 定义拟合函数：使用传入的数据集重新拟合模型
f.fun <- function(data) {
  # 从原模型提取系数作为起始值，提升拟合稳定性
  start_vals <- list(beta_0 = coef(obj)[1], beta_1 = coef(obj)[2])
  nls(y ~ nonlinear_func(x, beta_0, beta_1), 
      start = start_vals, 
      data = data)
}

library(hnp)
# 调用hnp函数生成带模拟包络的半正态图
hnp(nonlinear_model, newclass = TRUE, 
    diagfun = d.fun, simfun = s.fun, 
    fitfun = f.fun, data = data.set, n.sim = 200)

关键细节说明 #

• s.fun函数：基于模型的拟合值和残差标准差，生成服从正态分布的模拟残差，进而得到模拟响应值，这是模拟包络的数据源；
• f.fun函数：改用传入的data参数，并且从原模型提取系数作为起始值，避免固定起始值导致的拟合失败问题；
• n.sim参数：可调整模拟次数，默认200次，次数越多包络线越平滑，但计算时间会相应增加。

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