第10章 練習問題
以下は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("ch10-12ex.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' からマスクされています:
##
## %+%, alpha１０－１
幸福度のデータにどのような値があるかを確認し、9は欠損値にリコードし、うまくリコードできているか確認します。
table(data1$q04_1)##
## 1 2 3 4 5 9
## 310 1289 577 352 187 32data1$hap=recode(data1$q04_1,'1=5;2=4;3=3;4=2;5=1;9=NA')
table(data1$q04_1, data1$hap, useNA = "ifany")##
## 1 2 3 4 5 <NA>
## 1 0 0 0 0 310 0
## 2 0 0 0 1289 0 0
## 3 0 0 577 0 0 0
## 4 0 352 0 0 0 0
## 5 187 0 0 0 0 0
## 9 0 0 0 0 0 32幸福度のヒストグラムを作成します。
ggplot(data1, aes(x=hap))+geom_bar(fill="black", colour="white")+xlab("幸福度")+ylab("回答者数")+theme_classic()## Warning: Removed 32 rows containing non-finite values (stat_count).
１０－２
男女別の幸福度の分布を見ます。
table1 <- table(data1$gender, data1$hap)
table1##
## 1 2 3 4 5
## 1 132 227 381 660 137
## 2 55 125 196 629 173prop.table(table1,1)##
## 1 2 3 4 5
## 1 0.08588159 0.14769031 0.24788549 0.42940794 0.08913468
## 2 0.04668930 0.10611205 0.16638370 0.53395586 0.14685908男女で幸福度の分布に差があるか検定します。
chisq.test(table1, correct=FALSE)##
## Pearson's Chi-squared test
##
## data: table1
## X-squared = 79.423, df = 4, p-value = 2.308e-16p値は0.01より小さく、1％水準で帰無仮説を棄却します。したがって、男女で幸福度の分布は異なることがわかります。
１０－３
価値観に関する変数について99を欠損値にリコードし、クロス表で確認します。
data1 <- data1 %>% mutate(
q651=recode(q06_5_1,'99=NA'),
q652=recode(q06_5_2,'99=NA'),
q653=recode(q06_5_3,'99=NA'),
q654=recode(q06_5_4,'99=NA'),
q655=recode(q06_5_5,'99=NA'),
q656=recode(q06_5_6,'99=NA'),
q657=recode(q06_5_7,'99=NA'))table(data1$q06_5_1, data1$q651, useNA = "ifany")##
## 1 2 3 4 5 6 7 8 9 10 <NA>
## 1 408 0 0 0 0 0 0 0 0 0 0
## 2 0 170 0 0 0 0 0 0 0 0 0
## 3 0 0 254 0 0 0 0 0 0 0 0
## 4 0 0 0 236 0 0 0 0 0 0 0
## 5 0 0 0 0 428 0 0 0 0 0 0
## 6 0 0 0 0 0 282 0 0 0 0 0
## 7 0 0 0 0 0 0 139 0 0 0 0
## 8 0 0 0 0 0 0 0 168 0 0 0
## 9 0 0 0 0 0 0 0 0 95 0 0
## 10 0 0 0 0 0 0 0 0 0 294 0
## 99 0 0 0 0 0 0 0 0 0 0 273table(data1$q06_5_2, data1$q652, useNA = "ifany")##
## 1 2 3 4 5 6 7 8 9 10 <NA>
## 1 107 0 0 0 0 0 0 0 0 0 0
## 2 0 33 0 0 0 0 0 0 0 0 0
## 3 0 0 85 0 0 0 0 0 0 0 0
## 4 0 0 0 127 0 0 0 0 0 0 0
## 5 0 0 0 0 291 0 0 0 0 0 0
## 6 0 0 0 0 0 197 0 0 0 0 0
## 7 0 0 0 0 0 0 163 0 0 0 0
## 8 0 0 0 0 0 0 0 285 0 0 0
## 9 0 0 0 0 0 0 0 0 226 0 0
## 10 0 0 0 0 0 0 0 0 0 1078 0
## 99 0 0 0 0 0 0 0 0 0 0 155table(data1$q06_5_3, data1$q653, useNA = "ifany")##
## 1 2 3 4 5 6 7 8 9 10 <NA>
## 1 189 0 0 0 0 0 0 0 0 0 0
## 2 0 113 0 0 0 0 0 0 0 0 0
## 3 0 0 211 0 0 0 0 0 0 0 0
## 4 0 0 0 228 0 0 0 0 0 0 0
## 5 0 0 0 0 537 0 0 0 0 0 0
## 6 0 0 0 0 0 389 0 0 0 0 0
## 7 0 0 0 0 0 0 240 0 0 0 0
## 8 0 0 0 0 0 0 0 285 0 0 0
## 9 0 0 0 0 0 0 0 0 117 0 0
## 10 0 0 0 0 0 0 0 0 0 237 0
## 99 0 0 0 0 0 0 0 0 0 0 201table(data1$q06_5_4, data1$q654, useNA = "ifany")##
## 1 2 3 4 5 6 7 8 9 10 <NA>
## 1 457 0 0 0 0 0 0 0 0 0 0
## 2 0 200 0 0 0 0 0 0 0 0 0
## 3 0 0 311 0 0 0 0 0 0 0 0
## 4 0 0 0 246 0 0 0 0 0 0 0
## 5 0 0 0 0 534 0 0 0 0 0 0
## 6 0 0 0 0 0 349 0 0 0 0 0
## 7 0 0 0 0 0 0 199 0 0 0 0
## 8 0 0 0 0 0 0 0 157 0 0 0
## 9 0 0 0 0 0 0 0 0 50 0 0
## 10 0 0 0 0 0 0 0 0 0 81 0
## 99 0 0 0 0 0 0 0 0 0 0 163table(data1$q06_5_5, data1$q655, useNA = "ifany")##
## 1 2 3 4 5 6 7 8 9 10 <NA>
## 1 600 0 0 0 0 0 0 0 0 0 0
## 2 0 344 0 0 0 0 0 0 0 0 0
## 3 0 0 426 0 0 0 0 0 0 0 0
## 4 0 0 0 312 0 0 0 0 0 0 0
## 5 0 0 0 0 409 0 0 0 0 0 0
## 6 0 0 0 0 0 193 0 0 0 0 0
## 7 0 0 0 0 0 0 67 0 0 0 0
## 8 0 0 0 0 0 0 0 83 0 0 0
## 9 0 0 0 0 0 0 0 0 47 0 0
## 10 0 0 0 0 0 0 0 0 0 85 0
## 99 0 0 0 0 0 0 0 0 0 0 181table(data1$q06_5_6, data1$q656, useNA = "ifany")##
## 1 2 3 4 5 6 7 8 9 10 <NA>
## 1 87 0 0 0 0 0 0 0 0 0 0
## 2 0 46 0 0 0 0 0 0 0 0 0
## 3 0 0 89 0 0 0 0 0 0 0 0
## 4 0 0 0 97 0 0 0 0 0 0 0
## 5 0 0 0 0 303 0 0 0 0 0 0
## 6 0 0 0 0 0 226 0 0 0 0 0
## 7 0 0 0 0 0 0 157 0 0 0 0
## 8 0 0 0 0 0 0 0 278 0 0 0
## 9 0 0 0 0 0 0 0 0 272 0 0
## 10 0 0 0 0 0 0 0 0 0 999 0
## 99 0 0 0 0 0 0 0 0 0 0 193table(data1$q06_5_7, data1$q657, useNA = "ifany")##
## 1 2 3 4 5 6 7 8 9 10 <NA>
## 1 270 0 0 0 0 0 0 0 0 0 0
## 2 0 133 0 0 0 0 0 0 0 0 0
## 3 0 0 249 0 0 0 0 0 0 0 0
## 4 0 0 0 262 0 0 0 0 0 0 0
## 5 0 0 0 0 487 0 0 0 0 0 0
## 6 0 0 0 0 0 316 0 0 0 0 0
## 7 0 0 0 0 0 0 226 0 0 0 0
## 8 0 0 0 0 0 0 0 248 0 0 0
## 9 0 0 0 0 0 0 0 0 122 0 0
## 10 0 0 0 0 0 0 0 0 0 257 0
## 99 0 0 0 0 0 0 0 0 0 0 177クロンバックのα係数で信頼性の確認をします。
attach(data1)
data2 <- data.frame(q651,q652,q653,q654,q655,q656,q657)
alpha(data2)## Number of categories should be increased in order to count frequencies.##
## Reliability analysis
## Call: alpha(x = data2)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.76 0.77 0.76 0.32 3.3 0.0069 5.7 1.8 0.32
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.75 0.76 0.78
## Duhachek 0.75 0.76 0.78
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## q651 0.74 0.74 0.73 0.32 2.8 0.0078 0.0112 0.33
## q652 0.75 0.75 0.73 0.34 3.1 0.0073 0.0107 0.33
## q653 0.72 0.72 0.71 0.30 2.5 0.0084 0.0143 0.25
## q654 0.72 0.72 0.70 0.30 2.5 0.0084 0.0093 0.26
## q655 0.75 0.76 0.75 0.34 3.1 0.0073 0.0134 0.35
## q656 0.74 0.75 0.73 0.33 2.9 0.0075 0.0133 0.33
## q657 0.72 0.73 0.71 0.31 2.7 0.0082 0.0106 0.32
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## q651 2474 0.66 0.64 0.56 0.48 5.0 2.9
## q652 2592 0.59 0.58 0.48 0.40 7.6 2.7
## q653 2546 0.71 0.71 0.65 0.57 5.6 2.5
## q654 2584 0.71 0.72 0.67 0.58 4.4 2.4
## q655 2566 0.56 0.57 0.44 0.38 3.7 2.4
## q656 2554 0.62 0.61 0.52 0.44 7.6 2.6
## q657 2570 0.69 0.68 0.62 0.53 5.4 2.7α係数は0.76となっており、これら価値観の変数を合計しても良いことがわかります。
data1$value=q651+q652+q653+q654+q655+q656+q657合計した指標でヒストグラムを作成します。
ggplot(data1, aes(x=value))+geom_bar(fill="black", colour="white")+xlab("価値観")+ylab("回答者数")+theme_classic()## Warning: Removed 479 rows containing non-finite values (stat_count).
１０－４
最初に等分散性の検定を行います。
var.test(value~gender, data=data1)##
## F test to compare two variances
##
## data: value by gender
## F = 1.2364, num df = 1311, denom df = 955, p-value = 0.0004648
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 1.098181 1.390440
## sample estimates:
## ratio of variances
## 1.236385p値は0.01より小さいので分散は等しいとは言えないことがわかります。
そこで、分散が等しくないことを前提とした平均値の差の検定を行います。
t.test(value~gender, data=data1, var.equal=FALSE)##
## Welch Two Sample t-test
##
## data: value by gender
## t = -1.8042, df = 2168.7, p-value = 0.07133
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -1.8328557 0.0763308
## sample estimates:
## mean in group 1 mean in group 2
## 38.95960 39.83787p値は0.076で0.1は下回るので、10％水準では有意と判断することができますが、0.05は下回っていないので、5％水準で有意と判断することはできません。 したがって、10％水準では、男女で価値観のスコアに差があると言えますが、5％水準では差があると言えません。 本文と同様に、どの程度の信頼性を基準とするかで、少し結論が変わってきます。