第10章 個票データに親しむ
本章で使用する東大社研・若年パネル調査（JLPS-Y）wave1のオープンデータ（u001）は以下のURLで公開されていますので、取得してください。
https://csrda.iss.u-tokyo.ac.jp/infrastructure/urd/
分析結果の解釈については、書籍をご参照ください。以下はWindows10、R version 4.2.1で実行しています。
データを読み込むにあたって、作業場所の指定とそこにデータが置いてあることが前提になります。 R本体を操作している場合は、「ファイル」→「ディレクトリの変更」でデータの置いてある場所を指定するのが簡単です。 RStudioを操作している場合は、例えばデスクトップの「R」という名前のフォルダにデータがあるとすると
setwd("C:/Users/ユーザー名/Desktop/R")を最初に実行するのがよいでしょう。ユーザー名は各自異なるので注意です。
使用するデータを読み込み、表示します。なお、R version 4.1以前で読み込む際には
data1 <- read.csv("ファイル名.csv", fileEncoding = "UTF-8-BOM")のように、エンコードのオプションを指定する必要があります。version 4.2以降では必要ありませんので、以下ではオプションを指定していません。
data1 <- read.csv("u001.csv")本章で使用するパッケージを読み込みますが、パッケージはインストールしておく必要があります。 以下のようにコマンドでインストールするか、RやRStudioからクリックでインストールすることもできます。
install.packages(“stargazer”, dependencies = TRUE)パッケージの機能が使えるように読み込みます。
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.0 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(car)## 要求されたパッケージ carData をロード中です
##
## 次のパッケージを付け加えます: 'car'
##
## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## recode
##
## 以下のオブジェクトは 'package:purrr' からマスクされています:
##
## somelibrary(psych)##
## 次のパッケージを付け加えます: 'psych'
##
## 以下のオブジェクトは 'package:car' からマスクされています:
##
## logit
##
## 以下のオブジェクトは 'package:ggplot2' からマスクされています:
##
## %+%, alphalibrary(e1071)個人年収（ZQ47A）のデータを確認します。
table(data1$ZQ47A)##
## 1 2 3 4 5 6 7 8 9 10 11 12 14 99
## 113 61 60 133 163 176 118 76 16 7 1 1 42 33無回答（99）がいます。
配偶状態（ZQ50）と配偶者の生年-元号（ZQ52A）のクロス表を作成します。
table(data1$ZQ50,data1$ZQ52A)##
## 1 2 8 9
## 1 0 0 637 0
## 2 83 255 0 7
## 4 4 10 0 4未婚の回答者では配偶者の生年が非該当（8）になっています。
現在の働き方（JC_1）のデータを確認します。
table(data1$JC_1)##
## 1 2 3 4 5 6 7 8 10 11 12 99
## 7 489 150 31 5 18 24 6 126 48 92 4調査票の選択肢にはないデータ（10～12）があります。
最終学歴（ZQ23A）と最終学校の卒業（ZQ24）のクロス表を作成します。
table(data1$ZQ23A,data1$ZQ24)##
## 1 2 3 9
## 1 12 0 0 2
## 2 208 26 2 5
## 3 159 12 16 4
## 4 119 3 8 1
## 5 246 18 103 1
## 6 25 4 15 0
## 7 2 0 0 0
## 9 1 0 0 8それぞれの最終学歴に対して中退（2）、在学中（3）がいることがわかります。
表１０－３ 生活満足度（ZQ30D）の度数分布表
table(data1$ZQ30D)##
## 1 2 3 4 5 9
## 163 431 257 103 37 9無回答（9）がいます。
生活満足度をリコードします。
無回答をNAにリコードして、クロス表で確認します。
表10-5 クロス表でリコードの確認
data1$lsat=recode(data1$ZQ30D,'9=NA')
table(data1$ZQ30D,data1$lsat,useNA='ifany')##
## 1 2 3 4 5 <NA>
## 1 163 0 0 0 0 0
## 2 0 431 0 0 0 0
## 3 0 0 257 0 0 0
## 4 0 0 0 103 0 0
## 5 0 0 0 0 37 0
## 9 0 0 0 0 0 9個人年収をリコードします。
階級値をあてはめて、クロス表で確認します。
data1$resinc=recode(data1$ZQ47A, '1=0;2=12.5;3=50;4=100;5=200;6=300;7=400;8=500;9=700;10=1000;11=1500;12=2000;14:99=NA')
table(data1$ZQ47A, data1$resinc, useNA = "ifany")##
## 0 12.5 50 100 200 300 400 500 700 1000 1500 2000 <NA>
## 1 113 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 61 0 0 0 0 0 0 0 0 0 0 0
## 3 0 0 60 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 133 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 163 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 176 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 118 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 76 0 0 0 0 0
## 9 0 0 0 0 0 0 0 0 16 0 0 0 0
## 10 0 0 0 0 0 0 0 0 0 7 0 0 0
## 11 0 0 0 0 0 0 0 0 0 0 1 0 0
## 12 0 0 0 0 0 0 0 0 0 0 0 1 0
## 14 0 0 0 0 0 0 0 0 0 0 0 0 42
## 99 0 0 0 0 0 0 0 0 0 0 0 0 33精神的健康（ZA26A～ZQ26E）をリコードします。
リコードして、それぞれについてクロス表で確認します。ここではmutateを使っていますが、上で行ったように1つ1つリコードしてデータフレームに加えても問題ありません。
data1 <- data1 %>% mutate(
ZQ26Ar=recode(ZQ26A,'9=NA'),
ZQ26Br=recode(ZQ26B,'9=NA'),
ZQ26Cr=recode(ZQ26C,'1=5;2=4;3=3;4=2;5=1;9=NA'),
ZQ26Dr=recode(ZQ26D,'9=NA'),
ZQ26Er=recode(ZQ26E,'1=5;2=4;3=3;4=2;5=1;9=NA'))table(data1$ZQ26A, data1$ZQ26Ar, useNA="ifany")##
## 1 2 3 4 5 <NA>
## 1 30 0 0 0 0 0
## 2 0 114 0 0 0 0
## 3 0 0 333 0 0 0
## 4 0 0 0 264 0 0
## 5 0 0 0 0 247 0
## 9 0 0 0 0 0 12table(data1$ZQ26B, data1$ZQ26Br, useNA="ifany")##
## 1 2 3 4 5 <NA>
## 1 29 0 0 0 0 0
## 2 0 63 0 0 0 0
## 3 0 0 266 0 0 0
## 4 0 0 0 291 0 0
## 5 0 0 0 0 340 0
## 9 0 0 0 0 0 11table(data1$ZQ26C, data1$ZQ26Cr, useNA="ifany")##
## 1 2 3 4 5 <NA>
## 1 0 0 0 0 47 0
## 2 0 0 0 407 0 0
## 3 0 0 345 0 0 0
## 4 0 144 0 0 0 0
## 5 46 0 0 0 0 0
## 9 0 0 0 0 0 11table(data1$ZQ26D, data1$ZQ26Dr, useNA="ifany")##
## 1 2 3 4 5 <NA>
## 1 27 0 0 0 0 0
## 2 0 74 0 0 0 0
## 3 0 0 318 0 0 0
## 4 0 0 0 359 0 0
## 5 0 0 0 0 210 0
## 9 0 0 0 0 0 12table(data1$ZQ26E, data1$ZQ26Er, useNA="ifany")##
## 1 2 3 4 5 <NA>
## 1 0 0 0 0 78 0
## 2 0 0 0 321 0 0
## 3 0 0 453 0 0 0
## 4 0 106 0 0 0 0
## 5 30 0 0 0 0 0
## 9 0 0 0 0 0 125つの指標を合計して新しい変数（mhi_5）をつくるには、以下の①と②のどちらでもできます。
①基本的なコマンドを使う場合
data1$mhi_5 <- data1$ZQ26Ar+data1$ZQ26Br+data1$ZQ26Cr+data1$ZQ26Dr+data1$ZQ26Er②mutateを使う場合
data1 <- data1 %>% mutate(mhi_5=ZQ26Ar+ZQ26Br+ZQ26Cr+ZQ26Dr+ZQ26Er)5つの指標の合計点（mhi_5）を100点満点に換算して新しい変数（mhi5）を作成します。
data1 <- data1 %>% mutate(mhi5=(mhi_5-5)/20*100)記述統計を見てデータがうまくできていそうか確認します。
summary(data1$mhi_5)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 6.0 15.0 18.0 17.7 20.0 25.0 18summary(data1$mhi5)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 5.00 50.00 65.00 63.49 75.00 100.00 18MHI-5のクロンバックのα係数を計算します。
attach(data1)
data2 <- data.frame(ZQ26Ar,ZQ26Br,ZQ26Cr,ZQ26Dr,ZQ26Er)
alpha(data2)##
## Reliability analysis
## Call: alpha(x = data2)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.79 0.79 0.8 0.42 3.7 0.011 3.5 0.73 0.34
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.76 0.79 0.81
## Duhachek 0.77 0.79 0.81
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## ZQ26Ar 0.76 0.76 0.78 0.45 3.2 0.013 0.034 0.34
## ZQ26Br 0.70 0.71 0.69 0.38 2.4 0.016 0.022 0.32
## ZQ26Cr 0.77 0.76 0.74 0.44 3.2 0.012 0.033 0.39
## ZQ26Dr 0.73 0.73 0.72 0.40 2.7 0.014 0.026 0.34
## ZQ26Er 0.77 0.77 0.75 0.45 3.3 0.012 0.030 0.41
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## ZQ26Ar 988 0.72 0.70 0.58 0.52 3.6 1.08
## ZQ26Br 989 0.83 0.81 0.79 0.69 3.9 1.05
## ZQ26Cr 989 0.68 0.70 0.61 0.50 3.3 0.93
## ZQ26Dr 988 0.78 0.77 0.72 0.63 3.7 0.98
## ZQ26Er 988 0.66 0.69 0.59 0.49 3.3 0.88
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## ZQ26Ar 0.03 0.12 0.34 0.27 0.25 0.01
## ZQ26Br 0.03 0.06 0.27 0.29 0.34 0.01
## ZQ26Cr 0.05 0.15 0.35 0.41 0.05 0.01
## ZQ26Dr 0.03 0.07 0.32 0.36 0.21 0.01
## ZQ26Er 0.03 0.11 0.46 0.32 0.08 0.01図１０－３ 生活満足度のヒストグラム
ggplot(data1, aes(x=lsat))+geom_bar(fill="black", colour="white")+xlab("生活満足度")+ylab("回答者数（人）")+theme_classic()## Warning: Removed 9 rows containing non-finite values (stat_count).
表１０－６ 生活満足度の分布
table1 <- table(data1$lsat)
table1##
## 1 2 3 4 5
## 163 431 257 103 37prop.table(table1)##
## 1 2 3 4 5
## 0.16448032 0.43491423 0.25933401 0.10393542 0.03733602表１０－７ 性別と生活満足度のクロス表
table2 <- table(data1$lsat,data1$sex)
table2##
## 1 2
## 1 69 94
## 2 202 229
## 3 138 119
## 4 64 39
## 5 19 18prop.table(table2, margin=2)##
## 1 2
## 1 0.14024390 0.18837675
## 2 0.41056911 0.45891784
## 3 0.28048780 0.23847695
## 4 0.13008130 0.07815631
## 5 0.03861789 0.03607214独立性の検定を行います。
chisq.test(table2, correct=FALSE)##
## Pearson's Chi-squared test
##
## data: table2
## X-squared = 12.977, df = 4, p-value = 0.01139表１０－８ MHI-5の記述統計
summary(data1$mhi5)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 5.00 50.00 65.00 63.49 75.00 100.00 18var(data1$mhi5, na.rm=TRUE)## [1] 326.2686sd(data1$mhi5, na.rm=TRUE)## [1] 18.06291パッケージpsychのdescribe関数を利用する場合は以下のとおりです。fastオプションを使用しなければ次に計算する歪度も表示されます。
describe(data1$mhi5, fast=TRUE)## vars n mean sd min max range se
## X1 1 982 63.49 18.06 5 100 95 0.58歪度
skewness(data1$mhi5,na.rm=TRUE)## [1] -0.5237318図１０－４ 精神的健康度の分布
ggplot(data1, aes(x=mhi5))+geom_bar(fill="black", colour="white")+xlab("MHI-5")+ylab("回答者数（人）")+theme_classic()## Warning: Removed 18 rows containing non-finite values (stat_count).
精神的健康度の男女（sex）の平均値の差について検定します。
等分散性の検定を行います。
var.test(mhi5~sex, data=data1)##
## F test to compare two variances
##
## data: mhi5 by sex
## F = 0.82136, num df = 487, denom df = 493, p-value = 0.02966
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.6879528 0.9807405
## sample estimates:
## ratio of variances
## 0.8213632分散が等しくない場合の、平均値の差の検定を行います。
t.test(mhi5~sex, data=data1, var.equal=FALSE)##
## Welch Two Sample t-test
##
## data: mhi5 by sex
## t = 1.9312, df = 972.82, p-value = 0.05374
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -0.03586198 4.47984552
## sample estimates:
## mean in group 1 mean in group 2
## 64.61066 62.38866