使用TMB拟合含变量误差的线性回归时的斜率偏差问题

算力阿峰

2026-7-16

变量误差线性回归的TMB模型斜率估计偏差问题

我是Template Model Builder（TMB）的新手，正在通过拟合含变量误差的简单线性回归完成练习，目标后续将该方法应用到复杂模型。真实斜率为2，普通最小二乘法（OLS）得到的斜率估计值为1.81，TMB模型的估计值虽有改进，但仍低于真实值（95%置信区间为1.92至1.98），怀疑代码存在错误，寻求社区建议。

模拟与分析代码

#
# The following code examines the error in variables case of simple linear
# regression using maximum likelihood.
#
#-------------- Prepare Workspace: clear and load packages --------------------
#
rm(list=ls())
library(MASS)         # for mvnorm()
library(TMB)          # for Template Model Builder (TMB)
library(tidyverse)

#-------------------- Define model parameters ---------------------------------
#
model_name <-"eiv_regression"
set.seed(42)       # random number seed to fix the sequence for replication
b0 <- 10.0         # regression intercept
b1 <- 2.0          # regression slope
sigma_Y <- 5.0     # observation error on Y
lambda <- 1        # observation error ratio: sigma_X/sigma_Y
                   #   sigma_X = lambda * sigma_Y

N <- 1000     # number of observations (use lots to test)

#-------------------- Create Simulated Data Set -------------------------------
#
# simulate observation errors in x and Y
epsilon <- rnorm(N,0,sigma_Y)          # errors in Y
delta  <- rnorm(N,0,(lambda*sigma_Y))  # errors in X
# Simulate observed values of predictor variable
X_true <- runif(N,20,80)
X_obs <- X_true + delta
# Simulate observations as "true" value + Y error
Y_true <- b0 + b1*X_true
Y_obs <- Y_true + epsilon
# Make sure our estimates look like we think they should
plot(X_obs,Y_obs,pch=19,col="grey")
# Add the line estimated if we ignore errors in variables
lm_model <- lm(Y_obs~X_obs)
abline(b0,b1,lwd=3)
abline(coef(lm_model)[1],coef(lm_model)[2], lwd=3, lty=2)
# Print "naive" parameter estimates
iestTable <- as.data.frame(summary(lm_model)$coefficients)
rownames(estTable)[2] <- "X"
iestTable$True <- c(b0,b1)
iestTable <- estTable %>%
  select(True,Estimate,`Std. Error`)
iestTable
# autocorrelation in residuals?
acf(resid(lm_model))

#------------- Define Template Model as .CPP file -----------------------------

cat("
/*---------------------- eiv_regression.cpp --  -----------------------------*/
#include <TMB.hpp>
using Eigen::SparseMatrix;
using namespace density;

template<class Type>
Type objective_function<Type>::operator() ()
{
    // Data
    // Note: Transformations of raw data are performed in R
    DATA_VECTOR(Y);      // Response
    DATA_VECTOR(X);      // Predictor
    DATA_SCALAR(lambda)  // ratio sigma_X/sigma_Y as known

    // Book keeping 
    size_t N = Y.size();       // number of observations
                               // NOTE: use size_t for array indices

    // Parameters
    PARAMETER_VECTOR(beta);    // coefficients (intercept, slope)
    PARAMETER(log_sigma_Y);    // log sd for measurement errors of Y
                               //   logged to force positive values
    PARAMETER_VECTOR(del);     // vector of random observation errors for 
                               //   predictor sigma_X

    // Transform parameters to desired values for use and reporting
    Type sigma_Y = exp(log_sigma_Y);  // inverse log
    Type sigma_X = sigma_Y * lambda;

    // Initialize negative log likelihood (to minimize as optimization goal)
    Type nll = 0.;
    
    // X error: likelihood contribution of observation error in X
    nll -= sum(dnorm(del, Type(0), sigma_X, true));
    
    // PROCESS: predictions calculated as model (with predictor error)
    // Note: TMB array indices start at 0, R arrays start at 1
    vector<Type> Y_hat(N);
    for (size_t n=0; n<N; n++){                   // loop through all obs
      Y_hat(n) = beta(0) + beta(1)*(X(n) - del(n));
    }            
    
    // OBSERVATION: likelihood contribution of observed data calculated as
    //   predictions + observation errors. e_i ~ N(0,sigma_Y)
    for (size_t n=0; n<N; n++){                     // loop through all obs
      nll -= dnorm(Y(n), Y_hat(n), sigma_Y, true);
    }
    
    // Parameters to report to R, including Std. Error estimates based
    // on the Hessian and Delta Method
    ADREPORT(beta);        // regression coefficients
    ADREPORT(sigma_Y);     // sd estimate Y
    ADREPORT(sigma_X);     // sd estimate X
    ADREPORT(Y_hat);       // predicted values and their st.dev.
    
    return nll;            // return the negative log likelihood for the
                           // specific set of parameter values.
}
/*---------------------------------------------------------------------------*/
",file="eiv_regression.cpp")

#------------- Estimate model parameters via TMB ------------------------------
#
# Remove previous object file and DLL if present
if (file.exists(paste0(model_name, ".o"))) {
  file.remove(paste0(model_name,".o"))
}
if (file.exists(paste0(model_name, ".dll"))) {
  file.remove(paste0(model_name,".dll"))
}
#
# Compile Model
compile(paste0(model_name, ".cpp"))
#
# Load compiled model
dyn.load(dynlib(model_name))
#
# Define data list
Data <- list(Y=as.vector(Y_obs),       # Observed Response
             X=as.vector(X_obs),       # Observed Predictor
             lambda=lambda)            # sigma_X/sigma_Y assumption

# Define parameter list with initial values
Params <- list(beta=c(0,0),        # Regression coefficients
               log_sigma_Y=0,      # Y Observation Error (log)
               del=rep(0,N))       # errors in X

# Construct TMB function for optimization
Obj <- MakeADFun(data=Data, 
                 parameters=Params, 
                 DLL=model_name,
                 silent=T)
Obj$env$tracemgc <- FALSE
Obj$env$inner.control$trace <- FALSE

# Perform Optimization
Opt <- nlminb(start=Obj$par, objective=Obj$fn, gradient=Obj$gr,
              control=list(iter.max=1000, eval.max=500))
if(Opt$convergence != 0){
  print("Model had partial or non_convergence!")
}
# 
# Likelihood profiles to assess estimation performance
par(mfrow=c(3,1))
nEstPar <- 3
for (i in 1:nEstPar) {
  plot(tmbprofile(Obj,i, trace=0))
}
#
# Produce Report
Report <- sdreport(Obj)
ReportTable <- as.data.frame(summary(Report,"report")[1:4,])
rownames(ReportTable)[1:2] <- c("intercept","slope")
# Add CIs
ReportTable$CI95low <- ReportTable[,1] - ReportTable[,2]*1.96
ReportTable$CI95high <- ReportTable[,1] + ReportTable[,2]*1.96
# Add true values
ReportTable$True <- c(b0,b1,sigma_Y,lambda*sigma_Y)
# Add OLS values
OLS <- summary(lm(Y_obs~X_obs))
ReportTable$OLS <- c(OLS$coefficients[1:2,1],OLS$sigma,NA)
options(digits=3)
ReportTable
#
# Housekeeping: Unload model *.dll when done
dyn.unload(dynlib(model_name))

关键问题排查与修复建议

1. 随机效应的错误处理

代码中把del（X的观测误差）标记为PARAMETER_VECTOR，这是核心错误。del属于潜在随机变量，应该作为随机效应处理，TMB需要对其进行积分以得到边际似然，而非将其当作固定参数估计。

修复步骤：

在CPP文件中，将PARAMETER_VECTOR(del);替换为RANDOM_VARIABLE(del);

在R代码创建Obj对象时，指定随机变量：

Obj <- MakeADFun(data=Data, 
                 parameters=Params, 
                 DLL=model_name,
                 silent=T,
                 random="del")

这样TMB会使用Laplace近似积分掉del，得到正确的边际似然，消除过参数化导致的估计偏差。

2. 初始值优化

当前beta的初始值设为c(0,0)，log_sigma_Y设为0，这可能导致优化收敛到局部最优。建议用OLS的结果作为初始值：

OLS <- summary(lm(Y_obs~X_obs))
Params <- list(beta=coef(lm_model),        # 用OLS系数作为初始值
               log_sigma_Y=log(OLS$sigma), # 用OLS的sigma估计取对数
               del=rep(0,N))

更合理的初始值能提升优化的稳定性和准确性。

3. 代码笔误修复

模拟数据部分存在变量名笔误：

iestTable <- as.data.frame(summary(lm_model)$coefficients)
rownames(iestTable)[2] <- "X" # 原代码中estTable未定义，改为iestTable
iestTable$True <- c(b0,b1)
iestTable <- iestTable %>% # 同样改为iestTable
  select(True,Estimate,`Std. Error`)

修正后才能正确生成OLS估计结果表格。

4. 收敛验证

优化后除了检查Opt$convergence，还可以通过Obj$gr(Opt$par)查看梯度是否接近0（所有元素绝对值应小于1e-6），确认优化是否收敛到全局最优。

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

使用TMB拟合含变量误差的线性回归时的斜率偏差问题

模拟与分析代码 #

关键问题排查与修复建议 #

1. 随机效应的错误处理 #

2. 初始值优化 #

3. 代码笔误修复 #