1 数量化Ⅰ類とは

数量化Ⅰ類は,ダミー変数による重回帰分析と等価である.(青木繁久「数量化Ⅰ類はダミー変数による重解分析である」2005)


2 データ変換

Rのlm()関数はFactor型データの場合,自動的にダミー変数化して解析をする.そこで,数値データで入力されているのであればFactor型に変換する,あるいはデータをFactor型で入力されているデータを読み込み解析することで数量化Ⅰ類(ダミー変数による重回帰分析)を実行することができる.

2.1 Factor型にデータ変換

数値でデータが入力されている場合,以下のcodeでFactor型に変換できる.Rになれていない場合,Excelで置換した方が確実で早いかもしれない.

d1 <- c(4,2,3,5,3,3,4,1,4,3)
d2 <- c(2,1,2,1,2,1,2,2,1,1)
d3 <- c(3,1,2,3,2,1,3,1,3,1)
d <- data.frame(d1 , d2 , d3)
colnames(d) <- c("mitoshi","sex","manzoku")

sexのカテゴリーをMとF,manzokuのカテゴリーをH,M,Lに変換

d$sex <- factor(d$sex , label = c("M" , "F"))  # sex:1=M,2=F
d$manzoku <- factor(d$manzoku , label = c("H" , "M" , "L")) # manzoku:1=H,2=M,3=L

summary()関数で5数要約と度数分布を確認

summary(d)  # summaryで内容を確認,sex,manzokuの度数分布
##     mitoshi    sex   manzoku
##  Min.   :1.0   M:5   H:4    
##  1st Qu.:3.0   F:5   M:2    
##  Median :3.0         L:4    
##  Mean   :3.2                
##  3rd Qu.:4.0                
##  Max.   :5.0


2.2 Factor型データで入力

性別がf,mで,満足度がh,m,lで入力されている場合,このまま解析できる.

dy <- c(4,2,3,5,3,3,4,1,4,3)
ds <- c("f","m","f","m","f","m","f","f","m","m")
dm <- c("l","h","m","l","m","h","l","h","l","h")
dd <- data.frame(dy,ds,dm)
colnames(dd) <- c("見通し","性別","満足度")
table(dd$性別);table(dd$満足度)
## 
## f m 
## 5 5
## 
## h l m 
## 4 4 2


3 解析

3.1 Factor型データで入力したデータで分析

out1 <- lm(見通し ~ 満足度 + 性別 , data = dd)
summary(out1)
## 
## Call:
## lm(formula = 見通し ~ 満足度 + 性別, data = dd)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.750 -0.375  0.125  0.250  0.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   1.5000     0.4082   3.674  0.01040 * 
## 満足度l       2.2500     0.3953   5.692  0.00127 **
## 満足度m       1.5000     0.5590   2.683  0.03638 * 
## 性別m         1.0000     0.4082   2.449  0.04983 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5401 on 6 degrees of freedom
## Multiple R-squared:  0.8491, Adjusted R-squared:  0.7737 
## F-statistic: 11.26 on 3 and 6 DF,  p-value: 0.007072

解析結果をみると,切片は1.5,満足度の係数は,lは2.25,mは1.5,hは0である.性別はmが1.0,fは0である.どのカテゴリーを基準にするかは自動的に決まる.

factor型変数の場合並び順がアルファベット順になるので,そのままだと満足度の基準カテゴリーはhになる.基準カテゴリをmにしたいので,m,l,hの順に,levels()でカテゴリの順番を変える.

dd$満足度 <- factor(dd$満足度, levels = c("m", "l", "h"))
out2 <- lm(見通し ~ 満足度 + 性別 , data = dd)
summary(out2)
## 
## Call:
## lm(formula = 見通し ~ 満足度 + 性別, data = dd)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.750 -0.375  0.125  0.250  0.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.0000     0.3819   7.856 0.000225 ***
## 満足度l       0.7500     0.5103   1.470 0.192039    
## 満足度h      -1.5000     0.5590  -2.683 0.036377 *  
## 性別m         1.0000     0.4082   2.449 0.049825 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5401 on 6 degrees of freedom
## Multiple R-squared:  0.8491, Adjusted R-squared:  0.7737 
## F-statistic: 11.26 on 3 and 6 DF,  p-value: 0.007072

解析結果をみると,切片は3.0,満足度の係数は,lは0.75,hは-1.5,mは0である.性別はmが1.0,fは0である.


予測式に係数を代入してみる.

library(equatiomatic)  ### 以下はRのコードではない ###
extract_eq(out1, use_coefs = TRUE) # 

\[ \operatorname{\widehat{見通し}} = 1.5 + 2.25(\operatorname{満足度}_{\operatorname{l}}) + 1.5(\operatorname{満足度}_{\operatorname{m}}) + 1(\operatorname{性別}_{\operatorname{m}}) \]

extract_eq(out2, use_coefs = TRUE) # 

\[ \operatorname{\widehat{見通し}} = 3 + 0.75(\operatorname{満足度}_{\operatorname{l}}) - 1.5(\operatorname{満足度}_{\operatorname{h}}) + 1(\operatorname{性別}_{\operatorname{m}}) \]


仮に,男性で満足度がM(2)とすると,

\(\hat{y}=2.5+1.5\times1+2.25\times0-1\times0=4\)

\(\hat{y}=1.5+2.25\times0+1.5\times1+1\times1=4\)

予測値の値は同じになることがわかる.


3.2 多重共線性診断:gvif(一般化されたVIF)

満足度(manzoku):3カテゴリー,性別(sex):2カテゴリー

library(car)
vif(out1)
##            GVIF Df GVIF^(1/(2*Df))
## 満足度 1.428571  2        1.093265
## 性別   1.428571  1        1.195229
vif(out2)
##            GVIF Df GVIF^(1/(2*Df))
## 満足度 1.428571  2        1.093265
## 性別   1.428571  1        1.195229

分析結果は,当然であるがVIFの値は同じ.

\(GVIF\):一般化されたVIF

\(GVIF^{(1/(2*Df))}\):調整された一般化標準誤差増加係数 (\(aGVIF\))

ただし,レベルが2を超えるカテゴリカル変数の場合,レベル数を調整して他の変数との比較が可能になる.\(aGVIF\)を使用する場合,判断は経験値の平方根を取る必要がある.\(\sqrt{2.5}\)は懸念され,\(\sqrt{5}\)または\(\sqrt{10}\)(2.2または3.2)は多重共線性あり.\(aGVIF\)乗して元の経験値の判断と比較することもできる.

4 残差評価

round(residuals(out1),2)
##     1     2     3     4     5     6     7     8     9    10 
##  0.25 -0.50  0.00  0.25  0.00  0.50  0.25 -0.50 -0.75  0.50
predict(out1)
##    1    2    3    4    5    6    7    8    9   10 
## 3.75 2.50 3.00 4.75 3.00 2.50 3.75 1.50 4.75 2.50
qqnorm(residuals(out1))
qqline(residuals(out1))


5 変数選択:step()

方法:“both”(ステップワイズ法:変数減増法), “backward”(変数減少法,) “forward”(変数増加法)

情報規準:“AIC”, “AICc”, “AICb1”, “AICb2”, “BIC”, “KIC”, “KICc”, “KICb1”, “KICb2”

デフォルト:backward,AIC

step(out1, direction = "both")
## Start:  AIC=-9.43
## 見通し ~ 満足度 + 性別
## 
##          Df Sum of Sq   RSS     AIC
## <none>                 1.75 -9.4297
## - 性別    1      1.75  3.50 -4.4982
## - 満足度  2      9.45 11.20  5.1333
## 
## Call:
## lm(formula = 見通し ~ 満足度 + 性別, data = dd)
## 
## Coefficients:
## (Intercept)      満足度l      満足度m        性別m  
##        1.50         2.25         1.50         1.00
step(out1, direction = "forward")
## Start:  AIC=-9.43
## 見通し ~ 満足度 + 性別
## 
## Call:
## lm(formula = 見通し ~ 満足度 + 性別, data = dd)
## 
## Coefficients:
## (Intercept)      満足度l      満足度m        性別m  
##        1.50         2.25         1.50         1.00

step()を使った変数選択では,アイテムのうちひとつでもカテゴリーが除かれた場合,アイテム全体を削除する必要がある.


6 正規化スコア(カテゴリースコア)

解析結果は基準になるカテゴリーは0となっている.そこで各カテゴリーに数量を合計が0という条件で配分しすべてのカテゴリーのスコアを計算する.例えば標本1のカテゴリー・ウエイトは,女性で生活満足度は不満足であるから,0と2.25である.各標本のカテゴリー・ウエイトの合計は性別が5,生活満足度が12.00である.標本サイズは10であるから,平均は性別が0.5,生活満足度は1.20である.次に,正規化していない偏回帰係数から平均を引いた値を求める.

<性別>

男性:\(1-0.5=0.5\),女性:\(0-0.5=-0.5\)  

<生活満足度>

満足:\(0-1.20=-1.20\),どちらともいえない:\(1.50-1.20=0.3\),不満足:\(2.25-1.20=1.05\)

<定数項>:定数項は従属変数の単純な平均値,定数項:3.2

上記計算ですべてのカテゴリーの数量を求めることができる.


予測式: \(\hat{y}=3.2+\{(-1.2)\times(満足度H)+0.3\times(満足度M)+1.05\times(満足度L)\}+\{0.5\times(性別男)+(-0.5)\times(性別女)\}\)


7 カテゴリースコアの計算

7.1 Factor型変換

Intercept.h <- mean(d$mitoshi)
Intercept.h
## [1] 3.2
out1$coefficients
## (Intercept)     満足度l     満足度m       性別m 
##        1.50        2.25        1.50        1.00
run.hm <- c(0, out1$coefficients[2:3]) # H:0,M:1.5,L:2.25というベクトルを作成
catesco.hm <- run.hm - c(run.hm%*%(table(d$manzoku)/length(d$manzoku)))# manzokuの平均をrun.hmの各要素から減算
catesco.hm
##         満足度l 満足度m 
##   -1.05    1.20    0.45
run.hs <- c(0 , out1$coefficients[4]) # m:0,f:-1というベクトルを作成
catesco.hs <- run.hs - c(run.hs%*%(table(d$sex)/length(d$sex)))# sexの平均をrun.hsの各要素から減算
catesco.hs
##       性別m 
##  -0.5   0.5

7.2 Factor型で入力

Intercept.n <- mean(dd$見通し)
Intercept.n
## [1] 3.2
out2$coefficients
## (Intercept)     満足度l     満足度h       性別m 
##        3.00        0.75       -1.50        1.00
run.nm <- c(0, out2$coefficients[2:3]) # h:0,l:2.25,m:1.5というベクトルを作成
catesco.nm <- run.nm - c(run.nm%*%(table(dd$満足度)/length(dd$満足度)))# 満足度の平均をrun.nmの各要素から減算
catesco.nm
##         満足度l 満足度h 
##    0.30    1.05   -1.20
run.ns <- c(0 , out2$coefficients[4]) # f:0,m:1というベクトルを作成
catesco.ns <- run.ns - c(run.ns%*%(table(dd$性別)/length(dd$性別)))# 性別の平均をrun.nsの各要素から減算
catesco.ns
##       性別m 
##  -0.5   0.5

Factor型変換,Factor型で入力ともにカテゴリースコアは同じである.


カテゴリースコア/範囲

満足度 H M L 範囲
-1.2 0.3 1.05 2.25
性別
0.5 -0.5 1


8 偏相関係数

yd <- c(4,2,3,5,3,3,4,1,4,3)
ss <- c(-.5,0.5,-0.5,0.5,-0.5,0.5,-0.5,-0.5,0.5,0.5)
ms <- c(1.05 , -1.2 , 0.3 ,1.05 , 0.3 , -1.2 , 1.05 , -1.2 , 1.05 ,-1.2)
ds <- data.frame(yd , ss , ms)
colnames(ds) <- c("見通し","性別","満足度")
ds
##    見通し 性別 満足度
## 1       4 -0.5   1.05
## 2       2  0.5  -1.20
## 3       3 -0.5   0.30
## 4       5  0.5   1.05
## 5       3 -0.5   0.30
## 6       3  0.5  -1.20
## 7       4 -0.5   1.05
## 8       1 -0.5  -1.20
## 9       4  0.5   1.05
## 10      3  0.5  -1.20


8.1 偏相関係数行列の計算:cor2pcor()

r <- cor(ds)
r
##           見通し       性別     満足度
## 見通し 1.0000000  0.1856953  0.8076889
## 性別   0.1856953  1.0000000 -0.2948839
## 満足度 0.8076889 -0.2948839  1.0000000
library(corpcor)
pcor <- cor2pcor(r)
pcor   # 1行目 or 1列目がアイテムの偏相関係数
##           [,1]       [,2]       [,3]
## [1,] 1.0000000  0.7523548  0.9185587
## [2,] 0.7523548  1.0000000 -0.7678689
## [3,] 0.9185587 -0.7678689  1.0000000
#pcor[,1]

カテゴリースコア/範囲/偏相関係数

満足度 H M L 範囲 偏相関係数
-1.2 0.3 1.05 2.25 0.75
性別 範囲 偏相関係数
0.5 -0.5 1 0.92

範囲と偏相関係数をみると,範囲は満足度,偏相関係数は性別の方が大きい.