VBA回归代码与Excel内置回归结果存小数差异致p-value错误，如何对齐？

我编写的数组VBA回归代码，输出结果与Excel内置回归函数在第12位小数后出现差异，进而导致p-value错误。注：R语言的回归输出与Excel回归函数结果一致。

测试数据 #

现有VBA代码 #

Sub regression1()

Dim y, x, rsq2, rsq3, SeCo, tstats, Pvalues As Variant
Dim Ubr As Long, Ubc As Long
Dim ii, jj, P, p1, p2, intp2 As Integer
Dim m, rsq1, rvar1 As Variant
Dim Coeff, rvar, rsq, rvar2, rvar3 As Variant
Dim results As Variant

yi = Range("B1:B13")

xi = Range("D1:F13")

Ubr = UBound(xi, 1)
Ubc = UBound(xi, 2)

ReDim x(1 To Ubr - 1, 1 To Ubc + 1)
ReDim y(1 To Ubr - 1, 1 To 1)
ReDim results(1 To 2, 1 To Ubc + Ubc + 4)


    For ii = 1 To Ubr - 1
        x(ii, 1) = 1
    Next ii
    
    For ii = 2 To Ubr
        For jj = 2 To Ubc + 1
            x(ii - 1, jj) = xi(ii, jj - 1)
        Next jj
    Next ii
    
    For ii = 2 To Ubr
        y(ii - 1, 1) = yi(ii, 1)
    Next ii
    
    
    With Application
     m = .MInverse(.MMult(.Transpose(x), x))
     Coeff = .MMult(m, .MMult(.Transpose(x), y))
     rvar1 = Evaluate(.MMult(.Transpose(y), .MMult(x, Coeff)))
     rvar2 = (UBound(x, 1) - UBound(x, 2))
     rvar3 = .SumSq(y)
     rvar = (rvar3 - rvar1) / rvar2
     rsq1 = .MMult(x, Coeff)
     
     ReDim rsq2(1 To Ubr - 1, 1)
     For ii = 1 To Ubr - 1
      rsq2(ii, 1) = y(ii, 1) - rsq1(ii, 1)
     Next ii
     
     ReDim rsq3(1 To Ubr - 1, 1)
     For ii = 1 To Ubr - 1
      rsq3(ii, 1) = y(ii, 1) - .Average(y)
     Next ii
     
     rsq = 1 - (.SumSq(rsq2) / .SumSq(rsq3))
     
     ReDim SeCo(1 To UBound(m, 1), 1 To 1)
     For ii = 1 To UBound(m, 1)
      SeCo(ii, 1) = (rvar * m(ii, ii)) ^ 0.5
     Next ii
     
     ReDim tstats(1 To UBound(m, 1), 1 To 1)
     For ii = 1 To UBound(m, 1)
      tstats(ii, 1) = Coeff(ii, 1) / SeCo(ii, 1)
     Next ii
    End With
            
    With WorksheetFunction
    
    ReDim Pvalues(1 To UBound(m, 1), 1 To 1)
    For ii = 1 To UBound(m, 1)
      Pvalues(ii, 1) = .T_Dist_2T(Abs(tstats(ii, 1)), UBound(y, 1) - 3)
     Next ii
     
    results(1, 1) = "Variables"
    results(1, 2) = "R-Square"
    results(1, 3) = "P-value (intercept)"
    
    For iii = 1 To Ubc
        results(1, iii + 3) = "P-Value" & iii
    Next iii
    
    results(1, 3 + Ubc + 1) = "B0 (Intercept)"
    
    For iii = 1 To Ubc
        results(1, iii + 4 + Ubc) = "B" & iii
    Next iii


Dim Name As Variant
Dim Fname As Variant

Fname = "Inter"
    For iii = 1 To Ubc
    
        Name = xi(1, iii)
    
        Fname = .Concat(Fname, " + ", Name)
    
    Next iii
    
    results(2, 1) = Fname
    
    results(2, 2) = rsq
    
    For iii = 1 To Ubc + 1
    
        results(2, iii + 2) = Pvalues(iii, 1)
    
    Next iii
    
        For iii = 1 To Ubc + 1
    
        results(2, iii + 3 + Ubc) = Coeff(iii, 1)
    
    Next iii
    
    
    End With

Range("I2", Range("I2").Offset(UBound(results, 1) - 1, UBound(results, 2) - 1)) = results

Range("S22", Range("S22").Offset(UBound(m, 1) - 1, 0)) = tstats

Range("V22", Range("V22").Offset(UBound(m, 1) - 1, 0)) = SeCo

End Sub

问题根源 #

1. 矩阵运算精度差异：手动用MInverse+MMult实现回归的方式，和Excel内置回归采用的QR分解/Cholesky分解算法精度不同，微小的浮点误差会在后续计算中被放大，最终导致p-value偏差
2. 自由度硬编码：p-value计算中直接写死UBound(y, 1) - 3，没有根据自变量数量动态计算，当自变量数量变化时会出错
3. 冗余循环操作：手动循环计算残差平方和的过程，额外引入了精度损耗

调整方案 #

1. 改用Excel内置LINEST函数（最可靠） #

LINEST是Excel原生回归实现，直接调用它可以完美匹配Excel数据分析工具的结果，彻底避免手动矩阵运算的精度问题。核心修改如下：

' 替换原有的矩阵计算代码段
With Application
    ' 调用LINEST获取完整回归结果，参数：y值数组、x值数组（含截距项）、是否包含截距、是否返回统计量
    Dim linestResults As Variant
    linestResults = .LinEst(y, x, True, True)
    
    ' 提取关键统计量
    Coeff = .Index(linestResults, 1) ' 系数（注意LINEST返回顺序是从最后一个自变量到截距）
    SeCo = .Index(linestResults, 2) ' 标准误
    rsq = .Index(linestResults, 3, 1) ' R²
    df = .Index(linestResults, 4, 2) ' 自由度（LINEST直接返回，无需手动计算）
    
    ' 计算t统计量
    ReDim tstats(1 To UBound(Coeff), 1 To 1)
    For ii = 1 To UBound(Coeff)
        tstats(ii, 1) = Coeff(ii, 1) / SeCo(ii, 1)
    Next ii
End With

2. 修正p-value的自由度计算 #

将硬编码的自由度改为动态获取的值：

With WorksheetFunction
    ReDim Pvalues(1 To UBound(Coeff), 1 To 1)
    For ii = 1 To UBound(Coeff)
        Pvalues(ii, 1) = .T_Dist_2T(Abs(tstats(ii, 1)), df)
    Next ii
    ' 其余代码保持不变
End With

3. 调整系数顺序（重要） #

LINEST返回的系数顺序是从最后一个自变量到截距，需要反转顺序才能和Excel回归结果一致：

' 反转Coeff、SeCo、tstats、Pvalues的顺序，匹配Excel输出
Dim tempCoeff As Variant
ReDim tempCoeff(1 To UBound(Coeff), 1 To 1)
For ii = 1 To UBound(Coeff)
    tempCoeff(ii, 1) = Coeff(UBound(Coeff) - ii + 1, 1)
Next ii
Coeff = tempCoeff

' 同理反转SeCo、tstats、Pvalues，代码略

完整修改后的代码 #

Sub regression1()

Dim y, x, SeCo, tstats, Pvalues As Variant
Dim Ubr As Long, Ubc As Long
Dim ii, iii As Integer
Dim Coeff, rsq, df As Variant
Dim results As Variant
Dim yi, xi As Variant
Dim Name, Fname As Variant

yi = Range("B1:B13")
xi = Range("D1:F13")

Ubr = UBound(xi, 1)
Ubc = UBound(xi, 2)

ReDim x(1 To Ubr - 1, 1 To Ubc + 1)
ReDim y(1 To Ubr - 1, 1 To 1)
ReDim results(1 To 2, 1 To Ubc + Ubc + 4)

' 构造含截距项的X矩阵
For ii = 1 To Ubr - 1
    x(ii, 1) = 1
Next ii
For ii = 2 To Ubr
    For jj = 2 To Ubc + 1
        x(ii - 1, jj) = xi(ii, jj - 1)
    Next jj
Next ii
For ii = 2 To Ubr
    y(ii - 1, 1) = yi(ii, 1)
Next ii

' 改用LINEST计算回归结果
With Application
    Dim linestResults As Variant
    linestResults = .LinEst(y, x, True, True)
    
    ' 提取统计量
    Coeff = .Index(linestResults, 1)
    SeCo = .Index(linestResults, 2)
    rsq = .Index(linestResults, 3, 1)
    df = .Index(linestResults, 4, 2)
    
    ' 计算t统计量
    ReDim tstats(1 To UBound(Coeff), 1 To 1)
    For ii = 1 To UBound(Coeff)
        tstats(ii, 1) = Coeff(ii, 1) / SeCo(ii, 1)
    Next ii
    
    ' 反转系数顺序，匹配Excel输出（LINEST返回顺序是x3、x2、x1、截距）
    Dim tempCoeff, tempSeCo, tempTstats As Variant
    ReDim tempCoeff(1 To UBound(Coeff), 1 To 1)
    ReDim tempSeCo(1 To UBound(SeCo), 1 To 1)
    ReDim tempTstats(1 To UBound(tstats), 1 To 1)
    For ii = 1 To UBound(Coeff)
        tempCoeff(ii, 1) = Coeff(UBound(Coeff) - ii + 1, 1)
        tempSeCo(ii, 1) = SeCo(UBound(SeCo) - ii + 1, 1)
        tempTstats(ii, 1) = tstats(UBound(tstats) - ii + 1, 1)
    Next ii
    Coeff = tempCoeff
    SeCo = tempSeCo
    tstats = tempTstats
End With

' 计算p-value并整理结果
With WorksheetFunction
    ReDim Pvalues(1 To UBound(Coeff), 1 To 1)
    For ii = 1 To UBound(Coeff)
        Pvalues(ii, 1) = .T_Dist_2T(Abs(tstats(ii, 1)), df)
    Next ii
    
    ' 构造结果表头
    results(1, 1) = "Variables"
    results(1, 2) = "R-Square"
    results(1, 3) = "P-value (intercept)"
    For iii = 1 To Ubc
        results(1, iii + 3) = "P-Value" & iii
    Next iii
    results(1, 3 + Ubc + 1) = "B0 (Intercept)"
    For iii = 1 To Ubc
        results(1, iii + 4 + Ubc) = "B" & iii
    Next iii
    
    ' 构造变量名称
    Fname = "Inter"
    For iii = 1 To Ubc
        Name = xi(1, iii)
        Fname = .Concat(Fname, " + ", Name)
    Next iii
    results(2, 1) = Fname
    
    ' 填充结果
    results(2, 2) = rsq
    For iii = 1 To Ubc + 1
        results(2, iii + 2) = Pvalues(iii, 1)
    Next iii
    For iii = 1 To Ubc + 1
        results(2, iii + 3 + Ubc) = Coeff(iii, 1)
    Next iii
End With

' 输出结果到工作表
Range("I2", Range("I2").Offset(UBound(results, 1) - 1, UBound(results, 2) - 1)) = results
Range("S22", Range("S22").Offset(UBound(tstats, 1) - 1, 0)) = tstats
Range("V22", Range("V22").Offset(UBound(SeCo, 1) - 1, 0)) = SeCo

End Sub

内容的提问来源于stack exchange，提问作者Mah e fatima