Rではデータフレームでデータセットを表現する
ここにいろいろな説明を書く 手順などをかく
mtcars
#データを表示する(データフレーム名を打ち込む)
mtcars
NA
? mtcars
#データの桁数と列数を表示する
dim(mtcars)
[1] 32 11
32行・・・32個の観察がある サンプルサイズが32である
11列・・・属性や特徴をあらわす変数が11個ある
#データの概要を表示する
summary(mtcars)
mpg cyl disp hp drat wt qsec
Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0 Min. :2.760 Min. :1.513 Min. :14.50
1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5 1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89
Median :19.20 Median :6.000 Median :196.3 Median :123.0 Median :3.695 Median :3.325 Median :17.71
Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7 Mean :3.597 Mean :3.217 Mean :17.85
3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0 3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90
Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0 Max. :4.930 Max. :5.424 Max. :22.90
vs am gear carb
Min. :0.0000 Min. :0.0000 Min. :3.000 Min. :1.000
1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
Median :0.0000 Median :0.0000 Median :4.000 Median :2.000
Mean :0.4375 Mean :0.4062 Mean :3.688 Mean :2.812
3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
Max. :1.0000 Max. :1.0000 Max. :5.000 Max. :8.000
燃費mpgの最小値は10.40で、最大値は33.90である 平均meanは20.09で、中央値medianは19.20である
四分位数:データを大きさ順に並べたときに25%、50%、75%に相当する位置の値。 第2四分位は中央値と同じになる。上の例では19.20 第1四分位数は25%に相当する位置の値。上の例では15.43 第3四分位数は75%に相当する位置の値。上の例では22.80
条件に応じて、データを選択抽出する
抽出・選択にはいくつかの形式がある。 基本的には条件式と、抽出対象との組み合わせで、抽出データを表現し、データの抽出を行う。
$は、アクセス演算子。データフレーム内の変数にアクセスします。
amはオートマチックギア・マニュアルギア表す
am = 0 のとき、オートマチックギア am = 1 のとき、マニュアルギア
mtcars$am = 0もしやると、mtcarsの中のam変数がすべて0になります。 = は代入演算子です。
以下にam = 0の時の条件式を示す
mtcars$am == 0
[1] FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE
[21] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
実行すると、FALSEとTRUEが並んだ結果が表示される。 FALSEはam=0ではないデータ、TRUEはam=0であるデータである。
Exercise 1の「データフレームの表示」で表示したデータと、上記のFALSEとTRUEを比較してください。きちんと条件式どおりになっていることが確認できます。
今度は逆にam = 1の条件を示す
mtcars$am == 1
[1] TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE
[21] FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
結果はam==0のときと反転しています。amは0 or 1のデータなので、am == 1は am == 0の反転した結果となります。
否定の条件を示すときは !=を使います。下記の条件式は「am == 0ではない」を示します。結果はam == 1と同じになります。
mtcars$am != 0
[1] TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[18] TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
まずデータフレームの任意のデータを選択抽出します $はアクセス演算子です [,]も、アクセス演算子です 以下のように、[,]演算子を使います。
(1行目・1列目)のデータにアクセスするには、[1,1]を指定します
mtcars[1,1]
[1] 21
(3行目・2列目)のデータにアクセスするには、[3,2]を指定します
mtcars[3,2]
[1] 4
3行目の全てのデータにアクセスするときには、[3,]を指定します
mtcars[3,]
1列目の全てのデータにアクセスするときには、[,1]を指定します
mtcars[,1]
[1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2 10.4 10.4 14.7 32.4 30.4 33.9
[21] 21.5 15.5 15.2 13.3 19.2 27.3 26.0 30.4 15.8 19.7 15.0 21.4
列の指定のときには、列名を使うこともできます
mtcars[,"mpg"]
[1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2 10.4 10.4 14.7 32.4 30.4 33.9
[21] 21.5 15.5 15.2 13.3 19.2 27.3 26.0 30.4 15.8 19.7 15.0 21.4
#オートマチックギアのデータだけ表示する
mtcars[mtcars$am == 0, ]
NA
#マニュアルギアのデータだけ表示する
mtcars[mtcars$am == 1, ]
NA
条件式には不等式も使うことができます
以下では、mpgが20未満のデータを表示しています
mtcars[mtcars$mpg < 20, ]
オートマティック・ギア(am=0)の群と、マニュアル・ギアの群で燃費mpgの平均値を比較します。
オートマティック・ギアの群の燃費mpgの平均値を算出します
#オートマティックギアの群のmpgのデータ
mtcars[mtcars$am==0, "mpg"]
[1] 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2 10.4 10.4 14.7 21.5
[16] 15.5 15.2 13.3 19.2
#オートマティックギアの群のmpgのデータの個数
length(mtcars[mtcars$am==0, "mpg"])
[1] 19
#オートマティックギアの群のmpgの平均値
mean(mtcars[mtcars$am==0, "mpg"])
[1] 17.14737
#オートマティックの群のmpgの標準偏差
sd(mtcars[mtcars$am==0, "mpg"])
[1] 3.833966
マニュアルギアの軍の燃費mpgの平均値を算出します
#マニュアルギアの群のmpgのデータ
mtcars[mtcars$am==1, "mpg"]
[1] 21.0 21.0 22.8 32.4 30.4 33.9 27.3 26.0 30.4 15.8 19.7 15.0 21.4
#マニュアルギアの群のmpgのデータの個数
length(mtcars[mtcars$am==1, "mpg"])
[1] 13
#マニュアルギアの群のmpgの平均値
mean(mtcars[mtcars$am==1, "mpg"])
[1] 24.39231
#マニュアルギアの群のmpgの標準偏差
sd(mtcars[mtcars$am==1, "mpg"])
[1] 6.166504
平均値の差の検定(t検定)をします
t.test(mtcars[mtcars$am==0, "mpg"], mtcars[mtcars$am==1, "mpg"])
Welch Two Sample t-test
data: mtcars[mtcars$am == 0, "mpg"] and mtcars[mtcars$am == 1, "mpg"]
t = -3.7671, df = 18.332, p-value = 0.001374
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-11.280194 -3.209684
sample estimates:
mean of x mean of y
17.14737 24.39231
(参考)Welchの方法のt検定は、以下のようにt値を算出して、t分布からの乖離をみます https://ja.wikipedia.org/wiki/ウェルチのt検定
統計量tの算出は通常のt統計量の算出方法と同じです。 自由度νの算出がやや複雑ですが、手計算でも検算できます
一応、統計量tの算出だけ手計算で書いておくと、以下のようになります
(17.14737 - 24.39231)/ sqrt(3.833966^2/19 + 6.166504^2/13)
[1] -3.767124
上記の手計算した統計量tは、t.test関数の出力の2行目にあるように、t=-3.7671とほぼ同じであることがわかります(違いは丸め誤差)。
Rでは機能を追加するためにパッケージを追加する
#パッケージの追加
install.packages("ggplot2")
URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.6/ggplot2_3.3.3.tgz' を試しています
Content type 'application/x-gzip' length 4072099 bytes (3.9 MB)
==================================================
downloaded 3.9 MB
The downloaded binary packages are in
/var/folders/2b/yr0kxxc57712kvx0j2h0m7h00000gp/T//RtmpMJPSQY/downloaded_packages
インストールしたパッケージは、一度、実行環境にロードする必要がある
#パッケージのロード
library(ggplot2)
#2変数のグラフ(散布図・プロット図)を作る
ggplot(data=mtcars) + geom_point(aes(x=wt, y=mpg))
2変数のプロット図に、第3変数の情報(am)を追加する
第3変数の情報は、カラー(プロット点の色)やシェイプ(プロット点の形状)を用いる
ggplot(data=mtcars) + geom_point(aes(x=wt, y=mpg, colour= am))
amは0 or 1の離散値であるが、連続量として処理されている。
一応、上のグラフでも意味はわかるが、本来は以下のように、as.factor関数を使うとよい
ggplot(data=mtcars) + geom_point(aes(x=wt, y=mpg, colour= as.factor(am)) )
凡例のタイトルが気持ち悪い人は、labs関数を使うとラベル(labels)のタイトルを変更できる
ggplot(data=mtcars) + geom_point(aes(x=wt, y=mpg, colour= as.factor(am)) ) + labs(color="am")
同様に、形状に情報を対応させることで、第4変数の情報も追加することが出来る
ggplot(data=mtcars) + geom_point(aes(x=wt, y=mpg, colour= as.factor(am), shape=as.factor(cyl))) + labs(color="am", shape="cyl")
もし、プロット点が小さい場合(見えにくい)は、sizeを指定してやると、プロット点のサイズを調節できる
ggplot(data=mtcars) + geom_point(aes(x=wt, y=mpg, colour= as.factor(am), shape=as.factor(cyl)), size=3) + labs(color="am", shape="cyl")
2変数間のプロット図を、複数の変数の組み合わせ散布図行列と呼ぶ
1.散布図行列のパッケージの追加
散布図行列を作成するためには、パッケージを追加し、ロードする。
install.packages("GGally")
URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.6/GGally_2.1.0.tgz' を試しています
Content type 'application/x-gzip' length 1636919 bytes (1.6 MB)
==================================================
downloaded 1.6 MB
The downloaded binary packages are in
/var/folders/2b/yr0kxxc57712kvx0j2h0m7h00000gp/T//RtmpMJPSQY/downloaded_packages
library(GGally)
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
散布図行列を作図する
・各変数の分布 ・相関係数 ・プロット図 が表示されていることを確認すること
ggpairs(mtcars)
たくさんの変数で散布図行列を書くと、 ・処理が遅くなる ・変数間の関係をとらえるのが難しいくなる ので、変数を絞って散布図行列を作図する
変数を絞る際に、 ・何が目的変数になるのか ・何が説明変数になるのか を考えながら変数を絞り込んでいくこと
ggpairs(mtcars[, c("mpg","wt","am","cyl")])
回帰分析に先立って、必ず 目的変数と説明変数の関係をプロット図で確認すること。
・変数間の大まかな関連性の把握 ・そもそも線形の関係なのか(2次や3次の関係ではないのか?) ・天井効果や床効果はあるのか ・外れ値の影響はあるのか
などは、プロット図を見ることですぐにわかる
今回の分析では 目的変数:燃費mpg 説明変数:それ以外の変数
library(ggplot2)
#グラフの描画
ggplot(data=mtcars) + geom_point(aes(x=wt, y=mpg))
NA
NA
より慎重にプロット図を作成するなら、散布図行列を作成しても良い
library(GGally)
ggpairs(mtcars[, c("mpg", "wt")])
散布図行列では、まず目的変数の分布に注目すること!
◯目的変数 目的変数の分布が正規分布に近ければ、通常の回帰分析を行えば良い
目的変数が特殊な分布をしてたときには、特殊な分析方法が必要になる
◯説明変数 説明変数に関しては、分布の仮定は特にない。 説明変数間の関係に注目すること。
説明変数の範囲についても注目すること。 もしグループごとに説明変数の範囲が異なる、ということがあれば注意すること。
もしも欠損値がある場合は欠損値のパターンについても考慮すること
lm関数は最小二乗法によって回帰分析をおこなう関数である 最小二乗法による推定をOLS推定と呼んでいる
#回帰式の推定
m1 = lm(data=mtcars, mpg ~ wt)
#結果を表示する
m1
Call:
lm(formula = mpg ~ wt, data = mtcars)
Coefficients:
(Intercept) wt
37.285 -5.344
分析結果について、もうすこし多くの情報を表示する
summary(m1)
Call:
lm(formula = mpg ~ wt, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-4.5432 -2.3647 -0.1252 1.4096 6.8727
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
wt -5.3445 0.5591 -9.559 1.29e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10
観測値、予測値、残差を表示する
cbind(observed=mtcars$mpg, predited=predict(m1), residual=resid(m1))
observed predited residual
Mazda RX4 21.0 23.282611 -2.2826106
Mazda RX4 Wag 21.0 21.919770 -0.9197704
Datsun 710 22.8 24.885952 -2.0859521
Hornet 4 Drive 21.4 20.102650 1.2973499
Hornet Sportabout 18.7 18.900144 -0.2001440
Valiant 18.1 18.793255 -0.6932545
Duster 360 14.3 18.205363 -3.9053627
Merc 240D 24.4 20.236262 4.1637381
Merc 230 22.8 20.450041 2.3499593
Merc 280 19.2 18.900144 0.2998560
Merc 280C 17.8 18.900144 -1.1001440
Merc 450SE 16.4 15.533127 0.8668731
Merc 450SL 17.3 17.350247 -0.0502472
Merc 450SLC 15.2 17.083024 -1.8830236
Cadillac Fleetwood 10.4 9.226650 1.1733496
Lincoln Continental 10.4 8.296712 2.1032876
Chrysler Imperial 14.7 8.718926 5.9810744
Fiat 128 32.4 25.527289 6.8727113
Honda Civic 30.4 28.653805 1.7461954
Toyota Corolla 33.9 27.478021 6.4219792
Toyota Corona 21.5 24.111004 -2.6110037
Dodge Challenger 15.5 18.472586 -2.9725862
AMC Javelin 15.2 18.926866 -3.7268663
Camaro Z28 13.3 16.762355 -3.4623553
Pontiac Firebird 19.2 16.735633 2.4643670
Fiat X1-9 27.3 26.943574 0.3564263
Porsche 914-2 26.0 25.847957 0.1520430
Lotus Europa 30.4 29.198941 1.2010593
Ford Pantera L 15.8 20.343151 -4.5431513
Ferrari Dino 19.7 22.480940 -2.7809399
Maserati Bora 15.0 18.205363 -3.2053627
Volvo 142E 21.4 22.427495 -1.0274952
残差のヒストグラムを見る
残差が正規分布に従うように、回帰直線が推定される
hist(resid(m1))
回帰分析の結果を見やすくするパッケージを追加する
install.packages("stargazer")
URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.6/stargazer_5.2.2.tgz' を試しています
Content type 'application/x-gzip' length 621898 bytes (607 KB)
==================================================
downloaded 607 KB
The downloaded binary packages are in
/var/folders/2b/yr0kxxc57712kvx0j2h0m7h00000gp/T//RtmpMJPSQY/downloaded_packages
library(stargazer)
Please cite as:
Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
#回帰分析の結果を表(回帰テーブル)にして出力する
stargazer(m1, type="text")
length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changed
===============================================
Dependent variable:
---------------------------
mpg
-----------------------------------------------
wt -5.344***
(0.559)
Constant 37.285***
(1.878)
-----------------------------------------------
Observations 32
R2 0.753
Adjusted R2 0.745
Residual Std. Error 3.046 (df = 30)
F Statistic 91.375*** (df = 1; 30)
===============================================
Note: *p<0.1; **p<0.05; ***p<0.01
第3変数の影響を調整するために重回帰分析を行う
第3変数が共通原因となっているとき、交絡変数と呼ぶ
交絡の影響を調整した上で、説明変数が目的変数に影響しているのかを推定する
#回帰式の推定
m1 = lm(data=mtcars, mpg ~ wt)
m2 = lm(data=mtcars, mpg ~ wt+ hp)
m1
Call:
lm(formula = mpg ~ wt, data = mtcars)
Coefficients:
(Intercept) wt
37.285 -5.344
m2
Call:
lm(formula = mpg ~ wt + hp, data = mtcars)
Coefficients:
(Intercept) wt hp
37.22727 -3.87783 -0.03177
分析結果を見やすく表示する。 重回帰分析の分析結果をみやすくまとめたテーブルを回帰テーブルと呼ぶ。
library(stargazer)
stargazer(m1,m2, type="text")
length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changednumber of rows of result is not a multiple of vector length (arg 2)
=================================================================
Dependent variable:
---------------------------------------------
mpg
(1) (2)
-----------------------------------------------------------------
wt -5.344*** -3.878***
(0.559) (0.633)
hp -0.032***
(0.009)
Constant 37.285*** 37.227***
(1.878) (1.599)
-----------------------------------------------------------------
Observations 32 32
R2 0.753 0.827
Adjusted R2 0.745 0.815
Residual Std. Error 3.046 (df = 30) 2.593 (df = 29)
F Statistic 91.375*** (df = 1; 30) 69.211*** (df = 2; 29)
=================================================================
Note: *p<0.1; **p<0.05; ***p<0.01
(参考)階段的回帰(hierarchical regression)
注意:階段的回帰(hierarchical regression)は、階層線形モデル(hierarchical linear model)とは全くの別物です
m1とm2では、m2の方が決定係数が大きくなった。 しかし、hpという変数を1つ加えて、複雑なモデルになっている。
なるべく簡潔なモデルで説明したほうがいいという意味ではm1の方がよい。
m1に対して、m2は良いモデルといえるだろうか?
これに答えるのが階段的回帰である。 (注意:階段的回帰(hierarchical regression)と階層的線形モデル(hierachical linear model)とは全く別物である)
anova(m1,m2)
Analysis of Variance Table
Model 1: mpg ~ wt
Model 2: mpg ~ wt + hp
Res.Df RSS Df Sum of Sq F Pr(>F)
1 30 278.32
2 29 195.05 1 83.274 12.381 0.001451 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
分散テーブルは、統計的有意(0.001451)にm2がm1よりもデータを説明している、と示している
参考までにF値とその確率のの算出方法も付記する
F値の計算
{(278.32-195.05)/1} / (195.05/29)
[1] 12.38057
FよりもF値が大きい確率(Pr(>F))
1- pf( 12.38057, 1, 29)
[1] 0.001451643
注意: データにフィットしたモデルが必ずしも良いかどうかはわからない。
データにフィットしていなくても、理論的に正しいモデルであれば、「なぜデータにフィットしなかったのか?」を考える事。
データにフィットしたから正しいとかんがえると、真のモデルからまったく遠いモデルになってしまう。
これは、対象としているデータが観察データであることが大きな原因である
偏回帰係数が有意かどうかを見ながら、理論をつくるのは危険である(慎重に行わなくてはいけない)
# . は目的変数外のすべての変数を説明変数に投入する、という意味
m10 = lm(data=mtcars, mpg ~ .)
summary(m10)
Call:
lm(formula = mpg ~ ., data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-3.4506 -1.6044 -0.1196 1.2193 4.6271
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.30337 18.71788 0.657 0.5181
cyl -0.11144 1.04502 -0.107 0.9161
disp 0.01334 0.01786 0.747 0.4635
hp -0.02148 0.02177 -0.987 0.3350
drat 0.78711 1.63537 0.481 0.6353
wt -3.71530 1.89441 -1.961 0.0633 .
qsec 0.82104 0.73084 1.123 0.2739
vs 0.31776 2.10451 0.151 0.8814
am 2.52023 2.05665 1.225 0.2340
gear 0.65541 1.49326 0.439 0.6652
carb -0.19942 0.82875 -0.241 0.8122
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.65 on 21 degrees of freedom
Multiple R-squared: 0.869, Adjusted R-squared: 0.8066
F-statistic: 13.93 on 10 and 21 DF, p-value: 3.793e-07
僅かにwtが10%水準で有意である。
もしもmtcarsとまったく同じデータを1000回くりかえし取得したら、 すべての変数が統計的に優位になる
#データサイズを100倍にする
mtcars1000 = mtcars[rep(seq_len(nrow(mtcars)), 1000), ]
#データサイズを表示する
dim(mtcars1000)
[1] 32000 11
m10.1000 = lm(data=mtcars1000, mpg ~ .)
stargazer(m10, m10.1000, type="text")
length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changed
==========================================================================
Dependent variable:
------------------------------------------------------
mpg
(1) (2)
--------------------------------------------------------------------------
cyl -0.111 -0.111***
(1.045) (0.027)
disp 0.013 0.013***
(0.018) (0.0005)
hp -0.021 -0.021***
(0.022) (0.001)
drat 0.787 0.787***
(1.635) (0.042)
wt -3.715* -3.715***
(1.894) (0.049)
qsec 0.821 0.821***
(0.731) (0.019)
vs 0.318 0.318***
(2.105) (0.054)
am 2.520 2.520***
(2.057) (0.053)
gear 0.655 0.655***
(1.493) (0.038)
carb -0.199 -0.199***
(0.829) (0.021)
Constant 12.303 12.303***
(18.718) (0.480)
--------------------------------------------------------------------------
Observations 32 32,000
R2 0.869 0.869
Adjusted R2 0.807 0.869
Residual Std. Error 2.650 (df = 21) 2.147 (df = 31989)
F Statistic 13.932*** (df = 10; 21) 21,223.120*** (df = 10; 31989)
==========================================================================
Note: *p<0.1; **p<0.05; ***p<0.01
上記の回帰テーブルで、2つのモデルはまったく同じ係数を推定している。 しかし、(2)はすべての係数が有意である。 その理由は、サンプルサイズが非常に大きいため、標準誤差が小さくなっているためである
統計的な有意検定は、サンプルサイズに大きく依存する。 だから、統計的な有意性だけをもとに因果性の有無を判断すると、大きく間違うことになる。 (間違う理由は他にももっとたくさんある)
戦略分析を目的とした場合、交互作用モデルと媒介作用モデルは、とても重要なモデルである
交互作用モデルは、第三変数(moderator, 調整変数Z)の水準によって、変数Xの効果が変わることを意味している
変数Xの効果を、増大するためには、変数Zの水準をどのように設定すればよいか、は戦略上、とても有用である。
媒介作用モデルは、変数Xの効果が第3変数Z(媒介変数)を介して、Yに影響を及ぼすパスがあることを意味している (変数Xには間接効果と直接効果が存在することを意味している)
媒介モデルは、間接効果・直接効果を考慮することで、変数Xの総合効果を推定することが可能となる。
ここでは、交互作用モデルを推定してみます
m3 = lm(data=mtcars, mpg ~ wt + hp + wt:hp)
stargazer(m1,m2,m3, type="text")
length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changednumber of rows of result is not a multiple of vector length (arg 2)number of rows of result is not a multiple of vector length (arg 2)
========================================================================================
Dependent variable:
--------------------------------------------------------------------
mpg
(1) (2) (3)
----------------------------------------------------------------------------------------
wt -5.344*** -3.878*** -8.217***
(0.559) (0.633) (1.270)
hp -0.032*** -0.120***
(0.009) (0.025)
wt:hp 0.028***
(0.007)
Constant 37.285*** 37.227*** 49.808***
(1.878) (1.599) (3.605)
----------------------------------------------------------------------------------------
Observations 32 32 32
R2 0.753 0.827 0.885
Adjusted R2 0.745 0.815 0.872
Residual Std. Error 3.046 (df = 30) 2.593 (df = 29) 2.153 (df = 28)
F Statistic 91.375*** (df = 1; 30) 69.211*** (df = 2; 29) 71.660*** (df = 3; 28)
========================================================================================
Note: *p<0.1; **p<0.05; ***p<0.01
交互作用モデルのマージナル効果図を作図するパッケージを追加する
install.packages("interplot")
URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.6/interplot_0.2.2.tgz' を試しています
Content type 'application/x-gzip' length 223136 bytes (217 KB)
==================================================
downloaded 217 KB
The downloaded binary packages are in
/var/folders/2b/yr0kxxc57712kvx0j2h0m7h00000gp/T//RtmpMJPSQY/downloaded_packages
library(interplot)
要求されたパッケージ abind をロード中です
要求されたパッケージ arm をロード中です
要求されたパッケージ MASS をロード中です
要求されたパッケージ Matrix をロード中です
要求されたパッケージ lme4 をロード中です
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
method from
print.data.table
arm (Version 1.11-2, built: 2020-7-27)
Working directory is /Users/tatsu/Dropbox/2020 技術経営論/第3回 経営学でつかう回帰分析/2020-01-19_RSTUDIOでの実行例
interplot(m=m3, var1="wt", var2="hp", hist=TRUE) +
labs(x="hp", y="ME of wt") +
geom_hline(yintercept=0, colour="red")+
geom_vline(xintercept=237, colour="red")
| ### Exercise 9 媒介モデル |
| 媒介モデルの推定方法を説明します。 |
| 媒介モデルは、構造方程式モデル(SEM)の1つであると考えることができます。 |
| 構造方程式モデルは、複数の回帰モデル(構造方程式)から成り立っている統計モデルです。 |
| 媒介モデルの推定は、2つの方法が頻繁に使われます。 |
| 1つめは、複数の重回帰モデルの推定を行い、そこからSEMを推論する方法です。 この方法は、従来から利用されてきた重回帰モデルと連続性があるので、理解しやすい面がありますが、すべてのパラメータを推定しているわけではない点に注意が必要です。 |
| 2つめは、SEMとして推定する方法です。 SEMとして推定すれば、当然、SEMを構成するすべてのパラメータを推定することになります。 解釈が明確になる、という利点もあります。 |
| RでSEMを推定する際には、lavaanを使います |
| lavaanでは、以下のようにモデル式(model_00)を記述します。 lavaanのモデル式については、 http://lavaan.ugent.be/tutorial/ を参照してください。 |
| 概要は以下のページにまとまっています http://lavaan.ugent.be/tutorial/syntax1.html |
| モデル式 model_00 = " y ~ b1x + b2z z ~ b3x ind := b2b3 total := b1 + (b2*b3) " |
| その後、モデル式をsem関数で推定します。 fit_sem_00 = sem(model_00, data=test_data_00) |
| 推定した結果をsummary関数で表示します summary(fit_sem_00, fit.measure=TRUE) |
| ### Exercise 9 媒介モデル(数値例1) |
| もし、lavaanをインストールしていなければ、 install.packages(“lavaan”) を実行します |
r install.packages("lavaan") |
| ### データの作成 |
| 媒介モデルのデータを作成する |
| ```r library(lavaan) |
| ``` |
This is lavaan 0.6-7 lavaan is BETA software! Please report any bugs. |
| ###作成したデータを図示してみる |
| ### 媒介モデルのデータを重回帰モデルで推定するとどうなるか |
| ```r stargazer(fit_00_01,fit_00_02, type=“text”) |
| ``` |
length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changednumber of rows of result is not a multiple of vector length (arg 2) |
| ``` |
| ===================================================================== Dependent variable: ————————————————- y (1) (2) |
x 0.581*** 0.301***
(0.033) (0.036)
z 0.437***
(0.032)
Constant 0.029 0.025
(0.034) (0.031)
| Observations 1,000 1,000 R2 0.236 0.358 Adjusted R2 0.235 0.357 Residual Std. Error 1.068 (df = 998) 0.979 (df = 997) F Statistic 308.370*** (df = 1; 998) 278.312*** (df = 2; 997) ===================================================================== Note: p<0.1; p<0.05; p<0.01 ``` |
| 本来、媒介モデルで生成されたデータを、重回帰モデルで推定すると、model(2)のように推定される |
| データ生成に使った構造方程式は以下のものである y = 0.3x + 0.5z …eq(1) z = 0.6x …eq(2) |
| model (2)では、eq1だけが推定されている。推定された値は y = 0.3x + 0.5z …eq(1) に対して y = 0.301x + 0.437z というように推定されている。 |
| 媒介モデルにおけるxの効果は、 (1) xからzを経由してyに影響するもの(間接効果) (2) xからyに直接的に影響するもの(直接効果) に分かれる。 |
| xのyに対する効果は、(1)+(2)を合計した総合効果である |
| 重回帰モデル(2)で推定されているxの効果は、x→yの間接効果だけであり、 x→yの効果(総合効果)は過小評価されていることがわかる。 |
| データ生成で設定したパラメータは以下の通りである。 x→yの間接効果は、0.60.5=0.3である x→yの総合効果は、0.3 + 0.60.5 = 0.6である。 |
| SEMでは、ind=0.280, total=0.581のように正しい値に近い値が推定されている。 |
| x→yの効果(総合効果)を推定するために、重回帰モデル(2)を用いると、 直接効果のみを推定してしまい(0.3)、総合効果(0.6)と大きな違いが出てしまう。 |
| x→yの総合効果を推定するためには、重回帰モデル(2)を使うのは間違いであり、 重回帰モデル(1)を使うのが正しい。重回帰モデルでは、x→yの総合効果が推定されている(0.581)。 |
| 本来x→yの効果(総合効果)を推定するためには、重回帰モデル(1)を用いるべきであるが、 その代わりに重回帰モデル(2)を用いると誤った値を推定してしまう。 |
| このように、正しいモデルに対して、間違った統計モデルを推定することを 「モデルの特定化の誤り(model misspecification)」と呼ぶ。 |
| 特定化の誤りが生じていると、x→yの効果の大きさを正しく推定できない |
| 媒介モデルでは、x→yの効果(総合効果)は、model(1)のような統計モデルを構築しないと、正しく推定できない |
| model (1)で推定された値は y = 0.581x である。 |
| データ生成に使ったパラメータをもとにすると、 x→yの総合効果は、0.3 + 0.6*0.5 = 0.6 である |
| model(1)では、x→yの総合効果に近い値が推定されている |
サンプルサイズが小さい場合に、重回帰モデル(1)と重回帰モデル(2)を比較する。
重回帰モデル(1)では、変数xの回帰係数は、統計的有意である しかし、重回帰モデル(2)では、変数xの回帰係数は、統計的有ではない
model_mod_00 = "
y ~ 0.3*x + 0.5*z
z ~ 0.6*x
"
test_data_00b = simulateData(model_mod_00,sample.nobs = 50, seed=123)
stargazer(fit_00_01b,fit_00_02b, type="text")
length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changednumber of rows of result is not a multiple of vector length (arg 2)
=================================================================
Dependent variable:
---------------------------------------------
y
(1) (2)
-----------------------------------------------------------------
x 0.525*** 0.229
(0.161) (0.150)
z 0.608***
(0.133)
Constant 0.134 0.038
(0.156) (0.133)
-----------------------------------------------------------------
Observations 50 50
R2 0.181 0.432
Adjusted R2 0.164 0.408
Residual Std. Error 1.075 (df = 48) 0.904 (df = 47)
F Statistic 10.618*** (df = 1; 48) 17.894*** (df = 2; 47)
=================================================================
Note: *p<0.1; **p<0.05; ***p<0.01
model_00 = "
y ~ b1*x + b2*z
z ~ b3*x
ind := b2*b3
total := b1 + (b2*b3)
"
fit_sem_00 = sem(model_00, data=test_data_00)
summary(fit_sem_00,fit.measure=TRUE)
lavaan 0.6-7 ended normally after 11 iterations
Estimator ML
Optimization method NLMINB
Number of free parameters 5
Number of observations 1000
Model Test User Model:
Test statistic 0.000
Degrees of freedom 0
Model Test Baseline Model:
Test statistic 813.275
Degrees of freedom 3
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 1.000
Tucker-Lewis Index (TLI) 1.000
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -2790.937
Loglikelihood unrestricted model (H1) -2790.937
Akaike (AIC) 5591.874
Bayesian (BIC) 5616.413
Sample-size adjusted Bayesian (BIC) 5600.532
Root Mean Square Error of Approximation:
RMSEA 0.000
90 Percent confidence interval - lower 0.000
90 Percent confidence interval - upper 0.000
P-value RMSEA <= 0.05 NA
Standardized Root Mean Square Residual:
SRMR 0.000
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Regressions:
Estimate Std.Err z-value P(>|z|)
y ~
x (b1) 0.301 0.036 8.270 0.000
z (b2) 0.437 0.032 13.801 0.000
z ~
x (b3) 0.639 0.030 21.149 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.y 0.956 0.043 22.361 0.000
.z 0.953 0.043 22.361 0.000
Defined Parameters:
Estimate Std.Err z-value P(>|z|)
ind 0.280 0.024 11.558 0.000
total 0.581 0.033 17.578 0.000
設定したb1,b2,b3はそれぞれ以下のように推定された
推定値b1=0.301 設定したパラメータの値はb1=0.3
推定値b2=0.437 設定したパラメータの値はb2=0.5
推定値b3=0.639 設定したパラメータの値はb3=0.6
それぞれ、設定したパラメータに近い値が推定されている。 (誤差変数の影響があるので、パラメータはずれて推定される)
間接効果(ind=b2*b3)は0.280と推定されている
総合効果(total=b1+b2*b3)は0.581と推定されている
媒介変数モデルがより重要なのは、b2(もしくはb3)が、マイナスのとき(直接効果b1と異符号のとき)である。このとき、意図せざる行為(間接効果が意図と逆の方向に作用するモデル)を表すモデルとなる
具体的には y = b1x + b2z z = b3*x のような媒介モデルで、b1,b3がプラス、b2がマイナスである媒介モデルを考える。
このとき、xは直接的にはyを増加させる(直接効果はプラス)であるが、 xの総合効果はゼロもしくはマイナスである可能性がある。
以下に数値例を掲げる
model_mod_01 = "
y ~ 0.3*x + (-0.5)*z
z ~ 0.6*x
"
test_data_01 = simulateData(model_mod_01,sample.nobs = 1000, seed=123)
ggpairs(test_data_01)
fit_01_01 = lm(data=test_data_01, y ~ x )
fit_01_02 = lm(data=test_data_01, y ~ x + z)
stargazer(fit_01_01, fit_01_02, type="text")
length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changednumber of rows of result is not a multiple of vector length (arg 2)
================================================================
Dependent variable:
--------------------------------------------
y
(1) (2)
----------------------------------------------------------------
x 0.013 0.337***
(0.038) (0.037)
z -0.582***
(0.032)
Constant 0.052 0.048
(0.037) (0.032)
----------------------------------------------------------------
Observations 1,000 1,000
R2 0.0001 0.248
Adjusted R2 -0.001 0.246
Residual Std. Error 1.162 (df = 998) 1.008 (df = 997)
F Statistic 0.112 (df = 1; 998) 164.324*** (df = 2; 997)
================================================================
Note: *p<0.1; **p<0.05; ***p<0.01
データ生成に用いたパラメータを元にすると、 x→yの直接効果は0.3 x→yの総合効果は0.3 + 0.6*(-0.5) = 0.0 である
重回帰モデル(2)をではx→yの効果は0.337と推定されている。 この値を、x→yの効果(総合効果)だと考えると、間違いを起こす。
重回帰モデル(2)で推定された0.337は、x→yの直接効果(x→yの総合効果の一部分) である。
x→yの効果(総合効果)は、重回帰モデル(1)で正しく推定されている。
さらに留意する点がある。 重回帰モデル(1)の決定係数R2が0.0001であり、著しく低い、という点である。
決定係数だけを頼りにして、モデルを選択すると、重回帰モデル(2)を選択してしまう(R2=0.248)。
しかし、いままで見てきたように、x→yの総合効果を表すモデルとして、重回帰モデル(2)は間違っている。 重回帰モデル(2)は、x→yの直接効果(0.3)を推定しているだけで、x→yの総合効果(0.0)を示していない。
つまり、決定係数の大小では、正しいモデルを選択できない。
さらに、誤ったモデルである重回帰モデル(2)の回帰係数が統計的有意である点も注意してほしい。
誤ったモデルであっても、回帰係数は統計的有意になる。
つまり、回帰係数が「有意である」「有意でない」といういことと、x→yの効果が正しく推定されているかどうか、とは関係がない。
「モデル特定化の誤り」は、理論構築上の問題である。 どのようなモデル式を作るか?の段階で起こる問題で、「モデル特定化の誤り」は起こる。
モデル特定化の誤りが起こっている時でも、統計的有意になることはある。
だから、「統計的有意になった」→「モデル(因果関係を考察する)」というのは、大変危険であり、(統計的な手法によほど熟達していない限り)、やめたほうがいい。
model_01 = "
y ~ b1*x + b2*z
z ~ b3*x
ind := b2*b3
total := b1 + (b2*b3)
"
fit_sem_01 = sem(model_01, data=test_data_01)
summary(fit_sem_01,fit.measure=TRUE)
lavaan 0.6-7 ended normally after 11 iterations
Estimator ML
Optimization method NLMINB
Number of free parameters 5
Number of observations 1000
Model Test User Model:
Test statistic 0.000
Degrees of freedom 0
Model Test Baseline Model:
Test statistic 545.526
Degrees of freedom 3
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 1.000
Tucker-Lewis Index (TLI) 1.000
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -2838.466
Loglikelihood unrestricted model (H1) -2838.466
Akaike (AIC) 5686.932
Bayesian (BIC) 5711.471
Sample-size adjusted Bayesian (BIC) 5695.590
Root Mean Square Error of Approximation:
RMSEA 0.000
90 Percent confidence interval - lower 0.000
90 Percent confidence interval - upper 0.000
P-value RMSEA <= 0.05 NA
Standardized Root Mean Square Residual:
SRMR 0.000
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Regressions:
Estimate Std.Err z-value P(>|z|)
y ~
x (b1) 0.337 0.037 9.034 0.000
z (b2) -0.582 0.032 -18.152 0.000
z ~
x (b3) 0.557 0.032 17.255 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.y 1.014 0.045 22.361 0.000
.z 0.988 0.044 22.361 0.000
Defined Parameters:
Estimate Std.Err z-value P(>|z|)
ind -0.324 0.026 -12.506 0.000
total 0.013 0.038 0.335 0.738
参考
重回帰モデルでは、R2(適合度)やp値(統計的有意性)から、正しいモデルを選択することができなかった。
しかし、SEMならば適合度指標を使って、正しいモデルを選択することができるだろうか?
SEMの適合度指標には、CFIやRMSEAがある。
結論から言うと、適合度指標を使ったとしても、SEMでも正しいモデルを選択することはできない
以下は、データ生成とは異なる分析モデルで、SEMを適用した事例。 生成モデルとは異なるモデルであるにもかかわらず、適合度指標は目安をクリアしている。
model_01b = "
y ~ b1*x + b2*z
"
fit_sem_00c = sem(model_01b, data=test_data_01)
summary(fit_sem_00c,fit.measure=TRUE)
lavaan 0.6-7 ended normally after 11 iterations
Estimator ML
Optimization method NLMINB
Number of free parameters 3
Number of observations 1000
Model Test User Model:
Test statistic 0.000
Degrees of freedom 0
Model Test Baseline Model:
Test statistic 284.906
Degrees of freedom 2
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 1.000
Tucker-Lewis Index (TLI) 1.000
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -1425.764
Loglikelihood unrestricted model (H1) -1425.764
Akaike (AIC) 2857.528
Bayesian (BIC) 2872.251
Sample-size adjusted Bayesian (BIC) 2862.723
Root Mean Square Error of Approximation:
RMSEA 0.000
90 Percent confidence interval - lower 0.000
90 Percent confidence interval - upper 0.000
P-value RMSEA <= 0.05 NA
Standardized Root Mean Square Residual:
SRMR 0.000
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Regressions:
Estimate Std.Err z-value P(>|z|)
y ~
x (b1) 0.337 0.037 9.034 0.000
z (b2) -0.582 0.032 -18.152 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.y 1.014 0.045 22.361 0.000