# 必要なパッケージの読み込み
library(ggplot2) # グラフ作成用
## Warning: パッケージ 'ggplot2' はバージョン 4.4.2 の R の下で造られました
library(RColorBrewer) # カラーパレット作成用
library(dplyr) # データ処理用
## Warning: パッケージ 'dplyr' はバージョン 4.4.2 の R の下で造られました
##
## 次のパッケージを付け加えます: 'dplyr'
## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## filter, lag
## 以下のオブジェクトは 'package:base' からマスクされています:
##
## intersect, setdiff, setequal, union
library(scales) # スケール調整用
library(ggpubr) # 公開用のプロット補助
library(readxl) # Excelファイルの読み込み用
library(car) # 検定（Levene検定）用
## 要求されたパッケージ carData をロード中です
##
## 次のパッケージを付け加えます: 'car'
## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## recode
library(rstudioapi) # RStudio環境でのファイル操作用
library(tools) # ファイル名操作用
# Rスクリプトのファイル名を取得
if (rstudioapi::isAvailable()) {
script_name <- basename(rstudioapi::getActiveDocumentContext()$path) # 現在のスクリプトのファイル名
} else {
script_name <- "interactive_session" # RStudio以外の場合はデフォルト名
}
# データの読み込み
file_path <- "C:/Users/ndsho/Desktop/241113_cocktailPaper/241116_NodaCocktailData_2D.xlsx"
data <- read_excel(file_path, sheet = "2D_MYH+") # 指定したシートからデータを読み込み
# 'Type' を因子型に変換し、レベルの順序を指定
data$Type <- factor(data$Type, levels = c("HS", "B27", "NODA")) # グループの順序を指定
# グループごとの平均と標準誤差の計算
summary_data <- data %>%
group_by(Type) %>%
summarise(
AVE = mean(`MYHarea%`, na.rm = TRUE), # 平均値
SE = sd(`MYHarea%`, na.rm = TRUE) / sqrt(n()) # 標準誤差
)
# 分散分析 (ANOVA) を実施
anova_result <- aov(`MYHarea%` ~ Type, data = data) # 分散分析モデルを構築
anova_summary <- summary(anova_result) # 分散分析の結果を取得
# Tukey HSD 検定を実施
tukey_result <- TukeyHSD(anova_result) # Tukeyの多重比較検定を実行
# Tukeyの結果をデータフレームに変換
tukey_df <- as.data.frame(tukey_result$Type) %>%
mutate(
comparison = rownames(.), # 比較の名前を取得
x1 = sapply(strsplit(comparison, "-"), function(x) which(levels(data$Type) == x[1])), # 比較対象1
x2 = sapply(strsplit(comparison, "-"), function(x) which(levels(data$Type) == x[2])), # 比較対象2
length = abs(x2 - x1) # 水平線の長さを計算
) %>%
filter(`p adj` < 1) %>% # 有意な比較のみフィルタリング
arrange(desc(length)) %>% # 水平線の長さで降順に並び替え
group_by(length) %>%
mutate(height_group = cur_group_id()) %>% # 同じ長さの比較を同じ高さに配置
ungroup() %>%
mutate(offset = (2 - height_group) * -0.02) # 高いグループほどオフセットを大きく設定
# カスタムテーマを定義
my_custom_theme <- function() {
theme_minimal(base_family = "sans") + # 最小限のテーマを基にカスタマイズ
theme(
text = element_text(size = 12), # テキストサイズ
plot.title = element_text(size = 16, face = "bold", hjust = 0.5), # タイトルの書式
axis.title = element_text(size = 14, face = "bold"), # 軸タイトルの書式
axis.text = element_text(size = 12), # 軸ラベルの書式
axis.line = element_line(size = 0.8, color = "black"), # 軸線の設定
axis.ticks = element_line(color = "black", size = 0.6), # 軸目盛りの設定
axis.ticks.length = unit(-0.25, "cm"), # 軸目盛りの長さ
panel.grid = element_blank(), # グリッド線を非表示
panel.background = element_rect(fill = "transparent", color = NA), # 背景を透明に
plot.background = element_rect(fill = "transparent", color = NA), # プロット背景を透明に
legend.position = "none" # 凡例を非表示
)
}
# プロット作成
my_colors <- c(
"HS" = brewer.pal(9, "Greys")[4], # HSの色設定
"B27" = brewer.pal(9, "Blues")[6], # B27の色設定
"NODA" = brewer.pal(9, "Reds")[8] # Nodaの色設定
)
# プロット作成
p <- ggplot(data, aes(x = Type, y = `MYHarea%`, fill = Type)) +
geom_bar(data = summary_data, aes(x = Type, y = AVE, fill = Type),
stat = "identity", width = 0.6, color = "black", size = 1, inherit.aes = FALSE) + # 平均値のバー
geom_errorbar(data = summary_data, aes(x = Type, ymin = AVE - SE, ymax = AVE + SE),
width = 0.2, color = "black", size = 0.8, inherit.aes = FALSE) + # エラーバー
geom_jitter(shape = 21, color = "black", width = 0.1, size = 3, alpha = 0.8) + # 個別データポイント
labs(x = "", y = "MYH + area [%]", title = script_name) + # ラベルとタイトル
scale_fill_manual(values = my_colors) + # カラー設定
scale_y_continuous(expand = c(0, 0), limits = c(0, 80)) + # y軸を0から始め、上限はデータに自動調整
my_custom_theme() # カスタムテーマ適用
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# 線とp値を追加
y_max <- max(summary_data$AVE + summary_data$SE) * 1.05 # 棒グラフ上部の空間を確保
height_unit <- 6 # 各比較の高さ間隔
mergin <- 5 # バーと線の間隔
for (i in seq_len(nrow(tukey_df))) {
group1 <- unlist(strsplit(tukey_df$comparison[i], "-"))[1] # 比較の1つ目のグループ
group2 <- unlist(strsplit(tukey_df$comparison[i], "-"))[2] # 比較の2つ目のグループ
offset <- tukey_df$offset[i] # 個別オフセット
height <- y_max + tukey_df$height_group[i] * height_unit # 各比較の高さ計算
x1 <- tukey_df$x1[i] # 左側のx座標
x2 <- tukey_df$x2[i] # 右側のx座標
# 垂直線と水平線を描画
p <- p +
annotate("segment", x = x1 + offset, xend = x1 + offset,
y = summary_data$AVE[x1] + summary_data$SE[x1] + mergin, yend = height, size = 0.8) +
annotate("segment", x = x2 - offset, xend = x2 - offset,
y = summary_data$AVE[x2] + summary_data$SE[x2] + mergin, yend = height, size = 0.8) +
annotate("segment", x = x1 + offset, xend = x2 - offset,
y = height, yend = height, size = 0.8) +
annotate("text", x = mean(c(x1 + offset, x2 - offset)),
y = height + 2, label = paste0("p = ", format(tukey_df$`p adj`[i], digits = 2)), size = 4)
}
# プロットの表示
print(p)
# 検定結果の出力
cat("ANOVA結果:\n")
## ANOVA結果:
print(anova_summary)
## Df Sum Sq Mean Sq F value Pr(>F)
## Type 2 4349 2174.5 38.95 1.15e-06 ***
## Residuals 15 837 55.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nTukey HSD検定結果:\n")
##
## Tukey HSD検定結果:
print(tukey_result)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = `MYHarea%` ~ Type, data = data)
##
## $Type
## diff lwr upr p adj
## B27-HS 29.463667 18.258738 40.66860 0.0000161
## NODA-HS 35.616333 24.411404 46.82126 0.0000017
## NODA-B27 6.152667 -5.052262 17.35760 0.3530325
# 正規性の確認
shapiro_results <- data %>%
group_by(Type) %>%
summarise(p_value = shapiro.test(`MYHarea%`)$p.value)
cat("\nShapiro-Wilk検定結果:\n")
##
## Shapiro-Wilk検定結果:
print(shapiro_results)
## # A tibble: 3 × 2
## Type p_value
## <fct> <dbl>
## 1 HS 0.0783
## 2 B27 0.104
## 3 NODA 0.918
# 等分散性の確認
levene_result <- leveneTest(`MYHarea%` ~ Type, data = data)
bartlett_result <- bartlett.test(`MYHarea%` ~ Type, data = data)
cat("\nLevene検定結果:\n")
##
## Levene検定結果:
print(levene_result)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.9761 0.3995
## 15
cat("\nBartlett検定結果:\n")
##
## Bartlett検定結果:
print(bartlett_result)
##
## Bartlett test of homogeneity of variances
##
## data: MYHarea% by Type
## Bartlett's K-squared = 2.2819, df = 2, p-value = 0.3195
# 現在のタイムスタンプを取得
timestamp <- format(Sys.time(), "%Y%m%d_%H%M%S") # フォーマット: YYYYMMDD_HHMMSS
# ファイル名を生成（拡張子を除去し、タイムスタンプを追加）
output_path <- paste0(
"C:/Users/ndsho/Documents/Rstudio/Figures/",
file_path_sans_ext(script_name), # ファイル名の拡張子を除去
"_", timestamp, ".png"
)
# プロットを画像として保存
ggsave(
filename = output_path, # 動的に生成したファイル名
plot = p, # 保存するプロット
width = 5, height = 4, # サイズ (インチ単位)
dpi = 300 # 解像度
)
cat("プロットが保存されました:", output_path)
## プロットが保存されました: C:/Users/ndsho/Documents/Rstudio/Figures/interactive_session_20241117_165614.png