記述統計からクロス集計まで
Koyo Yamamori
5/7/2022
最初に入れておく（RstudioでRmdファイルを新規作成したらついてくる）
knitr::opts_chunk$set(echo = TRUE)今回の内容
- csvデータをRで読み込む
- 記述統計（平均，標準偏差など）を求める
- 適切なヒストグラムを作成する
- クロス集計表を作成する
- クロス修正表に対応したグラフを作成する
- 出力を画像ファイルで得て保存する
準備
疑似データの作成
学力検査データ
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionkokuA.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)
setwd("/Users/Shared/R")
#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")
setwd("/Users/Shared/R")
#write.csv(ques, "ques.csv", row.names = F)本題
データの読み込み
csvファイルのあるディレクトリを指定し，R上での任意のデータ名を付けて，csvファイルを読み込む。
setwd("/Users/Shared/R")
# Windowsの場合には，おそらくsetwd("c:\Users\Shared\R")のような感じと思われる
ach <- read.csv("ach.csv")
# データ名，変数名は数字では始めない ""で囲み，拡張子も付ける
ques <- read.csv("ques.csv") 日本語が入っているデータの場合，文字化けが生じることがある。その場合には文字コードを指定する。データに日本語は入れない方が無難
ach <- read.csv("ach.csv", encoding="Shift_JIS")このようにデータが読み込まれている。headで最初の数行を出力する。
head(ach)## id kokuA kokuB sansA sansB
## 1 1 10 6 8 7
## 2 2 7 9 11 13
## 3 3 8 7 12 12
## 4 4 9 3 14 10
## 5 5 9 9 14 10
## 6 6 11 5 12 8head(ques)## id schl age gender q01 q02 q03 q04
## 1 1 10001 24 1 1 3 2 1
## 2 2 10001 44 1 1 3 2 1
## 3 3 10001 25 2 3 3 1 1
## 4 4 10001 27 1 3 1 1 2
## 5 5 10001 34 1 1 3 2 1
## 6 6 10001 48 1 3 1 2 3学力検査データの処理
通過率をデータに加える
素点ではなく通過率を扱いたいことがよくある（と思われる）。そうしたければ，変数を計算して追加する。たとえば，以下では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 6 8 7 0.7142857 0.6666667 0.5000 0.5384615
## 2 2 7 9 11 13 0.5000000 1.0000000 0.6875 1.0000000
## 3 3 8 7 12 12 0.5714286 0.7777778 0.7500 0.9230769
## 4 4 9 3 14 10 0.6428571 0.3333333 0.8750 0.7692308
## 5 5 9 9 14 10 0.6428571 1.0000000 0.8750 0.7692308
## 6 6 11 5 12 8 0.7857143 0.5555556 0.7500 0.6153846とりあえず記述統計をしたい
データ全体に対する要約(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 :5.500 Median :11.000
## Mean :1000.5 Mean : 8.329 Mean :5.212 Mean : 9.822
## 3rd Qu.:1500.2 3rd Qu.:11.000 3rd Qu.:7.000 3rd Qu.:12.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.6111 Median :0.6875
## Mean : 7.436 Mean :0.5950 Mean :0.5791 Mean :0.6139
## 3rd Qu.:10.000 3rd Qu.:0.7857 3rd Qu.:0.7778 3rd Qu.:0.7500
## 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.5720
## 3rd Qu.:0.7692
## Max. :1.0000標準偏差などの記述(descirptive statistics)を知りたい場合には，psychパッケージを使う。ない場合にはinstall.packaes("psych")としてインストールする。
library(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.33 3.29 9.00 8.63 2.97 0 14 14 -0.73
## kokuB 3 2000 5.21 2.41 5.50 5.38 2.22 0 9 9 -0.49
## sansA 4 2000 9.82 3.79 11.00 10.25 2.97 0 16 16 -0.94
## sansB 5 2000 7.44 3.11 8.00 7.71 2.97 0 13 13 -0.68
## kokuA_prop 6 2000 0.59 0.24 0.64 0.62 0.21 0 1 1 -0.73
## kokuB_prop 7 2000 0.58 0.27 0.61 0.60 0.25 0 1 1 -0.49
## sansA_prop 8 2000 0.61 0.24 0.69 0.64 0.19 0 1 1 -0.94
## sansB_prop 9 2000 0.57 0.24 0.62 0.59 0.23 0 1 1 -0.68
## kurtosis se
## id -1.20 12.91
## kokuA -0.04 0.07
## kokuB -0.55 0.05
## sansA 0.23 0.08
## sansB -0.15 0.07
## kokuA_prop -0.04 0.01
## kokuB_prop -0.55 0.01
## sansA_prop 0.23 0.01
## sansB_prop -0.15 0.01論文に載せられる表を出力する
上のように表示されたデータを，いわゆる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$m <- as.numeric(ach.table$m)
ach.table$sd <- as.numeric(ach.table$sd)
library(tidyverse)
library(gt)
ach.table %>% gt()| exam | m | sd |
|---|---|---|
| 国語A | 0.5949643 | 0.2350786 |
| 国語B | 0.5790556 | 0.2680215 |
| 算数A | 0.6139062 | 0.2366693 |
| 算数B | 0.5720385 | 0.2393384 |
わりときれいな表になっている。これを，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.61 | 0.24 |
| 算数B | 0.57 | 0.24 |
さらに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.61 | 0.24 |
| 算数B | 0.57 | 0.24 |
もっと工夫すると，次のようなことができる。
# 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.33 | 3.29 | 0.59 | 0.24 |
| 国語B | 5.21 | 2.41 | 0.58 | 0.27 |
| 算数A | 9.82 | 3.79 | 0.61 | 0.24 |
| 算数B | 7.44 | 3.11 | 0.57 | 0.24 |
タイトルも付けてみる
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.33 | 3.29 | 0.59 | 0.24 |
| 国語B | 5.21 | 2.41 | 0.58 | 0.27 |
| 算数A | 9.82 | 3.79 | 0.61 | 0.24 |
| 算数B | 7.44 | 3.11 | 0.57 | 0.24 |
注釈も付けてみる
gt.rem <- gt.title %>%
# ここは表の内容の脚注
tab_footnote(
footnote = "正答数",
location = cells_column_spanners(spanners = "素点")) %>%
tab_footnote(footnote = "平均",
locations = cells_column_labels(columns = c(m, m_prop))) %>%
tab_footnote(footnote = "標準偏差",
locations = cells_column_labels(columns = c(sd, sd_prop))) %>%
tab_options(table.border.bottom.width = 0) %>% #注の下の線は無し
tab_source_note(source_note = md("このページの「疑似データの作成」セクションで作ったデータから作表したもので*実際のデータではなく*__疑似データを集計したものである。__"))
gt.rem| Table 1 教科別の平均と標準偏差 |
||||
|---|---|---|---|---|
| 教科 | 素点1 | 通過率 | ||
| M2 | SD3 | M2 | SD3 | |
| 国語A | 8.33 | 3.29 | 0.59 | 0.24 |
| 国語B | 5.21 | 2.41 | 0.58 | 0.27 |
| 算数A | 9.82 | 3.79 | 0.61 | 0.24 |
| 算数B | 7.44 | 3.11 | 0.57 | 0.24 |
| このページの「疑似データの作成」セクションで作ったデータから作表したもので実際のデータではなく疑似データを集計したものである。 | ||||
| 1 正答数 | ||||
| 2 平均 | ||||
| 3 標準偏差 | ||||
出力を画像として保存する
あとから微修正ができなくなるが，信頼を得るために必要なことである。コードの下に出力されるのは生成されたpngファイルである。任意のディレクトリに保存されるので，あとは貼り付けるだけでよい。
setwd("/Users/Shared/R/Fig") # 保存先のディレクトリ
gtsave(gt.rem, "表の例_01.png")
setwd("/Users/Shared/R") # ディレクトリを元に戻しておくヒストグラムを作る
今一度，学力検査の疑似データの内容を見てみる。
head(ach)## id kokuA kokuB sansA sansB kokuA_prop kokuB_prop sansA_prop sansB_prop
## 1 1 10 6 8 7 0.7142857 0.6666667 0.5000 0.5384615
## 2 2 7 9 11 13 0.5000000 1.0000000 0.6875 1.0000000
## 3 3 8 7 12 12 0.5714286 0.7777778 0.7500 0.9230769
## 4 4 9 3 14 10 0.6428571 0.3333333 0.8750 0.7692308
## 5 5 9 9 14 10 0.6428571 1.0000000 0.8750 0.7692308
## 6 6 11 5 12 8 0.7857143 0.5555556 0.7500 0.6153846例えば，国語Aの正答数について，単純にヒストグラムを出力することは，たいしたことはない。
hist(ach$kokuA)
データの様子を眺める程度なら，これで問題ない。しかし，国語Aの場合の得点の範囲は0から14（14段階）なのに，棒の数は13しかない。
hist(ach$kokuA, breaks=20)
とかすれば，それなりの出力は得られるが，0点の人が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
## 48 58 55 51 69 97 104 166 246 276 283 220 180 129 18単純に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) # 棒の間を詰めてヒストグラム風にする
)
ただし，本来，barplotは，barの間が空いているのが普通なので，本当に大丈夫かというエラーメッセージが出る。
通過率でbarplotを作る
通過率のテーブルを作る
freq.kokuA <- table(ach$kokuA) # これは先ほどやったことと同じ
prop.kokuA <- freq.kokuA / sum(freq.kokuA) # 各セルの数を全体の数で割った表
prop.kokuA##
## 0 1 2 3 4 5 6 7 8 9 10
## 0.0240 0.0290 0.0275 0.0255 0.0345 0.0485 0.0520 0.0830 0.1230 0.1380 0.1415
## 11 12 13 14
## 0.1100 0.0900 0.0645 0.0090最大値が0.15あたりなので，これに合わせたbarplotを，先ほどと同様の方法で作る。
barplot(prop.kokuA,
xlab = "正答数", ylab = "割合", # 縦，横軸に名前を付ける
ylim = c(0, 0.2), # 縦軸を正答数ごとの人数の最大値よりも大きく
cex.names = 0.8, # 軸ラベルをやや小さくする
cex.axis = 0.8, # 縦軸のラベル（数字）をやや小さめにする
axis.lty = 1, # 横軸に線を引く
space = c(0, 0) # 棒の間を詰めてヒストグラム風にする
)
正答数の刻みを変えたい
数十項目あるテストの場合，1点刻みだとグラフが見づらい。上記のデータは0点から14点の範囲なので，3点刻みにしてみる。
# データフレームにして集計できるようにする
df.fq.kokuA <- data.frame(freq.kokuA)
colnames(df.fq.kokuA) <- c("score", "freq")
# 得点カテゴリの変数を追加する
cat.fq.kokuA <- df.fq.kokuA %>%
dplyr::mutate(
cat = case_when(
score == "0" | score == "1" | score == "2" ~ 1,
score == "3" | score == "4" | score == "5" ~ 2,
score == "6" | score == "7" | score == "8" ~ 3,
score == "9" | score == "10" | score == "11" ~ 4,
score == "12" | score == "13" | score == "14" ~ 5,
TRUE ~ 0)
)
cat.fq.kokuA## score freq cat
## 1 0 48 1
## 2 1 58 1
## 3 2 55 1
## 4 3 51 2
## 5 4 69 2
## 6 5 97 2
## 7 6 104 3
## 8 7 166 3
## 9 8 246 3
## 10 9 276 4
## 11 10 283 4
## 12 11 220 4
## 13 12 180 5
## 14 13 129 5
## 15 14 18 5こういうデータになる。さらに続ける。
# 得点カテゴリごとの人数の和をデータ化する
cat.kokuA <- data.frame(
cat.fq.kokuA %>%
group_by(cat) %>%
dplyr::summarise(sum(freq))
)
cat.kokuA## cat sum.freq.
## 1 1 161
## 2 2 217
## 3 3 516
## 4 4 779
## 5 5 327このデータのsum.freq.をヒストグラムにする
barplot(cat.kokuA$sum.freq.,
xlab = "正答数", ylab = "人数", # 縦，横軸に名前を付ける
ylim = c(0, 1000), # 縦軸を正答数ごとの人数の最大値よりも大きく
cex.names = 1.0, # 軸ラベルをやや小さくする
cex.axis = 0.8, # 縦軸のラベル（数字）をやや小さめにする
axis.lty = 1, # 横軸に線を引く
names.arg=c("0-2", "3-5", "6-8", "9-11", "12-14"))
画像として保存する
画質がいいのはpdfである
setwd("/Users/Shared/R/Fig") # 保存先のディレクトリ
#pdf("barplot.pdf", width=6, height=6) # Windowsの場合
#par(family="Japan1GothicBBB")# Windowsの場合
quartz(type="pdf", file="barplot.pdf", width=6, height=6) #Macintoshの場合
par(family="HiraKakuProN-W3")#Macintoshの場合
barplot(cat.kokuA$sum.freq.,
xlab = "正答数", ylab = "人数", # 縦，横軸に名前を付ける
ylim = c(0, 1000), # 縦軸を正答数ごとの人数の最大値よりも大きく
cex.names = 1.0, # 軸ラベルをやや小さくする
cex.axis = 0.8, # 縦軸のラベル（数字）をやや小さめにする
axis.lty = 1, # 横軸に線を引く
names.arg=c("0-2", "3-5", "6-8", "9-11", "12-14"))
dev.off()## quartz_off_screen
## 2setwd("/Users/Shared/R")質問紙調査データの例
データはこのような形
head(ques)## id schl age gender q01 q02 q03 q04
## 1 1 10001 24 1 1 3 2 1
## 2 2 10001 44 1 1 3 2 1
## 3 3 10001 25 2 3 3 1 1
## 4 4 10001 27 1 3 1 1 2
## 5 5 10001 34 1 1 3 2 1
## 6 6 10001 48 1 3 1 2 3回答者の属性
性別
table(ques$gender)##
## 1 2
## 162 138年齢
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
## 13 9 5 11 10 10 11 4 9 17 8 5 10 8 12 5 12 10 11 9 7 5 8 9 12 10
## 48 49 50 51 52 53
## 5 7 12 13 11 12年齢の平均や標準偏差
library(psych)
describe(ques$age)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 300 37.65 9.53 38 37.67 11.86 22 53 31 0.01 -1.25 0.55グラフにする
割合を積み上げ棒グラフで表現するとよい。
barplot(matrix(prop.table(table(ques$gender))), horiz =TRUE)
これだと見栄えが悪いので，工夫する。#以降改行までは処理されない。 グラフには数値は載せない（役所は違うようだが）
#setwd("/Users/Shared/R/Fig")
#bmp(filename="g01.bmp", width = 1000, height = 300)
par(mar=c(4,4,4,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()
#setwd("/Users/Shared/R")回答の年齢層は，層別するとよさそうである。 層別にする方法は，学力検査データの例と同じである。
# 年齢カテゴリの変数を追加する
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 73
## 2 2 96
## 3 3 83
## 4 4 48これをグラフにする
# 上記のデータを因子型にする
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 # 色をグラフと揃える
)
項目の集計
psychパッケージのdescribe関数で分かる。
library(psych) # 一度読み込めばよい
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.65 9.53 38.0 37.67 11.86 22 53 31
## gender 4 300 1.46 0.50 1.0 1.45 0.00 1 2 1
## q01 5 300 2.04 0.80 2.0 2.05 1.48 1 3 2
## q02 6 300 1.95 0.82 2.0 1.93 1.48 1 3 2
## q03 7 300 2.09 0.79 2.0 2.11 1.48 1 3 2
## q04 8 300 2.03 0.83 2.0 2.04 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.25 0.55
## gender 0.16 -1.98 0.03
## q01 -0.07 -1.45 0.05
## q02 0.10 -1.53 0.05
## q03 -0.16 -1.40 0.05
## q04 -0.06 -1.56 0.05クロス集計
例えば，年齢と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 24 1 全く 3 2 1 20代
## 2 2 10001 44 1 全く 3 2 1 40代
## 3 3 10001 25 2 頻繁 3 1 1 20代
## 4 4 10001 27 1 頻繁 1 1 2 20代
## 5 5 10001 34 1 全く 3 2 1 30代
## 6 6 10001 48 1 頻繁 1 2 3 40代データが変わっていることが分かる。
クロス集計自体は，簡単である。
tab01 <- table(ques2$cat, ques2$q01)
tab01##
## 全く 時々 頻繁
## 20代 24 19 30
## 30代 31 39 26
## 40代 23 31 29
## 50代以上 13 18 17用途によって様々な形式がある。
# 行と列の合計欄を付ける
tab02 <- addmargins(tab01)
tab02##
## 全く 時々 頻繁 Sum
## 20代 24 19 30 73
## 30代 31 39 26 96
## 40代 23 31 29 83
## 50代以上 13 18 17 48
## Sum 91 107 102 300# 行方向の割合
tab03 <- addmargins(prop.table(tab01, margin = 1))
tab03##
## 全く 時々 頻繁 Sum
## 20代 0.3287671 0.2602740 0.4109589 1.0000000
## 30代 0.3229167 0.4062500 0.2708333 1.0000000
## 40代 0.2771084 0.3734940 0.3493976 1.0000000
## 50代以上 0.2708333 0.3750000 0.3541667 1.0000000
## Sum 1.1996256 1.4150179 1.3853565 4.0000000# 有効数字を少なくする
round(tab03, 2)##
## 全く 時々 頻繁 Sum
## 20代 0.33 0.26 0.41 1.00
## 30代 0.32 0.41 0.27 1.00
## 40代 0.28 0.37 0.35 1.00
## 50代以上 0.27 0.38 0.35 1.00
## Sum 1.20 1.42 1.39 4.00# 列方向の割合
tab04 <- addmargins(prop.table(tab01, margin = 2))
round(tab04,2) # 和が100にならない場合がある##
## 全く 時々 頻繁 Sum
## 20代 0.26 0.18 0.29 0.74
## 30代 0.34 0.36 0.25 0.96
## 40代 0.25 0.29 0.28 0.83
## 50代以上 0.14 0.17 0.17 0.48
## 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 24.0000000 19.0000000 30.0000000
## tab.r.20 0.3287671 0.2602740 0.4109589
## tab.c.20 0.2637363 0.1775701 0.2941176
## tab.n.30 31.0000000 39.0000000 26.0000000
## tab.r.30 0.3229167 0.4062500 0.2708333
## tab.c.30 0.3406593 0.3644860 0.2549020
## tab.n.40 23.0000000 31.0000000 29.0000000
## tab.r.40 0.2771084 0.3734940 0.3493976
## tab.c.40 0.2527473 0.2897196 0.2843137
## tab.n.50 13.0000000 18.0000000 17.0000000
## tab.r.50 0.2708333 0.3750000 0.3541667
## tab.c.50 0.1428571 0.1682243 0.1666667これを，報告できるような表に変換する
library(tidyverse) # 必要に応じて読み込む
library(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代 度数 24.0000000 19.0000000 30.0000000
## tab.r.20 20代 行の割合 0.3287671 0.2602740 0.4109589
## tab.c.20 20代 列の割合 0.2637363 0.1775701 0.2941176
## tab.n.30 30代 度数 31.0000000 39.0000000 26.0000000
## tab.r.30 30代 行の割合 0.3229167 0.4062500 0.2708333
## tab.c.30 30代 列の割合 0.3406593 0.3644860 0.2549020
## tab.n.40 40代 度数 23.0000000 31.0000000 29.0000000
## tab.r.40 40代 行の割合 0.2771084 0.3734940 0.3493976
## tab.c.40 40代 列の割合 0.2527473 0.2897196 0.2843137
## tab.n.50 50代以上 度数 13.0000000 18.0000000 17.0000000
## tab.r.50 50代以上 行の割合 0.2708333 0.3750000 0.3541667
## tab.c.50 50代以上 列の割合 0.1428571 0.1682243 0.1666667cross.tab <- tab.40.df.2 %>% gt()
cross.tab| agegroup | out | 全く | 時々 | 頻繁 |
|---|---|---|---|---|
| 20代 | 度数 | 24.0000000 | 19.0000000 | 30.0000000 |
| 20代 | 行の割合 | 0.3287671 | 0.2602740 | 0.4109589 |
| 20代 | 列の割合 | 0.2637363 | 0.1775701 | 0.2941176 |
| 30代 | 度数 | 31.0000000 | 39.0000000 | 26.0000000 |
| 30代 | 行の割合 | 0.3229167 | 0.4062500 | 0.2708333 |
| 30代 | 列の割合 | 0.3406593 | 0.3644860 | 0.2549020 |
| 40代 | 度数 | 23.0000000 | 31.0000000 | 29.0000000 |
| 40代 | 行の割合 | 0.2771084 | 0.3734940 | 0.3493976 |
| 40代 | 列の割合 | 0.2527473 | 0.2897196 | 0.2843137 |
| 50代以上 | 度数 | 13.0000000 | 18.0000000 | 17.0000000 |
| 50代以上 | 行の割合 | 0.2708333 | 0.3750000 | 0.3541667 |
| 50代以上 | 列の割合 | 0.1428571 | 0.1682243 | 0.1666667 |
さらに，加工する。
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代 | |||
| 度数 | 24 | 19 | 30 |
| 行の割合 | 0.33 | 0.26 | 0.41 |
| 列の割合 | 0.26 | 0.18 | 0.29 |
| 30代 | |||
| 度数 | 31 | 39 | 26 |
| 行の割合 | 0.32 | 0.41 | 0.27 |
| 列の割合 | 0.34 | 0.36 | 0.25 |
| 40代 | |||
| 度数 | 23 | 31 | 29 |
| 行の割合 | 0.28 | 0.37 | 0.35 |
| 列の割合 | 0.25 | 0.29 | 0.28 |
| 50代以上 | |||
| 度数 | 13 | 18 | 17 |
| 行の割合 | 0.27 | 0.38 | 0.35 |
| 列の割合 | 0.14 | 0.17 | 0.17 |
グラフにする
集計表を作る際に作ったデータを流用すればよい。
# グラフ用データ
closs.graph.tab <- t(prop.table(tab01, margin = 1))
closs.graph.tab##
## 20代 30代 40代 50代以上
## 全く 0.3287671 0.3229167 0.2771084 0.2708333
## 時々 0.2602740 0.4062500 0.3734940 0.3750000
## 頻繁 0.4109589 0.2708333 0.3493976 0.3541667# 使う色を指定しておく
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 # 色をグラフと揃える
)