library("tidyverse")
participant <- read_csv("data/participant-info.csv")
data <- read_csv("data/ahi-cesd.csv") Tidyverseを用いたデータハンドリングと可視化
国里愛彦(専修大学)
データの説明
- 以下の論文のデータを使う。
- データは以下のfigshareで公開されている。
Woodworth et al.(2018)のデータ
対象者：一般参加者295名
介入：３つのポジティブ心理学的介入と統制条件に無作為割付。(1)“Using Signature Strengths（ウェブサイトで性格の強みを評価し，上位の強みのうちの１つを次の週に新しい別のやり方で実施する）”，(2)“Three Good Things(その日に起こった良いことを３つあげて，なぜそれが起こったのか原因説明とともに記録する)”, (3)“Gratitude Visit(親切してくれたけど感謝を示してない人に手紙を書く介入)”,(4)“Recording early memories(統制条件，幼いことを記憶を書く)”
アウトカム：Authentic Happiness Inventory (AHI) の24項目，Center for Epidemiological Studies Depression(CES-D) scaleの20項目
測定時期：介入前(0), 介入後(1), 1週間後(2), １ヶ月後(3)，３ヶ月後(4)，６ヶ月後(5)
Woodworth et al.(2018)のデータにアクセスして，ファイルをダウンロードして開きます。
データの読み込み
データの読み込み，データの前処理，可視化などはtidyverseパッケージ使うと便利です。
dataフォルダに保存したWoodworth et al.(2018)のデータを読み込みます。データはcsv形式だったので，read_csv()で読み込みます。ファイルはparticipant-info.csvとahi-cesd.csvの２種類があります。dataフォルダにいれているので，ファイル名の前に”data/“をいれています。こうすると，dataフォルダ内の特定のファイルを指定できます。
- qmdだと以下のような感じです。Rチャンク右端の緑の三角形をクリックするとコードが実行されますので，クリックしてください。
- コードを実行するとEnvironmentにオブジェクトが保存されます。右端のアイコンをクリックするとオブジェクト（今回はデータ）を見ることができます。
- participant-info.csvには，id, intervention, sex, age, educ, incomeのデータが入っています。
- ahi-cesd.csvには，id, occasion, elapsed.days, intervention, ah01…などのデータが入っています。
- 個人的にはエクセルファイルはcsv形式に変換して処理するのが良いですが，エクセルから直接データを読みたい場合は，readxlパッケージのread_excel()を使って読み込むことができます。
- foreignパッケージを使うと，‘Epi Info’, ‘Minitab’, ‘S’, ‘SAS’, ‘SPSS’, ‘Stata’, ‘Systat’, ’Weka’などのファイルも読み込めます。
ネイティブパイプの使用
Rではパイプ演算子という便利なものを使います。パイプ演算子でつなぐと前の処理の内容が次の処理に引き継がれるので，可読性高くコードが書けます。以前は
magrittrパッケージを用いた%>%だったのですが，R4.1からは標準で|>が使えます。|>が使えるように設定します。Tools -> Global options -> Codeで，Use native pipe operator…にチェックをいれてApplyする。
Windowsの場合Ctrl + Shift + m Macの場合Command + Shift + mのショートカットでパイプ演算子が入力できます。
データフレーム
Rではデータフレームという構造でデータを扱う。データフレームは，a行b列からなる表形式のデータであり，列には変数が含まれる。
read_csv()で読み込んだデータはデータフレームになっており，オブジェクト名を打つとデータフレームが表示される。
Rでの代表的なデータ型としては以下がある（read_csv()で読み込むと自動で判定をする）。
| データ型 | 説明 | 省略形 |
|---|---|---|
| logical | 論理型（TRUE, FALSE） | lgl |
| integer | 整数型（1,2） | int |
| double | 倍精度小数点型 (1.0,2.1) | dbl |
| character | 文字型 | chr |
| factor | 因子型 | fct |
| ordered | 順序型 | ord |
| Date | 日付型 | date |
データの前処理
前処理には
tidyverseパッケージに含まれているdplyrパッケージを使う。色々と機能はあるが，ミニマムには以下の８つを覚える。filter()行の絞り込みarrange()行の並び替えselect()指定した列の抽出mutate()列の追加group_by()グループ化summarize()要約pivot_longer()横データから縦データへ変換（時系列解析などは縦データが便利）pivot_wider()縦データから横データへ変換
filter()
filter()で行の絞り込みができます。性別を女性(1,男性は2)だけに絞り込む場合は以下のようにできます。パイプ演算子は，左にあるオブジェクトを右（改行していることが多い）の関数にいれる機能を持っています。パイプ演算子を使うと処理を重ねていけるので見やすいし便利です。今回は，participantデータをfilter()にいれて，sexが1のものだけを抽出しています。
パイプ演算子の挿入はWindowsの場合Ctrl + Shift + m Macの場合Command + Shift + mです。
participant |>
filter(sex == 1)# A tibble: 251 × 6
id intervention sex age educ income
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2 1 1 59 1 1
2 3 4 1 51 4 3
3 4 3 1 50 5 2
4 6 1 1 31 5 1
5 7 3 1 44 5 2
6 8 2 1 57 4 2
7 9 1 1 36 4 3
8 10 2 1 45 4 3
9 11 2 1 56 5 1
10 12 2 1 46 4 3
# ℹ 241 more rows- filter()では，XX以上でしぼりこみもできる。以下は，年齢が40以上で絞り込んでいる。
summarize()を使うと，変数名=関数（変数）で要約ができます。最小40歳，最大85歳ですね。
participant |>
filter(age >= 40) |>
summarize(min_age = min(age), max_age = max(age))# A tibble: 1 × 2
min_age max_age
<dbl> <dbl>
1 40 83以下のように４０歳以下かつ女性で絞り込むこともできます。
participant |>
filter(age <= 40 & sex == 2)# A tibble: 17 × 6
id intervention sex age educ income
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 4 2 35 5 3
2 15 4 2 27 2 2
3 44 2 2 36 4 3
4 130 4 2 21 1 1
5 140 3 2 31 3 1
6 193 2 2 25 5 2
7 204 4 2 32 5 1
8 210 1 2 28 4 2
9 234 3 2 33 5 3
10 235 2 2 18 2 1
11 237 4 2 33 3 3
12 243 4 2 26 4 2
13 251 4 2 31 5 2
14 263 1 2 32 5 2
15 269 3 2 23 4 1
16 284 3 2 19 4 1
17 288 3 2 26 5 1arrange()
arrange()で行の並び替えができます。指定した変数で昇順で並び替えます。
participant |>
arrange(age)# A tibble: 295 × 6
id intervention sex age educ income
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 235 2 2 18 2 1
2 267 4 1 19 2 1
3 284 3 2 19 4 1
4 20 3 1 21 4 1
5 130 4 2 21 1 1
6 230 2 1 21 2 1
7 23 1 1 22 4 1
8 122 1 1 23 5 2
9 138 1 1 23 2 1
10 247 1 1 23 4 1
# ℹ 285 more rows- arrange(-age)のように-を使うことで降順にすることもできます。
participant |>
arrange(-age)# A tibble: 295 × 6
id intervention sex age educ income
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 51 2 1 83 2 2
2 244 2 1 77 3 2
3 215 4 1 75 3 2
4 83 4 1 71 2 1
5 114 1 1 71 4 1
6 155 3 2 71 4 1
7 93 3 1 68 1 2
8 141 2 2 67 5 2
9 281 4 2 67 4 2
10 67 1 1 65 4 2
# ℹ 285 more rowsselect()
select()で指定した列の抽出ができます。例えば，id,age,sexだけ抽出したい場合は，以下のように選択できます。
participant |>
select(id,age,sex)# A tibble: 295 × 3
id age sex
<dbl> <dbl> <dbl>
1 1 35 2
2 2 59 1
3 3 51 1
4 4 50 1
5 5 58 2
6 6 31 1
7 7 44 1
8 8 57 1
9 9 36 1
10 10 45 1
# ℹ 285 more rowsmutate()
mutate()で列の追加ができます。例えば，dataにはCES-Dの得点が入っていますが，その中でポジティブ感情の項目（4,8,12,16）だけを合計して得点化するとします。確認のためにsummarizeもしています。
data |>
mutate(positive_score = cesd04 + cesd08 + cesd12 + cesd16) |>
summarize(mean = mean(positive_score),
sd = sd(positive_score),
min = min(positive_score),
max = max(positive_score))# A tibble: 1 × 4
mean sd min max
<dbl> <dbl> <dbl> <dbl>
1 12.2 3.28 4 16group_by()
group_by()でグループ化します。dataのCES-Dの合計得点(cesdTotal)の介入前の得点(occasionが0)に対して，interventionでグループ分けして，summarizeします。
data |>
filter(occasion == 0) |>
select(intervention, cesdTotal) |>
group_by(intervention) |>
summarize(mean = mean(cesdTotal),
sd = sd(cesdTotal),
n = n())# A tibble: 4 × 4
intervention mean sd n
<dbl> <dbl> <dbl> <int>
1 1 15.1 11.7 72
2 2 16.2 9.89 76
3 3 16.1 11.8 74
4 4 12.8 9.51 73pivot_longer()とpivot_wider()
pivot_longer()で横データから縦データへ変換，pivot_wider()で縦データから横データへ変換します。今回のdataはすでにlong形式になっています（occasionごとに横にデータが並んでいるのではなく，１つの参加者でoccationごとに複数行が並んでいます）。どうもこのデータは，同じ時点で複数回答している人がいるようです。idとoccasionでグループ化して，回数が１以上のデータを探すとid8のoccasion2,id64のoccasion4が二回回答されています。
data |>
group_by(id,occasion) |>
filter(n() > 1) |>
ungroup() # A tibble: 4 × 50
id occasion elapsed.days intervention ahi01 ahi02 ahi03 ahi04 ahi05 ahi06
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 8 2 21.9 2 2 1 1 1 1 1
2 8 2 38.9 2 2 2 1 1 2 1
3 64 4 103. 1 2 3 2 3 2 3
4 64 4 140. 1 3 3 2 3 2 3
# ℹ 40 more variables: ahi07 <dbl>, ahi08 <dbl>, ahi09 <dbl>, ahi10 <dbl>,
# ahi11 <dbl>, ahi12 <dbl>, ahi13 <dbl>, ahi14 <dbl>, ahi15 <dbl>,
# ahi16 <dbl>, ahi17 <dbl>, ahi18 <dbl>, ahi19 <dbl>, ahi20 <dbl>,
# ahi21 <dbl>, ahi22 <dbl>, ahi23 <dbl>, ahi24 <dbl>, cesd01 <dbl>,
# cesd02 <dbl>, cesd03 <dbl>, cesd04 <dbl>, cesd05 <dbl>, cesd06 <dbl>,
# cesd07 <dbl>, cesd08 <dbl>, cesd09 <dbl>, cesd10 <dbl>, cesd11 <dbl>,
# cesd12 <dbl>, cesd13 <dbl>, cesd14 <dbl>, cesd15 <dbl>, cesd16 <dbl>, …filter()でid8のoccasion2,id64のoccasion4の二回目の回答は除外します。select(id,occasion,intervention,cesdTotal)で少し変数をしぼって，pivot_wider(names_from = occasion, values_from = cesdTotal)で横データにします。最後に，一番前をdata_wide <-とするとdata_wideというオブジェクト名として保存されます。
data_wide <- data |>
filter(!(id == 8 & cesdTotal == 49)) |>
filter(!(id == 64 & cesdTotal == 16)) |>
select(id,occasion,intervention,cesdTotal) |>
pivot_wider(names_from = occasion, values_from = cesdTotal) - 今度は，
pivot_longer()で縦データを横データに変換します。pivot_longer(cols = c("移動する列"), names_to = "分ける変数名", values_to = "値の名前")で横データになります。ただ，元に戻るのではなく，横→縦変換時に欠測になっているので，横データに戻した時に欠測が含まれています。
data_wide |>
pivot_longer(cols = c(-id,-intervention), names_to = "occasion", values_to = "cesdTotal")# A tibble: 1,770 × 4
id intervention occasion cesdTotal
<dbl> <dbl> <chr> <dbl>
1 1 4 0 14
2 1 4 1 6
3 1 4 2 NA
4 1 4 3 NA
5 1 4 4 NA
6 1 4 5 NA
7 2 1 0 7
8 2 1 1 10
9 2 1 2 13
10 2 1 3 8
# ℹ 1,760 more rows記述統計と効果量
glimpse()
- Environmentでも確認できますが，glimpse()でデータの変数を確認できる。
participant |>
glimpse()Rows: 295
Columns: 6
$ id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
$ intervention <dbl> 4, 1, 4, 3, 2, 1, 3, 2, 1, 2, 2, 2, 4, 4, 4, 4, 3, 2, 1, …
$ sex <dbl> 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, …
$ age <dbl> 35, 59, 51, 50, 58, 31, 44, 57, 36, 45, 56, 46, 34, 41, 2…
$ educ <dbl> 5, 1, 4, 5, 5, 5, 5, 4, 4, 4, 5, 4, 5, 1, 2, 1, 4, 5, 3, …
$ income <dbl> 3, 1, 3, 2, 2, 1, 2, 2, 3, 3, 1, 3, 3, 2, 2, 1, 2, 2, 1, …summary()
- 連続変数に関しては，summary()で必要な記述統計を出すことができる（selectを組み合わせると便利）。
participant |>
select(age, educ, income) |>
summary() age educ income
Min. :18.00 Min. :1.00 Min. :1.000
1st Qu.:34.00 1st Qu.:3.50 1st Qu.:2.000
Median :44.00 Median :4.00 Median :2.000
Mean :43.76 Mean :3.98 Mean :2.044
3rd Qu.:52.00 3rd Qu.:5.00 3rd Qu.:3.000
Max. :83.00 Max. :5.00 Max. :3.000 table()
- 質的変数の場合は，table()でクロス集計を出すことができる。このデータは女性が多いことが分かる。
participant |>
select(sex) |>
table()sex
1 2
251 44 gtsummary()
gtsummaryパッケージのtbl_summary()を使うと簡単に変数の要約ができます。
library(gtsummary)
participant |>
select(-id) |>
tbl_summary()| Characteristic | N = 2951 |
|---|---|
| intervention | |
| 1 | 72 (24%) |
| 2 | 76 (26%) |
| 3 | 74 (25%) |
| 4 | 73 (25%) |
| sex | |
| 1 | 251 (85%) |
| 2 | 44 (15%) |
| age | 44 (34, 52) |
| educ | |
| 1 | 14 (4.7%) |
| 2 | 21 (7.1%) |
| 3 | 39 (13%) |
| 4 | 104 (35%) |
| 5 | 117 (40%) |
| income | |
| 1 | 73 (25%) |
| 2 | 136 (46%) |
| 3 | 86 (29%) |
| 1 n (%); Median (IQR) | |
byでグループ分けをして表を作ることもできる。
participant |>
select(-id) |>
tbl_summary(by = intervention)| Characteristic | 1, N = 721 | 2, N = 761 | 3, N = 741 | 4, N = 731 |
|---|---|---|---|---|
| sex | ||||
| 1 | 62 (86%) | 66 (87%) | 62 (84%) | 61 (84%) |
| 2 | 10 (14%) | 10 (13%) | 12 (16%) | 12 (16%) |
| age | 45 (36, 54) | 46 (36, 53) | 44 (34, 51) | 40 (32, 51) |
| educ | ||||
| 1 | 3 (4.2%) | 3 (3.9%) | 2 (2.7%) | 6 (8.2%) |
| 2 | 3 (4.2%) | 8 (11%) | 3 (4.1%) | 7 (9.6%) |
| 3 | 10 (14%) | 7 (9.2%) | 7 (9.5%) | 15 (21%) |
| 4 | 26 (36%) | 32 (42%) | 32 (43%) | 14 (19%) |
| 5 | 30 (42%) | 26 (34%) | 30 (41%) | 31 (42%) |
| income | ||||
| 1 | 21 (29%) | 20 (26%) | 17 (23%) | 15 (21%) |
| 2 | 31 (43%) | 34 (45%) | 36 (49%) | 35 (48%) |
| 3 | 20 (28%) | 22 (29%) | 21 (28%) | 23 (32%) |
| 1 n (%); Median (IQR) | ||||
gtsummary()の調整法としては，こちらが参考になります。
より柔軟に表を作成するパッケージとしては
gtパッケージがありますが，説明を省略します。
データの可視化
- データに対してすぐにでも高度な統計手法を適用したくなるが，まずは記述統計とともにデータの可視化を必ずしましょう。
- データの可視化は，(1)データに対する前処理が適切にできているかチェック，(2)データの測定法が適切にできているかチェック，(3)データの傾向をチェック，(4)データに適用する統計手法が適切かチェックするために行う。
- データの可視化を行ってないとエラーに気づかない可能性があり，結果を歪めてしまう可能性がある。
ggplot2
Rには可視化のための関数があるが，ggplot2パッケージがきれいに図が作成できて，コードも読みやすく書ける。
ggplot2は，グラフィックの文法 (grammar of graphics)という思想の元に作成されている。これは，有る決まった構造で作成されたレイヤーを重ねていって図を作成するという考え方です。
棒グラフ
- 棒グラフを題材にggplot2の使い方を学びましょう。ggplot2はtidyverseパッケージを読み込むと使えます。まずは，性別ごとの人数の棒グラフを描きます。女性か男性か図に書き込みたいので，sexの1,2の数値を女性と男性に変換した変数をsex2として追加します。
library(tidyverse)
participant <- participant |>
mutate(sex2 = ifelse(sex ==1, '女性','男性'))- まずggplotに使用するデータを指定します。この段階では図を書くキャンバスが出力されるだけです。
ggplot(data = participant) 
- キャンバスにレイヤーを追加します。追加する際にはパイプ演算子ではなくて
+を使います。aes()では，x軸, y軸の変数や色などを指定します。ここでは，１変数のsex2だけ指定します。macの場合は文字化けしているのですが，ちょっと置いておきます。
ggplot(data = participant) +
aes(x=sex2)
- 続けて，
geom_bar()で，棒グラフを重ねます。
ggplot(data = participant) +
aes(x=sex2) +
geom_bar()
- 文字化けが気になってきたので，plotのテーマはtheme_grayのまま
base_family = "HiraKakuPro-W3"を指定します。
ggplot(data = participant) +
aes(x=sex2) +
geom_bar() +
theme_gray(base_family = "HiraKakuPro-W3")
- ちょっと灰色の背景だと論文に載せにくいなと思ったかもしれません。その場合は，
theme_classic()を使います。あれ？グラフが浮いている感じして気持ち悪いですね。
ggplot(data = participant) +
aes(x=sex2) +
geom_bar() +
theme_classic(base_family = "HiraKakuPro-W3") 
scale_y_continuous(expand=c(0,0))を追加すると浮かなくなります。
ggplot(data = participant) +
aes(x=sex2) +
geom_bar() +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0))
- せっかくなので，男女でグラフの色を分けたいと思います。
aes()内のfillにsex2を指定します。
ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar() +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0))
- ちょっと色が派手で灰色と白くらいにしたいです。
ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar() +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0)) +
scale_fill_manual(values=c("white","gray"))
あ！線がない！geom_bar()内に色を追加します。
ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar(color="black") +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0)) +
scale_fill_manual(values=c("white","gray"))
- 図はいいですが，x軸とy軸のラベルを変えたいですね。
labs()を追加します。
ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar(color="black") +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0)) +
scale_fill_manual(values=c("white","gray")) +
labs(x = "性別", y = "人数")
- この図の場合は，凡例は不要ですね。
theme(legend.position = "none")を追加します。
ggplot(data = participant) +
aes(x=sex2,fill=sex2) +
geom_bar(color="black") +
theme_classic(base_family = "HiraKakuPro-W3") +
scale_y_continuous(expand=c(0,0)) +
scale_fill_manual(values=c("white","gray")) +
labs(x = "性別", y = "人数") +
theme(legend.position = "none")
使えそうな図が作れました。ggplot2は面倒な感じもしなくもないですが，キャンバスにちょっとずつレイヤーを足していくことで図を作る感じをイメージするといいです。その都度，「ggplot 凡例 消す」とかで検索すると簡単に情報が手に入ると思います（chatGPTも便利です）。
ヒストグラム
量的な変数についてはヒストグラムを書いて分布を確認します。
dataからoccasion=0のCES-D得点のヒストグラムを描きます。filterでoccasion=0に絞り込んで，それをggplotに投げます。パイプ演算子を使っているので，filterされたdataがggplotに入るので，aesだけ指定します。geom_histogram()でヒストグラムが描けます。
data |>
filter(occasion == 0) |>
ggplot(aes(x=cesdTotal)) +
geom_histogram()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
- theme_classic()にしたり，グラフが浮かないようにしたり，x軸のラベルを変えたりして，調整します。
data |>
filter(occasion == 0) |>
ggplot(aes(x=cesdTotal)) +
geom_histogram() +
theme_classic() +
scale_y_continuous(expand=c(0,0)) +
labs(x = "CES-D") `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
- ヒストグラムも良いですが，密度プロット（density plot）もあると便利です。ヒストグラムと重ねて描くために，
geom_histogram(aes(y=after_stat(density)),fill='white', color='black')という指定をしています。また，geom_density(fill='black', alpha = 0.5)のalphaは重ねた際の透過度になります。
data |>
filter(occasion == 0) |>
ggplot(aes(x=cesdTotal)) +
geom_histogram(aes(y=after_stat(density)),fill='white', color='black') +
geom_density(fill='black', alpha = 0.5) +
theme_classic() +
scale_y_continuous(expand=c(0,0)) +
labs(x = "CES-D", y="") `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
散布図
- dataからahiTotalとcesdTotalの散布図を描いてみます。相関や回帰の前には散布図を描くようにしましょう。散布図は
geom_point()で描けます。
data |>
filter(occasion == 0) |>
ggplot(aes(x = cesdTotal, y = ahiTotal)) +
geom_point()
- 論文用に調整します。きれいに負の相関をしています。
data |>
filter(occasion == 0) |>
ggplot(aes(x = cesdTotal, y = ahiTotal)) +
geom_point() +
theme_classic() +
labs(x = "CES-D", y = "AHI")
- 散布図内にサブグループがある場合に，色を塗り分けたいこともあるかもしれない。fillとcolorに群の情報をいれると塗り分けができる。なお，interventionは数値型になっているので因子型に変換している。
data |>
filter(occasion == 0) |>
mutate(intervention = as.factor(intervention)) |>
ggplot(aes(x = cesdTotal, y = ahiTotal,
fill = intervention, color = intervention)) +
geom_point() +
theme_classic() +
labs(x = "CES-D", y = "AHI")
- ３つの変数の関係をみる上では，バブルプロットも便利です。sizeに３つ目の変数を指定すれば描けます。今回はあまり意味のあるものではないので，そこから意味は読み取れないですね。
data |>
filter(occasion == 0) |>
mutate(intervention = as.factor(intervention)) |>
ggplot(aes(x = cesdTotal, y = ahiTotal, size = intervention)) +
geom_point() +
theme_classic() +
labs(x = "CES-D", y = "AHI")Warning: Using size for a discrete variable is not advised.
折れ線
時間的な変化や回帰直線などをプロットする際には，折れ線が使える。
今回は，occasionが0~5まで５時点あるので，id=1〜10までプロットしてみる。工夫としては，idで10以下のデータにしぼったこと，idは因子型にしたことになる。折れ線は，
geom_line()を使う。
data |>
filter(id <= 10) |>
mutate(id = as.factor(id)) |>
ggplot(aes(x = occasion, y = cesdTotal, color = id)) +
geom_line() 
theme_classic()などの見栄えの調整をした。geom_point()も重ねると点がわかりやすいので追加している。
data |>
filter(id <= 10) |>
mutate(id = as.factor(id)) |>
ggplot(aes(x = occasion, y = cesdTotal, color = id)) +
geom_line() +
geom_point() +
theme_classic() +
labs(x = "Time", y = "CES-D")
箱ひげ図
- 群間の比較では箱ひげ図が便利と思われる。箱ひげ図は，最小，第１四分位(下から25%)，中央値，第３四分位(下から75%)，最大値，外れ値をプロットする。
- dataの介入前の(3)“Gratitude Visit(親切してくれたけど感謝を示してない人に手紙を書く介入)”と(4)“Recording early memories(統制条件，幼いことを記憶を書く)”を比較してみる。箱ひげ図は
geom_boxplot()でプロットできる（stat_boxplot()はエラーバーを出すため）。
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention = as.factor(intervention)) |>
ggplot(aes(x = intervention, y = cesdTotal)) +
stat_boxplot(geom = "errorbar", width = 0.3) +
geom_boxplot()
- 介入名を入れたり，見栄えの調整をする。
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_boxplot(geom = "errorbar", width = 0.3) +
geom_boxplot() +
theme_classic() +
labs(x = "Intervention", y = "CES-D")
- 箱ひげ図よりもデータの分布をみたい場合は，バイオリンプロットが便利です。
geom_violin(trim = FALSE)で描けます。
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
geom_violin(trim = FALSE) +
theme_classic() +
labs(x = "Intervention", y = "CES-D")
- バイオリンプロットにさらにデータポイントを描く（ドットプロット）とよりイメージしやすくなります。
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
geom_violin(trim = FALSE) +
geom_jitter(position=position_jitter(0.1))+
theme_classic() +
labs(x = "Intervention", y = "CES-D")
レインクラウドプロット(Raincloud plot)
- 箱ひげ図，バイオリンプロット，ドットプロットを組み合わせたものとして，レインクラウドプロットがある。ggdistパッケージで描けます。まず，
stat_halfeye()でバイオリンプロットの半分みたいなものを描きます。
library(ggdist)
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye()
- ドットプロット追加します。
stat_halfeye()の信頼区間っぽいやつを削除するために，point_color = NA, .width = 0を指定します。
library(ggdist)
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye(width = .5, point_color = NA, .width = 0) +
stat_dots(side = "left", dotsize = 1, justification = 1.05, binwidth = 1)
- 箱ひげ図を追加します。
library(ggdist)
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
geom_boxplot(width = .1,outlier.shape = NA) 
- 最終調整をします。なぜか左に寄っちゃうので，
coord_cartesian(clip = 'off', xlim = c(1, 2))で調整しました。
library(ggdist)
data |>
filter(intervention == 3 | intervention == 4) |>
filter(occasion == 0) |>
mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |>
ggplot(aes(x = intervention2, y = cesdTotal)) +
stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
geom_boxplot(width = .1,outlier.shape = NA) +
theme_classic() +
labs(x = "Intervention", y = "CES-D") +
coord_cartesian(clip = 'off', xlim = c(1, 2))
- ４群の比較もしてみます。４つなのでifelse()ではなく，
case_when()を使っています。