tidyplotsの使い方　（Part3　実験データのプロットと理論曲線への当てはめ・減衰定数計算）

library(tidyverse) # データ整形・可視化のための総合パッケージ群を読み込み
library(tidyplots) # グラフ作成を簡単にするパッケージを読み込み
library(readr) # CSVファイルなどを読み込むパッケージを読み込み

#データ読み込み
sampledata <- read_csv("sampledata.csv") # CSVファイルを読み込んでsampledataという名前のオブジェクトに保存。read_csv()関数はreadrパッケージの関数

# 生データから変化量(F/F0)への標準化
# 各行の値を1~50行までの平均値で除算してF/F0を計算し、さらに−１することで変化の振幅を計算しやすいようにした。

sample_normalized <- sampledata %>% # %>%はパイプ演算子。データを次の処理に渡す
  mutate(across( # mutate()は列を変更・追加する関数。across()は複数の列に同じ処理を適用
    .cols = -1,   # .cols = -1 は1列目(Time列)以外の全ての列を処理対象とする。マイナスは「除外」の意味
    .fns = ~ .x / mean(.x[1:50]) -1  # .fns は適用する関数。~ は関数定義の記号。.x は処理対象の列。.x[1:50]は1~50行目。mean()で平均値を計算し、各値をその平均で割ってから1を引く
    ))

# tidydataに変換
sample_tidy <- sample_normalized |> # |>は新しいパイプ演算子(%>%と同じ意味)
  pivot_longer( # 横長データを縦長データに変換する関数
    cols = starts_with("Puff"), # starts_with("Puff")で"Puff"で始まる列名を全て選択
    names_to = "Puff",     # 元の列名を格納する新しい列の名前を"Puff"とする
    values_to = "FdivF0"  # 元の値を格納する新しい列の名前を"FdivF0"とする
    )

# 試行ごとにフィッティングを実行するために入れ物を作成
puff_list <- unique(sample_tidy$Puff) # unique()で重複を除いたユニークな値を取得。sample_tidy$PuffでPuff列のデータを取り出し、Puff実験の種類リストを作成
pred_table_list <- list() # 空のリストを作成。後で予測値を保存する
fit_models <- list()  # 空のリストを作成。フィッティングモデルを保存する

for (puff in puff_list) { # for文で繰り返し処理。puff_listの各要素について以下の処理を実行
  # 試行ごとのデータを抽出
  puff_data <- sample_tidy |> # パイプでデータを渡す
    filter(Puff == puff, Time > 16) # filter()で条件に合うデータだけを抽出。Puff == puffは現在処理中のPuff実験、Time > 16は16秒より後のデータのみ使用(サンプルデータをみると16秒以降のデータを使うと良さそうだったので)
  
  # 非線形最小二乗法で指数関数モデルを当てはめる
  # モデル: y = a * exp(-t/tau) + c
  
  # 初期値をデータから推定
    y_min <- min(puff_data$FdivF0, na.rm = TRUE) # min()で最小値を取得。na.rm = TRUEは欠損値(NA)を除外して計算
    y_end <- mean(tail(puff_data$FdivF0, 10), na.rm = TRUE)  # tail()でデータの最後の10行を取得し、mean()でその平均値を計算。最後の10点の平均
    a_init <- y_min - y_end  # 振幅の初期推定値(最小値 - 終値)
    c_init <- y_end  # オフセットの初期推定値(終値)
    
    # カーブフィッティング本体    
    fit_model <- nls( # nls()は非線形最小二乗法(Non-Linear least Squares)の関数。データを理論曲線に当てはめる
      FdivF0 ~ a * exp(-(Time - 16)/tau) + c, # ~は「左辺は右辺で説明される」という意味。exp()は指数関数(eの累乗)。指数関数減衰モデル: y = a × e^(-(t-16)/τ) + c
      data = puff_data, # 使用するデータを指定
      start = list(a = a_init, tau = 10, c = c_init) # startはパラメータの初期値を指定。tau = 10は時定数の初期値を10秒と仮定
      )
    
    # モデルを保存
    fit_models[[puff]] <- fit_model # フィッティング結果をリストに保存。[[puff]]はリストの要素名を指定
    
    # 予測値テーブルの作成
    time_seq <- seq(16, 90, by = 0.1) # seq()で等間隔の数列を作成。16から90まで0.1刻みの時間点を作成
    pred_values <- predict(fit_model, newdata = data.frame(Time = time_seq)) # predict()でモデルから予測値を計算。newdataで予測する時間点を指定。data.frame()でデータフレーム形式にする
    pred_table_list[[puff]] <- tibble( # tibble()でデータフレームを作成(tidyverse版)。予測値をテーブル形式でリストに保存
      Time = time_seq, # 時間列
      FdivF0 = pred_values, # 予測値列
      Puff = puff # Puff実験名列
    )
}

# 全試行の予測値を結合
pred_table <- bind_rows(pred_table_list) # bind_rows()で複数のデータフレームを縦に結合。全Puff実験の予測値を1つのテーブルにまとめる

sample_tidy |> # パイプでデータを渡す
  tidyplot(x = Time, y = FdivF0, color=Puff) |> # tidyplot()でグラフ作成開始。x, yで横軸と縦軸を指定。colorで色分けする変数を指定
  
  # 実データを点で表示
  add_data_points(alpha=0.5) |> # add_data_points()で散布図として点を追加。alpha=0.5で透明度を50%に設定

  # X軸Y軸の範囲と目盛の設定
  adjust_x_axis(limits = c(10, 60), breaks = seq(10, 60, 5)) |> # adjust_x_axis()でX軸の設定を調整。limitsで軸の範囲を10~60に設定。breaksで目盛りを10, 15, 20...60に設定
  adjust_y_axis(limits = c(-0.7, 0.1)) |> # Y軸の範囲を-0.7~0.1に設定
  add_reference_lines(x=15, y=0) |> # add_reference_lines()で参照線を追加。x=15とy=0の位置に線を引く
  
  # タイトルや軸名
  adjust_title(fontsize=14, "実測値と指数関数減衰曲線") |> # グラフのタイトルを設定。フォントサイズを14ポイントに
  adjust_x_axis_title("時間(秒)") |> # X軸ラベルを設定
  adjust_y_axis_title("変化量\n(投与前を０に標準化)") |> # Y軸ラベルを設定。\nは改行を表す
  
  # グラフの大きさやデザインテーマの設定
  adjust_size(width = 14, height = 7, unit = "cm") |> # グラフのサイズを幅14cm、高さ7cmに設定
  # theme_ggplot2() |> # コメントアウト(実行されない)。必要に応じてテーマを変更可能
  
  # 各試行のフィッティング曲線を赤い破線で重ね書き
  # この機能はtidyplotsにないので大本のggplot2の機能を利用
  add(ggplot2::geom_line( # add()でggplot2の機能を追加。geom_line()で線グラフを追加
    data = pred_table, # 使用するデータ(予測値テーブル)を指定
    aes(x = Time, y = FdivF0, color = Puff), # aes()は美的マッピング(どのデータを何に対応させるか)
    linewidth = 0.5, # 線の太さを0.5に設定
    alpha = 1, # 透明度を1(不透明)に設定
    inherit.aes = FALSE)) |> # inherit.aes = FALSEで親のaes設定を継承しない
  
  # 減衰曲線の注釈を記載(tidyplotsの機能に戻って)
  add_annotation_text( # add_annotation_text()でグラフに注釈テキストを追加
    x = 40, y = -0.25, # テキストの位置(x座標40、y座標-0.25)
    size=12, # サイズ指定(この引数は実際には使われていない可能性あり)
    "減衰曲線\n y = a * exp(-t-16/tau) + c", # 表示するテキスト内容。\nで改行
    fontsize = 12  ) # フォントサイズを12ポイントに設定

# フィッティング係数を格納するデータフレームを作成
coefficients_table <- tibble( # tibble()で空のデータフレームを作成
  Puff = character(), # character()で文字列型の空ベクトルを作成
  a = numeric(), # numeric()で数値型の空ベクトルを作成(振幅)
  tau = numeric(), # 時定数
  c = numeric(), # オフセット
  a_se = numeric(), # aの標準誤差
  tau_se = numeric(), # tauの標準誤差
  c_se = numeric(), # cの標準誤差
  residual_se = numeric() # 残差の標準誤差
)

for (puff in names(fit_models)) { # names()でリストの要素名を取得。保存した全てのモデルについてループ処理
  # 係数と標準誤差を取得
  coefs <- coef(fit_models[[puff]]) # coef()でモデルの係数(パラメータ推定値)を取得
  summ <- summary(fit_models[[puff]]) # summary()でモデルの詳細な統計情報を取得
  se <- summ$coefficients[, "Std. Error"] # 標準誤差を取り出す。[, "Std. Error"]で"Std. Error"列を取得
  res_se <- summ$sigma # sigmaは残差の標準誤差(モデルの当てはまりの良さを表す)
  
  # テーブルに追加
  coefficients_table <- coefficients_table |> # パイプでデータを渡す
    add_row( # add_row()でデータフレームに新しい行を追加
      Puff = puff, # Puff実験名
      a = coefs["a"], # 振幅の推定値
      tau = coefs["tau"], # 時定数の推定値
      c = coefs["c"], # オフセットの推定値
      a_se = se["a"], # 振幅の標準誤差
      tau_se = se["tau"], # 時定数の標準誤差
      c_se = se["c"], # オフセットの標準誤差
      residual_se = res_se # 残差の標準誤差
    )
}

# 係数一覧表を表示
cat("\n=== フィッティング係数一覧表 ===\n") # cat()でテキストを出力。\nは改行
print(coefficients_table) # print()でオブジェクトを表示

# 平均値と標準偏差を計算(小数２桁)
cat("\n=== 係数の平均値 ± 標準偏差 ===\n", # cat()で整形したテキストを出力
    "a (振幅):", # ラベル
    round(mean(coefficients_table$a),2),"±",round(sd(coefficients_table$a),2),"\n", # mean()で平均値、sd()で標準偏差を計算。round(..., 2)で小数点以下2桁に丸める
    
    "tau (減衰率):", 
    round(mean(coefficients_table$tau),2),"±",round(sd(coefficients_table$tau),2),"\n", # tauの平均±標準偏差
    
    "c (オフセット):", # ラベル
    round(mean(coefficients_table$c),2), "±",round(sd(coefficients_table$c),2), "\n") # cの平均±標準偏差

## === フィッティング係数一覧表 ===
## # A tibble: 6 × 8
##   Puff       a   tau       c    a_se tau_se     c_se residual_se
##   <chr>  <dbl> <dbl>   <dbl>   <dbl>  <dbl>    <dbl>       <dbl>
## 1 Puff1 -0.584  2.71 -0.0448 0.00315 0.0211 0.000313     0.00454
## 2 Puff2 -0.570  4.03 -0.0557 0.00298 0.0313 0.000383     0.00531
## 3 Puff3 -0.431  4.58 -0.0745 0.00165 0.0264 0.000231     0.00314
## 4 Puff4 -0.453  6.79 -0.0683 0.00178 0.0426 0.000334     0.00417
## 5 Puff5 -0.413  7.64 -0.0995 0.00121 0.0365 0.000250     0.00301
## 6 Puff6 -0.370  2.38 -0.0631 0.0106  0.0976 0.000972     0.0142

## 
## === 係数の平均値 ± 標準偏差 ===
##  a (振幅): -0.47 ± 0.09 
##  tau (減衰率): 4.69 ± 2.14 
##  c (オフセット): -0.07 ± 0.02