使用R中simex结合nls+幂均值方法时遇测量误差报错求助

错误提示 #

Error: measurement.error is constant 0 in column(s) 1

场景说明 #

用R的nls()拟合了带异方差校正（拟合值幂函数加权）的非线性模型，尝试用simex包校正自变量的测量误差，明明指定了非零测量误差，却触发上述错误。

复现代码 #

library(simex)
set.seed(123456789)
x = runif(n = 1000, min = 1, max = 3.6)
x_err = x + rnorm(n = 1000, mean = 0, sd = 0.1)
y_mean = 100/(1+10^(log10(100)-x)*0.75)
y_het = y_mean + rnorm(n = 1000, mean = 0, sd = 10*x^-2)
y_het = ifelse(y_het > 0, y_het, 0)
w = (100/(1+10^(log10(100)-x_err)*0.75))^-2

nls_fit = nls(y_het ~ 100/(1+10^((log10(k)-x_err)*h)), start = list("k" = 100, "h" = 0.75), weights = 1/w)

simex(nls_fit, SIMEXvariable = "x_err", measurement.error = 0.1, asymptotic = F)

原因与解决方法 #

错误原因 #

simex无法从带权重的nls模型中正确识别自变量的测量误差参数，加权操作改变了模型内部的变量结构，导致工具误判测量误差为0。

解决方法 #

方法1：分离加权与测量误差校正步骤

先拟合无权重的nls模型，用simex完成测量误差校正后，再手动应用异方差加权：

# 拟合无权重的基础模型
nls_fit_unweighted = nls(y_het ~ 100/(1+10^((log10(k)-x_err)*h)), start = list("k" = 100, "h" = 0.75))
# 用simex校正测量误差
simex_result = simex(nls_fit_unweighted, SIMEXvariable = "x_err", measurement.error = 0.1, asymptotic = F)

# 基于校正后的参数重新计算权重并拟合
corrected_params = coef(simex_result)
y_fit_corrected = 100/(1+10^((log10(corrected_params["k"])-x_err)*corrected_params["h"]))
w_corrected = y_fit_corrected^-2
final_fit = nls(y_het ~ 100/(1+10^((log10(k)-x_err)*h)), start = corrected_params, weights = 1/w_corrected)

方法2：用nlme处理异方差（更兼容）

nlme包的gnls()函数可以原生建模异方差结构，且simex对gnls模型的支持更好，无需拆分步骤：

library(nlme)
# 拟合带异方差的gnls模型（用拟合值的幂函数定义方差）
gnls_fit = gnls(y_het ~ 100/(1+10^((log10(k)-x_err)*h)), 
                start = list("k" = 100, "h" = 0.75),
                weights = varPower(form = ~fitted(.)))
# 直接应用simex校正测量误差
simex_gnls = simex(gnls_fit, SIMEXvariable = "x_err", measurement.error = 0.1, asymptotic = F)

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