データ解析：第1回授業　　多変量解析とは何か？&Rの初歩

I.多変量データの分析

1.多変量データとは何か？

  • 多変量データ：複数の変数・変量（variables・variates）からなるデータ。以下の表のように、複数の変数から構成されるデータセットである。
  • 多変量解析：複数の変数・変量からなるデータに対して、データの特性を要約したり、関心のある変数間の関係性を回帰分析等により明らかにすること

2.変数の呼称

  • 原因と結果に対応する呼称

  • 変数データ解析の性質に応じた呼称

3.モデルとは？

  • 多変量解析の中でも、多変量回帰分析（重回帰分析）とは何か？
  • 推定のためのモデル（estimation model）：従属変数＝独立変数 (lhs=rhs)で表される

実際のモデルとはどのようなものか？

  • 検証する仮説群

仮説1：学歴が高いほど、賃金が上がる

仮説2：職業訓練を受けるほど、賃金が上がる

仮説3：専門的な就業年数が長いほど、賃金が上がる

仮説4：仕事に対する動機づけが高いほど、賃金が上がる

これを数式によって表現すると・・・ \[ y_{wage}=\beta_{0}+\beta_{1}x_{education}+\beta_{2}x_{training}+\beta_{3}x_{expertise}+\beta_{4}x_{motivation}+\epsilon. \]

各要素の説明

\(y_{wage}:賃金\)

\(\beta_{0}：定数項（切片）\)

\(x_{education}：学歴\)

\(x_{training}：職業訓練受講歴\)

\(x_{expertise}：専門的な就業年数\)

\(x_{motivation}：仕事に対する動機づけ\)

本授業では、こうした重回帰分析の多様な形態を学ぶことを目指すが、それ以外の重回帰分析の手法も学ぶ。

4．データの種類と多変量解析

データの種類によって、適した分析方法を見極めることがデータ解析には求められる。データの種類と適した分析手法について確認する。

  • 観察データ(observation data)：観察から得られたデータ
    • 横断面データ：特に個人に関する意識調査のデータはサーヴェイ・データ(survey data)と呼ばれる。実験的処置が含まれない1時点からなるサーヴェイデータは、観察データであり、横断面データである。

- パネル・データ：観察データのパネル・データにおいて、個人に関するサーヴェイ・データは、横断面情報\(i\)の方が時系列情報\(t\)より多い、横断面＝時系列データ(cross sectional-time series data: CSTS)データである。これに対して、OECD諸国に関するパネルデータ等の場合、横断面情報\(i\)の方が時系列情報\(t\)より少ないこともあり、その場合は時系列＝横断面データ(time series-cross section data: TSCS data)と呼ばれる。

  • 時系列データ(time-series data):ある特定の集団、地域、国に関する長期間のデータが集積されたもの。経済成長率、日経平均株価、支持率などのデータがある

  • 実験データ(experimental data)：実験によって得られたデータ
    • 横断面データ(cross-section data)：ある1時点における個人、グループ、集団、場所・領域などに関する情報を集めたデータであり、実験データにおける横断面データには、必ず介入(intervention)・処置(treatment)と統制(control)・対照(contrast)を分ける変数が含まれる

  • パネル・データ(panel data)/横断面＝時系列データ(cross section-time-series data: CSTS data)：個人、グループ、場所に関する複数時点の観察が含まれたデータ。実験データの場合、多くの場合に横断面\(i\)に関する情報が、時系列\(t\)に関する情報より多く含まれる。

II.因果関係の問題

1. 相関関係(correlation)があるということ

相関係数\(r\)が0.6の場合：

相関係数\(r\)が0.6の場合：

相関係数\(r\)が0.02の場合(つまり\(r=0\)を否定できない場合)：

ちなみに相関に関する図は、次のようなコードで出力していた。こうしてコードを書けるようになっていこう！

set.seed(123) # 同じ結果を得るためにseedを設定
n <- 100 # サンプルサイズ

# 相関が0.6になるようなデータセット生成
x <- rnorm(n)
y <- -0.02 * x + sqrt(1 - (-0.02)^2) * rnorm(n)
df <- data.frame(x, y)

#必要なパッケージの導入
library(ggplot2)
library(ggalt)

#描画
p <- ggplot(df, aes(x=x, y=y)) +
  geom_point() + # 散布図
  geom_smooth(method='lm', se=FALSE, linetype='dashed', color='blue') + # 回帰直線
  annotate("text", x = mean(df$x), y = mean(df$y), label = "r=0.02", color = "red") + # 相関係数のテキストを中央に配置
  coord_fixed(ratio = 1) + # アスペクト比を1に設定
  geom_encircle()+# プロット全体を囲む
  theme_bw()

print(p)

2. 相関係数の計算：まずはざっとR(RStudio)に触れてみよう！

相関関係について考えるために、まずはRを使ってピアソン積率相関係数を計算してみよう。

df_wel<-read.csv("welfare_state.csv")
attach(df_wel)
head(df_wel)
##    id year     pop  po65 socx_pub    per_aged leftseat
## 1 AUL 1960 10275.0 874.9       NA 0.085148418     36.9
## 2 AUL 1961 10508.2 894.6       NA 0.085133515     37.9
## 3 AUL 1962 10700.5 913.6       NA 0.085379188     49.2
## 4 AUL 1963 10906.9 933.0       NA  0.08554218     48.5
## 5 AUL 1964 11121.6 948.1       NA 0.085248525     41.0
## 6 AUL 1965 11340.9 966.3       NA 0.085204878     41.0

分析結果を解釈する。結果の見方がわかった上で、結果の解釈の仕方を覚えているだろうか？

per_aged<-as.numeric(per_aged)
cor.test(socx_pub,per_aged)
## 
##  Pearson's product-moment correlation
## 
## data:  socx_pub and per_aged
## t = 18.372, df = 649, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5319956 0.6333041
## sample estimates:
##       cor 
## 0.5849266

では、相関関係があるということは、因果関係があることを意味するだろうか？

  • 高齢者人口割合が増える→社会保障費が拡充される

  • 社会保障費が拡充され福祉が充実する→人々の寿命が伸びた

  • あるいはその両方の相乗効果・シナジーなのかも？

  • いや、人々の寿命を伸ばし、社会保障費も上げる第3の要因があるのかも？！（例：平和で戦争のない民主主義体制の継続）

2つの要素間に相関関係があることが、ただちに両者の因果関係を必ずしも意味しない。これはどういう意味なのだろうか・・・

3.因果関係の根本問題：もしも実験ができるなら・・・

  • 因果関係の根本問題：Nyman-Rubinの因果モデル

各個体それぞれにおいて処置（treaetment　他に処理・介入等とも表現）を受けた時と受けなかった時の両方が観察可能な時：その処置による結果（outocome）の差を定量的に観測することで、処置による因果効果を確かめることができる

⇒その根本問題とは？：「処置を受けた個体」と「処置を受けなかった個体」、どちらかは存在せず観察できないので、個体に関する処置の差（unit level causal effect）を比較することはできない

  • 因果効果の測定をめぐる代替的手法

以下のStabe Unit Treatment Value Assumption (SUTVA)のもとで、集団レヴェルでの因果効果(population level causal effect)を評価する；

（1）観察される個体の反応は、他の個体がどのような処理を受けてもその影響を受けない→他の個体とは独立

（2）処理水準はそれぞれ1つの水準に定められる

集団間の結果の差を平均処置効果（average causal effect: ATE）として処置の効果、すなわち因果効果とみなす

  • 実験的手法

    • ランダム化比較試験（randomized controll treatment: RCT）：ランダム・サンプリングにより統制群(control group)と処置群（treatment group）に被験者を無作為割り当て（random assaignment）することによって、処置の集合的な因果効果を見ることによって、因果効果を特定する

    • 自然実験（natural experiment）：あるイヴェントが生じたことにより、「偶然に」処置を受けたグループが生じた際、処置をを受けたグループと受けていない統制群との間で結果の差を観察する（例は口頭で説明します）

⇒こうした実験的手法によって集められたデータを「実験データ」と呼ぶ。因果効果を直接的に測定しようとする設計に基づていることから、より因果効果を適切に測りやすいと考えられている

しかし、考えてみよう。いつでもこうした実験的な手法が採用できるだろうか・・・？

4.相関関係と因果関係：多変量解析の技法

先ほどの相関と因果の問題に戻ろう。

  • 要因間に相関が認められるとしてもその相関関係が、

（1）2つの要因間の因果関係によって生じたものなのか

（2） 2つの要因の共通原因となる第三の要因があって、その要因の変動によって生じたものなのか

を区別する必要がある（黒木学『構造的因果モデルの基礎』共立出版、2017年、3頁）

  • 第3の要因：交絡(covariates)の問題　

    • 【交絡(covariates)】民主主義が人々の寿命にも、民主主義が社会保障支出の増大もともに促し、そのもとで「高齢者の人口割合→社会保障支出」という因果が観察されているのが観察されているのかもしれない

- 【因果媒介効果(causal mediation effect)】高齢者人口の増加が成熟した民主主義を促し、それが民主主義的福祉国家のもとでの社会保障支出をもたらしているかもしれない

関心のある2変数間の因果関係に影響を与える、第3番目以降の変数、すなわち多変量を扱うのが多変量解析であり、何を明らかにするかに応じて利用する分析手法が変わる

5. 多変量解析の技法

多変量解析とひとくちに言っても、その中には多くの技法が含まれる。「I-3．モデル」の説明で取り上げた「多変量回帰」はその一例だが、他にも次のようなデータの要約、因果性の検証ための技法がある。

  • 重回帰分析 (multivariates regrression analysis)：ある1つの目的変数に対する因果性を仮定し、その要因と考えられる複数の変数からの因果性を同定する（identify）ための方法
result<-lm(socx_pub~per_aged+leftseat)
summary(result)
## 
## Call:
## lm(formula = socx_pub ~ per_aged + leftseat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.1099 -2.9080 -0.0915  2.8955 12.9471 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.87238    0.99570   1.880   0.0605 .  
## per_aged    117.18442    7.02275  16.686  < 2e-16 ***
## leftseat      0.06025    0.01050   5.737 1.48e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.074 on 648 degrees of freedom
##   ( 559 個の観測値が欠損のため削除されました )
## Multiple R-squared:  0.3739, Adjusted R-squared:  0.372 
## F-statistic: 193.5 on 2 and 648 DF,  p-value: < 2.2e-16
  • 構造方程式モデリング (Structural Equation Modeling (SEM))：複数の変数間に複数の因果を仮定し、共分散行 列や相関行列をもとに因果性を同定するための方法
library(lavaan)　#構造方程式モデリングのためのパッケージをロード
library(semPlot)

DF2 <- data.frame(per_aged, socx_pub , leftseat)　#分析に利用するデータフレームの作成

model <- '          #モデルの設定
 
  # 回帰
  socx_pub~leftseat
  per_aged~leftseat
  socx_pub~ per_aged
 '
fit2 <- sem(model, data=DF2)　#推定
summary(fit2)
## lavaan 0.6-12 ended normally after 1 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         5
## 
##                                                   Used       Total
##   Number of observations                           651        1210
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   socx_pub ~                                          
##     leftseat          0.060    0.010    5.750    0.000
##   per_aged ~                                          
##     leftseat          0.000    0.000    6.852    0.000
##   socx_pub ~                                          
##     per_aged        117.184    7.007   16.725    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .socx_pub         16.523    0.916   18.042    0.000
##    .per_aged          0.001    0.000   18.042    0.000
semPaths(fit2, "model", "est", intercepts = FALSE)　#推定結果の描画

  • 因子分析 (Factor Analysis)
library(caret) #ダミー変数作成のために利用
library(psych)　#因子分析のための```fa()```コマンドのために利用

politician<-read.csv("politician2017.csv", fileEncoding="SJIS")　#データを読み込み
attach(politician)

tmp<-data.frame(as.factor(Q1_1)) #データフレームを作成

tmp1 <- dummyVars(~.,data=tmp) #ダミー変数を効率良く作成
policy<- as.data.frame(predict(tmp1, tmp))

result <- fa(policy,factors=2,rotation="varimax")　#因子分析
result
## Factor Analysis using method =  minres
## Call: fa(r = policy, factors = 2, rotation = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
##                      MR1      h2     u2 com
## as.factor.Q1_1..1  -1.04 1.07334 -0.073   1
## as.factor.Q1_1..2   0.06 0.00411  0.996   1
## as.factor.Q1_1..3   0.09 0.00839  0.992   1
## as.factor.Q1_1..4   0.05 0.00265  0.997   1
## as.factor.Q1_1..5   0.25 0.06499  0.935   1
## as.factor.Q1_1..6   0.10 0.01025  0.990   1
## as.factor.Q1_1..7   0.07 0.00472  0.995   1
## as.factor.Q1_1..9   0.02 0.00025  1.000   1
## as.factor.Q1_1..10  0.05 0.00265  0.997   1
## as.factor.Q1_1..11  0.04 0.00129  0.999   1
## as.factor.Q1_1..12  0.08 0.00567  0.994   1
## as.factor.Q1_1..13  0.06 0.00381  0.996   1
## as.factor.Q1_1..14  0.02 0.00050  0.999   1
## as.factor.Q1_1..15  0.04 0.00182  0.998   1
## as.factor.Q1_1..16  0.10 0.00949  0.991   1
## as.factor.Q1_1..99  0.07 0.00535  0.995   1
## 
##                 MR1
## SS loadings    1.20
## Proportion Var 0.07
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## The degrees of freedom for the null model are  120  and the objective function was  22.11 with Chi Square of  10122.98
## The degrees of freedom for the model are 104  and the objective function was  21.79 
## 
## The root mean square of the residuals (RMSR) is  0.06 
## The df corrected root mean square of the residuals is  0.06 
## 
## The harmonic number of observations is  465 with the empirical chi square  350.31  with prob <  2e-28 
## The total number of observations was  465  with Likelihood Chi Square =  9962.43  with prob <  0 
## 
## Tucker Lewis Index of factoring reliability =  -0.139
## RMSEA index =  0.451  and the 90 % confidence intervals are  0.445 0.46
## BIC =  9323.66
## Fit based upon off diagonal values = 0.37
  • 主成分分析 (Prinsipal Component Analysis:PCA)
  • クラスター分析 (Clauster Analysis)
data<-iris

distance <- dist(data) #ユークリッド距離
hc <- hclust(distance, "ward.D2")　#　樹形図

plot(hc) #デンドログラムの描画

上記いずれの方法もデータの特性を要約し、関心のある要素間の関係性を分析するものである。分散分析法（Analysis of Variance: ANOVA）が、どちらかというと2変数間（1つの独立変数と1つ従属変数）の平均値の差や分散の差に焦点を当てるものである一方で、複数の変数の要約・因果性の特定を試みていることを理解しておこう。

　

III.Rを使ったデータ分析

1．統計分析用のソフトウェア

データ解析を行うに際して、利用可能なソフトウェアは以下のように多くある。

  • Excel：Microsoft Office自体は有料だが、PCに標準装備のことも多い

  • HAD：Excelがあれば、それにアドインすることで無料で利用可能。こちらよりダウンロード可能。以下のSPSSに対応する分析が、無料で入手できるツールによって可能

  • SPSS：有料。Graphic User Interface (GUI)により簡易な操作で各種の統計分析が可能

  • STATA：有料。優れたGUIを有しつつ、コマンド入力により高度な分析も可能。また高度な分析を行う際のアドインも充実

2．Rを利用するメリット

  • 無料

  • コマンド入力とpackageにより高度な分析が可能

  • Rユーザーの層が厚いことから、かなり高度な分析であっても、solutionをオンライン上で検索して発見することができる

3．実践：RStudioのインストール

①まずはじめにRstudioを各自のパソコンにインストールし利用できるようにしよう。

②次にR本体をインストールしよう。 ※Macの場合は、「Mac R download」で検索し、ダウンロードしてください。

4．実践：Projectを作ろう

①RStudioを起動したら、まずは「Project」を作るとよい。Projectによって、1つの分析プロセスを構成し管理していくことがスムーズな作業を進める鍵になる。以下の赤く囲んだところから「New Project」を選択。

②次にProjectを作るフォルダを指定する。新しくフォルダを作成するならNew Directoryを選択、既存のフォルダを利用するならExisting Directoryを選択する。以下では、New Directoryを選択するものとして説明するが、どちらでも難易度が異なるということはない。

③続いて、New Projectを選択したら・・・

④BrowseボタンからProjectを作成したいフォルダを選択し、Directory NameでProjectの名前を入力する。可能な限りローマ字表記で。

⑤ここでは「practice」という名前のプロジェクトを作成したが、下の画面のように、画面右上に自分の入力した名前が表示されていればprojectの作成は成功。

⑥次に画面の右下、Filesというタブを見てみよう。Projectが保存されているフォルダ内のファイルを確認できるが、そこに拡張子.Rprojのプロジェクトファイルも格納されているのがわかる。

5.実践：データの読み込み

最初に利用するデータを読み込むことから始めよう。

①Rでは、他のソフトウェアで作られたデータセットを読み込むことが多い。最も多いパターンは、Excelによる.csv(comma separated variables)データである。ここで注意するべきことは、データファイルの保存場所である。利用するデータを必ず、作成したProjextファイルと同じフォルダ内に保存しよう。保存が成功しているなら、下の画面のように、Filesタブ内に、データファイルの存在を確認することができる。

②いよいよコマンドを入力していくが、コマンドを入力するのは「R Script」という画面である。ここにコマンドを書き、スクリプトとして保存することで、毎回同じ作業を繰り返して行うことができる。そしてスクリプトを共有することで、他の分析者とも「再現可能性」をもって、作業を追試しあうことができる。スクリプト画面を立ち上げてみよう。

画面左上の「＋」マークから、「R Script」を選択するとスクリプトを記入する画面が立ち上がる。スクリプトに書いたコマンドを実行した結果は、下の「Console」画面に表示される。それによって生成されたオブジェクトは、左上の「Environment」タブに表示されていく。

③では実際に、データを読み込むコマンドを入力してみよう。次のようなコマンドを打つとデータを読み込むことができる。attach()コマンドによって、データファイル内の変数をスムーズに分析に利用することができるので、データの読み込みと同時にこのコマンド処理もしておくようにしよう。

df_cws<-read.csv("cws.csv")
attach(df_cws)

上記のコマンドについて説明する。最初に書いたcws」というのは、データを格納するオブジェクトであり、任意の名前をつけてもらって構わない。data frameを略して、dfで始まる名前を付けることも多い。「<-」はオブジェクトにread.csvのコマンドによって、()内の作業をしたオブジェクトを格納するという記号である。「小なり（<）」と「ハイフン（-）」からなっている。次のattach()はデータフレーム内に含まれる各変数を分析に利用できるようにするためのコマンドである。

スクリプトにはこのように書いていくことになる。最初の赤く囲んだ部分は、このスクリプトタイトルである。ポイントは#マーク。#で始まる行、あるいは文言はメモとして認識され、コマンドとは認識されない。だから7行目の###-------load dataもコマンドとは認識されないメモ行である。そしてすべての行を選択し、左上の青で囲んだ「Run」をクリックするとコマンドが実行される。

ここで実行された結果を「Console」で確認しよう。赤で囲んだのが1回目の試行だが、attach()コマンドをattacをミスタイプしている。するとエラーが赤文字で表示される。多くの場合、赤文字で表示さているときには、その作業がコマンドのエラーによりうまくなされていないことを表している。2回目はの試行では、赤文字は出ておらず、次のコマンドを待機する>が表示されている。

そして右上の「Environment」には、新たに生成されたデータオブジェクトとして「df_cws」が表示されている。

6.実践：パッケージのインストール

①Rでは、様々な形式のデータファイルを読み込むことができるが、そのためにはパッケージ（package）のインストールが必要になる。パッケージのインストールは次のコマンドで実行できる。

install.packages("openxlsx")

ここでは、エクセルファイルを幅広く読み込んだり書きだしたりするための、openxlsxパッケージを例に書いたが、簡単な分析から高度な分析まで多様なパッケージをインストールし分析に利用できるのがRの強みである。

②他にも以下のように、画面右下の「Packages」タブから、「install」を選択。

ウィンドウ内に「openxlsx」と入力してパッケージをインストールすることも可能。

③そしてExcelの「.xlsx」ファイルを読み込む場合には、以下のlibrary(openxlsx)によってまずはパッケージを読み込む。続いて、read.xlsx()によって、次のようにコードを書けば、エクセルファイルの”cws”シートのところにあるデータを読み込むことができる。

library(openxlsx)
df_cws<-read.xlsx("cws.xlsx", sheet = "cws")
attach(df_cws)

7.保存

①スクリプトの保存が最も大事。下の保存のアイコンをクリック。

②任意の名前を付けて「保存」すると、スクリプト画面の名称が変更になり（赤）、右下の「Files」画面に拡張子「.R」のファイルが追加されているのがわかる。

これ以降は、「保存」のアイコンをクリックすることで上書きすることが可能。

IV. R Markdownを使ってみよう

R Markdownの意義

あらゆる分野の分析において、分析結果の追試可能性、再現性はとても重要である。他の研究者が、同じデータや同じ条件下でのデータを使って同じ分析をしたときに、同様の分析結果が得られるのかを確かめることが科学においては重視される。その一連の作業を効率よく進めるために、RにはR上でコマンドを入力し分析をするだけではなく、分析結果の報告といった文章情報も統合的に入力する機能として「R Markdown」機能が備わっている。データ解析Bの授業では、この機能を積極的に利用し、レポートを効率的に、再現性を保ったうえで書いていく技術を身に着けよう。

R Markdownの実践

①再び、画面左上の「＋」マークから、「R Markdown」を選択。

②「Title」欄に任意の名前を入力し、現時点では、「HTML」出力を選択して「OK」。

③デフォルトで、次のような例示画面が出力される。それぞれの行の意味は説明とともに、以下図を参照してもらいたい。

④この段階で早速、以下図の赤で囲んだ「knit」ボタンをクリックしてみよう。markdownを保存することが求められるので、任意の名前を入力し「save」。

⑤以下のようなHTML画面が表示される。

そしてこの説明のためのHTML画面も、実は・・・