Rの関数:sens.slope {trend}

Rの関数から sens.slope {trend} を確認します。

関数 `sens.slope` とは

sens.slope は、時系列データのトレンドの傾き（勾配）を推定するための Sen’s Slope（センの傾き） を計算する関数です。

通常、トレンドの傾きを求める際には最小二乗法（OLS）による線形回帰が用いられますが、OLSは外れ値（異常値）の影響を受けやすいという弱点があります。

一方、Sen’s Slope はノンパラメトリックな手法であり、データ対の傾きの「中央値」を採用するため、外れ値に対して非常に頑健（ロバスト）であるという特徴を持ちます。また、データが正規分布に従う必要もありません。

内部アルゴリズムの概要

全てのデータの組み合わせ \((x_i, x_j)\) (\(i < j\)) について、2点間の傾き \(d_k = (x_j - x_i) / (j - i)\) を計算します。
計算された全ての傾き \(d_k\) の 中央値（median） を Sen’s Slope (\(b_{sen}\)) とします。
信頼区間は、Mann-Kendall検定の統計量 \(S\) の分散に基づき、傾きの順位を利用して算出されます。

関数 `sens.slope` の引数

library(trend)
args(sens.slope)

function (x, conf.level = 0.95) 
NULL

x
- 解析対象となる数値ベクトルです。通常は等間隔で測定された時系列データを指定します。
- 数値型である必要があります。また、コード内の na.fail(x) により、欠損値（NA）が含まれている場合はエラーとなります。
conf.level
- 傾きの信頼区間（Confidence Interval）を計算する際の信頼水準を指定します。
- デフォルト: 0.95（95%信頼区間）。

シミュレーションコード

以下に、sens.slope の特徴である「外れ値への強さ」を確認するためのサンプルデータを用いたシミュレーションコードを示します。

このシミュレーションでは、明確な上昇トレンドを持つデータに対し、意図的に極端な外れ値を混入させます。

その際、通常の線形回帰（lm）と sens.slope がそれぞれどのような傾きを推定するかを比較し、本関数の有用性を検証します。

なお、有意水準は5%とします。

シミュレーションデータの生成

40, 45時点目に極端な負の外れ値（スパイク）を混入させます。

# パッケージの読み込み
library(ggplot2)

# 乱数シードの固定
seed <- 20251227
set.seed(seed)

# 設定: 50時点の時系列データ
n <- 50
time_index <- 1:n

# 真のトレンド: 傾き 2.0
true_slope <- 2.0
true_intercept <- 10.0

# データの生成 (基本トレンド + ノイズ)
y <- true_intercept + true_slope * time_index + rnorm(n, mean = 0, sd = 5)

# 外れ値の混入
# データの40番目と45番目に、トレンドから大きく外れた異常値を設定する
# 通常の回帰分析なら、この値に引っ張られて傾きが小さくなるはずです
outlier_indices <- c(40, 45)
y[outlier_indices] <- y[outlier_indices] - 150 # 本来の値より大幅に小さくする

# データフレーム化
df_sim <- data.frame(
  Time = time_index,
  Value = y
)

cat("--- データ概要 ---\n")
cat(sprintf("観測数: %d\n", n))
cat(sprintf("真の傾き: %.1f\n", true_slope))

# 生成されたサンプルデータの可視化
p1 <- ggplot(data = df_sim, mapping = aes(x = Time, y = Value)) +
  geom_point() +
  theme_minimal() +
  labs(
    title = "生成されたサンプルデータ", x = "時間 (Time)",
    y = "観測値 (Value)"
  )

print(p1)

--- データ概要 ---
観測数: 50
真の傾き: 2.0

推定の実行

# 1. 通常の線形回帰 (OLS)
model_lm <- lm(Value ~ Time, data = df_sim)
slope_lm <- coef(model_lm)["Time"]
intercept_lm <- coef(model_lm)["(Intercept)"]

# 2. Sen's Slope (sens.slope)
res_sens <- sens.slope(df_sim$Value, conf.level = 0.95)
slope_sens <- res_sens$estimates[["Sen's slope"]]

# Sen's Slopeの切片計算
# sens.slope関数自体は「傾き」のみを返すため、
# ここでは、切片を median(y - slope * t) として求めます
intercept_sens <- median(df_sim$Value - slope_sens * df_sim$Time)

# 結果の表示
cat("--- 通常の線形回帰の結果 ---\n")
print(summary(model_lm))

cat("--- Sen's Slopeの結果 ---\n")
print(res_sens)

--- 通常の線形回帰の結果 ---

Call:
lm(formula = Value ~ Time, data = df_sim)

Residuals:
     Min       1Q   Median       3Q      Max 
-132.267   -1.786    5.079   13.055   20.192 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  17.8604     8.2964   2.153   0.0364 *  
Time          1.4740     0.2832   5.206 3.99e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 28.89 on 48 degrees of freedom
Multiple R-squared:  0.3608,    Adjusted R-squared:  0.3475 
F-statistic:  27.1 on 1 and 48 DF,  p-value: 3.989e-06

--- Sen's Slopeの結果 ---

    Sen's slope

data:  df_sim$Value
z = 7.7459, n = 50, p-value = 9.494e-15
alternative hypothesis: true z is not equal to 0
95 percent confidence interval:
 1.779092 2.042651
sample estimates:
Sen's slope 
   1.921817

線形回帰（OLS）の結果

Coefficients: Time Estimate 1.4740
- 真の値（2.0）と比較すると、約 26% も過小評価されています。これは、シミュレーションで混入させた「2つの極端に小さな値（外れ値）」に回帰直線が引っ張られ、傾きが寝てしまった（緩やかになってしまった）ことを意味します。
Residuals (Min -132.267):
- 残差の最小値が大きなマイナス値を示しています。これはモデルが捉えきれない異常値の存在を示唆しています。
Multiple R-squared (0.3608):
- 決定係数が約 0.36 と低く、データの変動の多くを説明できていません。外れ値によるノイズが全体のフィッティングを阻害している結果です。

Sen’s Slopeの結果

sample estimates: Sen’s slope 1.921817
- sens.slope の推定結果は 1.921817 となりました。真の値（2.0）との誤差は 0.08 程度であり、高い精度で本来のトレンドを捉えています。OLSが外れ値に引きずられて 1.47 まで数値を落としたのに対し、Sen’s Slope は中央値ベースの計算を行うことで外れ値を無視し、大多数のデータが示す「真の傾き」を維持しました。
信頼区間 (1.779 ～ 2.042):
- 95%信頼区間の中に、真の値である 2.0 が含まれています。

可視化による比較

# 可視化による比較
# 予測線を引くためのデータ作成
df_sim$Pred_LM <- intercept_lm + slope_lm * df_sim$Time
df_sim$Pred_Sens <- intercept_sens + slope_sens * df_sim$Time
# True Trend = 10.0 + 2.0 * Time
df_sim$True_Trend <- true_intercept + true_slope * df_sim$Time

p2 <- ggplot(df_sim, aes(x = Time, y = Value)) +
  # 1. 実データ（背景の線と点）
  geom_line(color = "gray85") +
  geom_point(size = 2, alpha = 0.5) +

  # 2. 外れ値を強調表示（赤丸）
  geom_point(
    data = df_sim[outlier_indices, ], aes(x = Time, y = Value),
    color = "red", size = 4, shape = 1
  ) +

  # 3. 真のトレンド（黒色・実線・やや透明）
  geom_line(aes(y = True_Trend, color = "真のトレンド (正解)"),
    linewidth = 1.5, alpha = 0.3
  ) +

  # 4. 通常の線形回帰（青色・破線）
  geom_line(aes(y = Pred_LM, color = "線形回帰 (OLS)"),
    linewidth = 1, linetype = "dashed"
  ) +

  # 5. Sen's Slope（緑色・実線）
  geom_line(aes(y = Pred_Sens, color = "Sen's Slope (頑健)"),
    linewidth = 1
  ) +

  # 色と凡例の設定
  scale_color_manual(
    name = "推定手法 / 基準",
    values = c(
      "真のトレンド (正解)" = "black",
      "線形回帰 (OLS)" = "blue",
      "Sen's Slope (頑健)" = "darkgreen"
    )
  ) +
  labs(
    title = "外れ値が存在する場合のトレンド推定比較",
    x = "時間 (Time)",
    y = "観測値 (Value)"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(p2)