餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

コーシー分布も扱えるランク基準線形モデル推定(Rank-based Estimation for Linear Models)

Kloke and McKean (2012) "Rfit: Rank-based Estimation for Linear Models," The R Journal, Vol.4(2)にアルゴリズムの概要と利用方法の詳細が書いてあるので説明する必要はないのですが、U検定やBrunner-Munzel検定を使ってきた外れ値がある観測値がコーシー分布らしいデータに線形回帰や分散分析をかけるのに便利なパッケージがありました。

ランク基準線形モデル推定は、ランク*1を用いた線形回帰で一致推定量を計算する手法で、OLSでは一致推定量が得られない誤差項がコーシー分布である場合なども、頑強な結果を示してくれます。ノンパラ手法の紹介ページで存在を知ったのですが、カーネル回帰のように大きなサンプルサイズが必要になるわけではなく、U検定やBrunner-Munzel検定を使う状況で、統制変数を加えてみたいときなどに気軽に使えそうです。

誤差項がコーシー分布

誤差項がコーシー分布の場合は、OLSよりも真のパラメーター(i.e. 切片項が1,xの係数が2,zの係数が-3)に近い推定量を出してきます。

# データセットを作成
df01 <- with(list(n = 90), {
	set.seed(520)
	x <- rnorm(n, sd = 2)
	z <- runif(n, min = -1, max = 1)
	g <- 1:n %% 3
	gcoef <- c(0, 0.05, -0.05)
	y_c <- 1 + 2*x - 3*z  + rcauchy(n, scale = 1)
	y_n <- 1 + 2*x - 3*z  + rnorm(n)
	data.frame(y_c, y_n, x, z, g = as.factor(g))
})

# 誤差項がコーシー分布のときの推定
library("Rfit")
r_rfit_c <- rfit(y_c ~ x + z + g, data = df01) # 推定
summary(r_rfit_c) # 推定結果の確認
summary(lm(y_c ~ x + z + g, data = df01)) # OLSの場合


誤差項が正規分布

誤差項が正規分布の場合は、OLSとだいたい同じ結果です。

r_rfit_n <- rfit(y_n ~ x + z + g, data = df01) # 推定
summary(r_rfit_n) # 推定結果の確認
summary(lm(y_n ~ x + z + g, data = df01)) # OLSの場合


モデル選択

モデル選択はlmとanovと同様にできます。

drop.test(r_rfit_c, rfit(y_c ~ x + z, data = df01)) # モデル選択


分散分析

一元配置分散分析はoneway.rfitで、n元配置分散分析はraovでできます。

r_oneway_anova <- with(df01, oneway.rfit(y, g)) # 一元配置分散分析
summary(r_oneway_anova, method = "tukey") # 多重比較の補正をして結果表示

ぽちっとな感がある実用的なソリューションですね。

*1:誤差項のベクトルをRのrank関数で計算した値。