ムーア・ペンローズ型一般化逆行列を使った重回帰分析 Last modified: Jun 28, 2004
目的
重回帰分析を行う
（R では lm 関数を使うのがよいが，個々の解析結果を取り出すのが面倒なので）
（mreg 関数では，ランク落ちの場合に解が得られないので）
使用法
library(MASS) # R のセッションで，一度だけ行う
mreg2(dat)
引数
dat データ行列。従属変数は最終列におき，それ以外は独立変数と見なす。
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/mreg2.R", encoding="euc-jp")
# 重回帰分析（mreg 関数を 3 カ所書き換えた）
mreg2 <- function(dat) # データ行列
{
dat <- subset(dat, complete.cases(dat)) # 欠損値を持つケースを除く
n <- nrow(dat) # ケース数
nc <- ncol(dat) # 列数
nv <- nc-1 # 独立変数の個数（最後の一つは従属変数）
if (is.null(colnames(dat))) { # 変数名が付いていないときには仮の名前を付ける
colnames(dat) <- paste("Var", 1:nc, sep="")
}
r <- cor(dat) # 相関係数行列
m <- colMeans(dat)[-nc] # 独立変数の平均値ベクトル
# betas <- solve(r[-nc, -nc], r[,nc][-nc]) # 標準化偏回帰係数
betas <- ginv(r[-nc, -nc]) %*% r[,nc][-nc] # 上の 1 行をこれに置き換え
variance <- var(dat)*(n-1) # 変動共変動行列
prop <- diag(variance) # 対角成分
prop <- (prop/prop[nc])[-nc] # 偏回帰係数に変換するための係数（独立変数の変動/従属変数の変動）
b <- betas/sqrt(prop) # 偏回帰係数
Sr <- variance[,nc][-nc]%*%b # 回帰による変動
St <- variance[nc, nc] # 全変動
Se <- St-Sr # 誤差変動
SS <- c(Sr, Se, St) # 平方和（変動）ベクトル
dfr <- nv # 回帰による変動の自由度
dfe <- n-nv-1 # 誤差変動の自由度
dft <- n-1 # 全変動の自由度
df <- c(dfr, dfe, dft) # 自由度ベクトル
Ms <- SS/df # 平均平方ベクトル
f <- Ms[1]/Ms[2] # F 値
fvalue <- c(f, NA, NA)
p <- c(pf(f, dfr, dfe, lower.tail=FALSE), NA, NA) # P 値
b0 <- mean(dat[,nc])-sum(b*m) # 定数項
b <- c(b, b0) # 偏回帰係数ベクトル
# inv <- solve((n-1)*cov(dat)[-nc, -nc])
inv <- ginv((n-1)*cov(dat)[-nc, -nc]) # 上の 1 行をこれに置き換え
SEb <- c(sapply(1:nv, function(i) sqrt(inv[i, i]*Ms[2])), sqrt((1/n+m%*%inv%*%m)*Ms[2]))
tval <- b/SEb # 偏回帰係数の有意性検定
pval <- pt(abs(tval), n-nv-1, lower.tail=FALSE)*2 # P 値
# tolerance <- 1/diag(solve(cor(dat)[-nc, -nc])) # トレランス
tolerance <- 1/diag(ginv(cor(dat)[-nc, -nc])) # 上の 1 行をこれに置き換え
result <- data.frame(b, SEb, tval, pval, c(betas, NA), c(tolerance, NA))
rownames(result) <- c(colnames(dat)[1:nv], "定数項")
colnames(result) <- c("偏回帰係数", "標準誤差", "t 値", "P 値", "標準化偏回帰係数", "トレランス")
R2 <- 1-Se/St # 重相関係数の二乗
R <- sqrt(R2) # 重相関係数
R2s <- 1-Ms[2]/Ms[3] # 自由度調整済み重相関係数の二乗
loglik <- 0.5*(sum(-n*(log(2*pi)+1-log(n)+log(Se)))) # 対数尤度
AIC <- 2*nc+2-2*loglik # AIC
Rs <- c("重相関係数"=R, "重相関係数の二乗"=R2, "自由度調整済重相関係数の二乗"=R2s,
"対数尤度"=loglik, "AIC"=AIC)
anova <- data.frame(SS, df, Ms, fvalue, p) # 分散分析表
rownames(anova) <- c("回帰", "残差", "全体")
colnames(anova) <- c("平方和", "自由度", "平均平方", "F 値", "P 値")
return(list(result=result, anova=anova, Rs=Rs))
}
使用例
dat <- matrix(c( # 3 列目が従属変数
1.2, 1.9, 0.9,
1.6, 2.7, 1.3,
3.5, 3.7, 2.0,
4.0, 3.1, 1.8,
5.6, 3.5, 2.2,
5.7, 7.5, 3.5,
6.7, 1.2, 1.9,
7.5, 3.7, 2.7,
8.5, 0.6, 2.1,
9.7, 5.1, 3.6
), byrow=TRUE, nc=3)
dat2 <- cbind(dat[,1], dat)
library(MASS) # R のセッションで，一度だけ行う（何回も行ってもエラーにはならないが，無駄）
mreg2(dat2)
出力結果例
dat2 は，以下のようになる。つまり，左の 2 列が全く同じ値を取る。
> dat2
[,1] [,2] [,3] [,4]
[1,] 1.2 1.2 1.9 0.9
[2,] 1.6 1.6 2.7 1.3
[3,] 3.5 3.5 3.7 2.0
[4,] 4.0 4.0 3.1 1.8
[5,] 5.6 5.6 3.5 2.2
[6,] 5.7 5.7 7.5 3.5
[7,] 6.7 6.7 1.2 1.9
[8,] 7.5 7.5 3.7 2.7
[9,] 8.5 8.5 0.6 2.1
[10,] 9.7 9.7 5.1 3.6
このような場合，通常のやり方では逆行列が求まらないので解が求まらない。
> mreg(dat2)
Error in drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")) :
Lapack routine dgesv: system is exactly singular
ムーア・ペンローズ型一般逆行列を使用すると解は求まる。
> mreg2(dat2)
$result # Var1 と Var2 に対する偏回帰係数が同じである。
B SE(B) t value P value beta tolerance
Var1 0.1023086 0.004085162 25.04395 2.668281e-07 0.3353364 3.934304
Var2 0.1023086 0.004085162 25.04395 2.668281e-07 0.3353364 3.934304
Var3 0.2866338 0.011666965 24.56798 2.990948e-07 0.6579264 0.983576
Constant 0.1491756 0.058873576 2.53383 4.444936e-02 NA NA
$anova
SS d.f. MS F value P value
Regression 6.67164369 3 2.223881231 470.558 1.655702e-07
Residual 0.02835631 6 0.004726051 NA NA
Total 6.70000000 9 0.744444444 NA NA
$Rs
R R square adj. R sq.
0.9978816 0.9957677 0.9936516
元のデータ dat を使った場合と比較すると，B は半分になっていることがわかる。
> mreg(dat)
$result
B SE(B) t value P value beta tolerance
Var1 0.2046172 0.007564251 27.050556 2.419227e-08 0.6706727 0.983576
Var2 0.2866338 0.010801511 26.536454 2.763814e-08 0.6579264 0.983576
Constant 0.1491756 0.054506340 2.736849 2.904999e-02 NA NA
$anova # 残差の自由度が異なるので，分散分析表の P 値も異なる
SS d.f. MS F value P value
Regression 6.67164369 2 3.335821847 823.4765 4.931875e-09
Residual 0.02835631 7 0.004050901 NA NA
Total 6.70000000 9 0.744444444 NA NA
$Rs
R R square adj. R sq. # R と R square は同じであるが，adjusted R square は異なる
0.9978816 0.9957677 0.9945585
解説ページ
参考文献
上田太一郎著「データマイニング実践集」共立出版株式会社，162 〜 171 ページ
ランク落ちしているデータ行列において重回帰分析を行うときにムーア・ペンローズ型一般逆行列を使う方法が記載されている。上田氏の 著書では，国立筑波技術短期大学の小池将貴教授が作ったムーア・ペンローズ型一般逆行列を求める関数が紹介されている。
ムーア・ペンローズ型一般逆行列を求める関数は，R にもあるので，ここにこの方法を紹介することにした。
数量化I類や数量化II類を重回帰分析や判別分析を用いて行う場合において，ダミー変数の作成時にランク落ちを避けるために冗長なダミー変数を 1 個削除するという，余分な（？）ことをしなくてもすむなどのメリットもある。
直前のページへ戻る E-mail to Shigenobu AOKI
