データ解析：第2回授業　Rを用いた記述統計・相関分析

I.記述統計

今日使うデータと目的：日本の内閣支持率に関するデータを使って、内閣支持率はどのような特徴を持ち、何と関係していると考えられるのか考えてみよう。

1.研究の設計

  • 研究の問い：内閣支持率や与党支持率はどのような特徴を持ち、どのように変化してきたのか。経済変数との関係はどのようなものか。
  • 利用するデータ：時事通信社による毎月の内閣支持率に関する時系列形式の世論調査データ（1960年6月～現在）
  • 利用する変数：内閣支持率、自民党支持率、景気評価、消費者物価指数、失業率

2. 記述統計

データの読み込み

今日使うデータを読み込む。

cabinet<-read.csv("cabinet.csv")
attach(cabinet)

シンプルな記述統計

もっともシンプルな記述統計は以下のコマンドによって出力できる。

summary(cabinet)
##    yearMon             time_id         approve           dis       
##  Length:679         Min.   :  1.0   Min.   : 4.40   Min.   : 6.00  
##  Class :character   1st Qu.:170.5   1st Qu.:28.75   1st Qu.:28.75  
##  Mode  :character   Median :340.0   Median :37.60   Median :34.40  
##                     Mean   :340.0   Mean   :36.44   Mean   :36.56  
##                     3rd Qu.:509.5   3rd Qu.:43.30   3rd Qu.:42.55  
##                     Max.   :679.0   Max.   :78.40   Max.   :83.70  
##                                                                    
##    ruling_ldp           q3c1              q3c2             q3c3      
##  Min.   :-1.0000   Min.   : 0.0000   Min.   : 0.000   Min.   : 0.00  
##  1st Qu.: 1.0000   1st Qu.: 0.1000   1st Qu.: 3.200   1st Qu.:48.90  
##  Median : 1.0000   Median : 0.3000   Median : 4.000   Median :58.50  
##  Mean   : 0.8527   Mean   : 0.6076   Mean   : 4.597   Mean   :56.18  
##  3rd Qu.: 1.0000   3rd Qu.: 0.4000   3rd Qu.: 5.300   3rd Qu.:64.50  
##  Max.   : 1.0000   Max.   :11.5000   Max.   :14.800   Max.   :76.70  
##                    NA's   :44        NA's   :44       NA's   :44     
##       q3c4            q3c5             q3c6            q5c1       
##  Min.   : 0.00   Min.   : 0.000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:25.90   1st Qu.: 3.400   1st Qu.:0.400   1st Qu.:0.1000  
##  Median :29.90   Median : 4.900   Median :0.800   Median :0.2000  
##  Mean   :31.11   Mean   : 5.533   Mean   :1.196   Mean   :0.3546  
##  3rd Qu.:36.55   3rd Qu.: 6.800   3rd Qu.:1.700   3rd Qu.:0.4000  
##  Max.   :53.10   Max.   :29.400   Max.   :6.200   Max.   :3.3000  
##  NA's   :44      NA's   :44       NA's   :44      NA's   :33      
##       q5c2             q5c3            q5c4            q5c5      
##  Min.   : 0.000   Min.   : 0.00   Min.   : 0.00   Min.   : 0.00  
##  1st Qu.: 4.300   1st Qu.:42.70   1st Qu.:19.20   1st Qu.: 2.30  
##  Median : 7.450   Median :51.25   Median :25.95   Median : 4.20  
##  Mean   : 8.109   Mean   :50.24   Mean   :26.64   Mean   : 5.79  
##  3rd Qu.:11.500   3rd Qu.:59.40   3rd Qu.:33.60   3rd Qu.: 7.70  
##  Max.   :26.600   Max.   :69.30   Max.   :53.90   Max.   :35.70  
##  NA's   :33       NA's   :33      NA's   :33      NA's   :33     
##       q5c6             cpi             unemp         nav_high        
##  Min.   : 0.000   Min.   : 18.20   Min.   :1.000   Length:679        
##  1st Qu.: 4.300   1st Qu.: 49.75   1st Qu.:1.600   Class :character  
##  Median : 7.000   Median : 87.30   Median :2.600   Mode  :character  
##  Mean   : 8.405   Mean   : 74.19   Mean   :2.807                     
##  3rd Qu.:10.700   3rd Qu.: 96.70   3rd Qu.:4.000                     
##  Max.   :23.000   Max.   :100.70   Max.   :5.500                     
##  NA's   :33                                                          
##       ldp       
##  Min.   :11.70  
##  1st Qu.:23.12  
##  Median :28.20  
##  Mean   :28.29  
##  3rd Qu.:33.38  
##  Max.   :52.90  
##  NA's   :1

summary()コマンドで分かる統計量は以下の4種類のみ

  • 平均
  • 最大値：観察
  • 最小値
  • 最小値

より詳しい記述統計を算出するには、別のコマンドが必要。その前に・・・

変数の作成

「q3c1-6」、「Q5c1-6」の各変数をもとに経済評価変数を作成しておこう。

  • 「q3c1-6」は次の質問と回答をもとにしたもの(暮らし向き評価) 質問文：あなたの暮らし向きは、昨年の今ごろと比べてどうですか。楽になっていますか、苦しくなっていますか。

    • 大変楽になった(q3c1)
    • やや楽になった(q3c2)
    • 変わりない(q3c3)
    • やや苦しくなった(q3c4)
    • 大変苦しくなった(q3c5)
    • わからない(q3c6)

  • 「q5c1-6」は次の質問と回答をもとにしたもの(景気評価) 質問文：あなたの暮らし向きは、昨年の今ごろと比べてどうですか。楽になっていますか、苦しくなっていますか。

    • 確かに良くなってきたと思う(q5c1)
    • やや良くなってきたと思う(q5c2)
    • 変わらないと思うい(q5c3)
    • やや悪くなってきたと思う(q5c4)
    • 確かに悪くなってきたと思う(q5c5)
    • わからない(q5c6)

暮らし向き評価の「好評価」と「悪評価」を表す変数を作るためにはどうしたらいいだろうか？

#暮らし向き好評価
life_good<-q3c1+q3c2

#暮らし向き悪評価
life_bad<-q3c4+q3c5

#景気好評価
soc_good<-q5c1+q5c2
  
#景気悪評価
soc_bad<-q5c4+q5c5

上記の変数を既存のデータセットに組み込むためには、「$」記号を利用する。

#暮らし向き好評価
cabinet$life_good<-q3c1+q3c2

#暮らし向き悪評価
cabinet$life_bad<-q3c4+q3c5

#景気好評価
cabinet$soc_good<-q5c1+q5c2
  
#景気悪評価
cabinet$soc_bad<-q5c4+q5c5

経済評価変数を組み込んだうえで、以下の分析を進めよう。

より詳しい記述統計

summary()コマンドでは、ばらつき、分布の特徴についてわかる情報は限られている。より詳しい記述統計を得るために、psychパッケージ内のコマンドを利用する。

library(psych)
describe(cabinet)
##            vars   n   mean     sd median trimmed    mad  min   max range  skew
## yearMon*      1 679 340.00 196.15 340.00  340.00 252.04  1.0 679.0 678.0  0.00
## time_id       2 679 340.00 196.15 340.00  340.00 252.04  1.0 679.0 678.0  0.00
## approve       3 679  36.44  11.24  37.60   36.36  10.38  4.4  78.4  74.0  0.15
## dis           4 679  36.56  12.40  34.40   35.68   9.64  6.0  83.7  77.7  0.72
## ruling_ldp    5 679   0.85   0.52   1.00    1.00   0.00 -1.0   1.0   2.0 -3.26
## q3c1          6 635   0.61   1.64   0.30    0.26   0.15  0.0  11.5  11.5  4.85
## q3c2          7 635   4.60   2.67   4.00    4.27   1.48  0.0  14.8  14.8  1.25
## q3c3          8 635  56.18  11.24  58.50   57.20  11.12  0.0  76.7  76.7 -1.39
## q3c4          9 635  31.11   7.57  29.90   30.87   7.41  0.0  53.1  53.1 -0.07
## q3c5         10 635   5.53   3.78   4.90    5.08   2.37  0.0  29.4  29.4  2.44
## q3c6         11 635   1.20   0.99   0.80    1.06   0.74  0.0   6.2   6.2  1.13
## q5c1         12 646   0.35   0.45   0.20    0.26   0.15  0.0   3.3   3.3  2.25
## q5c2         13 646   8.11   4.87   7.45    7.77   5.19  0.0  26.6  26.6  0.60
## q5c3         14 646  50.24  11.25  51.25   50.99  12.53  0.0  69.3  69.3 -0.74
## q5c4         15 646  26.64   9.39  25.95   26.30  10.67  0.0  53.9  53.9  0.27
## q5c5         16 646   5.79   4.97   4.20    4.94   3.26  0.0  35.7  35.7  2.09
## q5c6         17 646   8.40   5.04   7.00    7.77   4.30  0.0  23.0  23.0  0.94
## cpi          18 679  74.19  28.85  87.30   77.52  15.86 18.2 100.7  82.5 -0.85
## unemp        19 679   2.81   1.31   2.60    2.72   1.78  1.0   5.5   4.5  0.42
## nav_high*    20 679 340.00 196.15 340.00  340.00 252.04  1.0 679.0 678.0  0.00
## ldp          21 678  28.29   6.70  28.20   28.36   7.56 11.7  52.9  41.2 -0.06
## life_good    22 635   5.20   2.67   4.40    4.80   1.63  0.0  15.2  15.2  1.39
## life_bad     23 635  36.64  10.13  34.90   36.05   9.19  0.0  74.1  74.1  0.46
## soc_good     24 646   8.46   5.19   7.60    8.07   5.41  0.0  28.4  28.4  0.66
## soc_bad      25 646  32.43  13.62  30.15   31.36  14.31  0.0  77.1  77.1  0.64
##            kurtosis   se
## yearMon*      -1.21 7.53
## time_id       -1.21 7.53
## approve        0.36 0.43
## dis            0.60 0.48
## ruling_ldp     8.63 0.02
## q3c1          23.45 0.06
## q3c2           1.92 0.11
## q3c3           3.94 0.45
## q3c4           1.25 0.30
## q3c5           9.77 0.15
## q3c6           0.77 0.04
## q5c1           5.90 0.02
## q5c2          -0.19 0.19
## q5c3           0.92 0.44
## q5c4          -0.56 0.37
## q5c5           6.35 0.20
## q5c6          -0.18 0.20
## cpi           -0.95 1.11
## unemp         -1.05 0.05
## nav_high*     -1.21 7.53
## ldp           -0.39 0.26
## life_good      1.66 0.11
## life_bad       1.66 0.40
## soc_good      -0.08 0.20
## soc_bad       -0.08 0.54
  • vars：変数の列番号
  • n: 観察数
  • mean：平均値\(\mu=\frac{1}{n}\sum^{n}_{i-1}x_i\)
  • sd: 標準偏差 \(sd=\sqrt{\frac{1}{n-1}\sum^{n}_{i=1}(x_{i}-\mu)^2}\)

  • median：中央値
  • timmed：刈り込み平均
  • mad: 中央絶対偏差値
  • min: 最小値
  • max: 最大値
  • range: 範囲
  • skew: 歪度 \(skewness=\frac{n}{(n-1)(n-2)}\sum^{n}_{i=1}(\frac{x_{i}-\mu}{s})^3\)
##   Unsuccessful convergence.
##   Unsuccessful convergence.
##   Successful convergence.

  • kurtosis: 尖度とても複雑な算出式！
##   Successful convergence.
##   Unsuccessful convergence.
##   Successful convergence.

  • se: 標準誤差

記述統計からどのようなことが明らかになるだろうか？

3．変数の分布の確認

内閣支持率、与党支持率、経済変数、それぞれにはどのような分布の特徴があるだろうか？

上に挙げたヒストグラムと同じものを、各自で描画してみよう。

①まずは単純な頻度を表すヒストグラムを描画する。

hist(approve)

②次に、密度をもとにヒストグラムを描画してみよう。

hist(approve, freq=F) #freq=Fが「頻度で描画しない」という指示

③ここに密度関数の曲線を併記しよう。

hist(approve, freq=F)
lines(density(approve,lwd = 2))  #density()が密度関数。lines()によって密度関数の曲線をoverlayしている

④ヒストグラムのタイトルやX軸、y軸のラベルを変更してみよう。

hist(approve, freq=F, main="内閣支持率のヒストグラム",xlab="値", ylab="密度")
lines(density(approve,lwd = 2))

内閣支持率の分布について、どのような特徴が見てとれるだろうか。そして分布の形状を、記述統計と照らし合わせて考えてみよう。

そして他の変数についても、各自でヒストグラムを作成し、結果を解釈していこう。以下と同じような図形が描画できただろうか？

II.相関分析

2変数間の関係を確かめるために相関を調べる。ここでは、内閣支持率と自民党支持率の関係について検討することから始めよう。

1.散布図による確認

散布図は次のようなplot()関数によって、シンプルに描画することが可能。

plot(ldp,approve)

これもヒストグラムの時と同様に、各軸とメインタイトルの名前を変えることができる。

plot(ldp,approve, xlab="自民党支持率",ylab="内閣支持率",main="内閣支持率と自民党支持率の関係に関する散布図")

ここでもう1つのアドバンスな図の描画にチャレンジしよう。「ggplot2」というパッケージを使うことで、より見やすく、効率的な図の作成ができる。

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

g <- ggplot(data=cabinet, aes(x = ldp, y = approve)) #gというオブジェクトにグラフの土台となる情報を格納
g <- g + geom_point() # gというオブジェクトを、geom_pointコマンドによって散布図で出力
plot(g) #完成しているオブジェクトgを描画

軸ラベルとメインタイトルを変える時には、デフォルトのプロットとは異なり以下のように設定する。

g <- ggplot(data=cabinet, aes(x = ldp, y = approve)) 
g <- g + geom_point()+xlab("自民党支持率")+ylab("内閣支持率")+
  ggtitle("内閣支持率と自民党支持率の関係に関する散布図") #xlab()がx軸、ylab()がy軸それぞれの軸ラベルで、ggtitle()がメインタイトル
plot(g)

また、ggplot2で図を描画する際に、目盛りの幅(scale of axis)を変更することもできる。

g <- ggplot(data=cabinet, aes(x = ldp, y = approve)) 
g <- g + geom_point()+xlab("自民党支持率")+ylab("内閣支持率")+
  ggtitle("内閣支持率と自民党支持率の関係に関する散布図")+ #xlab()がx軸、ylab()がy軸それぞれの軸ラベルで、ggtitle()がメインタイトル
  scale_x_continuous(breaks=seq(0,80,5)) #x軸の目盛りを「0から80までにして、目盛りの区切りを5に変更する」という意味
plot(g)

このように見てくると、内閣支持率と自民党支持率の間には何らかの関係がありそうに見える。これを相関分析によって確かめていこう。

2.ピアソンの積率相関分析

基本統計量と相関係数の関係

もう1度基本的な統計量を確認しておこう。

  • 1変数の平均：\(\bar{x}=\frac{1}{n}\sum^{n}_{i=1}x_{i}\).
  • 1変数の不偏分散：\(s^{2}=\frac{1}{n-1}\sum^{n}_{i-1}(x_{i}-\bar{x})^2\).
  • 2変数間の共分散：\(S_{xy}=\frac{1}{n-1}\sum^{n}_{i-1}(x_{i}-\bar{x})(y_{i}-\bar{x})\).

このもとで相関係数\(r_{xy}\)は、

\[r_{xy}=\frac{s_{xy}}{s_{x}s_{y}}=\frac{\sum^n_{i=1}(x_{i}-\bar{x})(y_i-\bar{y})}{\sqrt{\sum^{n}_{i=1}(x_{i}-\bar{x})^2\sum^{n}_{i=1}(y_{i}-\bar{y})^2}}\].

内閣支持率と自民党支持率の相関分析

内閣支持率と自民党支持率の間のピアソン積率相関分析を計算しよう。

cor.test(approve,ldp)
## 
##  Pearson's product-moment correlation
## 
## data:  approve and ldp
## t = 7.3349, df = 676, p-value = 6.373e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2003112 0.3398605
## sample estimates:
##       cor 
## 0.2715123

分析結果は次のような構成になっている。

t検定

ここで統計的検定、t検定のことを復習しておこう。相関係数に関するt検定の帰無仮説は、

  • 帰無仮説(\(H_0\))：真の相関係数は0
  • 対抗仮説(\(H_A\))：真の相関係数は0ではない

t検定のための、\(r_{xy}\)のt統計量は以下の式で求まる。

\[t=r\sqrt{\frac{n-2}{1-r^2}}\]

母集団の中から無作為に抽出された標本の相関係数はt分布に従うことから、t分布をもとに上記のt統計量を評価する。t統計量の評価の方法を覚えているだろうか。

上記の図に、自分で1%水準で統計的に有意となる場合のt値のラインを描画してみよう。

内閣支持率と自民党支持率の相関分析の推定結果

  • 続いて、内閣支持率と自民党支持率の相関係数の検定結果について考えてみよう。

t = 7.3349, df = 676, p-value = 6.373e-13

という推定結果はどのように解釈できるだろうか？

  • 相関係数の値の解釈 \(r_{xy} > 0\)：正の相関 \(r_{xy} < 0\)：負の相関 \(r_{xy} = 0\)：無相関 \(r_{xy} = 1\)：正の完全相関 \(r_{xy} = -1\): 負の完全相関

内閣支持率と他の変数との相関

他の変数と内閣支持率との相関についても推定しよう。暮らし向き評価と景気評価についてはどうだろうか。

  • 内閣支持率と暮らし向き好評価の相関
## 
##  Pearson's product-moment correlation
## 
## data:  approve and life_good
## t = 2.5011, df = 633, p-value = 0.01263
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.0212820 0.1753795
## sample estimates:
##       cor 
## 0.0989238
  • 内閣不支持率と暮らし向き低評価の相関
## 
##  Pearson's product-moment correlation
## 
## data:  dis and life_bad
## t = 5.0256, df = 633, p-value = 0.0000006538
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1199037 0.2695789
## sample estimates:
##       cor 
## 0.1958819

  • 内閣支持率／不支持率、景気好評価値／悪評価の相関はどのように解釈できるだろうか？

  • 自民党支持率と各種経済評価の相関はどのように解釈できるだろうか？

3.質的変数の相関分析

データの読み込み

利用するデータ：東京大学朝日新聞共同調査・政治家調査　2017年衆議院議員候補者調査 ※コードブックも参照してください。

politician<-read.csv("politician2017.csv", fileEncoding="SJIS")
attach(politician)

ケンドール順位相関係数(Kendall’s rank correlation coefficient)

-安全保障に関する政策態度間の相関の検証

利用する質問文

  1. 日本の防衛力はもっと強化すべきだ (Q4_1)
  1. 賛成
  2. どちらかと言えば賛成
  3. どちらとも言えない
  4. どちらかと言えば反対
  5. 反対
  6. 無回答
  1. 他国からの攻撃が予想される場合には先制攻撃もためらうべきではない (Q4_2)
  1. 賛成
  2. どちらかと言えば賛成
  3. どちらとも言えない
  4. どちらかと言えば反対
  5. 反対
  6. 無回答

①「無回答」の「99」という値を欠損値に変換する。

describe(Q4_1)
##    vars   n mean    sd median trimmed  mad min max range skew kurtosis  se
## X1    1 465 5.03 17.21      2    1.78 1.48   1  99    98 5.24    25.67 0.8

最大値が「99」であるが、この数値に順序としての意味はないので数値を欠損値として扱う必要がある。よって変数の変換を行う。

defense<-Q4_1 #Q4_1変数をもとに、新たに「defense」という名前の変数を作成
defense[Q4_1==99]<-NA #「Q4_1=99」にあたる観察を「NA」に置き換える

describe(defense)
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 450  1.9 1.13      2    1.69 1.48   1   5     4 1.25     0.81 0.05

同様に先制攻撃についても変換。

premptive<-Q4_2 #Q4_2変数をもとに、新たに「premptive」という名前の変数を作成
premptive[Q4_2==99]<-NA #「Q4_2=99」にあたる観察を「NA」に置き換える

describe(premptive)
##    vars   n mean  sd median trimmed  mad min max range skew kurtosis   se
## X1    1 440 3.33 1.1      3    3.35 1.48   1   5     4 0.06    -0.58 0.05

②質的変数であることから、先ほどのピアソン積率相関係数を利用することは適切ではない。よって質的な変数のための、ケンドール順位相関係数によって、質的変数の相関を計算する。

2つの変量、\(x\)と\(y\)がある時、対応する2つの観察\((x_s,y_s)\)と\((x_t,y_t)\)を取り出したときに、

\(P\):\(x_s\)と\(x_t\)、\(y_s\)と\(y_t\)の大小関係が同じ向きである組の数

\(Q\) =\(x_s\)と\(x_t\)と\(y_s\)と\(y_t\)、との大小関係が異なる向きである組の数

\(T_x\) = \(x_s=x_t\)である組の数

\(T_y\) =\(y_s=y_t\)である組の数

\(N\)=\(\frac{n(n-1)}{2}\) （組の総数）

であると設定すると、ケンドールの順位相関係数\(\tau\)は、次式によって定義されてる。 \[ \tau = \frac{P - Q} {\sqrt{N - T_x}\sqrt{N - T_y}} \]

cor.test(defense,premptive,method="k")
## 
##  Kendall's rank correlation tau
## 
## data:  defense and premptive
## z = 12.112, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.4934934

ケンドールの順位相関検定では、z検定により有意性を評価する。政策態度、防衛強化と先制攻撃容認の間の相関関係はどのように解釈できるだろうか。

③では、上記の作業をもとにして、「小さな政府」志向の政策態度と防衛強化の相関を計算してみよう。

## 
##  Kendall's rank correlation tau
## 
## data:  defense and welfare
## z = 11.966, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.4948624

上記と同じ推定結果が得られただろうか。そしてこの結果はどのように解釈できるだろうか。