数量化Ⅰ類(ダミー変数による重回帰分析)
©2023 Nobuyuki Otake
最終作成日:2024-02-04
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=Lsummary()関数で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.02.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 23 解析
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.195229vif(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.50predict(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.50qqnorm(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.00step(out1, direction = "forward")## Start: AIC=-9.43
## 見通し ~ 満足度 + 性別##
## Call:
## lm(formula = 見通し ~ 満足度 + 性別, data = dd)
##
## Coefficients:
## (Intercept) 満足度l 満足度m 性別m
## 1.50 2.25 1.50 1.00step()を使った変数選択では，アイテムのうちひとつでもカテゴリーが除かれた場合，アイテム全体を削除する必要がある．
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.2out1$coefficients## (Intercept) 満足度l 満足度m 性別m
## 1.50 2.25 1.50 1.00run.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.45run.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.57.2 Factor型で入力
Intercept.n <- mean(dd$見通し)
Intercept.n## [1] 3.2out2$coefficients## (Intercept) 満足度l 満足度h 性別m
## 3.00 0.75 -1.50 1.00run.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.20run.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.5Factor型変換，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.208.1 偏相関係数行列の計算:cor2pcor()
r <- cor(ds)
r## 見通し 性別 満足度
## 見通し 1.0000000 0.1856953 0.8076889
## 性別 0.1856953 1.0000000 -0.2948839
## 満足度 0.8076889 -0.2948839 1.0000000library(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 |
範囲と偏相関係数をみると，範囲は満足度，偏相関係数は性別の方が大きい．