データの読み込み

今日使うデータを読み込む。まずは.csvデータを読み込んでみよう。.csv形式が、最も一般的なデータのファイル形式だというのはわかってもらえるだろう。その分、最もシンプルにデータを読み込める。Rを使って共同で分析をするときに、なるべくなら.csv形式でやり取りするように心がけよう。

df_cabinet<-read.csv("cabinet.csv")
#行っていること：cabinet.csvファイルを読み込み、df_cabinetにデータフレーム・オブジェクトとして格納

他にもさまざまな形式のデータが存在する。.xlsxもおなじみのデータの1つ。ではどのように読み込めばいいだろうか。まずは、\(\texttt{readxl}\)パッケージを使って、read_excel()関数を使う場合。

# パッケージの読み込み
library(readxl) #最初のインストールは、install.packages("readxl")で行う

# xlsxファイルなので、複数のシート名が存在。シート名を指定
df_cabinet <- read_excel("cabinet.xlsx", sheet = "Sheet1")

次は、\(\texttt{openxlsx}\)パッケージから、read.xlsx()関数を使う場合。

library(openxlsx)

df_cabinet<- read.xlsx("cabinet.xlsx", sheet = 1)

\(\texttt{xlsx}\)パッケージから`read.xlsx()関数を使う場合

library(xlsx)

df_cabinet <- read.xlsx("cabinet.xlsx")

最後に、\(\texttt{rio}\)パッケージから`import()関数を使う場合

library(rio)

df_cabinet <- import("cabinet.xlsx")

また、データが統計ソフトウェアStataを使って作成され、.dta形式で配布されることもしばしばある。その場合には、以下のように、\(\texttt{foreign}\)パッケージから`read.dta()関数を使う場合

library(foreign)

df_cabinet <- read.dta("cabinet.dta")

他にも、\(\texttt{haven}\)パッケージから`read_dta()関数を使う場合がある。

library(haven)

data_haven <- read_dta("cabinet.dta")

いずれの関数でも、_と.の使い分けがあるので、その点は注意しておこう。少しのスペルミス、大文字小文字の判別ミスもあいまいにできないのがRコマンドの入力だと理解しておいてほしい。

最も簡単な記述統計

まずは最も簡単な記述統計を出してもらう。

summary(df_cabinet)

##    yearMon             time_id         approve           dis       
##  Length:444         Min.   :236.0   Min.   : 4.40   Min.   : 6.00  
##  Class :character   1st Qu.:346.8   1st Qu.:29.07   1st Qu.:29.77  
##  Mode  :character   Median :457.5   Median :38.80   Median :35.40  
##                     Mean   :457.5   Mean   :37.72   Mean   :37.96  
##                     3rd Qu.:568.2   3rd Qu.:45.62   3rd Qu.:44.02  
##                     Max.   :679.0   Max.   :78.40   Max.   :83.70  
##                                                                    
##    ruling_ldp           q3c1            q3c2            q3c3      
##  Min.   :-1.0000   Min.   :0.000   Min.   :1.300   Min.   :32.90  
##  1st Qu.: 1.0000   1st Qu.:0.100   1st Qu.:3.100   1st Qu.:56.77  
##  Median : 1.0000   Median :0.200   Median :3.800   Median :62.45  
##  Mean   : 0.7748   Mean   :0.238   Mean   :3.911   Mean   :60.92  
##  3rd Qu.: 1.0000   3rd Qu.:0.300   3rd Qu.:4.600   3rd Qu.:66.50  
##  Max.   : 1.0000   Max.   :0.800   Max.   :7.700   Max.   :76.70  
##                    NA's   :2       NA's   :2       NA's   :2      
##       q3c4            q3c5             q3c6             q5c1       
##  Min.   :16.50   Min.   : 1.000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:24.90   1st Qu.: 3.325   1st Qu.:0.4000   1st Qu.:0.1000  
##  Median :28.55   Median : 4.600   Median :0.6000   Median :0.1000  
##  Mean   :29.32   Mean   : 4.941   Mean   :0.6771   Mean   :0.2527  
##  3rd Qu.:32.80   3rd Qu.: 6.175   3rd Qu.:0.9000   3rd Qu.:0.3000  
##  Max.   :47.10   Max.   :16.000   Max.   :6.2000   Max.   :2.0000  
##  NA's   :2       NA's   :2        NA's   :2        NA's   :2       
##       q5c2             q5c3            q5c4            q5c5       
##  Min.   : 0.400   Min.   :18.60   Min.   : 9.40   Min.   : 0.100  
##  1st Qu.: 3.600   1st Qu.:47.40   1st Qu.:19.60   1st Qu.: 2.300  
##  Median : 6.400   Median :56.60   Median :26.40   Median : 4.200  
##  Mean   : 7.016   Mean   :53.96   Mean   :27.27   Mean   : 5.888  
##  3rd Qu.: 9.775   3rd Qu.:61.50   3rd Qu.:34.20   3rd Qu.: 7.700  
##  Max.   :19.500   Max.   :69.30   Max.   :53.90   Max.   :35.700  
##  NA's   :2        NA's   :2       NA's   :2       NA's   :2       
##       q5c6             cpi             unemp         nav_high        
##  Min.   : 2.000   Min.   : 72.60   Min.   :1.900   Length:444        
##  1st Qu.: 3.925   1st Qu.: 88.78   1st Qu.:2.600   Class :character  
##  Median : 4.900   Median : 96.10   Median :3.400   Mode  :character  
##  Mean   : 5.611   Mean   : 93.52   Mean   :3.512                     
##  3rd Qu.: 7.475   3rd Qu.: 97.80   3rd Qu.:4.500                     
##  Max.   :11.600   Max.   :100.70   Max.   :5.500                     
##  NA's   :2                                                           
##       ldp       
##  Min.   :11.70  
##  1st Qu.:21.80  
##  Median :25.20  
##  Mean   :26.03  
##  3rd Qu.:31.30  
##  Max.   :40.20  
##  NA's   :1

ここで、summary()コマンドによって算出できる統計量は以下の4種類のみ

• Mean （平均）
• Median （中央値）
• Max. （最大値）
• Min. （最小値）
• 1st Qu（第1四分位、25％目の値）
• 3st Qu（第3四分位、75％目の値）

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

変数の作成

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

• 「q3c1-6」は次の質問と回答をもとにしたもの(暮らし向き評価)

質問文：あなたの暮らし向きは、昨年の今ごろと比べてどうですか。楽になっていますか、苦しくなっていますか。

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

• 「q5c1-6」は次の質問と回答をもとにしたもの(景気評価)

質問文：世間の景気をどう見ますか。先月と変わらないと思いますか、悪くなってきたと思いますか、良くなってきたと思いますか。

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

暮らし向き評価の「肯定的評価」と「否定的評価」を表す変数を作るためにはどうしたらいいだろうか？最も簡単に、これらの変数を既存のデータセットに組み込むためには、「$」（ドル記号・あるオブジェクトの中から特定の要素抜き出すための演算子）を利用する。

#暮らし向き肯定的評価
df_cabinet$life_good <- df_cabinet$q3c1 + df_cabinet$q3c2

#暮らし向き否定的評価
df_cabinet$life_bad <- df_cabinet$q3c4 + df_cabinet$q3c5

#景気肯定的評価
df_cabinet$soc_good <- df_cabinet$q5c1 + df_cabinet$q5c2
  
#景気否定的評価
df_cabinet$soc_bad <- df_cabinet$q5c4 + df_cabinet$q5c5  

しかしこれだと、繰り返し、\(\texttt{df_cabinet}\)を書き続けなければならない。効率性において劣るコードとなる。そこで、効率的なコードの作成、そしてtidyなデータの構築を、この授業でも習得していこう。

そのために不可欠なのが、\(\texttt{tidyverse}\) package。tidyなデータの作成、データの美しい可視化に必要なパッケージをさらにパッケージ化したものと理解するとよい。最初に、\(\texttt{tidyverse}\)パッケージをしっかりインストールしておこう。

なお、\(\texttt{tidyverse}\)の中に含まれるパッケージは次の通り。

tidyverseに含まれるパッケージの一覧

library(tidyverse) #tidyverseの読み込み

df_cabinet <- df_cabinet |>
  mutate(　#tidyverse内のmutate関数を使って、データを補完する作業を行うよう指示
    # 暮らし向き肯定的評価
    life_good = q3c1 + q3c2,
    # 暮らし向き否定的評価
    life_bad  = q3c4 + q3c5,
    # 景気肯定的評価
    soc_good  = q5c1 + q5c2,
    # 景気否定的評価
    soc_bad   = q5c4 + q5c5
  )

#df_cabinetデータの最初5行を確認しよう。新たな変数が格納されているだろうか？
head(df_cabinet)

##      yearMon time_id approve  dis ruling_ldp q3c1 q3c2 q3c3 q3c4 q3c5 q3c6 q5c1
## 1 1980-01-01     236    27.4 44.6          1  0.2  3.8 47.0 40.1  7.0  2.0  0.2
## 2 1980-02-01     237    27.6 44.0          1  0.0  3.0 40.4 45.1 10.3  1.2  0.4
## 3 1980-03-01     238    25.5 45.1          1  0.1  3.1 36.3 47.1 12.4  1.1  0.1
## 4 1980-04-01     239    25.8 47.2          1  0.4  3.5 32.9 45.5 16.0  1.6  0.2
## 5 1980-05-01     240    29.1 42.1          1  0.3  3.3 37.0 45.1 13.2  1.1  0.0
## 6 1980-06-01     241    27.3 43.9          1  0.1  3.8 40.6 44.0 10.1  1.4  0.1
##   q5c2 q5c3 q5c4 q5c5 q5c6  cpi unemp nav_high  ldp life_good life_bad soc_good
## 1  5.3 42.0 36.2  5.7 10.7 72.6   1.9 6,776.60 30.1       4.0     47.1      5.5
## 2  4.4 38.5 39.6  8.4  8.8 73.1   1.9 6,839.86 29.8       3.0     55.4      4.8
## 3  5.9 45.2 33.9  5.9  9.0 73.5   1.9 6,794.04 28.9       3.2     59.5      6.0
## 4  7.8 39.2 35.3  7.7  9.8 74.9   2.0 6,904.81 31.1       3.9     61.5      8.0
## 5  7.5 42.3 31.8  6.8 11.6 75.6   2.0 6,882.65 28.9       3.6     58.3      7.5
## 6  7.7 49.7 27.0  5.5 10.0 75.9   1.9 6,870.70 30.0       3.9     54.1      7.8
##   soc_bad
## 1    41.9
## 2    48.0
## 3    39.8
## 4    43.0
## 5    38.6
## 6    32.5

経済評価変数を組み込むことができたので、以下でより詳しい記述統計の分析を進めよう。

より詳しい記述統計

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

library(psych)
describe(df_cabinet)

##            vars   n   mean     sd median trimmed    mad   min   max range  skew
## yearMon*      1 444 222.50 128.32 222.50  222.50 164.57   1.0 444.0 443.0  0.00
## time_id       2 444 457.50 128.32 457.50  457.50 164.57 236.0 679.0 443.0  0.00
## approve       3 444  37.72  12.11  38.80   37.69  11.79   4.4  78.4  74.0  0.07
## dis           4 444  37.96  13.08  35.40   37.07  10.16   6.0  83.7  77.7  0.68
## ruling_ldp    5 444   0.77   0.63   1.00    0.97   0.00  -1.0   1.0   2.0 -2.44
## q3c1          6 442   0.24   0.14   0.20    0.22   0.15   0.0   0.8   0.8  0.83
## q3c2          7 442   3.91   1.20   3.80    3.85   1.19   1.3   7.7   6.4  0.48
## q3c3          8 442  60.92   7.54  62.45   61.62   6.60  32.9  76.7  43.8 -0.85
## q3c4          9 442  29.32   6.11  28.55   28.92   5.71  16.5  47.1  30.6  0.58
## q3c5         10 442   4.94   2.20   4.60    4.73   1.93   1.0  16.0  15.0  1.07
## q3c6         11 442   0.68   0.49   0.60    0.62   0.44   0.0   6.2   6.2  3.74
## q5c1         12 442   0.25   0.34   0.10    0.18   0.15   0.0   2.0   2.0  2.72
## q5c2         13 442   7.02   4.40   6.40    6.62   4.60   0.4  19.5  19.1  0.68
## q5c3         14 442  53.96   9.96  56.60   54.94   9.19  18.6  69.3  50.7 -0.87
## q5c4         15 442  27.27   9.58  26.40   26.94  10.97   9.4  53.9  44.5  0.31
## q5c5         16 442   5.89   5.22   4.20    4.98   3.41   0.1  35.7  35.6  2.17
## q5c6         17 442   5.61   2.13   4.90    5.43   1.78   2.0  11.6   9.6  0.72
## cpi          18 444  93.52   6.43  96.10   94.47   3.04  72.6 100.7  28.1 -1.21
## unemp        19 444   3.51   1.06   3.40    3.47   1.33   1.9   5.5   3.6  0.22
## nav_high*    20 444 222.50 128.32 222.50  222.50 164.57   1.0 444.0 443.0  0.00
## ldp          21 443  26.03   6.19  25.20   26.15   6.67  11.7  40.2  28.5 -0.04
## life_good    22 442   4.15   1.26   4.00    4.09   1.19   1.5   8.4   6.9  0.50
## life_bad     23 442  34.26   8.00  33.05   33.66   7.86  19.5  61.5  42.0  0.68
## soc_good     24 442   7.27   4.64   6.50    6.82   4.67   0.4  21.1  20.7  0.76
## soc_bad      25 442  33.16  14.01  30.45   32.03  14.53  10.2  77.1  66.9  0.66
##            kurtosis   se
## yearMon*      -1.21 6.09
## time_id       -1.21 6.09
## approve        0.18 0.57
## dis            0.42 0.62
## ruling_ldp     3.98 0.03
## q3c1           0.57 0.01
## q3c2           0.19 0.06
## q3c3           0.42 0.36
## q3c4          -0.25 0.29
## q3c5           1.75 0.10
## q3c6          35.72 0.02
## q5c1           8.28 0.02
## q5c2          -0.26 0.21
## q5c3           0.24 0.47
## q5c4          -0.80 0.46
## q5c5           6.77 0.25
## q5c6          -0.67 0.10
## cpi            0.44 0.31
## unemp         -1.26 0.05
## nav_high*     -1.21 6.09
## ldp           -0.62 0.29
## life_good      0.23 0.06
## life_bad       0.05 0.38
## soc_good      -0.05 0.22
## soc_bad       -0.22 0.67

この操作を、\(\texttt{tidyverse}\)に基づいて記述すると次のようになる。%>%(パイプ関数)を用いて、オブジェクト内への操作を連続的に指示できる。パイプ関数によって、「そのデータ（オブジェクト）内にパイプ関数以下の操作を加える」と伝えることになる。

df_cabinet %>% 　#df_cabinetデータ内に対して、describe()で記述統計を出すように指示する、という意味
  describe()

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

• kurtosis: 尖度、とても複雑な算出式！

ところで先の記述統計の出力結果を見せられても、読者は「？」となるばかり。分析結果は、表にして、一覧性がある状態で出力しなくてはならない。それに全部の変数の記述統計が必要ではないし、全部の統計量が必要なわけでもない。

library(flextable) #出力用の表作成に適したパッケージ
library(officer) #microsoftの機能と結びつけるためのパッケージ

#------ df_cabinetについての全部の記述統計を取得
desc_all <- psych::describe(df_cabinet[,-1])

#------ describe() の結果は行名に変数名が格納されている。ここに変数名という列名を付加
desc_tbl <- desc_all %>% 
  rownames_to_column(var = "variable")

#------ 出力対象となる変数を指定
needed_vars <- c("approve", "dis", "ldp", "cpi", "unemp", "life_good", "life_bad", "soc_good", "soc_bad")

#------ 必要な変数かつ、表示項目のみ抽出したデータフレームを作成（mean, sd, median, min, max, skew）
desc_subset <- desc_tbl %>% 
  filter(variable %in% needed_vars) %>%  #filter()変数でvariablesの中でneed_varsの行に当たるものを選ぶと指定
  dplyr::select(variable, mean, sd, median, min, max, skew) %>% #これらの統計量の列を選択すると指定
  mutate(across(where(is.numeric), ~ round(., 3))) #数字(numeric)に当たるものに関しては小数点3位以下に丸めると指定

#------ flextable を用いて見栄えの良い表を作成
ftdesc <- flextable(desc_subset)
ftdesc <- autofit(ftdesc)
ftdesc #右側の「Viewer」タブに表が表示される。

# officerパッケージでWord文書を作成し、flextableを挿入して保存
doc <- read_docx()
doc <- body_add_par(doc, "記述統計量", style = "heading 1")
doc <- body_add_flextable(doc, value = ftdesc)
print(doc, target = "記述統計.docx")

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

もう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}}\].

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

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

df_cabinet %>% 
  with(cor.test(approve, ldp))

## 
##  Pearson's product-moment correlation
## 
## data:  approve and ldp
## t = 6.2066, df = 441, p-value = 1.25e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1954265 0.3669104
## sample estimates:
##       cor 
## 0.2834326

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

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 = 6.2066,   df = 441,   p-value = 1.25e-09

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

• 相関係数の値の解釈

• 信頼区間(confidence intervals)

95%信頼区間とは？：ある試行（ここでは内閣支持率と自民党支持率を集め相関係数を求めるという試行）を繰り返したときに、各試行で得られる推定値の95%が収まる区間の推定

⇒無限回施行を繰り返せば、真の値は95%の確率でこの範囲内に位置する

どういうことを意味しているのか、実際にデータを1000回リサンプリングすることで確かめてみよう。

bootstrapで、「相関係数が信頼区間95%内に95%収まるように」と指示しているわけではないので、1000回中950回の試行で相関係数が信頼区間内に収まったのは偶然。しかし、まさに95%のサンプルにおいて、相関係数が信頼区間95%の中に収まることを実感してもらえたのではないだろうか。

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

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

• 内閣支持率と暮らし向き好評価の相関

## 
##  Pearson's product-moment correlation
## 
## data:  approve and life_good
## t = 4.9528, df = 440, p-value = 0.000001045
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1395139 0.3162886
## sample estimates:
##       cor 
## 0.2297957

• 内閣不支持率と暮らし向き低評価の相関

## 
##  Pearson's product-moment correlation
## 
## data:  dis and life_bad
## t = 4.5065, df = 440, p-value = 0.000008458
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1191078 0.2974902
## sample estimates:
##       cor 
## 0.2100464

データの読み込み

利用するデータ：東京大学朝日新聞共同調査・政治家調査　2017年衆議院議員候補者調査

※コードブックを参照。

まずはデータの読み込み。fileEncoding="SJIS"は何を指示しているのだろうか。いちど、fileEncoding="SJIS"をつけずにデータを読み込んでみよう。どのようなエラーが出るだろうか。

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

データの前処理と自作関数の作成

このデータをもとに、安全保障に関する態度と価値観との間の相関について、検討していこう。政治家に対するサーベイデータである「東京大学朝日新聞共同調査・政治家調査　2017年衆議院議員候補者調査」は、主に質的な変数（第1回授業資料参照）から構成されている。よって、質的な変数を分析に利用するためのデータの前処理にまずは取り組まなくてはならない。

ここでは、「防衛力に関する態度」、「先制攻撃に対する態度」、「首相の靖国神社参拝に対する態度」の間の関係・相関を分析する。それぞれの質問文について、まずは整理しておこう。

【利用する質問文】

(1) 日本の防衛力はもっと強化すべきだ (Q4_1)

1. 賛成

2. どちらかと言えば賛成

3. どちらとも言えない

4. どちらかと言えば反対

5. 反対

99. 無回答

(2) 他国からの攻撃が予想される場合には先制攻撃もためらうべきではない (Q4_2)

1. 賛成

2. どちらかと言えば賛成

3. どちらとも言えない

4. どちらかと言えば反対

5. 反対

99. 無回答

(5) 首相には靖国神社に参拝してほしい（Q4_5）

1. 賛成

2. どちらかと言えば賛成

3. どちらとも言えない

4. どちらかと言えば反対

5. 反対

99. 無回答

この状態でも、分析に用いることができそうだが、いくつか問題がある。まず「99」は欠損値に当てはめられた意味のない値なので、\(\texttt{NA}\)に変換する必要がある。さらに、賛成が1、、反対が5となっているが、賛成の強度を表す方が直感にも適合的なので、スコアを逆の順に付す必要がある。

df_pol %>%
with(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

\(\texttt{tidyverse}\)に含まれる\(\texttt{dplyr}\)内のcase_when()関数を使って、変数の再コード化(recode)を行う。

df_pol <- df_pol %>%
  mutate( #mutateで以下の操作の変数を埋めていく
    # Q4_1 から defense を作成
    defense = case_when(
      Q4_1 == 1  ~ 5, #逆順に指定していっているのが分かるだろう
      Q4_1 == 2  ~ 4,
      Q4_1 == 3  ~ 3,
      Q4_1 == 4  ~ 2,
      Q4_1 == 5  ~ 1,
      Q4_1 == 99 ~ NA_real_,　#欠損値の指定
      TRUE       ~ NA_real_
    ),
    # Q4_2 から preemptive を作成
    preemptive = case_when(
      Q4_2 == 1  ~ 5,
      Q4_2 == 2  ~ 4,
      Q4_2 == 3  ~ 3,
      Q4_2 == 4  ~ 2,
      Q4_2 == 5  ~ 1,
      Q4_2 == 99 ~ NA_real_,
      TRUE       ~ NA_real_
    ),
    # Q4_5 から yasukuni を作成
    yasukuni = case_when(
      Q4_5 == 1  ~ 5,
      Q4_5 == 2  ~ 4,
      Q4_5 == 3  ~ 3,
      Q4_5 == 4  ~ 2,
      Q4_5 == 5  ~ 1,
      Q4_5 == 99 ~ NA_real_,
      TRUE       ~ NA_real_
    )
  )

ここで、Rにおいて極めて便利な機能であるfunction()設定について紹介する。Rでは関数を自作できる。

例えば、上のコードは、3つの変数を作るためにかなり長いものになっている。3変数ならまだいいかもしれないが、政策評価に関わる変数は、10以上存在している（コードブック参照）。意識調査でよくみられるマトリクス形式の質問であれば、さらに多くの変数群の前処理を、同じコードで繰り返すことになる。その時コードは冗長となり、効率化が損なわれる。自作関数を作成し、その関数のもとに作業を繰り返す効率化された、洗練されたコードを考えてみたい。

まずは、上のcase_when()を使ったコードを関数化しよう。関数を作成するためのコードは、function(x){}である。xを関数に代入したときに、何を返すか（return）を指示する仕組みになっている。ここでは、recode_var()という関数を自作してみよう。

recode_var <- function(x) {
  case_when(
    x == 1  ~ 5,
    x == 2  ~ 4,
    x == 3  ~ 3,
    x == 4  ~ 2,
    x == 5  ~ 1,
    x == 99 ~ NA_real_,
    TRUE    ~ NA_real_
  )
}

この関数は短いものだが、長大な関数になることもしばしばある。その場合に、scriptに長い関数が書かれるよりも、簡便に関数を呼び出す方が追試の際の利便性も格段に上がる。

そこで、関数の設定後に関数を.Rファイルとして保存し、呼び出す方法を知っておこう。まずは以下のように、Scriptファイルを新たに立ち上げる。

次に、関数の部分を貼り付ける。それを関数と同じ名前のファイルとして保存する。

そして分析に使うScript内で以下のように記述すれば、自作関数を呼び出せる。但し、自作関数は同じ作業ディレクトリ内に収めておかなくてはならない。

source("recode_var.R")

この関数がどのように働くか、確かめてみたい。いま、コードブックを見ると、以下のような変数がある。これを\(\texttt{welfare}\)変数として定めるべく、上記の関数recode_var()を利用する。

(6) 社会福祉など政府のサービスが悪くなっても、お金のかからない小さな政府の方が良い (Q4_6)

1. 賛成

2. どちらかと言えば賛成

3. どちらとも言えない

4. どちらかと言えば反対

5. 反対

99. 無回答

df_pol <- df_pol %>%
  mutate( 
    welfare = recode_var(Q4_6)
  )

既述のcase_when()を使ったコードをより、はるかにコードが短くなり効率化されているのが分かる。これを複数の変数に繰り返し適用できるとよい。

こうした繰り返し作業を可能にするのが、for(i in ...)によるループ関数である。まずは操作を行いたい変数名のベクトル（リスト）を作成する。次に新しくつけたい名前のベクトルも指示する。その後にループ関数を適用する。

#------- 操作を行いたい変数の名前をまとめるベクトル（リスト）の作成
orig_vars <- c("Q4_1", "Q4_2", "Q4_5")

#------- 新しい変数名のベクトル（リスト）の作成
new_vars  <- c("defense", "preemptive", "yasukuni")

#------- forループでまとめて処理して、それらにrecode_var()関数を適用
for(i in length(orig_vars)) {
  df_pol[[ new_vars[i] ]] <- recode_var(df_pol[[ orig_vars[i] ]])
}

上記のコードでは、次のことを行っている。

• iは1から3なのだが、より一般的に、\(\texttt{orig_vars}\)の長さ分を繰り返すという指示で、length()関数（長さはいくつ？）によって定義

• length(orig_vars)分の回数、\(\texttt{df_pol}\)内の\(i\)番目の\(\texttt{orig_vars}\)に対して、recode_var()を実施する

• 上記の結果を\(\texttt{df_pol}\)内に、\(i\)番目の\(\texttt{new_vars}\)に対応させて格納

また、length()だと長さをうまく認識しない場合があるので、その時には、seq_along()も利用可能である（多くは、length()でエラーなく対応可能だが為念）。

for(i in seq_along(orig_vars)) {
  df_cabinet[[ new_vars[i] ]] <- reverse_code(df_cabinet[[ orig_vars[i] ]])
}

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

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

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}} \]

df_pol %>%
with(cor.test(defense,preemptive,method="k"))

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

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

さらに、できることなら複数の変数間の順位相関係数を一覧にして確かめられるとよい。以下のような操作が可能である。

#------- 分析に必要な変数名についてのベクトルの作成
var_names <- c("defense", "preemptive", "yasukuni", "welfare")

#------- 各組み合わせごとに相関検定を行い、結果を結合
corr_results <- var_names %>% 
  combn(2, simplify = FALSE) %>% #すべてのペアを形成するという指示
  map_dfr(function(pair) { #map_dfr() で各組に対して cor.test() を実行し、その結果（相関係数と p 値）を tibble 形式でまとめるという流れ
    test <- cor.test(df_pol[[pair[1]]], df_pol[[pair[2]]], method = "kendall") #相関係数とその検定結果をtestに格納
    tibble(
      変数1 = pair[1],
      変数2 = pair[2],
      相関係数 = round(test$estimate, 3),　#testのうちのestimateオブジェクトを抜き出して小数点3位でまとめる
      p値  = round(test$p.value, 3)　#testのうちのp.valueオブジェクトを抜き出して小数点3位でまとめる
    )
  })

#------- 結果を表示
print(corr_results)

## # A tibble: 6 × 4
##   変数1      変数2      相関係数   p値
##   <chr>      <chr>         <dbl> <dbl>
## 1 defense    preemptive    0.493     0
## 2 defense    yasukuni      0.592     0
## 3 defense    welfare       0.495     0
## 4 preemptive yasukuni      0.456     0
## 5 preemptive welfare       0.297     0
## 6 yasukuni   welfare       0.376     0

上記の結果を一覧表にして示しても良いのだが、

#------- ペアプロット作成のためのパッケージ
library(GGally)

# df_pol から分析に使う変数のみを抽出して列名変更
data_sel <- df_pol %>%
  dplyr::select(defense, preemptive, yasukuni, welfare) %>%
  rename(
    "防衛力" = defense,
    "先制攻撃" = preemptive,
    "靖国参拝" = yasukuni,
    "小さな政府" = welfare
  )

#------- 下段三角パネルを作成するための自作関数（詳細は下記）
lower_fn <- function(data, mapping, ...){
  ggplot(data = data, mapping = mapping) +
    geom_point(alpha = 0.6) +
    stat_density2d(aes(fill = ..level..), geom = "polygon", alpha = 0.3) +
    scale_fill_viridis_c() +
    theme_minimal()
}

#------- 散布図をもとにした2変数間の関係、単変量の分布、相関係数値を一覧できるペアプロットを作成（詳細は下記）
pairs_plot <- ggpairs(
  data_sel,
  lower = list(continuous = lower_fn),
  diag  = list(continuous = "densityDiag"),
  upper = list(continuous = wrap("cor", size = 4))
)

pairs_plot

下側三角パネルの散布図と密度の描画に関する説明

stat_density2d(aes(fill = ..level..), geom = "polygon", alpha = 0.3) +
scale_fill_viridis_c() +

• stat_density2d(): 二次元の密度推定を行い、その結果をポリゴン（塗りつぶされた領域）として表示

• ..level..: stat_density2d() 内で計算される「密度レベル」を示す内部変数。..level.. の前後を二重のドット（..）で囲む記法は、\(\texttt{ggplot2}\) の統計変換（stat）の内部で計算された変数にアクセスするため。

• geom = "polygon": 密度の輪郭を描かず、密度に関する図形の内側が塗りつぶされる多角形としてプロッするという指示

• alpha = 0.3: 多角形の透明度を 0.3（30% の不透明度）に設定

• scale_fill_viridis_c(): 色覚配慮のカラースケール。密度レベルに応じたグラデーションを展開

上側三角パネルの散布図と密度の描画に関する説明：階段状のプロットなので、散布図行列の下側三角（\(\texttt{lower}\)）、対角線上(\(\texttt{diag}\))、上側三角(\(\texttt{upper}\))に何を配置するのかを指示したうえで作図している。具体的には…

pairs_plot <- ggpairs(
  data_sel,
  lower = list(continuous = lower_fn),
  diag  = list(continuous = "densityDiag"),
  upper = list(continuous = wrap("cor", size = 5))
)

• lower = list(continuous = lower_fn): 散布図行列の下側のパネルで表示するプロットを定義

  • continuous = lower_fn : 連続変数の組み合わせについてあらかじめ定義したカスタム関数lower_fn()を適用

• diag = list(continuous = "densityDiag"): 対角線上のパネルに表示するプロットを定義

  • continuous = "densityDiag": 対角パネルには各変数の密度プロット（≒density histogram）を描画。単変量の分布を明示

• upper = list(continuous = wrap("cor", size = 5)): 散布図行列の上側（三角領域）のパネルで表示するプロットを定義

  • continuous = wrap("cor", size = 3): 各変数の組み合わせについて相関係数がテキスト（数値）として、size=3で表示。wrap("cor", size = 3) は、パッケージ\(\texttt{GGally}\)に用意されている関数で、相関係数の表示を指示