Rによる判別分析、プロビット、ロジット分析、ポアソン回帰
options(width=100,warn=-1,digits=5,messages=1,tidy=TRUE)
## 倒産企業と非倒産企業の仮想データを準備
# ya1:倒産企業の自己資本比率 ya2:倒産企業のインタレスト
#カバレッジ比率
# yb1:非倒産企業の自己資本比率 yb2:非倒産企業の
ya1<-c(15.5,5.6,28.3,12.7,-78.7)
#インタレストカバレッジ比率
# 便宜的に倒産企業は０ 非倒産企業は１ を識別標識zに
#入力する
ya2<-c(1.14,1.03,-0.42,-8.05,0.58)
yb1<-c(11.4,13.7,30.25,34.06,54.24,52.23,35.61,36.47,59.7,46.8,34.3,33.4,47.7,38.1)
#始めからデータフレーム型でデータを用意すればよいのだ
yb2<-c(2.3,1.93,2.43,5.5,7.76,6.05,6.31,5.2,22.36,18.92,3.21,3.53,3.91,-3.65)
#が、線形判別分析の計算練習で、通常の多変量解析による行
#列演算で得た結果とパッケージMASSのldaを使って得た結果を#比較するために、大変に冗長な手順になるが,わざとベクトル#でデータを準備している。
(za1<-rep(0,length(ya1)))
## [1] 0 0 0 0 0
ya<-cbind(ya1,ya2,za1)
length(ya1)
## [1] 5
colnames(ya)<-c("x1","x2","z")
ya
## x1 x2 z
## [1,] 15.5 1.14 0
## [2,] 5.6 1.03 0
## [3,] 28.3 -0.42 0
## [4,] 12.7 -8.05 0
## [5,] -78.7 0.58 0
(zb1<-rep(1,length(yb1)))
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
length(yb1)
## [1] 14yb<-cbind(yb1,yb2,zb1)
yb
## yb1 yb2 zb1
## [1,] 11.40 2.30 1
## [2,] 13.70 1.93 1
## [3,] 30.25 2.43 1
## [4,] 34.06 5.50 1
## [5,] 54.24 7.76 1
## [6,] 52.23 6.05 1
## [7,] 35.61 6.31 1
## [8,] 36.47 5.20 1
## [9,] 59.70 22.36 1
## [10,] 46.80 18.92 1
## [11,] 34.30 3.21 1
## [12,] 33.40 3.53 1
## [13,] 47.70 3.91 1
## [14,] 38.10 -3.65 1
colnames(yb)<-c("x1","x2","z")
yb
## x1 x2 z
## [1,] 11.40 2.30 1
## [2,] 13.70 1.93 1
## [3,] 30.25 2.43 1
## [4,] 34.06 5.50 1
## [5,] 54.24 7.76 1
## [6,] 52.23 6.05 1
## [7,] 35.61 6.31 1
## [8,] 36.47 5.20 1
## [9,] 59.70 22.36 1
## [10,] 46.80 18.92 1
## [11,] 34.30 3.21 1
## [12,] 33.40 3.53 1
## [13,] 47.70 3.91 1
## [14,] 38.10 -3.65 1
##以下で倒産企業と非倒産企業の説明変数全部の行列が作ら
#れる
yall
## x1 x2 z
## [1,] 15.50 1.14 0
## [2,] 5.60 1.03 0
## [3,] 28.30 -0.42 0
## [4,] 12.70 -8.05 0
## [5,] -78.70 0.58 0
## [6,] 11.40 2.30 1
## [7,] 13.70 1.93 1
## [8,] 30.25 2.43 1
## [9,] 34.06 5.50 1
## [10,] 54.24 7.76 1
## [11,] 52.23 6.05 1
## [12,] 35.61 6.31 1
## [13,] 36.47 5.20 1
## [14,] 59.70 22.36 1
## [15,] 46.80 18.92 1
## [16,] 34.30 3.21 1
## [17,] 33.40 3.53 1
## [18,] 47.70 3.91 1
## [19,] 38.10 -3.65 1
#####################################
## パッケージMASS ldaによる線形判別分析
suppressPackageStartupMessages(library(MASS))
dyall<-as.data.frame(yall)
dyall
## x1 x2 z
## 1 15.50 1.14 0
## 2 5.60 1.03 0
## 3 28.30 -0.42 0
## 4 12.70 -8.05 0
## 5 -78.70 0.58 0
## 6 11.40 2.30 1
## 7 13.70 1.93 1
## 8 30.25 2.43 1
## 9 34.06 5.50 1
## 10 54.24 7.76 1
## 11 52.23 6.05 1
## 12 35.61 6.31 1
## 13 36.47 5.20 1
## 14 59.70 22.36 1
## 15 46.80 18.92 1
## 16 34.30 3.21 1
## 17 33.40 3.53 1
## 18 47.70 3.91 1
## 19 38.10 -3.65 1
(ZZZ<- lda(z~ . ,data=dyall))
## Call:
## lda(z ~ ., data = dyall)
##
## Prior probabilities of groups:
## 0 1
## 0.26316 0.73684
##
## Group means:
## x1 x2
## 0 -3.320 -1.1440
## 1 37.711 6.1257
##
## Coefficients of linear discriminants:
## LD1
## x1 0.033300
## x2 0.071224
(CC<-apply(ZZZ$means%*%ZZZ$scaling,2,mean))
## LD1
## 0.75002#z=0.0333＊(自己資本比率)＋0.0712＊（インタレストカバレ
#ッジ）-0.75
#ゼロが判別点となり、上式で得た得点がゼロ以下では倒産、
#ゼロより大きけば非倒産と判別される。定数項は倒産企業の
#各変数の平均値と非倒産企業の各説明変数の平均値との中点
#を式0.0333＊(自己資本比率)＋0.0712＊（インタレストカバ
#レッジ）が通過するように計算されている。
##判別関数によるデータの判別結果
#左の1列目は当初の区分（倒産０ 非倒産１） 2列目は
#モデルによる判別結果
data.frame(dyall[,3],predict(ZZZ)$class)
## dyall...3. predict.ZZZ..class
## 1 0 1
## 2 0 1
## 3 0 1
## 4 0 0
## 5 0 0
## 6 1 1
## 7 1 1
## 8 1 1
## 9 1 1
## 10 1 1
## 11 1 1
## 12 1 1
## 13 1 1
## 14 1 1
## 15 1 1
## 16 1 1
## 17 1 1
## 18 1 1
## 19 1 1
#度数分布表を見ると、判別点はゼロになるが若干この近辺で
# いくつかの企業が混在している
plot(ZZZ,dimen=1)
# 予測的中率をみるために下記のようなテーブルを作ってみ
#ると
table(実際=dyall[,3],予測=predict(ZZZ)$class)
## 予測
## 実際 0 1
## 0 2 3
## 1 0 14
#実際は倒産企業である会社をモデル予測では非倒産と判定し
#ている例が３つあるが、非倒産企業についてはすべて的中さ
#せていることが分かる
###################
### Probit 分析の適用
##計算過程で警告が出るが計算手順の概略を示すことが目的な#のでそれを無視してすすめる。
suppressPackageStartupMessages(library(glm.predict))
## Error in library(glm.predict): there is no package
#called 'glm.predict'
tousanprobit <- glm( z ~ x1 + x2,data = dyall, family = binomial(link = "probit"))
summary(tousanprobit)
##
## Call:
## glm(formula = z ~ x1 + x2, family = binomial(link = "probit"),
## data = dyall)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5934 0.0000 0.0026 0.0669 1.1268
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1174 1.8853 -1.12 0.26
## x1 0.1020 0.0824 1.24 0.22
## x2 0.4480 0.4313 1.04 0.30
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 21.9007 on 18 degrees of freedom
## Residual deviance: 7.8578 on 16 degrees of freedom
## AIC: 13.86
##
## Number of Fisher Scoring iterations: 11
#プロビットモデルによる各社の確率
#ここでは0.5を境界として０に近いほど倒産確率は高く、
#逆に１に近いほど非倒産の確率が高いと判断する
round(predict(tousanprobit,type="response"),digits=2)
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 0.49 0.14 0.72 0.00 0.00 0.53 0.56 0.98 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.55
plot(fitted(tousanprobit))
table(実際 = dyall$z,
予測 = round(fitted(tousanprobit)))
## 予測
## 実際 0 1
## 0 4 1
## 1 0 14
#この例では倒産企業の的中率がたまたま悪いが、実際の分析
#では大量のデータを使うのでどの手法がよいかは分からな
#い。
#ロジットモデルの適用
summary(tousanlogit)
##
## Call:
## glm(formula = z ~ x1 + x2, family = binomial(link = "logit"),
## data = dyall)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6326 -0.0001 0.0284 0.1184 1.1052
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.626 3.347 -1.08 0.28
## x1 0.176 0.147 1.20 0.23
## x2 0.779 0.745 1.05 0.30
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 21.9007 on 18 degrees of freedom
## Residual deviance: 7.9328 on 16 degrees of freedom
## AIC: 13.93
##
## Number of Fisher Scoring iterations: 9
#ここでも0.5を境界として０に近いほど倒産確率は高く、
#逆に１に近いほど非倒産の確率が高いと判断する
round(predict(tousanlogit,type="response"),digits=2)
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 0.50 0.14 0.74 0.00 0.00 0.54 0.57 0.97 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.99 0.99 1.00 0.56
table(実際 = dyall$z,予測 = round(fitted(tousanlogit)))
## 予測
## 実際 0 1
## 0 4 1
## 1 0 14
#この例では倒産企業の的中率がよくなっているが、実際の分
#析では大量のデータを使うのでなんともいえない。
#ロジットモデルの結果を3次元グラフで示すと以下のようにな
#る。ライブラリーlibrary(rockchalk)を使うと素人でもグラ
#フを描きやすい。
suppressPackageStartupMessages(library(rockchalk))
plotPlane(tousanlogit, plotx1 = "x1", plotx2 = "x2")
plotPlane(tousanlogit, plotx1 = "x1", plotx2 = "x2", phi = 20,theta=60)
##ポアソン回帰の適用
#ポアソン回帰モデルは一定期間の事故発生件数、特許申請件
#数、犯罪発生件数など発生頻度の少ない事象をポアソン分布
#を仮定してモデル化されている。
#倒産についていえば一定期間に同一企業が何回も倒産するこ
#とはないのでカウントデータは０と１だけで２以上のケース
#は存在しないことになるが、幸にも,ポアソン分布の期待値
#（ここでは平均倒産回数）が小さい場合は回数が
#２以上の確率は無視できるくらい小さくなることが知られて
#おり、２値データ#モデルを十分近似できるといわれている。
tousanpo <- glm(z ~ x1 + x2,data = dyall, family = poisson)
summary(tousanpo)
##
## Call:
## glm(formula = z ~ x1 + x2, family = poisson, data = dyall)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.143 -0.331 -0.013 0.274 0.718
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.079844 0.702754 -1.54 0.12
## x1 0.023138 0.020934 1.11 0.27
## x2 0.000526 0.045798 0.01 0.99
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8.5507 on 18 degrees of freedom
## Residual deviance: 5.7137 on 16 degrees of freedom
## AIC: 39.71
##
## Number of Fisher Scoring iterations: 5
round(fitted(tousanpo))
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 0 0 1 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1
table(実際 = dyall$z, 予測 = round(fitted(tousanpo)))
## 予測
## 実際 0 1
## 0 4 1
## 1 2 12
#ポアソン回帰で倒産か非倒産かの二者択一をさせている。
###########################################################
##ここからは全くの蛇足になるが、行列演算による判別分析の#係数等の推定値とライブラリによる係数のアウトプットの違
#いについて検討を加えてみた。
#多変量解析の教科書には線形代数を使った計算過程が示され
#ているので自己学習のために、上記と同じサンプルについて
#判別関数の係数を計算してみた。
#各郡の平均値の差 のベクトル
d1
## [1] -41.031
d2<-mean(ya2)-mean(yb2)
d2
## [1] -7.2697
D<-c(d1,d2)
D
## [1] -41.0314 -7.2697
# 平均偏差ベクトル
Sa1<-(ya1-mean(ya1))
Sa1
## [1] 18.82 8.92 31.62 16.02 -75.38
Sa2<-(ya2-mean(ya2))
Sa2
## [1] 2.284 2.174 0.724 -6.906 1.724
# 郡Aの行列
SSa<-cbind(Sa1,Sa2)
SSa
## Sa1 Sa2
## [1,] 18.82 2.284
## [2,] 8.92 2.174
## [3,] 31.62 0.724
## [4,] 16.02 -6.906
## [5,] -75.38 1.724
SSa<-as.matrix(SSa,nrow(SSa),ncol(SSa))
SSa
## Sa1 Sa2
## [1,] 18.82 2.284
## [2,] 8.92 2.174
## [3,] 31.62 0.724
## [4,] 16.02 -6.906
## [5,] -75.38 1.724
# 郡１の偏差平方和 行列
# 郡１の共分散行列 サンプル数－１ 不偏推定量
acov<-SSasqd/(nrow(SSa)-1)
acov
## Sa1 Sa2
## Sa1 1843.09 -38.830
## Sa2 -38.83 15.283
#郡２の偏差ベクトル
Sb1<-(yb1-mean(yb1))
Sb1
## [1] -26.31143 -24.01143 -7.46143 -3.65143 16.52857 14.51857 -2.10143 -1.24143 21.98857
## [10] 9.08857 -3.41143 -4.31143 9.98857 0.38857
Sb2<-(yb2-mean(yb2))
Sb2
## [1] -3.825714 -4.195714 -3.695714 -0.625714 1.634286 -0.075714 0.184286 -0.925714 16.234286
## [10] 12.794286 -2.915714 -2.595714 -2.215714 -9.775714
SSb<-cbind(Sb1,Sb2)
SSb
## Sb1 Sb2
## [1,] -26.31143 -3.825714
## [2,] -24.01143 -4.195714
## [3,] -7.46143 -3.695714
## [4,] -3.65143 -0.625714
## [5,] 16.52857 1.634286
## [6,] 14.51857 -0.075714
## [7,] -2.10143 0.184286
## [8,] -1.24143 -0.925714
## [9,] 21.98857 16.234286
## [10,] 9.08857 12.794286
## [11,] -3.41143 -2.915714
## [12,] -4.31143 -2.595714
## [13,] 9.98857 -2.215714
## [14,] 0.38857 -9.775714
# 郡２の行列
SSb<-as.matrix(SSb,nrow(SSb),ncol(SSb))
SSb
## Sb1 Sb2
## [1,] -26.31143 -3.825714
## [2,] -24.01143 -4.195714
## [3,] -7.46143 -3.695714
## [4,] -3.65143 -0.625714
## [5,] 16.52857 1.634286
## [6,] 14.51857 -0.075714
## [7,] -2.10143 0.184286
## [8,] -1.24143 -0.925714
## [9,] 21.98857 16.234286
## [10,] 9.08857 12.794286
## [11,] -3.41143 -2.915714
## [12,] -4.31143 -2.595714
## [13,] 9.98857 -2.215714
## [14,] 0.38857 -9.775714
# 郡２の偏差平方和 行列
SSbsqd<-t(SSb)%*%SSb
SSbsqd
## Sb1 Sb2
## Sb1 2524.0 726.40
## Sb2 726.4 592.82
# 郡２の共分散行列 サンプル数－１ 不偏推定量
bcov<-SSbsqd/(nrow(SSb)-1)
bcov
## Sb1 Sb2
## Sb1 194.156 55.877
## Sb2 55.877 45.601
#郡１と郡２のプールされた分散行列 各郡の
#（サンプル数－１）での加重平均
Spl<-acov*(nrow(SSa)-1)/(nrow(SSa)+nrow(SSb)-2)+bcov*(nrow(SSb)-1)/(nrow(SSa)+nrow(SSb)-2)
Spl
## Sa1 Sa2
## Sa1 582.141 33.593
## Sa2 33.593 38.468
#プールされた分散行列の逆行列
invSpl<-solve(Spl)
invSpl
## Sa1 Sa2
## Sa1 0.0018090 -0.0015797
## Sa2 -0.0015797 0.0273755
#上記の逆行列に平均差ベクトルＤ を乗じ 係数を求める
(alpha<-invSpl%*%D)
## [,1]
## Sa1 -0.06274
## Sa2 -0.13419#この係数はライブラリーMASSのldaを使って求めた数値
#x1 0.03329955 x2 0.07122402 と一見すると全く異なって
#いる。
# しかし、両者の係数ベクトルを作り正規化して比べると同一となっている。
# ldaで計算した係数ベクトルを作り正規化する
V1norm<-sqrt(t(V1)%*%V1)
(VV1<-V1/V1norm[1])
## LD1
## x1 0.42353
## x2 0.90588
## 同様に行列計算で求めた係数を正規化すると
aa<-c(alpha[1:2])
(V<-aa)
## [1] -0.06274 -0.13419
Vnorm<-sqrt(t(V)%*%V)
(VV<-V/Vnorm[1])
## [1] -0.42353 -0.90588
##実質的に同じ結果を得る。正方行列の固有値λに対する固有
#ベクトルのｋ倍もλに対する固有ベクトルとなるということ示
#している。
#念のためldaで求めた係数ベクトルを行列計算による係数ベク
#トルで除してみると、以下のように、すべて同一倍率になっ
#ている。
ZZZ$scaling/aa
## LD1
## x1 -0.53076
## x2 -0.53076