林の数量化
Kosugitt
Friday, December 19, 2014
はじめに
今回のお題は林の数量化理論です。
非線形なモデリングというシリーズでやってきましたが，ここではそもそも数字ではないものに数字を割り当てる，という意味で「数量化」と呼ばれるモデルを一般に扱います。
ちなみに，これらのモデルは林知己夫という日本の偉大なメトリシャンがつくりあげたモデル群であり，その業績をたたえて名前が冠されています。林先生自身は統計パッケージもなかった時代に，必要に迫られてその都度データに沿ったモデリングを行っただけであり，最初から理論化を考えて作られたわけではありません。
数量化理論には様々なモデルがあり，これらを分類してI類，II類，III類，IV類，V類，VI類，と名付けられています（飽戸弘による）。今振り返って考えてみると，これらのモデルは 名義尺度変数 に対するアプローチであったと言えます。
今回はこれらのうち，I～III類について解説します。
ついでに
今回は，群馬大学の青木先生のWeb site(http://aoki2.si.gunma-u.ac.jp/R/) を参考にしながら進めます。インターネット環境を整えた上で読み進めてください。
数量化I類
数量化I類は，回帰分析の一種です。その中でも特に，独立変数が名義尺度水準になっているものを指します。 早速やってみましょう。
# データの準備
sample <- read.csv("sample_utf8.csv",head=T,na.strings="*")
sample$gender <- factor(sample$gender,label=c("male","female"))
qt1.dat <- subset(sample,select=c("class","gender"))
# 青木先生の数量化I類のサイトからソースコードを取ってくる。
source("http://aoki2.si.gunma-u.ac.jp/R/src/qt1.R", encoding="euc-jp")
result.qt1 <- qt1(qt1.dat,y=sample$height)## [1] "class" "gender"print(result.qt1)## カテゴリースコア
## class.A 1.11415
## class.B -1.66564
## class.C 0.51773
## gender.male 6.60073
## gender.female -6.60073
## 定数項 151.38350plot(result.qt1)
plot(result.qt1,which="fitness")
実はこの方法は，lm関数をつかった普通の回帰分析でもできます。 Rの回帰分析は独立変数がカテゴリカルな場合，自動的に，数量化のような数値の割り当てという考え方を応用してくれるのです。
result.lm <- lm(height~gender+class,data=sample)この時の結果は，例えば性別の場合，最初のラベル(male)を0としたとき，２つ目のラベル(female)が-13.2015である，という表現になっています。数量化のgender.male,gender.femaleのスコアも13.20146(=6.60073-(-6.60073))離れていることを確認して下さい。
数量化II類
数量化の２つ目のモデルも回帰分析の一種ですが，説明される従属変数のほうが名義尺度である場合です。
# データの準備
sample$height.c <- ifelse(sample$height>157,1,
ifelse(sample$height>150,2,3))
sample$height.c <- factor(sample$height.c,label=c("high","mid","short"))
sample$weight.c <- ifelse(sample$weight>70,1,
ifelse(sample$weight>55,2,3))
sample$weight.c <- factor(sample$weight.c,label=c("heavy","fat","light"))
qt2.dat <- subset(sample,select=c("gender","height.c","weight.c"))
# 青木先生の数量化II類のサイトからソースコードを取ってくる。
source("http://aoki2.si.gunma-u.ac.jp/R/src/qt2.R", encoding="euc-jp")
result.qt2 <- qt2(qt2.dat[,2:3],group=qt2.dat[,1])
print(result.qt2)##
## カテゴリー・スコア
##
## 解1
## height.c.high -1.48357
## height.c.mid -0.28644
## height.c.short 0.88501
## weight.c.heavy -0.02271
## weight.c.fat -0.04466
## weight.c.light 0.04861plot(result.qt2,which="category.score")
このモデルは判別分析Descriminant Analysisとも呼ばれているものと同じです。
library(MASS)
result.lda <- lda(gender~height+weight,data=sample)
print(result.lda)## Call:
## lda(gender ~ height + weight, data = sample)
##
## Prior probabilities of groups:
## male female
## 0.5 0.5
##
## Group means:
## height weight
## male 157.9624 61.3864
## female 144.8046 52.2994
##
## Coefficients of linear discriminants:
## LD1
## height -0.18456560
## weight 0.02973659prediction.lda <- predict(result.lda)
plot(prediction.lda$x,col=sample$gender)
plot(sample$gender,prediction.lda$x)
数量化III類
数量化III類は，回帰分析のような外的な変数を持たず，因子分析のように変数間関係からその構造を明らかにする方法です。 早速やってみましょう。
#データの準備。前回のタイタニックデータ。
data(Titanic)
Titanic## , , Age = Child, Survived = No
##
## Sex
## Class Male Female
## 1st 0 0
## 2nd 0 0
## 3rd 35 17
## Crew 0 0
##
## , , Age = Adult, Survived = No
##
## Sex
## Class Male Female
## 1st 118 4
## 2nd 154 13
## 3rd 387 89
## Crew 670 3
##
## , , Age = Child, Survived = Yes
##
## Sex
## Class Male Female
## 1st 5 1
## 2nd 11 13
## 3rd 13 14
## Crew 0 0
##
## , , Age = Adult, Survived = Yes
##
## Sex
## Class Male Female
## 1st 57 140
## 2nd 14 80
## 3rd 75 76
## Crew 192 20tmp <- data.frame(Titanic)
Titanic.df <- data.frame(
Class = rep(tmp$Class, tmp$Freq),
Sex = rep(tmp$Sex, tmp$Freq),
Age = rep(tmp$Age, tmp$Freq),
Survived = rep(tmp$Survived, tmp$Freq)
)
summary(Titanic.df)## Class Sex Age Survived
## 1st :325 Male :1731 Child: 109 No :1490
## 2nd :285 Female: 470 Adult:2092 Yes: 711
## 3rd :706
## Crew:885head(Titanic.df)## Class Sex Age Survived
## 1 3rd Male Child No
## 2 3rd Male Child No
## 3 3rd Male Child No
## 4 3rd Male Child No
## 5 3rd Male Child No
## 6 3rd Male Child No# 青木先生の数量化III類のサイトからソースコードを取ってくる。
source("http://aoki2.si.gunma-u.ac.jp/R/src/qt3.R", encoding="euc-jp")
# ダミー変数を作るコードも使う
source("http://aoki2.si.gunma-u.ac.jp/R/src/make.dummy.R", encoding="euc-jp")
result.qt3 <- qt3(Titanic.df)
print(result.qt3)## 解1 解2 解3 解4 解5 解6
## 固有値 0.44508 0.30504 0.25001 0.20504 0.17852 0.11632
## 相関係数 0.66714 0.55231 0.50001 0.45281 0.42251 0.34105
##
## カテゴリースコア
##
## 解1 解2 解3 解4 解5 解6
## Class.1st -1.72668 -2.22959 -1.78000 0.10939 3.45776 -0.01730
## Class.2nd -0.97619 0.45721 4.78409 1.11605 1.07466 -0.66424
## Class.3rd -0.19576 1.93742 -1.31812 1.28287 -0.72882 -0.87405
## Class.Crew 1.10462 -0.87402 0.16455 -1.42298 -1.03447 0.91753
## Sex.Male 0.64092 -0.00439 0.00504 -0.20310 0.32260 -0.72796
## Sex.Female -2.36051 0.01616 -0.01856 0.74801 -1.18813 2.68107
## Age.Child -1.95131 5.32791 0.04765 -6.01163 2.68136 1.11569
## Age.Adult 0.10167 -0.27760 -0.00248 0.31323 -0.13971 -0.05813
## Survived.No 0.76367 0.34444 -0.00750 0.46849 0.58877 0.80044
## Survived.Yes -1.60038 -0.72183 0.01572 -0.98178 -1.23385 -1.67743plot(result.qt3)
ここでのプロトに見られるように，関係の深い変数同士は近くに，関係の浅い変数は遠くになるように数値が割り振られます。
興味深いことに，これは変数のプロットだけでなくケースのプロットも行います。
plot(result.qt3,which="sample.score")
この方法は，双対尺度法や対応分析と呼ばれる方法と同じものです。
library(MASS)
Titanic.dmy <- make.dummy(Titanic.df)
result.corresp <- corresp(Titanic.dmy,nf=3)
biplot(result.corresp)
まとめ
これらの手法は，数値化されていないものに数値を割り振るという考え方で分析するものです。 逆に言えば，測定されている(名義尺度水準以上，つまり言語化されている)ものであれば，その構造を見たり，説明に使ったりできるということです。
- 数量化I類は名義尺度で予測するモデル
- 数量化II類は予測される変数が名義尺度の分類モデル
- 数量化III類は名義尺度の主成分分析（因子分析）モデル
課題
Merry Christmas!