LEC07: グラフの実践3（ggplot2）
Kenri Kodaka
2025-11-11
[準備1]
データフレームをインポート
◯◯◯ 使用するライブラリ
ggpplot2ライブラリ（グラフ描画パッケージ）およびdplyr（データ操作）を使用するために、tidyverseパッケージを、また、複数グラフ配置のためにpatchworkライブラリを
install.packages("tidyverse")
install.packages("patchwork")library(tidyverse)
library(patchwork)
# 以下はtidyverseに含まれます。
# library(ggplot2)
# library(dplyr)◯◯◯ ブッダの耳錯覚
今回の演習では、 「ブッダの耳錯覚」のデータを扱います。
「ブッダの耳錯覚」については以下の映像等を参考にしてください。
【論文発表】耳がびよ〜んと伸びる皮膚変形錯覚『ブッダの耳錯覚』の実験が、国際論文誌「i-Perception」に掲載されました（7.23）。この錯覚は誰か一人の相手がいれば、あとは何の道具もいりません。絵にあるように、一方の手で耳たぶを引っ張るのと同時に、もう一方の手で見えない皮膚の伸びをパント… pic.twitter.com/3L3mz9a39w
— 「」kenrikodaka.as (@kenrikodaka) August 3, 2024
今回扱う実験では、ブッダの耳錯覚において、耳の上を自分で掴むこと（セルフアンカリング）による効果を検証しています。
ブッダの耳錯覚におけるセルフアンカリング
- [論文PDF]｜耳介上部をつまむとブッダの耳錯覚の効果が強化される, 2025年度日本認知科学会第42回大会
◯◯◯ データフレームをインポート
まずは以下のCSVを読み込みましょう。
url = "https://lab.kenrikodaka.com/_download/csv/buddhaexp2024.csv"
source = read.csv(url)最初の6行を確認します。
head(source)## SBJ Q1T0B0 Q1T0B1 Q1T1B0 Q1T1B1 Q2T0B0 Q2T0B1 Q2T1B0 Q2T1B1 Q3T0B0 Q3T0B1 ## 1 1 1.0 3.5 2.5 5.5 0.0 2.5 1.5 3.0 0.0 2.5 ## 2 2 1.5 4.0 1.5 4.5 2.5 3.0 1.5 4.5 0.5 1.5 ## 3 3 0.0 5.0 0.0 6.0 0.0 3.5 0.0 4.5 0.0 0.0 ## 4 4 1.0 3.0 3.0 3.5 1.5 0.0 1.5 3.0 0.5 1.5 ## 5 5 2.0 3.0 4.0 5.5 2.0 4.0 4.0 5.5 0.0 0.0 ## 6 6 0.0 1.0 0.0 1.5 0.0 0.0 0.0 0.0 0.0 0.0 ## Q3T1B0 Q3T1B1 Q4T0B0 Q4T0B1 Q4T1B0 Q4T1B1 E0T0B0 E0T0B1 E0T1B0 E0T1B1 E1T0B0 ## 1 1.5 3 0.0 0.0 0.0 0.0 4.25 15.45 20.10 1.70 36.305 ## 2 0.0 0 0.0 0.0 0.0 0.0 13.15 14.15 11.10 11.10 16.550 ## 3 0.0 0 0.0 0.0 0.0 0.0 -0.25 3.05 1.10 0.20 4.550 ## 4 3.0 3 2.5 2.0 2.5 2.0 3.60 8.20 -1.50 -0.20 5.600 ## 5 0.5 0 0.0 0.5 1.5 0.5 -2.10 -1.30 1.40 3.30 -1.950 ## 6 0.0 0 0.0 0.0 0.0 0.0 -2.00 -0.60 4.15 0.75 3.450 ## E1T0B1 E1T1B0 E1T1B1 ORDER ## 1 49.005 32.80 50.20 0 ## 2 54.950 41.35 46.35 1 ## 3 55.450 4.80 77.30 0 ## 4 9.100 5.20 6.80 0 ## 5 40.850 30.00 75.00 1 ## 6 1.950 0.40 2.80 0
◯◯◯ 変数の内容
以下で、各変数について解説します。
- SBJ: 被験者のID（1-16）｜16人分の実験データが含まれています。
nrow(source)## [1] 16
変数名[XX][TaBb]の説明
[XX] Q1-Q4はアンケートの設問に対する主観評価値（0-6の7段階）
（主観評価実験）
- Q1:「自分の耳たぶが通常よりも伸びている感じがした。」
- Q2:「何も無い空間に「見えない」皮膚が存在しているように感じた。」
- Q3:「自分の耳全体が下に移動しているような感じがした。」
- Q4:「自分の耳全体の面積が小さくなる感じがした。」
（行動実験）
- E0:耳上端の主観位置の錯覚前後の降下量（cm）
- E1:耳下端の主観位置の錯覚前後の降下量（cm）
[TaBb] Taは耳上部、Bbは耳たぶ側の条件
- T0B0:（体験者が）耳上端に触れない x （実験者が）耳下端をつまむ（だけ）
- T0B1:（体験者が）耳上端に触れない x （実験者が）耳下端の下を引っ張る
- T1B0:（体験者が）耳上端の上をつまむ x （実験者が）耳下端をつまむ（だけ）
- T1B1:（体験者が）耳上端の上をつまむ x （実験者が）耳下端の下を引っ張る
ORDER 行動実験の順序（今回の演習では、ORDERは無視します）
- 0：T1B1→T0B1（錯覚条件に関してT1が先、その後にT0）
- 1：T0B1→T1B1（錯覚条件に関してT0が先、その後にT1）
◯◯◯ データフレームの整形
ひとまず
sourceをdat_にコピーします。
dat_ = sourceまずはアンケート実験の結果を収納するデータフレーム
dat_Qを作成します。 以下の変数を追加することを想定しています。
- rate（評価値）
- stmt（質問項目）｜Q1〜Q4
- etop（耳上端条件）｜T0：触れない、T1：つまむ
- ebot（耳下端条件）｜B0：つまむ（だけ）、B1：引っ張る
dat_Q = dat_ %>%
pivot_longer(
cols = 2:17, #2行目から17行までが対象
names_to = "type",
values_to = "rate"
) %>%
select(SBJ,type,rate)
dat_Q #途中経過（typeは後に削除）## # A tibble: 256 × 3 ## SBJ type rate ## <int> <chr> <dbl> ## 1 1 Q1T0B0 1 ## 2 1 Q1T0B1 3.5 ## 3 1 Q1T1B0 2.5 ## 4 1 Q1T1B1 5.5 ## 5 1 Q2T0B0 0 ## 6 1 Q2T0B1 2.5 ## 7 1 Q2T1B0 1.5 ## 8 1 Q2T1B1 3 ## 9 1 Q3T0B0 0 ## 10 1 Q3T0B1 2.5 ## # ℹ 246 more rows
#typeに含まれる文字に応じて、変数の値を決めます。
## str_detect(A,"str")は、文字列Aに"str"が含まれる時TRUEを返す。
dat_Q = dat_Q %>%
mutate(
stmt = case_when(
str_detect(type,"Q1") ~ "Q1",
str_detect(type,"Q2") ~ "Q2",
str_detect(type,"Q3") ~ "Q3",
str_detect(type,"Q4") ~ "Q4"
),
etop = case_when(
str_detect(type,"T0") ~ "T0",
str_detect(type,"T1") ~ "T1"
),
ebot = case_when(
str_detect(type,"B0") ~ "B0",
str_detect(type,"B1") ~ "B1"
)
) %>%
select(-type)
dat_Q## # A tibble: 256 × 5 ## SBJ rate stmt etop ebot ## <int> <dbl> <chr> <chr> <chr> ## 1 1 1 Q1 T0 B0 ## 2 1 3.5 Q1 T0 B1 ## 3 1 2.5 Q1 T1 B0 ## 4 1 5.5 Q1 T1 B1 ## 5 1 0 Q2 T0 B0 ## 6 1 2.5 Q2 T0 B1 ## 7 1 1.5 Q2 T1 B0 ## 8 1 3 Q2 T1 B1 ## 9 1 0 Q3 T0 B0 ## 10 1 2.5 Q3 T0 B1 ## # ℹ 246 more rows
ひきつづき、行動実験の結果を収納するデータフレーム
dat_Dを作成します。 以下の変数を追加することを想定しています。
- drift（位置変化量）｜錯覚前後で位置がどれだけ変化したか（cm）
- epos（変化部位）｜TOP：耳上端、BOT：耳たぶ
- etop（耳上端条件）｜T0：触れない、T1：つまむ
- ebot（耳下端条件）｜B0：つまむ（だけ）、B1：引っ張る
dat_D = dat_ %>%
pivot_longer(
cols = 18:25, #2行目から17行までが対象
names_to = "type",
values_to = "drift"
) %>%
select(SBJ,type,drift)
dat_D #途中経過（typeは後に削除）## # A tibble: 128 × 3 ## SBJ type drift ## <int> <chr> <dbl> ## 1 1 E0T0B0 4.25 ## 2 1 E0T0B1 15.4 ## 3 1 E0T1B0 20.1 ## 4 1 E0T1B1 1.7 ## 5 1 E1T0B0 36.3 ## 6 1 E1T0B1 49.0 ## 7 1 E1T1B0 32.8 ## 8 1 E1T1B1 50.2 ## 9 2 E0T0B0 13.2 ## 10 2 E0T0B1 14.2 ## # ℹ 118 more rows
#typeに含まれる文字に応じて、変数の値を決めます。
dat_D = dat_D %>%
mutate(
epos = case_when(
str_detect(type,"E0") ~ "TOP",
str_detect(type,"E1") ~ "BOT"
),
etop = case_when(
str_detect(type,"T0") ~ "T0",
str_detect(type,"T1") ~ "T1"
),
ebot = case_when(
str_detect(type,"B0") ~ "B0",
str_detect(type,"B1") ~ "B1"
)
) %>%
select(-type)
dat_D## # A tibble: 128 × 5 ## SBJ drift epos etop ebot ## <int> <dbl> <chr> <chr> <chr> ## 1 1 4.25 TOP T0 B0 ## 2 1 15.4 TOP T0 B1 ## 3 1 20.1 TOP T1 B0 ## 4 1 1.7 TOP T1 B1 ## 5 1 36.3 BOT T0 B0 ## 6 1 49.0 BOT T0 B1 ## 7 1 32.8 BOT T1 B0 ## 8 1 50.2 BOT T1 B1 ## 9 2 13.2 TOP T0 B0 ## 10 2 14.2 TOP T0 B1 ## # ℹ 118 more rows
[Graph1_QuestionnaireBoxplot]
アンケート結果をボックスプロットでまとめる
# [ggplot]
## interactionは、複数の変数の組み合わせをグラフの属性にマップできる
## interaction(etop,ebot)は「T0.B0」「T1.B0」「T0.B1」「T1.B1」の4つの要素を作る
gp = ggplot(dat_Q,aes(x=interaction(etop,ebot),
y=rate,fill=interaction(etop,ebot)))
# [geom_boxplot]ボックスプロットによる描画
gp = gp + geom_boxplot()
gp = gp + scale_x_discrete(
labels=c("T0B0","T1B0","T0B1","T1B1"))
gp = gp + scale_y_continuous(limits=c(0,6),breaks=0:6)
gp = gp + scale_fill_discrete(
labels=c("T0 x B0","T1 x B0","T0 x B1","T1 x B1"))
gp = gp + labs(x = "CONDITION", y = "RATE", fill="CONDITION")
gp = gp + facet_grid(. ~ stmt)gp + theme_minimal() +
ggtitle("Graph1_QuestionnaireBoxplot")
[Graph2_QuestionnaireBar]
アンケート結果を棒グラフでまとめる
まずは、
summarizeを使って、個々の条件の平均値・標準誤差を算出します。 エラーバーに使用される「標準誤差」は、標準偏差をサンプル数の平方根で割ることで計算されます。
dat_stat = dat_Q %>%
group_by(stmt,etop,ebot) %>%
summarize(
avg = mean(rate), #平均
sd = sd(rate), #標準偏差
n = n(),#サンプル数
se = sd / sqrt(n) #標準誤差
) %>%
select(stmt,etop,ebot,avg,se)
dat_stat## # A tibble: 16 × 5 ## # Groups: stmt, etop [8] ## stmt etop ebot avg se ## <chr> <chr> <chr> <dbl> <dbl> ## 1 Q1 T0 B0 1.28 0.341 ## 2 Q1 T0 B1 3.28 0.382 ## 3 Q1 T1 B0 1.69 0.433 ## 4 Q1 T1 B1 4.59 0.317 ## 5 Q2 T0 B0 1.25 0.351 ## 6 Q2 T0 B1 2.47 0.402 ## 7 Q2 T1 B0 1.53 0.378 ## 8 Q2 T1 B1 3.53 0.402 ## 9 Q3 T0 B0 0.344 0.156 ## 10 Q3 T0 B1 1.09 0.286 ## 11 Q3 T1 B0 0.469 0.216 ## 12 Q3 T1 B1 0.844 0.312 ## 13 Q4 T0 B0 0.25 0.177 ## 14 Q4 T0 B1 0.281 0.151 ## 15 Q4 T1 B0 0.281 0.177 ## 16 Q4 T1 B1 0.219 0.137
◯◯◯ Graph2A（エラーバーなし）
ここでは、エラーバーのない棒グラフを作成します。
gp_a = ggplot(dat_stat,aes(x=interaction(etop,ebot), y=avg,
fill=interaction(etop,ebot)))
# geom_barは数え上げ、geom_colは数値をそのまま表示することに注意
## position_dodgeは横に並べるの意味
gp_a = gp_a + geom_col(position=position_dodge(width=0.9))
gp_a = gp_a + scale_x_discrete(labels=c("T0B0","T1B0","T0B1","T1B1"))
gp_a = gp_a + scale_y_continuous(limits=c(0,6),breaks=0:6)
gp_a = gp_a + scale_fill_discrete(
labels=c("T0 x B0","T1 x B0","T0 x B1","T1 x B1"))
gp_a = gp_a + facet_grid(. ~ stmt)
gp_a = gp_a + labs(x="CONDITION",y="RATE",fill="CONDITION")gp_a + theme_minimal() +
ggtitle("Graph2A_QuestionnaireBar")
◯◯◯ Graph2B（エラーバーあり）
geom_errorbarで、エラーバーを追加します。
gp_b = ggplot(dat_stat,aes(x=interaction(etop,ebot), y=avg,
fill=interaction(etop,ebot)))
gp_b = gp_b + geom_col(position=position_dodge(width=0.9))
#[geom_errobar]
## aes内のyminはエラーバーの下限ymaxは上限に対応
## yminは「平均値-標準誤差」、ymaxは「平均値+標準誤差」に対応
## widthは横幅（1が棒グラフの横幅いっぱい）
gp_b = gp_b + geom_errorbar(aes(ymin = avg-se, ymax = avg+se), width=0.3)
gp_b = gp_b + scale_x_discrete(labels=c("T0B0","T1B0","T0B1","T1B1"))
gp_b = gp_b + scale_y_continuous(limits=c(0,6),breaks=0:6)
gp_b = gp_b + scale_fill_discrete(
labels=c("T0 x B0","T1 x B0","T0 x B1","T1 x B1"))
gp_b = gp_b + facet_grid(. ~ stmt)
gp_b = gp_b + labs(x="CONDITION",y="RATE",fill="CONDITION")gp_b + theme_minimal() +
ggtitle("Graph2B_QuestionnaireBar")
[Graph3_EarPositionDrift]
錯覚前後の耳（上端・下端）の主観位置の変化
まずはdat_Dの構造を復習します。
- drift：錯覚前後の位置の変化（cm）
- epos：耳が上端（TOP）か、下端（BOT）か
- etop：耳の上を被験者がつままない (T0) or つまむ（T1）
- ebot：耳の下を実験者がつまむ（B0） or 引っ張る（B1）
dat_D## # A tibble: 128 × 5 ## SBJ drift epos etop ebot ## <int> <dbl> <chr> <chr> <chr> ## 1 1 4.25 TOP T0 B0 ## 2 1 15.4 TOP T0 B1 ## 3 1 20.1 TOP T1 B0 ## 4 1 1.7 TOP T1 B1 ## 5 1 36.3 BOT T0 B0 ## 6 1 49.0 BOT T0 B1 ## 7 1 32.8 BOT T1 B0 ## 8 1 50.2 BOT T1 B1 ## 9 2 13.2 TOP T0 B0 ## 10 2 14.2 TOP T0 B1 ## # ℹ 118 more rows
グラフ化の前準備として、eposの順番を左からTOP→BOTとします。
dat_D$epos = factor(dat_D$epos,levels=c("TOP","BOT"))# [ggplot]
gp = ggplot(dat_D,aes(x=interaction(etop,ebot),
y=drift,fill=epos))
# [geom_hline] y=0とy=50に横線を引く。solidは実線、dottedは点線
gp = gp + geom_hline(yintercept = 0, linetype="solid",linewidth=0.5)
gp = gp + geom_hline(yintercept = 50, linetype="dotted",linewidth=0.55)
# [geom_boxplot]ボックスプロットによる描画
gp = gp + geom_boxplot()
gp = gp + scale_x_discrete(
limits=c("T0.B0","T1.B0","T0.B1","T1.B1"),
labels=c("T0B0","T1B0","T0B1","T1B1"))
gp = gp + scale_y_continuous(limits=c(80,-5),
breaks=10*(8:-1),trans="reverse")
# 今回はbrewerパレットを使います（scale_fill_discreteと色の割り当てが異なる）。
# RColorBrewer::display.brewer.all() #色のパレットを見て、pairedを選ぶ。
gp = gp + scale_fill_brewer(labels=c("TOP","BOTTOM"),palette="Paired")
gp = gp + labs(x = "CONDITION",
y = "SUBJECTIVE EAR-LOCATION'S DRIFT(cm)",
fill = "EAR'S POSITION")gp + theme_gray() +
ggtitle("Graph3_EarPositionDrift")
[Graph4_EarSizeDeformation]
錯覚前後の耳のサイズの変化
最初に
dat_Dの構造を確認。
dat_D## # A tibble: 128 × 5 ## SBJ drift epos etop ebot ## <int> <dbl> <fct> <chr> <chr> ## 1 1 4.25 TOP T0 B0 ## 2 1 15.4 TOP T0 B1 ## 3 1 20.1 TOP T1 B0 ## 4 1 1.7 TOP T1 B1 ## 5 1 36.3 BOT T0 B0 ## 6 1 49.0 BOT T0 B1 ## 7 1 32.8 BOT T1 B0 ## 8 1 50.2 BOT T1 B1 ## 9 2 13.2 TOP T0 B0 ## 10 2 14.2 TOP T0 B1 ## # ℹ 118 more rows
耳のサイズの（下方向）変化（dsize）はebotからetopを引いたものです。 以下で、dsize = ebot - etopとして、新たな属性dsizeを追加しdat_D2とします。
この際、
mutateで計算しやすいように、あらかじめ、二つあるeposの測定値（drift）を単一の行に戻します。これはpivot_longerと逆の作用となるので、pivot_widerを使います。
dat_D2 = dat_D %>%
pivot_wider(
names_from = epos, #新たに列の変数にするもの
values_from = drift, #対象となる値
names_glue = "drift_{epos}" #列名の先頭にdrift_をつける。（TOP、BOTがデフォルト）
)
dat_D2 #途中経過## # A tibble: 64 × 5 ## SBJ etop ebot drift_TOP drift_BOT ## <int> <chr> <chr> <dbl> <dbl> ## 1 1 T0 B0 4.25 36.3 ## 2 1 T0 B1 15.4 49.0 ## 3 1 T1 B0 20.1 32.8 ## 4 1 T1 B1 1.7 50.2 ## 5 2 T0 B0 13.2 16.6 ## 6 2 T0 B1 14.2 55.0 ## 7 2 T1 B0 11.1 41.4 ## 8 2 T1 B1 11.1 46.4 ## 9 3 T0 B0 -0.25 4.55 ## 10 3 T0 B1 3.05 55.4 ## # ℹ 54 more rows
ここまでくれば、
mutateを用いた行ごとの変数間の計算で、変形量dsizeを簡単に算出することができます。
dat_D2 = dat_D2 %>%
mutate(dsize = drift_BOT - drift_TOP) %>%
select(-drift_TOP,-drift_BOT) #個々の変化量は捨てます。
dat_D2## # A tibble: 64 × 4 ## SBJ etop ebot dsize ## <int> <chr> <chr> <dbl> ## 1 1 T0 B0 32.1 ## 2 1 T0 B1 33.6 ## 3 1 T1 B0 12.7 ## 4 1 T1 B1 48.5 ## 5 2 T0 B0 3.4 ## 6 2 T0 B1 40.8 ## 7 2 T1 B0 30.2 ## 8 2 T1 B1 35.2 ## 9 3 T0 B0 4.8 ## 10 3 T0 B1 52.4 ## # ℹ 54 more rows
# [ggplot]
gp = ggplot(dat_D2,aes(x=interaction(etop,ebot),
y=dsize,fill=interaction(etop,ebot)))
# [geom_boxplot]ボックスプロットによる描画
gp = gp + geom_boxplot()
gp = gp + scale_x_discrete(
labels=c("T0B0","T1B0","T0B1","T1B1"))
gp = gp + scale_y_continuous(limits=c(0,100),breaks=10*(0:10))
gp = gp + scale_fill_discrete(
labels=c("T0 x B0","T1 x B0","T0 x B1","T1 x B1"))
gp = gp + labs(x = "CONDITION", y = "DEFORMATION SIZE (cm)",
fill = "CONDITION")gp + theme_minimal() +
ggtitle("Graph4_EarSizeDeformation")
[確認課題]
ApMedia07_Work
- Graph4_EarSizeDeformationをエラーバー付きの棒グラフにしてください。
- 締切は11月中とします. -ファイル名は「2350xx_work7.R」としてください。