forked from yukiyanai/quant-methods-R
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathch07.R
180 lines (153 loc) · 6.73 KB
/
ch07.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
## ch07.R
##
## 浅野正彦・矢内勇生. 2018. 『Rによる計量政治学』オーム社
## 第7章 統計的推定
##
## Created: 2018-11-22 Yuki Yanai
## Modified: 2018-11-24 YY
## 2021-05-28 YY
## tidyverse パッケージを読み込む
library("tidyverse")
if (capabilities("aqua")) { # Macかどうか判定し、Macの場合のみ実行
theme_set(theme_gray(base_size = 10,
base_family = "HiraginoSans-W3"))
}
####################################################
## 7.1 母集団と標本
####################################################
## Rコードなし
####################################################
## 7.2 標本分布
####################################################
## 10万人の母集団の身長を設定する
## 身長の母平均は170cm、母標準偏差は6とする
set.seed(2018-11-20) # 乱数の種を設定する。
population <- rnorm(100000, mean = 170, sd = 6)
mean(population)
sd(population)
## 図7.3 母集団の分布
hist_pop <- data_frame(pop = population) %>%
ggplot(aes(x = pop, y = ..density..)) + # 古い書き方
#ggplot(aes(x = pop, y = after_stat(density))) + # 新しい書き方
geom_histogram(binwidth = 1, color = "black") +
xlim(140, 200) + ylim(0, 0.6) +
labs(x = "身長 (cm)", y = "確率密度")
print(hist_pop)
## 上で定義した母集団から、サイズ10の標本を500個抽出する
## 結果を保存するために、10行x500列の空の行列を用意する
size10 <- matrix(NA, nrow = 10, ncol = 500)
for (i in 1:500) { # 標本抽出を500回繰り返す
size10[, i] <- sample(population, size = 10, replace = FALSE)
}
## それぞれの標本での平均値を計算し、データフレームの1変数として保存する
df <- data_frame(means_s10 = colMeans(size10))
## 図7.4
## 標本平均の標本分布をヒストグラムにする
hist_size10 <- df %>%
ggplot(aes(x = means_s10, y = ..density.. )) + # 古い書き方
#ggplot(aes(x = means_s10, y = after_stat(density))) + # 新しい書き方
geom_histogram(binwidth = 1, color = "black") +
xlim(140, 200) + ylim(0, 0.6) +
labs(x = expression(paste("身長の標本平均 ", bar(x), " (cm)")),
y = "確率密度")
print(hist_size10)
## 標本サイズ10の場合の標本平均のばらつき（標準誤差）を計算する
sd(df$means_s10)
## 上で定義した母集団から、サイズ90の標本を500個抽出する
## 結果を保存するために、90行x500列の空の行列を用意する
size90 <- matrix(NA, nrow = 90, ncol = 500)
for (i in 1:500) { # 標本抽出を500回繰り返す
size90[, i] <- sample(population, size = 90, replace = FALSE)
}
## それぞれの標本での平均値を計算し、データフレームの1変数として保存する
df <- mutate(df, means_s90 = colMeans(size90))
## 図7.5
## 標本平均の標本分布をヒストグラムにする
hist_size90 <- df %>%
ggplot(aes(x = means_s90, y = ..density.. )) + # 古い書き方
#ggplot(aes(x = means_s90, y = after_stat(density))) + # 新しい書き方
geom_histogram(binwidth = 1, color = "black") +
xlim(140, 200) + ylim(0, 0.6) +
labs(x = expression(paste("身長の標本平均 ", bar(x), " (cm)")),
y = "確率密度")
print(hist_size90)
## 標本サイズ90の場合の標本平均のばらつき（標準誤差）を計算する
sd(df$means_s90)
## 標本サイズ10の場合と90の場合の標本のばらつきの比
sd(df$means_s10) / sd(df$means_s90)
####################################################
## 7.3 母平均の推定と信頼区間
####################################################
## 図7.6 t分布
df_t <- data_frame(x = seq(-3, 3, length.out = 100)) %>%
mutate(`1` = dt(x, df = 1),
`5` = dt(x, df = 5),
`99` = dt(x, df = 99)) %>%
gather(key = "df", value = "density", `1`:`99`)
t_dist <- ggplot(df_t, aes(x = x,
y = density,
color = df,
linetype = df)) +
geom_line() +
labs(x = "", y = "確率密度") +
scale_x_continuous(breaks = -3:3) +
scale_color_discrete(name = "自由度") +
scale_linetype_discrete(name = "自由度") +
guides(shape = guide_legend(reverse = TRUE))
print(t_dist)
## 自由度99のt分布で、データの95%が収まる範囲の上限値を求める
qt(p = 0.025, df = 99, lower.tail = FALSE)
# または、
qt(p = 0.975, df = 99, lower.tail = TRUE)
## 同様に、下限値を求める
qt(p = 0.025, df = 99, lower.tail = TRUE)
## 上で定義した母集団（身長の母平均が約170）から
## 標本サイズ100の標本を1,000個抽出し、
## それぞれの標本について信頼区間を求める
n <- 100
res <- matrix(NA, nrow = n, ncol = 1000)
for (i in 1:1000) {
res[, i] <- sample(population, size = n, replace = FALSE)
}
## それぞれの標本の平均と標準偏差を求めてデータフレームの変数にする。
## その後、標準誤差 se と95%信頼区間の下限 lb、上限ub を計算するn
df_ci <- data_frame( # 古い書き方
#tibble( # 新しい書き方
id = 1:1000,
means = colMeans(res),
sds = apply(res, MARGIN = 2, FUN = sd)) %>%
mutate(se = sds / sqrt(n),
lb = means + qt(p = 0.025, df = 99, lower.tail = TRUE) * se,
ub = means + qt(p = 0.025, df = 99, lower.tail = FALSE) * se)
## 結果の一部を見てみる
head(df_ci)
## 信頼区間の中に真の母平均が含まれているかどうか判定する変数を作る
(pop_mean <- mean(population)) # 真の母平均
df_ci <- mutate(df_ci,
success = (lb <= pop_mean & ub >= pop_mean))
## 母平均を区間内に捉えた95%信頼区間の割合を計算する
mean(df_ci$success)
## 1000個の標本から50個をランダムに選び、95%信頼区間を図示する。
set.seed(2018-11-22)
ci_plt <- sample_n(df_ci, size = 50) %>%
ggplot(aes(x = reorder(id, means), y = means,
ymin = lb, ymax = ub, color = success)) +
geom_hline(yintercept = pop_mean, color = "dodgerblue") +
geom_linerange() +
geom_point() +
scale_color_discrete(name = "母数を含む?") +
labs(x = "", y = "標本平均と95%信頼区間") +
coord_flip()
print(ci_plt)
## 例7.4：衆院データを使って信頼区間を求める
HR <- read_rds("data/hr-data.Rds")
## Rdsファイルの読み込みがうまくいかない場合は以下を実行
#download.file(url = "https://git.io/fxhQU",
# destfile = "data/hr-data.csv")
#HR <- read_csv("data/hr-data.csv")
## 候補者の年齢ageの50%信頼区間を求める
age_ci50 <- t.test(HR$age, conf.level = 0.5)
age_ci50$conf.int
## 候補者の年齢ageの95%信頼区間を求める
age_ci95 <- t.test(HR$age)
age_ci95$conf.int