R Notebook

---

Exercise 1 データの簡単な操作

1. データの表示

Rではデータフレームでデータセットを表現する

ここにいろいろな説明を書く 手順などをかく

mtcars

#データを表示する(データフレーム名を打ち込む)
mtcars

NA

2. データのヘルプ(説明)の表示

? mtcars

3. データの次元の表示

#データの桁数と列数を表示する
dim(mtcars)

[1] 32 11

32行・・・32個の観察がある サンプルサイズが32である

11列・・・属性や特徴をあらわす変数が11個ある

4. データの概要の表示

#データの概要を表示する
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

---

Exercise 2 データの抽出

条件に応じて、データを選択抽出する

抽出・選択にはいくつかの形式がある。 基本的には条件式と、抽出対象との組み合わせで、抽出データを表現し、データの抽出を行う。

$は、アクセス演算子。データフレーム内の変数にアクセスします。

1. 条件式の表現

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

2. 任意のデータを選択抽出する

まずデータフレームの任意のデータを選択抽出します $はアクセス演算子です [,]も、アクセス演算子です 以下のように、[,]演算子を使います。

（1行目・１列目）のデータにアクセスするには、[1,1]を指定します

mtcars[1,1]

[1] 21

（3行目・2列目）のデータにアクセスするには、[3,2]を指定します

mtcars[3,2]

[1] 4

３行目の全てのデータにアクセスするときには、[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

3. 条件式と抽出対象を組み合わせて選択抽出をする

#オートマチックギアのデータだけ表示する
mtcars[mtcars$am == 0, ]

NA

#マニュアルギアのデータだけ表示する
mtcars[mtcars$am == 1, ]

NA

条件式には不等式も使うことができます

以下では、mpgが20未満のデータを表示しています

mtcars[mtcars$mpg < 20, ]

---

Exercise 3 平均値の差の検定

オートマティック・ギア(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関数の出力の２行目にあるように、t=-3.7671とほぼ同じであることがわかります（違いは丸め誤差）。

---

Exercise 4 グラフのプロット

1. パッケージのインストール

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

2. パッケージのロード

インストールしたパッケージは、一度、実行環境にロードする必要がある

#パッケージのロード
library(ggplot2)

3. プロット図を作成する

#2変数のグラフ(散布図・プロット図)を作る
ggplot(data=mtcars) + geom_point(aes(x=wt, y=mpg))

4. 第３変数の情報を追加する

2変数のプロット図に、第３変数の情報(am)を追加する

第３変数の情報は、カラー（プロット点の色）やシェイプ（プロット点の形状）を用いる

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")

同様に、形状に情報を対応させることで、第４変数の情報も追加することが出来る

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")

---

Exercise 5. 散布図行列を作成する

２変数間のプロット図を、複数の変数の組み合わせ散布図行列と呼ぶ

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

2. 散布図行列のパッケージをロードする

library(GGally)

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

3. 散布図行列を作成する

散布図行列を作図する

・各変数の分布 ・相関係数 ・プロット図 が表示されていることを確認すること

ggpairs(mtcars)

4. 適切な散布図行列を作成する

たくさんの変数で散布図行列を書くと、 ・処理が遅くなる ・変数間の関係をとらえるのが難しいくなる ので、変数を絞って散布図行列を作図する

変数を絞る際に、 ・何が目的変数になるのか ・何が説明変数になるのか を考えながら変数を絞り込んでいくこと

ggpairs(mtcars[, c("mpg","wt","am","cyl")])

---

Exercise 6 回帰分析をする（単回帰）

1. プロット図を作成する

回帰分析に先立って、必ず 目的変数と説明変数の関係をプロット図で確認すること。

・変数間の大まかな関連性の把握 ・そもそも線形の関係なのか（２次や３次の関係ではないのか？） ・天井効果や床効果はあるのか ・外れ値の影響はあるのか

などは、プロット図を見ることですぐにわかる

今回の分析では 目的変数：燃費mpg 説明変数：それ以外の変数

library(ggplot2)

#グラフの描画
ggplot(data=mtcars) + geom_point(aes(x=wt, y=mpg))

NA
NA

より慎重にプロット図を作成するなら、散布図行列を作成しても良い

library(GGally)
ggpairs(mtcars[, c("mpg", "wt")])

散布図行列では、まず目的変数の分布に注目すること!

◯目的変数 目的変数の分布が正規分布に近ければ、通常の回帰分析を行えば良い

目的変数が特殊な分布をしてたときには、特殊な分析方法が必要になる

◯説明変数 説明変数に関しては、分布の仮定は特にない。 説明変数間の関係に注目すること。

説明変数の範囲についても注目すること。 もしグループごとに説明変数の範囲が異なる、ということがあれば注意すること。

もしも欠損値がある場合は欠損値のパターンについても考慮すること

2. 回帰分析(単回帰)を行う

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))

3. 回帰分析結果の見方

回帰分析の結果を見やすくするパッケージを追加する

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

---

Exercise 7 重回帰分析を行う

1. なぜ重回帰分析をするか?

第３変数の影響を調整するために重回帰分析を行う

第３変数が共通原因となっているとき、交絡変数と呼ぶ

交絡の影響を調整した上で、説明変数が目的変数に影響しているのかを推定する

#回帰式の推定
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  

2. 複数の回帰モデルの比較

分析結果を見やすく表示する。 重回帰分析の分析結果をみやすくまとめたテーブルを回帰テーブルと呼ぶ。

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

3. (参考)階段的回帰

(参考)階段的回帰(hierarchical regression)

注意：階段的回帰(hierarchical regression)は、階層線形モデル(hierarchical linear model)とは全くの別物です

m1とm2では、m2の方が決定係数が大きくなった。 しかし、hpという変数を１つ加えて、複雑なモデルになっている。

なるべく簡潔なモデルで説明したほうがいいという意味では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

注意： データにフィットしたモデルが必ずしも良いかどうかはわからない。

データにフィットしていなくても、理論的に正しいモデルであれば、「なぜデータにフィットしなかったのか？」を考える事。

データにフィットしたから正しいとかんがえると、真のモデルからまったく遠いモデルになってしまう。

これは、対象としているデータが観察データであることが大きな原因である

4. 星占いはやってはだめ

偏回帰係数が有意かどうかを見ながら、理論をつくるのは危険である（慎重に行わなくてはいけない）

# . は目的変数外のすべての変数を説明変数に投入する、という意味
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)はすべての係数が有意である。 その理由は、サンプルサイズが非常に大きいため、標準誤差が小さくなっているためである

統計的な有意検定は、サンプルサイズに大きく依存する。 だから、統計的な有意性だけをもとに因果性の有無を判断すると、大きく間違うことになる。 （間違う理由は他にももっとたくさんある）

---

Exercise 8 交互作用モデル

戦略分析を目的とした場合、交互作用モデルと媒介作用モデルは、とても重要なモデルである

交互作用モデルは、第三変数（moderator, 調整変数Z）の水準によって、変数Xの効果が変わることを意味している

変数Xの効果を、増大するためには、変数Zの水準をどのように設定すればよいか、は戦略上、とても有用である。

媒介作用モデルは、変数Xの効果が第３変数Z(媒介変数)を介して、Yに影響を及ぼすパスがあることを意味している （変数Xには間接効果と直接効果が存在することを意味している）

媒介モデルは、間接効果・直接効果を考慮することで、変数Xの総合効果を推定することが可能となる。

ここでは、交互作用モデルを推定してみます

1. 交互作用モデルを推定する

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

2. マージナル効果図を作図する

交互作用モデルのマージナル効果図を作図するパッケージを追加する

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

2. マージナル効果図を作成する

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 技術経営論/第３回 経営学でつかう回帰分析/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)の１つであると考えることができます。

構造方程式モデルは、複数の回帰モデル（構造方程式）から成り立っている統計モデルです。

媒介モデルの推定は、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

SEMとして媒介モデルを推定する

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と推定されている

---

数値例(2)：意図せざる行為を表す媒介モデル

媒介変数モデルがより重要なのは、b2（もしくはb3）が、マイナスのとき(直接効果b1と異符号のとき)である。このとき、意図せざる行為(間接効果が意図と逆の方向に作用するモデル)を表すモデルとなる

具体的には y = b1x + b2z z = b3*x のような媒介モデルで、b1,b3がプラス、b2がマイナスである媒介モデルを考える。

このとき、xは直接的にはyを増加させる（直接効果はプラス)であるが、 xの総合効果はゼロもしくはマイナスである可能性がある。

以下に数値例を掲げる

SEMで推定する

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の効果が正しく推定されているかどうか、とは関係がない。

「モデル特定化の誤り」は、理論構築上の問題である。 どのようなモデル式を作るか？の段階で起こる問題で、「モデル特定化の誤り」は起こる。

モデル特定化の誤りが起こっている時でも、統計的有意になることはある。

だから、「統計的有意になった」→「モデル（因果関係を考察する）」というのは、大変危険であり、（統計的な手法によほど熟達していない限り）、やめたほうがいい。

SEMで推定する

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