library(dplyr)
library(psych)
library(tidyverse)
library(gt)
記述統計・作表・作図・クロス集計
広島大学大学院 教育における統計分析の応用
記述統計・作表・作図・クロス集計
広島大学大学院 教育における統計分析の応用
使うライブラリ
1 サンプルデータを取り込む
1.1 サンプルデータのダウンロード
- 学力検査を模したデータ
- https://lab.educpsychol.com/sample/hiroshima/ach.csv
- 質問紙調査を模したデータ
- https://lab.educpsychol.com/sample/hiroshima/ques.csv
取り込んだcsvファイルは，以後この講義で自身がRを用いた分析を行うファイルを置くディレクトリに格納する。
なお，これらのデータは以下のように作成した。
# 学力検査
kokuA.norm <-data.frame(as.integer(rnorm(1200, mean = 10, sd = 2)))
colnames(kokuA.norm) <- c("kokuA")
kokuA.norm$kokuA <- as.numeric(kokuA.norm$kokuA)
kokuA.norm.0.14 <- kokuA.norm %>%
dplyr::filter(kokuA >= 1 & kokuA <= 14)
kokuA.uni <- data.frame(
as.integer(runif(2000 - nrow(kokuA.norm.0.14), min = 0, max = 14))
)
colnames(kokuA.uni) <- c("kokuA")
kokuA.moto <- rbind(kokuA.norm.0.14, kokuA.uni)
kokuA.spl <- sample(kokuA.moto, length(kokuA.moto))
kokuB.norm <-data.frame(as.integer(rnorm(1200, mean = 7, sd = 2)))
colnames(kokuB.norm) <- c("kokuB")
kokuB.norm$kokuB <- as.numeric(kokuB.norm$kokuB)
kokuB.norm.0.9 <- kokuB.norm %>%
dplyr::filter(kokuB >= 1 & kokuB <= 9)
kokuB.uni <- data.frame(
as.integer(runif(2000 - nrow(kokuB.norm.0.9), min = 0, max = 9))
)
colnames(kokuB.uni) <- c("kokuB")
kokuB.moto <- rbind(kokuB.norm.0.9, kokuB.uni)
kokuB.spl <- sample(kokuB.moto, length(kokuB.moto))
sansA.norm <-data.frame(as.integer(rnorm(1200, mean = 12, sd = 2)))
colnames(sansA.norm) <- c("sansA")
sansA.norm$sansA <- as.numeric(sansA.norm$sansA)
sansA.norm.0.16 <- sansA.norm %>%
dplyr::filter(sansA >= 1 & sansA <= 16)
sansA.uni <- data.frame(
as.integer(runif(2000 - nrow(sansA.norm.0.16), min = 0, max = 16))
)
colnames(sansA.uni) <- c("sansA")
sansA.moto <- rbind(sansA.norm.0.16, sansA.uni)
sansA.spl <- sample(sansA.moto, length(sansA.moto))
sansB.norm <-data.frame(as.integer(rnorm(1200, mean = 9, sd = 2)))
colnames(sansB.norm) <- c("sansB")
sansB.norm$sansB <- as.numeric(sansB.norm$sansB)
sansB.norm.0.13 <- sansB.norm %>%
dplyr::filter(sansB >= 1 & sansB <= 13)
sansB.uni <- data.frame(
as.integer(runif(2000 - nrow(sansB.norm.0.13), min = 0, max = 13))
)
colnames(sansB.uni) <- c("sansB")
sansB.moto <- rbind(sansB.norm.0.13, sansB.uni)
sansB.spl <- sample(sansB.moto, length(sansB.moto))
id <- data.frame(c(1:2000)); colnames(id) <- c("id")
ach.spl <- cbind(id, kokuA.spl, kokuB.spl, sansA.spl, sansB.spl)
write.csv(ach.spl, "ach.csv", row.names = F)
# 質問紙調査
id <- data.frame(c(1:300))
schl <- data.frame(c(rep(c(10001:10015), each = 10),
rep(c(20001:20015), each = 10)))
age <- data.frame(as.integer(runif(300 , min = 22, max = 54)))
gender <- data.frame(as.integer(runif(300 , min = 1, max = 3)))
# 整数の場合1以上3未満として範囲を1から2にする
q01 <- data.frame(as.integer(runif(300 , min = 1, max = 4)))
q02 <- data.frame(as.integer(runif(300 , min = 1, max = 4)))
q03 <- data.frame(as.integer(runif(300 , min = 1, max = 4)))
q04 <- data.frame(as.integer(runif(300 , min = 1, max = 4)))
ques <- cbind(id, schl, age, gender, q01, q02, q03, q04)
colnames(ques) <- c("id", "schl", "age", "gender",
"q01", "q02", "q03", "q04")
write.csv(ques, "ques.csv", row.names = F)1.2 csvファイルを読み込む
- データ名，変数名は数字では始めないこと。
""で囲み，拡張子も付けること。
ach <- read.csv("ach.csv") #学力
ques <- read.csv("ques.csv") #質問紙- 日本語が入っているデータの場合，文字化けが生じることがある。その場合には文字コードを指定する。データに日本語は入れない方が無難
ach <- read.csv("ach.csv", encoding="Shift_JIS")- 読み込んだデータを確認する
head(ach) id kokuA kokuB sansA sansB
1 1 10 5 10 10
2 2 8 6 10 6
3 3 12 5 9 11
4 4 7 4 8 8
5 5 9 9 8 11
6 6 9 7 10 9head(ques) id schl age gender q01 q02 q03 q04
1 1 10001 41 2 3 2 2 1
2 2 10001 40 1 2 2 3 2
3 3 10001 42 2 1 2 2 2
4 4 10001 45 2 3 2 2 1
5 5 10001 27 2 2 3 1 3
6 6 10001 26 2 3 3 3 12 学力検査データの処理の演習
2.1 通過率をデータに加える
- 素点ではなく通過率を扱いたいことがよくある（と思われる）。そうしたければ，変数を計算して追加する。
- たとえば，以下のように
achというデータのkokuA_propに，achというデータのkokuAを14（テスト冊子の全項目数）で割ったものを代入する。
ach$kokuA_prop <- ach$kokuA / 14
ach$kokuB_prop <- ach$kokuB / 9
ach$sansA_prop <- ach$sansA / 16
ach$sansB_prop <- ach$sansB / 13- 以下のように列が追加されている
head(ach) id kokuA kokuB sansA sansB kokuA_prop kokuB_prop sansA_prop sansB_prop
1 1 10 5 10 10 0.7142857 0.5555556 0.6250 0.7692308
2 2 8 6 10 6 0.5714286 0.6666667 0.6250 0.4615385
3 3 12 5 9 11 0.8571429 0.5555556 0.5625 0.8461538
4 4 7 4 8 8 0.5000000 0.4444444 0.5000 0.6153846
5 5 9 9 8 11 0.6428571 1.0000000 0.5000 0.8461538
6 6 9 7 10 9 0.6428571 0.7777778 0.6250 0.69230772.2 とりあえず記述統計をしたい
データ全体に対する要約(summary)統計量は以下のようにして求められる
summary(ach) id kokuA kokuB sansA
Min. : 1.0 Min. : 0.000 Min. :0.000 Min. : 0.000
1st Qu.: 500.8 1st Qu.: 7.000 1st Qu.:4.000 1st Qu.: 8.000
Median :1000.5 Median : 9.000 Median :6.000 Median :11.000
Mean :1000.5 Mean : 8.259 Mean :5.208 Mean : 9.851
3rd Qu.:1500.2 3rd Qu.:11.000 3rd Qu.:7.000 3rd Qu.:13.000
Max. :2000.0 Max. :14.000 Max. :9.000 Max. :16.000
sansB kokuA_prop kokuB_prop sansA_prop
Min. : 0.000 Min. :0.0000 Min. :0.0000 Min. :0.0000
1st Qu.: 6.000 1st Qu.:0.5000 1st Qu.:0.4444 1st Qu.:0.5000
Median : 8.000 Median :0.6429 Median :0.6667 Median :0.6875
Mean : 7.497 Mean :0.5899 Mean :0.5787 Mean :0.6157
3rd Qu.:10.000 3rd Qu.:0.7857 3rd Qu.:0.7778 3rd Qu.:0.8125
Max. :13.000 Max. :1.0000 Max. :1.0000 Max. :1.0000
sansB_prop
Min. :0.0000
1st Qu.:0.4615
Median :0.6154
Mean :0.5767
3rd Qu.:0.7692
Max. :1.0000 - 標準偏差などの記述(descirptive statistics)を知りたい場合には，psychパッケージを使う。
- ない場合には
install.packaes("psych")としてインストールする。
describe(ach) vars n mean sd median trimmed mad min max range skew
id 1 2000 1000.50 577.49 1000.50 1000.50 741.30 1 2000 1999 0.00
kokuA 2 2000 8.26 3.35 9.00 8.58 2.97 0 14 14 -0.78
kokuB 3 2000 5.21 2.46 6.00 5.39 2.97 0 9 9 -0.55
sansA 4 2000 9.85 3.85 11.00 10.30 2.97 0 16 16 -0.98
sansB 5 2000 7.50 3.04 8.00 7.77 2.97 0 13 13 -0.71
kokuA_prop 6 2000 0.59 0.24 0.64 0.61 0.21 0 1 1 -0.78
kokuB_prop 7 2000 0.58 0.27 0.67 0.60 0.33 0 1 1 -0.55
sansA_prop 8 2000 0.62 0.24 0.69 0.64 0.19 0 1 1 -0.98
sansB_prop 9 2000 0.58 0.23 0.62 0.60 0.23 0 1 1 -0.71
kurtosis se
id -1.20 12.91
kokuA -0.03 0.07
kokuB -0.59 0.05
sansA 0.25 0.09
sansB -0.14 0.07
kokuA_prop -0.03 0.01
kokuB_prop -0.59 0.01
sansA_prop 0.25 0.01
sansB_prop -0.14 0.012.3 論文に載せられる表を出力する
- 上のように表示されたデータを，いわゆるcut and pasteすれば表を作ることは出来る。
- しかし，cut and paste は誤りが発生しやすいので，機械に作業させるのが望ましい。
- まず，表の元となるデータを作る。
exam <- c("国語A", "国語B", "算数A", "算数B") #変数名
m <- c(mean(ach$kokuA_prop), mean(ach$kokuB_prop),
mean(ach$sansA_prop), mean(ach$sansB_prop))
sd <- c(sd(ach$kokuA_prop), sd(ach$kokuB_prop),
sd(ach$sansA_prop), sd(ach$sansB_prop))
ach.table <- data.frame(cbind(exam, m, sd))
ach.table exam m sd
1 国語A 0.589928571428571 0.23924336793883
2 国語B 0.578666666666667 0.273043101958322
3 算数A 0.61565625 0.240446768610445
4 算数B 0.576692307692308 0.234200463024955- 表をきれいに作るパッケージ使って，とりあえず出力させてみる。
- 使うパッケージは
library(tidyverse),library(gt) - コードチャンクを以下のようにすると，1行おきの網掛けがかからない
::: {.cell .plain}
:::
::: {.cell .plain}
```{.r .cell-code}
ach.table$m <- as.numeric(ach.table$m)
ach.table$sd <- as.numeric(ach.table$sd)
ach.table %>% gt()| exam | m | sd |
|---|---|---|
| 国語A | 0.5899286 | 0.2392434 |
| 国語B | 0.5786667 | 0.2730431 |
| 算数A | 0.6156563 | 0.2404468 |
| 算数B | 0.5766923 | 0.2342005 |
:::
- わりときれいな表になっている。
- これを，APA styleに合わせてみる。まず，有効数字を揃える。
tab.1 <- ach.table %>% gt()
tab.2 <- tab.1 %>%
fmt_number(
columns = c("m", "sd"),
decimals = 2) - 列名を振り直して，統計量を示すアルファベットを斜体にする
tab.3 <- tab.2 %>%
cols_label(
exam = md("教科"), # md は markdown記法
m = md("*M*"),
sd = md("*SD*")
)
tab.3| 教科 | M | SD |
|---|---|---|
| 国語A | 0.59 | 0.24 |
| 国語B | 0.58 | 0.27 |
| 算数A | 0.62 | 0.24 |
| 算数B | 0.58 | 0.23 |
- さらにAPAスタイルに合わせる
tab.4 <- tab.3 %>%
tab_options(table.width = pct(25), # 表全体の幅をやや拡げる
table_body.hlines.width = 0, # tableの中の水平線は引かない
column_labels.border.top.width = 2,
column_labels.border.top.color = "black",
column_labels.border.bottom.width = 2,
column_labels.border.bottom.color = "black",
# 最上列上下の水平線は細く黒く
table_body.border.bottom.color = "black" #最下線は黒
) %>%
cols_align(align = "center", columns = c("m", "sd"))
# 2列目移行は中央寄せ(APA 7thから)
tab.4 | 教科 | M | SD |
|---|---|---|
| 国語A | 0.59 | 0.24 |
| 国語B | 0.58 | 0.27 |
| 算数A | 0.62 | 0.24 |
| 算数B | 0.58 | 0.23 |
- もっと工夫すると，次のようなことができる。
# 4教科のテストについて
exam <- c("国語A", "国語B", "算数A", "算数B")
# 素点の平均と標準偏差を求め
m <- c(mean(ach$kokuA), mean(ach$kokuB),
mean(ach$sansA), mean(ach$sansB))
sd <- c(sd(ach$kokuA), sd(ach$kokuB),
sd(ach$sansA), sd(ach$sansB))
# 通過率の平均と標準偏差を求め
m_prop <- c(mean(ach$kokuA_prop), mean(ach$kokuB_prop),
mean(ach$sansA_prop), mean(ach$sansB_prop))
sd_prop <- c(sd(ach$kokuA_prop), sd(ach$kokuB_prop),
sd(ach$sansA_prop), sd(ach$sansB_prop))
# これらをまとめた表にして
ach.tab <- data.frame(cbind(exam, m, sd, m_prop, sd_prop))
# MとSDを数値型にして
ach.tab$m <- as.numeric(ach.tab$m)
ach.tab$sd <- as.numeric(ach.tab$sd)
ach.tab$m_prop <- as.numeric(ach.tab$m_prop)
ach.tab$sd_prop <- as.numeric(ach.tab$sd_prop)
# APAスタイルにあわせた出力を得て
gt.tab <- ach.tab %>%
gt() %>%# gtパッケージで表
fmt_number(columns = c("m", "sd", "m_prop", "sd_prop"), decimals = 2) %>%
# 有効数字揃え
cols_label(exam = md("教科"), m = md("*M*"), sd = md("*SD*"),
m_prop = md("*M*"), sd_prop = md("*SD*")) %>%
# 列名を振り直してMとSDは斜体
tab_spanner(label = "素点", columns = c("m", "sd")) %>%
tab_spanner(label = "通過率", columns = c("m_prop", "sd_prop")) %>%
tab_options(table.width = pct(50), # 表全体の幅をやや拡げる
table_body.hlines.width = 0, # tableの中の水平線は引かない
# spannerの下線を適切にするための設定
column_labels.border.top.width = 2,
column_labels.border.top.color = "black",
column_labels.border.bottom.width = 1,
column_labels.border.bottom.color = "black",
table.border.top.width = 2,
table_body.border.top.color = "black", #最上線は黒
table.border.bottom.width = 2,
table_body.border.bottom.color = "black" #最下線は黒
) %>%
# 2列目以降は中央寄せ(APA 7thから)
cols_align(align = "center",
columns = c("m", "sd", "m_prop", "sd_prop"))
gt.tab| 教科 |
素点
|
通過率
|
||
|---|---|---|---|---|
| M | SD | M | SD | |
| 国語A | 8.26 | 3.35 | 0.59 | 0.24 |
| 国語B | 5.21 | 2.46 | 0.58 | 0.27 |
| 算数A | 9.85 | 3.85 | 0.62 | 0.24 |
| 算数B | 7.50 | 3.04 | 0.58 | 0.23 |
- タイトルも付けてみる
gt.title <- gt.tab %>%
tab_header(title = md("**Table 1**<br>教科別の平均と標準偏差")) %>%
tab_style(
style=cell_text(align="left"),
locations = cells_title("title")) %>%
tab_options(
heading.title.font.size = px(16),
table.border.top.width = 0, #タイトルの上の線はなし
column_labels.border.top.width = 3, # タイトル下の線はあり 調整のため3
column_labels.border.top.color = "black"
)
gt.title| Table 1 教科別の平均と標準偏差 |
||||
| 教科 |
素点
|
通過率
|
||
|---|---|---|---|---|
| M | SD | M | SD | |
| 国語A | 8.26 | 3.35 | 0.59 | 0.24 |
| 国語B | 5.21 | 2.46 | 0.58 | 0.27 |
| 算数A | 9.85 | 3.85 | 0.62 | 0.24 |
| 算数B | 7.50 | 3.04 | 0.58 | 0.23 |
- 出力を画像として保存する
- あとから微修正ができなくなるが，信頼を得るために必要なことである。
- コードの下に出力されるのは生成されたpngファイルである。任意のディレクトリに保存されるので，あとは貼り付けるだけでよい。
gtsave(gt.title, "表の例_01.png")
gtsave(gt.title, "表の例_01.pdf") # pdf の方が TeX で扱いやすい- 作ったpngファイルをここで出力させると以下の通りである。
2.4 ヒストグラムを作る
- 今一度，学力検査の疑似データの内容を見てみる。
head(ach) id kokuA kokuB sansA sansB kokuA_prop kokuB_prop sansA_prop sansB_prop
1 1 10 5 10 10 0.7142857 0.5555556 0.6250 0.7692308
2 2 8 6 10 6 0.5714286 0.6666667 0.6250 0.4615385
3 3 12 5 9 11 0.8571429 0.5555556 0.5625 0.8461538
4 4 7 4 8 8 0.5000000 0.4444444 0.5000 0.6153846
5 5 9 9 8 11 0.6428571 1.0000000 0.5000 0.8461538
6 6 9 7 10 9 0.6428571 0.7777778 0.6250 0.6923077- 例えば，国語Aの正答数について，単純にヒストグラムを出力することは，たいしたことはない。
hist(ach$kokuA)
- データの様子を眺める程度なら，これで問題ない。
- しかし，国語Aの場合の得点の範囲は0から14（14段階）なのに，棒の数は13しかない。
hist(ach$kokuA, breaks=20)
- これで，それなりの出力は得られるが，0点の人が1点の人と隣り合わせになってしまう。
- 仕方ないので，棒グラフを愚直に描く。
2.4.1 度数分布表を作ってbarplotにする
まずは，国語Aについて，度数分布表を作ってみる。
# 度数分布表をfreq.kokuAというデータにする
freq.kokuA <- table(ach$kokuA)
freq.kokuA
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
61 60 57 60 75 74 96 146 254 295 283 237 176 104 22 単純にbarplotを作ると，以下のようになる。
barplot(freq.kokuA)
このままだと，あまり見た目がよくないので，調整する。
barplot(freq.kokuA,
xlab = "正答数", ylab = "人数", # 縦，横軸に名前を付ける
ylim = c(0, 300), # 縦軸を正答数ごとの人数の最大値よりも大きく
cex.names = 0.8, # 軸ラベルをやや小さくする
cex.axis = 0.8, # 縦軸のラベル（数字）をやや小さめにする
axis.lty = 1, # 横軸に線を引く
space = c(0, 0) # 棒の間を詰めてヒストグラム風にする
)Warning in space + width: longer object length is not a multiple of shorter
object length
Warning in space + width: longer object length is not a multiple of shorter
object length
- 本来，barplotは，barの間が空いているのが普通なので，本当に大丈夫かというエラーメッセージが出るが気にしない。
3 質問紙調査データの例
データはこのような形
head(ques) id schl age gender q01 q02 q03 q04
1 1 10001 41 2 3 2 2 1
2 2 10001 40 1 2 2 3 2
3 3 10001 42 2 1 2 2 2
4 4 10001 45 2 3 2 2 1
5 5 10001 27 2 2 3 1 3
6 6 10001 26 2 3 3 3 13.1 回答者の属性を見てみる
- 性別
table(ques$gender)
1 2
149 151 - 年齢
table(ques$age)
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
10 13 13 10 13 14 6 4 7 9 6 9 10 7 11 12 9 6 9 9 10 9 9 14 11 5
48 49 50 51 52 53
4 14 5 9 15 8 - 年齢の平均や標準偏差
describe(ques$age) vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 300 37.18 9.56 37 37.11 12.6 22 53 31 0.01 -1.27 0.553.2 グラフにする
- 割合を積み上げ棒グラフで表現するとよい。
barplot(matrix(prop.table(table(ques$gender))), horiz =TRUE)
- これだと見栄えが悪いので，工夫する。
#以降改行までは処理されない。
# 凡例が入るようにする（上下の余白を詰める）
par(mar = c(3, 4, 1, 12))
# グラフ用に行列を渡す
barplot(matrix(prop.table(table(ques$gender))),
col = c("gray1", "gray50"),
horiz = TRUE)
## 凡例をプロットエリア外に描けるようにする
par(xpd = TRUE)
# 凡例
legend(x = par()$usr[2], y = par()$usr[4], # 凡例がグラフエリア右上
legend = c("女性", "男性"),
xjust = -.05, # プロットエリアから少し離して凡例を載せる
fill = c("gray1", "gray50") # 色をグラフと揃える
)
- 画像ファイルを出力する際には，以下のようにすると調整したものが保存される。
#画像ファイルで高さを低くして保存する
png(filename = "g01.png", width = 7, height = 1.5, units = "in", res = 150)
# 凡例が入るようにする（上下の余白を詰める）
par(mar = c(3, 4, 1, 12))
# グラフ用に行列を渡す
barplot(matrix(prop.table(table(ques$gender))),
col = c("gray1", "gray50"),
horiz = TRUE)
## 凡例をプロットエリア外に描けるようにする
par(xpd = TRUE)
# 凡例
legend(x = par()$usr[2], y = par()$usr[4], # 凡例がグラフエリア右上
legend = c("女性", "男性"),
xjust = -.05, # プロットエリアから少し離して凡例を載せる
fill = c("gray1", "gray50") # 色をグラフと揃える
)
dev.off()quartz_off_screen
2 - 作ったpngファイルをここで出力させると以下の通りである。
3.2.1 別の層分けをしてグラフを作る
- 回答の年齢層は，層別するとよさそうである。
- 層別にする方法は，学力検査データの例と同じである。
# 年齢カテゴリの変数を追加する
cat.age <- ques %>%
dplyr::mutate(
cat = case_when(
age < 30 ~ 1,
age < 40 ~ 2,
age < 50 ~ 3,
TRUE ~ 4)
)
# 年齢カテゴリごとの人数の和をデータ化する
cat.age.sum <- data.frame(table(cat.age$cat))
cat.age.sum Var1 Freq
1 1 83
2 2 86
3 3 94
4 4 37- これをグラフにする
# 上記のデータを因子型にする
cat.age.sum$fac <- factor(cat.age.sum$Var1,
levels = c(1,2,3,4),
labels = c("20代","30代","40代","50代以上"),
ordered=TRUE)
# 使う色を指定しておく
color=c("gray1","gray25","gray49","gray73")
# 凡例も
par(mar=c(4,4,4,12)) # 凡例が入るように
barplot(matrix(prop.table(cat.age.sum$Freq)), # グラフ用に行列を渡す
col = color,
horiz = TRUE)
##凡例をプロットエリア外に描けるようにする
par(xpd = TRUE)
# 凡例
legend(x = par()$usr[2], y = par()$usr[4], # 凡例がグラフエリア右下
legend = cat.age.sum$fac,
xjust = -.05, #プロットエリアから少し離して凡例を載せる
fill = color # 色をグラフと揃える
)
3.2.2 項目の集計
psychパッケージのdescribe関数で分かる。
describe(ques) vars n mean sd median trimmed mad min max range
id 1 300 150.50 86.75 150.5 150.50 111.19 1 300 299
schl 2 300 15008.00 5008.36 15008.0 15008.00 7413.00 10001 20015 10014
age 3 300 37.18 9.56 37.0 37.11 12.60 22 53 31
gender 4 300 1.50 0.50 2.0 1.50 0.00 1 2 1
q01 5 300 1.97 0.84 2.0 1.96 1.48 1 3 2
q02 6 300 1.99 0.81 2.0 1.99 1.48 1 3 2
q03 7 300 2.08 0.83 2.0 2.10 1.48 1 3 2
q04 8 300 2.13 0.80 2.0 2.16 1.48 1 3 2
skew kurtosis se
id 0.00 -1.21 5.01
schl 0.00 -2.01 289.16
age 0.01 -1.27 0.55
gender -0.01 -2.01 0.03
q01 0.06 -1.58 0.05
q02 0.01 -1.48 0.05
q03 -0.16 -1.53 0.05
q04 -0.24 -1.41 0.053.3 クロス集計
- 例えば，年齢とq01とのクロス集計をしてみる。
# 年齢して層別化してデータに追加する
ques2 <- ques # 元のデータを保持しておく
ques2 <- ques2 %>%
dplyr::mutate(
cat = case_when(
age < 30 ~ 1,
age < 40 ~ 2,
age < 50 ~ 3,
TRUE ~ 4)
)
# カテゴリを因子型にする
## 年齢
ques2$cat <- factor(ques2$cat,
levels = c(1,2,3,4),
labels = c("20代","30代","40代","50代以上"),
ordered=TRUE)
## q01
ques2$q01 <- factor(ques2$q01,
levels = c(1,2,3),
labels = c("全く","時々","頻繁"),
ordered=TRUE)
head(ques2) id schl age gender q01 q02 q03 q04 cat
1 1 10001 41 2 頻繁 2 2 1 40代
2 2 10001 40 1 時々 2 3 2 40代
3 3 10001 42 2 全く 2 2 2 40代
4 4 10001 45 2 頻繁 2 2 1 40代
5 5 10001 27 2 時々 3 1 3 20代
6 6 10001 26 2 頻繁 3 3 1 20代データが変わっていることが分かる。
クロス集計自体は，簡単である。
tab01 <- table(ques2$cat, ques2$q01)
tab01
全く 時々 頻繁
20代 32 21 30
30代 32 35 19
40代 32 24 38
50代以上 14 10 13- 用途によって様々な形式がある。
# 行と列の合計欄を付ける
tab02 <- addmargins(tab01)
tab02
全く 時々 頻繁 Sum
20代 32 21 30 83
30代 32 35 19 86
40代 32 24 38 94
50代以上 14 10 13 37
Sum 110 90 100 300# 行方向の割合
tab03 <- addmargins(prop.table(tab01, margin = 1))
tab03
全く 時々 頻繁 Sum
20代 0.3855422 0.2530120 0.3614458 1.0000000
30代 0.3720930 0.4069767 0.2209302 1.0000000
40代 0.3404255 0.2553191 0.4042553 1.0000000
50代以上 0.3783784 0.2702703 0.3513514 1.0000000
Sum 1.4764391 1.1855782 1.3379827 4.0000000# 有効数字を少なくする
round(tab03, 2)
全く 時々 頻繁 Sum
20代 0.39 0.25 0.36 1.00
30代 0.37 0.41 0.22 1.00
40代 0.34 0.26 0.40 1.00
50代以上 0.38 0.27 0.35 1.00
Sum 1.48 1.19 1.34 4.00# 列方向の割合
tab04 <- addmargins(prop.table(tab01, margin = 2))
round(tab04,2) # 和が100にならない場合がある
全く 時々 頻繁 Sum
20代 0.29 0.23 0.30 0.82
30代 0.29 0.39 0.19 0.87
40代 0.29 0.27 0.38 0.94
50代以上 0.13 0.11 0.13 0.37
Sum 1.00 1.00 1.00 3.00- ここまでのデータが揃うと，以下のような表を作ることができる。
# 人数の行列を分割する
tab.n.20 <- tab02[1, -4] # 1行目 合計列なし
tab.n.30 <- tab02[2, -4] # 2行目 合計列なし
tab.n.40 <- tab02[3, -4] # 3行目 合計列なし
tab.n.50 <- tab02[4, -4] # 4行目 合計列なし
# 行方向の割合を分割
tab.r.20 <- tab03[1, -4] # 1行目 合計列なし
tab.r.30 <- tab03[2, -4] # 2行目 合計列なし
tab.r.40 <- tab03[3, -4] # 3行目 合計列なし
tab.r.50 <- tab03[4, -4] # 4行目 合計列なし
# 列方向の割合を分割
tab.c.20 <- tab04[1, -4] # 1行目 合計列なし
tab.c.30 <- tab04[2, -4] # 2行目 合計列なし
tab.c.40 <- tab04[3, -4] # 3行目 合計列なし
tab.c.50 <- tab04[4, -4] # 4行目 合計列なし
# 分割したデータを積み上げ直す
tab.40 <- rbind(tab.n.20, tab.r.20, tab.c.20,
tab.n.30, tab.r.30, tab.c.30,
tab.n.40, tab.r.40, tab.c.40,
tab.n.50, tab.r.50, tab.c.50)
tab.40 全く 時々 頻繁
tab.n.20 32.0000000 21.0000000 30.0000000
tab.r.20 0.3855422 0.2530120 0.3614458
tab.c.20 0.2909091 0.2333333 0.3000000
tab.n.30 32.0000000 35.0000000 19.0000000
tab.r.30 0.3720930 0.4069767 0.2209302
tab.c.30 0.2909091 0.3888889 0.1900000
tab.n.40 32.0000000 24.0000000 38.0000000
tab.r.40 0.3404255 0.2553191 0.4042553
tab.c.40 0.2909091 0.2666667 0.3800000
tab.n.50 14.0000000 10.0000000 13.0000000
tab.r.50 0.3783784 0.2702703 0.3513514
tab.c.50 0.1272727 0.1111111 0.1300000- これを，報告できるような表に変換する
- 以下では
tidyverse,gtパッケージを使う。必要に応じて読み込むこと。
tab.40.df <- data.frame(tab.40) # データの形を変換する
# 年齢のグループを付ける
tab.40.df$agegroup <- c(1,1,1,2,2,2,3,3,3,4,4,4)
tab.40.df$agegroup <- factor(tab.40.df$agegroup,
levels = c(1,2,3,4),
labels = c("20代","30代","40代","50代以上"),
ordered = TRUE)
# 出力の内容の名前を付ける
tab.40.df$out <- c(1,2,3,1,2,3,1,2,3,1,2,3)
tab.40.df$out <- factor(tab.40.df$out,
levels = c(1,2,3),
labels = c("度数","行の割合","列の割合"),
ordered = FALSE)
# 列を並び替える
tab.40.df.2 <- tab.40.df[c(4,5,1:3)]
tab.40.df.2 agegroup out 全く 時々 頻繁
tab.n.20 20代 度数 32.0000000 21.0000000 30.0000000
tab.r.20 20代 行の割合 0.3855422 0.2530120 0.3614458
tab.c.20 20代 列の割合 0.2909091 0.2333333 0.3000000
tab.n.30 30代 度数 32.0000000 35.0000000 19.0000000
tab.r.30 30代 行の割合 0.3720930 0.4069767 0.2209302
tab.c.30 30代 列の割合 0.2909091 0.3888889 0.1900000
tab.n.40 40代 度数 32.0000000 24.0000000 38.0000000
tab.r.40 40代 行の割合 0.3404255 0.2553191 0.4042553
tab.c.40 40代 列の割合 0.2909091 0.2666667 0.3800000
tab.n.50 50代以上 度数 14.0000000 10.0000000 13.0000000
tab.r.50 50代以上 行の割合 0.3783784 0.2702703 0.3513514
tab.c.50 50代以上 列の割合 0.1272727 0.1111111 0.1300000cross.tab <- tab.40.df.2 %>% gt()
cross.tab| agegroup | out | 全く | 時々 | 頻繁 |
|---|---|---|---|---|
| 20代 | 度数 | 32.0000000 | 21.0000000 | 30.0000000 |
| 20代 | 行の割合 | 0.3855422 | 0.2530120 | 0.3614458 |
| 20代 | 列の割合 | 0.2909091 | 0.2333333 | 0.3000000 |
| 30代 | 度数 | 32.0000000 | 35.0000000 | 19.0000000 |
| 30代 | 行の割合 | 0.3720930 | 0.4069767 | 0.2209302 |
| 30代 | 列の割合 | 0.2909091 | 0.3888889 | 0.1900000 |
| 40代 | 度数 | 32.0000000 | 24.0000000 | 38.0000000 |
| 40代 | 行の割合 | 0.3404255 | 0.2553191 | 0.4042553 |
| 40代 | 列の割合 | 0.2909091 | 0.2666667 | 0.3800000 |
| 50代以上 | 度数 | 14.0000000 | 10.0000000 | 13.0000000 |
| 50代以上 | 行の割合 | 0.3783784 | 0.2702703 | 0.3513514 |
| 50代以上 | 列の割合 | 0.1272727 | 0.1111111 | 0.1300000 |
- さらに，加工する。
cross.tab <- tab.40.df.2 %>%
gt(groupname_col = "agegroup") %>%
fmt_number(
columns = c(3:5),
rows = c(1, 4, 7, 10), # 1,4,7,10行目は実数
decimals = 0) %>%
fmt_number(
columns = c(3:5),
rows = c(2:3, 5:6, 8:9, 11:12), # 1,4,7,10行目以外は割合
decimals = 2) %>%
cols_label(
agegroup = md("年齢"),
out = ("")) %>%
tab_spanner(label = "頻度", columns = c(3:5)) %>%
cols_align(align = "left", columns = c(2)) %>%
tab_options(table_body.hlines.width = 0, # tableの中の水平線は引かない
# table.border.top.width = 0, #タイトルの上の線を消す，
# table.border.bottom.width = 0, #注の下の線を消す
heading.title.font.size = px(16), #タイトルのフォントサイズ
# row_group.border.bottom.width = 0, #セクション名の下の線を消す
# stub.border.width = 0,
column_labels.border.top.width = 3, # 変数名の行の線
column_labels.border.top.color = "black",
column_labels.border.bottom.width = 2, # 変数名の行の線
column_labels.border.bottom.color = "black",
table_body.border.bottom.color = "black", #テーブルの下線を黒く
row_group.border.top.color = "black", #セクション名の上の線を黒
row_group.border.top.width = 1, #セクション名の上の線を細く
row_group.border.bottom.color = "white", #セクション名の上の線を消す
row_group.border.bottom.width = 1, #セクション名の上の線を細く
table.border.top.width = 0, #タイトルの上の線はなし
table.width = pct(30)) %>%
tab_header(title = md("**Table 2**<br>年齢となんとかの頻度の関係")) %>%
tab_style(
style=cell_text(align="left"),
locations = cells_title("title")
)
cross.tab| Table 2 年齢となんとかの頻度の関係 |
|||
頻度
|
|||
|---|---|---|---|
| 全く | 時々 | 頻繁 | |
| 20代 | |||
| 度数 | 32 | 21 | 30 |
| 行の割合 | 0.39 | 0.25 | 0.36 |
| 列の割合 | 0.29 | 0.23 | 0.30 |
| 30代 | |||
| 度数 | 32 | 35 | 19 |
| 行の割合 | 0.37 | 0.41 | 0.22 |
| 列の割合 | 0.29 | 0.39 | 0.19 |
| 40代 | |||
| 度数 | 32 | 24 | 38 |
| 行の割合 | 0.34 | 0.26 | 0.40 |
| 列の割合 | 0.29 | 0.27 | 0.38 |
| 50代以上 | |||
| 度数 | 14 | 10 | 13 |
| 行の割合 | 0.38 | 0.27 | 0.35 |
| 列の割合 | 0.13 | 0.11 | 0.13 |
3.3.1 グラフにする
- 集計表を作る際に作ったデータを流用すればよい。
# グラフ用データ
closs.graph.tab <- t(prop.table(tab01, margin = 1))
closs.graph.tab
20代 30代 40代 50代以上
全く 0.3855422 0.3720930 0.3404255 0.3783784
時々 0.2530120 0.4069767 0.2553191 0.2702703
頻繁 0.3614458 0.2209302 0.4042553 0.3513514# 使う色を指定しておく
color=c("gray25","gray49","gray73")
# 凡例も
par(mar=c(4,4,4,12)) # 凡例が入るように
barplot(closs.graph.tab, # グラフ用に行列を渡す
col = color,
horiz = TRUE,
main = "年齢となんとかの頻度の関係",
xlab = "割合",
ylab = "年齢")
##凡例をプロットエリア外に描けるようにする
par(xpd = TRUE)
# 凡例
legend(x = par()$usr[2], y = par()$usr[4], # 凡例がグラフエリア右下
legend = c("全く","時々","頻繁"),
xjust = -.05, #プロットエリアから少し離して凡例を載せる
fill = color # 色をグラフと揃える
)