Rの関数:nnls {nnls}

Rの関数から nnls {nnls} を確認します。

nnls関数とは

nnls関数は、非負最小二乗法（Non-Negative Least Squares）を実装した関数です。これは、制約付きの最小二乗問題の一種です。

通常の線形回帰問題は、与えられた行列Aとベクトルbに対して、残差二乗和 ||Ax - b||^2 を最小化するような係数ベクトルxを求めます。しかし、問題によっては、係数xの各要素が物理的な量（例えば、物質の濃度や構成比率など）を表すため、負の値を取ることが許されない場合があります。

nnls関数は、そのような全ての係数 x_i が x_i >= 0 を満たすという制約条件下で、残差二乗和 ||Ax - b||^2 を最小化するxを求めるために使用されます。

(注意)

単回帰式を表す場合、説明変数をx、係数をAやbとして、 y = Ax + b と表記しますが、ここでは回帰分析を b = Ax として表記します。

  • A:

    計画行列 (Design Matrix) と呼ばれる行列です。この行列Aが、説明変数のデータセットに相当します。行列の各列がひとつの説明変数（例: A1, A2…）、各行がひとつの観測サンプルを表します。

  • x:

    係数ベクトル (Coefficient Vector) です。このベクトルがが、モデルが推定しようとしている未知の係数の集まりです。nnls関数の目的は、このベクトルxの最適な値を求めることです。

  • b:

    応答変数ベクトル (Response Vector) です。予測したい対象である目的変数の観測値の集まりです。

nnls関数の引数

library(nnls)
args(nnls)
function (A, b) 
NULL

  • A:

    これは計画行列（Design Matrix）であり、説明変数を行ごとに、サンプルを列ごとに配置した数値行列です。回帰分析における説明変数のデータに相当します。m x n の次元を持ち、mは観測数（標本サイズ）、nは説明変数の数を表します。

  • b:

    これは応答変数（Response Vector）であり、予測したい対象の観測値を集めた数値ベクトルです。回帰分析における目的変数に相当します。このベクトルの長さ m は、行列 A の行数と一致しなければなりません。

シミュレーション

以下に、nnls関数を確認するためのシミュレーションコードを示します。

ここでは、二つの説明変数(A1とA2)がほぼ同一という、極端な多重共線性状況を設定したサンプルデータを作成します。

さらに、真の係数に0と正の値を混ぜて設定し、通常の最小二乗法（OLS）の結果と比較することで、nnlsの非負制約の有効性を確認します。

# シミュレーション結果の再現性を確保するために乱数の種を設定します
seed <- 20251115
set.seed(seed)

# 観測数（m）と説明変数の数（n）を設定します
m <- 80
n <- 5

# 基となる説明変数 A1 を作成します
A1 <- runif(m, 10, 20)

# A1 とほぼ同一の説明変数 A2 を作成します
# A2 は A1 に微小なノイズを加えたものです
A2 <- A1 + rnorm(m, mean = 0, sd = 1e-4)

# 他の独立した説明変数を作成します
A3 <- runif(m, 5, 15)
A_other <- matrix(runif(m * (n - 3)), nrow = m)

# 計画行列 A を結合して作成します
A <- cbind(A1, A2, A3, A_other)
colnames(A) <- paste0("A", 1:n)


# 真の係数ベクトル x_true を作成します
# ほぼ同一の変数 A1, A2 のうち、A1 のみに効果があるとします
x_true <- c(10, 0, 5, 0, 0)

# 応答変数ベクトル b を作成します
noise <- rnorm(m, sd = 1)
b <- A %*% x_true + noise

# nnls関数とOLSによるモデル推定
# nnls関数による推定
cat("--- nnls関数による推定 ---\n")
nnls_result <- nnls(A, b)
print(nnls_result)
--- nnls関数による推定 ---
Nonnegative least squares model
x estimates: 9.98637 0 5.013002 0 0 
residual sum-of-squares: 55.61
reason terminated: The solution has been computed sucessfully.
# 通常の最小二乗法（OLS）による推定
cat("--- lm関数による推定 ---\n")
ols_result <- lm(b ~ A - 1)
print(summary(ols_result))
--- lm関数による推定 ---

Call:
lm(formula = b ~ A - 1)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.95693 -0.45884 -0.00636  0.55686  1.86854 

Coefficients:
      Estimate Std. Error t value Pr(>|t|)    
AA1  148.77309  856.05741   0.174   0.8625    
AA2 -138.76130  856.05721  -0.162   0.8717    
AA3    5.02562    0.02688 186.984   <2e-16 ***
AA4   -0.43529    0.31357  -1.388   0.1692    
AA5   -0.64219    0.28001  -2.293   0.0246 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8218 on 75 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 9.825e+05 on 5 and 75 DF,  p-value: < 2.2e-16
cat("--- 推定結果の比較 ---\n")

# 結果をデータフレームにまとめます
results_df <- data.frame(
  True_Coefficients = x_true,
  NNLS_Coefficients = coef(nnls_result),
  OLS_Coefficients = coef(ols_result)
)

print(round(results_df, 3))
--- 推定結果の比較 ---
    True_Coefficients NNLS_Coefficients OLS_Coefficients
AA1                10             9.986          148.773
AA2                 0             0.000         -138.761
AA3                 5             5.013            5.026
AA4                 0             0.000           -0.435
AA5                 0             0.000           -0.642

1. 通常の最小二乗法（OLS）に見られる問題点

まず、OLS_Coefficientsの列、特にA1とA2の行を確認します。

  • A1の係数は 148.773 という大きな正の値、一方でA2の係数は -138.761 という大きな負の値になっています。
  • この現象は、ほぼ同一の変数がモデルに投入されたことで発生する典型的な多重共線性の弊害の一つです。OLSは残差の二乗和を最小化するために、互いに打ち消し合うような極端な係数を算出してしまっています。このような解は不安定であり、少しデータが変わるだけで係数が大きく変動する可能性があります。

2. nnlsによる問題の解決と、その利点

次に、NNLS_Coefficientsの列を確認します。

  • nnlsは、A1の係数を 9.986 と推定し、これは真の係数である 10 に近い値です。一方で、A2の係数は 0.000 となっています。
  • この結果は、「応答変数の変動は、主に変数A1によって説明され、A2は（A1と重複するため）寄与しない」と解釈できます。
  • 全ての係数が0以上でなければならないという非負制約により、先のOLSのように、大きな正の係数を大きな負の係数で打ち消すことは出来ず、nnlsは相関の高い変数群の中から、結果を最もよく説明する変数を探索し、他の冗長な変数の係数を0にするという、一種の変数選択のような挙動を示しています。

以上です。