gls {nlme}のスクラッチ再現

###############################################
# 一般化最小二乗法（GLS）のスクラッチ再現と
# nlme::gls() 関数の結果比較シミュレーション
###############################################

# ----------------------------------------------------------
# 【目的】
# 誤差項に相関構造（AR(1)型）を持つサンプルデータを生成し、
# 1. 既知の共分散行列 Σ を用いた「理論通りのGLS推定（スクラッチ実装）」
# 2. Rの関数 gls() による推定結果
# が一致することを確認する。
# ----------------------------------------------------------

# ----------------------------------------------------------
# 1. 必要なパッケージの読み込み
# ----------------------------------------------------------
library(MASS) # 多変量正規乱数生成 (mvrnorm)
library(nlme) # 一般化最小二乗法 gls() の関数を提供

# ----------------------------------------------------------
# 2. 乱数シード設定（再現性確保）
# ----------------------------------------------------------
seed <- 20251019
set.seed(seed)

# ----------------------------------------------------------
# 3. モデル設定：サンプルサイズと誤差構造
# ----------------------------------------------------------
n <- 100 # 標本サイズ
rho <- 0.6 # 誤差のAR(1)型相関係数（隣接誤差間の相関）
sigma2 <- 1 # 誤差の分散（全体スケール）

# AR(1)型共分散行列 Σ を構築
# Σ_ij = σ^2 * ρ^|i−j|
Sigma <- sigma2 * rho^abs(outer(1:n, 1:n, "-"))
# → outer(1:n, 1:n, "-") は i-j の行列を作る。
#   その絶対値を取りρのべき乗を計算し、σ^2を掛けて共分散行列を得る。

cat("--- 生成した共分散行列Σの一部を確認 ---\n")
print(Sigma[1:7, 1:7])

# ----------------------------------------------------------
# 4. 説明変数（X）と真のパラメータ設定
# ----------------------------------------------------------
x1 <- runif(n, 0, 10) # 一様乱数：説明変数1
x2 <- runif(n, 0, 5) # 一様乱数：説明変数2

# デザイン行列（切片を含む）
X <- cbind(1, x1, x2)

# 真の回帰係数（母数）
beta_true <- c(2, 0.5, -0.3)
# y = 2 + 0.5*x1 - 0.3*x2 + 誤差

# ----------------------------------------------------------
# 5. 相関構造を持つ誤差項 ε の生成
# ----------------------------------------------------------
eps <- MASS::mvrnorm(1, mu = rep(0, n), Sigma = Sigma)
# → 平均0、共分散Σに従う多変量正規乱数ベクトル

cat("\n\n--- 生成した相関構造を持つ誤差項 ε の一部を確認 ---\n")
print(str(eps))

# ----------------------------------------------------------
# 6. 応答変数 y の生成
# ----------------------------------------------------------
y <- X %*% beta_true + eps
# → 線形モデル y = Xβ + ε の生成

# データフレーム化
df <- data.frame(y, x1, x2, group = factor(rep(1, n)))
# 注意：
#   gls() の correlation 引数では「グループ単位で相関構造を推定」する。
#   全観測が同一系列に属するように group=1 を指定するのがポイント。
#   （id=1:n のようにしてしまうと、各観測が別系列と解釈され相関が無効化される）

# ----------------------------------------------------------
# 7. GLSのスクラッチ実装（既知Σを使用）
# ----------------------------------------------------------
# 一般化最小二乗法の理論式：
# β̂_GLS = (X' Σ^{-1} X)^{-1} X' Σ^{-1} y
Sigma_inv <- solve(Sigma)
beta_gls <- solve(t(X) %*% Sigma_inv %*% X) %*% t(X) %*% Sigma_inv %*% y

# このβ_glsは「Σを既知として正確にGLS推定を行った結果」であり、
# 理論上は最尤推定値と一致する。

# ----------------------------------------------------------
# 8. Rのgls()関数による推定
# ----------------------------------------------------------
model_gls <- gls(
  y ~ x1 + x2, # 回帰式
  data = df, # データフレーム
  correlation = corAR1( # 誤差構造：AR(1)モデルを指定
    value = rho, # 相関係数ρを既知として固定
    form = ~ 1 | group, # 1つのグループ（系列）内でAR(1)を仮定
    fixed = TRUE # ρを推定せず、指定値を固定
  ),
  method = "ML" # 最尤法（REMLでもほぼ同じだが、比較を厳密にするためML）
)

# gls() は誤差共分散構造をcorrelationオブジェクトから内部的に構築し、
# GLS推定を行う関数である。
# ここでは、理論上と同じρを固定しているため、
# スクラッチ実装の結果とほぼ一致するはずである。

# ----------------------------------------------------------
# 9. 結果の比較出力
# ----------------------------------------------------------
cat("\n\n--- スクラッチ実装と関数glsによる係数の比較 ---\n")
print(cbind(
  beta_true         = beta_true, # 真のパラメータ値
  beta_gls_scratch  = as.vector(beta_gls), # 手計算（スクラッチ）によるGLS推定
  beta_gls_function = coef(model_gls) # gls()関数による推定値
))