In ykunisato/somatori_kunisato_2019_replication_study: Code and data of somatori kunisato 2019 replication study

はじめに

本コードは, 杣取・国里(2019)の「アンへドニア(anhedonia)と遅延割引: Lempert & Pizzagalli (2010)の追試」で用いられている統計解析を行うためのものである。 匿名化処理後のデータ("data_anonymized")を読み込み, 以下の解析を行う。

• 割引率の推定

• Lempert & Pizzagalli (2010)と同じ相関分析

• AICに基づいたモデル比較および相関分析

• R2を用いたモデルの当てはまりの検討

• 可視化によるモデルの当てはまりの検討

• AUCを用いた相関分析

年齢は参加者を特定できてしまう可能性があるため, 本コードおよび匿名化処理後のデータでは層別化したものを用いる。 本文中に掲載されている図表は, 本コードによって作成したものを文字サイズ等加工して掲載されている。そのため, 本コードの出力と多少見え方が異なる可能性がある。

割引率の推定には，第一著者の環境(Macbook Pro, 3.1 GHz Intel Core i7, メモリ16GB)で5-6時間ほどかかる。そのため推定済みの結果をdata_estimaed.rdsとして保存しており， 以下では使用している。もし，割引率の推定も検討されたい場合は，計算負荷が高いのでBinderでの実行を避け，各自のローカル環境で実施を推奨する。

使用するRパッケージ

library(tidyverse)

set.seed(1234)

データの読み込み

dat <- read.csv("../data/data_anonymized.csv", header = T) %>%
  # IVがドル単位(100円 = 1ドル)になっているため円単位に変換
  mutate(d0_iv = d0_iv*100,
         d2_iv = d2_iv*100,
         d30_iv = d30_iv*100,
         d180_iv = d180_iv*100,
         d365_iv = d365_iv*100)

Lempert & Pizzagalli (2010)と同様の相関分析

Lempert & Pizzagalli (2010)と同様に, 双曲割引関数, 指数関数, q-指数関数を用いて割引率を最尤推定し, 双曲割引関数によって得られた割引率とアンヘドニアの相関を検討する。

割引率を推定する自作R関数を読み込む

source( "../R/estimate_discounting_rate.R" )

最尤推定

双曲割引関数, 指数関数, q-指数関数それぞれについて最尤推定を行い, 元データに推定値を結合するが(#したコード)，以下では推定済みの結果を読み込んでいる。

# 最尤推定の処理は計算負荷が高いので，既に推定済みの結果を読み込む。
dat <- readRDS(file = "../data/data_estimated.rds")

# 以下が最尤推定になる。計算負荷が高いので，Binderでは実行しない。
# dat <- dat %>%
#  split(.$id) %>%
#  map(mutate_hdf, 100000) %>%
#  map(mutate_exp, 100000) %>%
#  map_dfr(mutate_qexp, 100000)

# 推定結果を保存
#saveRDS(object = dat, file = "data_estimated.rds")

相関分析

# 割引率が0に推定されてしまった参加者がいた場合にそのデータを除外する 
temp <- dat %>%
  filter(k_hdf != 0)

# 対数をとった割引率とアンヘドニアの相関係数の検定
cor.test(log(temp$k_hdf), temp$shaps)

# 散布図を作成
temp %>%
  ggplot(aes(x = log(k_hdf), y = shaps)) +
  geom_point(colour = "dimgray", alpha = 0.7) +
  geom_smooth(method = "lm",
              colour = "dimgray", fill = "dimgray", alpha = 0.2) +
  # windowsの場合はbase_family引数は削除
  theme_gray(base_size = 14, base_family = "HiraKakuPro-W3") +
  xlab("log(k)") + ylab("SHAPS")

AICに基づいたモデル比較および相関分析

AICを計算

aic_hdf <- sum(2*dat$fvalmin_hdf + 2 * 2)
aic_exp <- sum(2*dat$fvalmin_exp + 2 * 2)
aic_qexp <- sum(2*dat$fvalmin_qexp + 2 * 3)

cat("HDF: ", aic_hdf, "\n", "EXP: ", aic_exp, "\n", "Qexp: ", aic_qexp)

モデル比較

aics <- c(aic_hdf, aic_exp, aic_qexp)

if(which(aics == min(aics)) == 1){
  model <- "k_hdf"
}
if(which(aics == min(aics)) == 2){
  model <- "k_exp"
}
if(which(aics == min(aics)) == 3){
  model <- "k_qexp"
}

cat("Best model is: ", model)

相関分析

# 割引率が0に推定されてしまった参加者がいた場合にそのデータを除外する
temp <- dat %>%
  select(model, shaps) %>%
  filter_(paste(model, "!=", 0))

cor.test(log(temp[,1]), temp[,2])

temp %>%
  ggplot(aes(x = log(temp[,1]), y = temp[,2])) +
  geom_point(colour = "dimgray", alpha = 0.7) +
  geom_smooth(method = "lm",
              colour = "dimgray", fill = "dimgray", alpha = 0.2) +
  # windowsの場合はbase_family引数は削除
  theme_gray(base_size = 14, base_family = "HiraKakuPro-W3") +
  xlab("log(k)") + ylab("SHAPS")

R2によるモデルの当てはまりの検討

モデルの予測値を追加

dat <- dat %>%
  mutate(d0_iv_hdf_pred = 1000/(1+k_hdf*0),
         d2_iv_hdf_pred = 1000/(1+k_hdf*2),
         d30_iv_hdf_pred = 1000/(1+k_hdf*30),
         d180_iv_hdf_pred = 1000/(1+k_hdf*180),
         d365_iv_hdf_pred = 1000/(1+k_hdf*365),
         d0_iv_exp_pred = 1000*exp(-1*k_exp*0),
         d2_iv_exp_pred = 1000*exp(-1*k_exp*2),
         d30_iv_exp_pred = 1000*exp(-1*k_exp*30),
         d180_iv_exp_pred = 1000*exp(-1*k_exp*180),
         d365_iv_exp_pred = 1000*exp(-1*k_exp*365),
         d0_iv_qexp_pred = 1000/((1+(1-q_qexp)*k_qexp*0))^(1/(1-q_qexp)),
         d2_iv_qexp_pred = 1000/((1+(1-q_qexp)*k_qexp*2))^(1/(1-q_qexp)),
         d30_iv_qexp_pred = 1000/((1+(1-q_qexp)*k_qexp*30))^(1/(1-q_qexp)),
         d180_iv_qexp_pred = 1000/((1+(1-q_qexp)*k_qexp*180))^(1/(1-q_qexp)),
         d365_iv_qexp_pred = 1000/((1+(1-q_qexp)*k_qexp*365))^(1/(1-q_qexp)))

R2を計算する自作R関数を読み込む

source( "../R/compute_r2.R" )

R2を比較

双曲割引関数, 指数関数, q-指数関数それぞれについて, 予測値をもとにR2値を算出し, 元データに結合する。

dat <- dat %>%
  split(.$id) %>%
  map_dfr(calcu_R2)

cat("HDF: ", mean(dat$R2_hdf), "\n",
    "EXP: ", mean(dat$R2_exp), "\n",
    "QEXP: ", mean(dat$R2_qexp))

dat %>%
  select(id, R2_hdf, R2_exp, R2_qexp) %>%
  gather(key = model, value = R2, -id) %>%
  mutate(model = recode(model,
                        "R2_hdf" = "双曲割引関数",
                        "R2_exp" = "指数関数",
                        "R2_qexp" = "q-指数関数") %>%
           as.factor() %>%
           fct_relevel("双曲割引関数", "指数関数", "q-指数関数")) %>%
  ggplot() +
  geom_histogram(aes(x = R2),
                 alpha = 0.5, colour = "dimgray", fill = "gray") +
  scale_x_continuous(breaks = seq(-20, 1, 2)) +
  xlab("R2") + ylab("度数") +
  # windowsの場合はbase_family引数は削除
  theme_gray(base_size = 16, base_family = "HiraKakuPro-W3") +
  facet_wrap(~model)

可視化によるモデルの当てはまりの検討

双曲割引関数, 指数関数, q-指数関数の3つのモデルにおけるR2値の平均が最大の参加者および最小の参加者それぞれについて, 縦軸をIV, 横軸を遅延時間とし, 観測値に対する予測値の当てはまり度合いを検討する。なお, 各点を観測値, 実線を予測値とする。

HDF

dat %>%
  filter(meanR2 %in% c(max(meanR2), min(meanR2))) %>%
  arrange(meanR2) %>%
  select(id, d0_iv, d2_iv, d30_iv, d180_iv, d365_iv, R2_hdf) %>%
  mutate(id = paste("ID: ", id, sep = "") %>%
           factor() %>%
           fct_inorder()) %>%
  gather(key = delay, value = iv, -id, -R2_hdf) %>%
  mutate(pred = dat %>%
           filter(meanR2 %in% c(max(meanR2), min(meanR2))) %>%
           arrange(meanR2) %>%
           select(id, R2_hdf,
                  d0_iv_hdf_pred, d2_iv_hdf_pred, d30_iv_hdf_pred,
                  d180_iv_hdf_pred, d365_iv_hdf_pred) %>%
           gather(key = delay, value = pred, -id, -R2_hdf) %>%
           pull(pred),
         delay = recode(delay,
                        "d0_iv" = 0,
                        "d2_iv" = 2,
                        "d30_iv" = 30,
                        "d180_iv" = 180,
                        "d365_iv" = 365),
         # expression表現にするために値を作成
         R2_hdf = paste("R^2", "==", round(R2_hdf, 3), sep = "")) %>%
  ggplot() + 
  geom_point(aes(x = delay, y = iv),
             size = 2) + 
  geom_line(aes(x = delay, y = pred)) +
  geom_text(aes(x = 200, y = 800, label = R2_hdf),
            size = 5, parse = TRUE) +
  # windowsの場合はbase_family引数は削除
  theme_gray(base_size = 16, base_family = "HiraKakuPro-W3") + 
  xlab("遅延時間") + ylab("IV/IV_predicted") + 
  facet_wrap(~id)

EXP

dat %>%
  filter(meanR2 %in% c(max(meanR2), min(meanR2))) %>%
  arrange(meanR2) %>%
  select(id, d0_iv, d2_iv, d30_iv, d180_iv, d365_iv, R2_exp) %>%
  mutate(id = paste("ID: ", id, sep = "") %>%
           factor() %>%
           fct_inorder()) %>%
  gather(key = delay, value = iv, -id, -R2_exp) %>%
  mutate(pred = dat %>%
           filter(meanR2 %in% c(max(meanR2), min(meanR2))) %>%
           arrange(meanR2) %>%
           select(id, R2_exp,
                  d0_iv_exp_pred, d2_iv_exp_pred, d30_iv_exp_pred,
                  d180_iv_exp_pred, d365_iv_exp_pred) %>%
           gather(key = delay, value = pred, -id, -R2_exp) %>%
           pull(pred),
         delay = recode(delay,
                        "d0_iv" = 0,
                        "d2_iv" = 2,
                        "d30_iv" = 30,
                        "d180_iv" = 180,
                        "d365_iv" = 365),
         # expression表現にするために値を作成
         R2_exp = paste("R^2", "==", round(R2_exp, 3), sep = "")) %>%
  ggplot() + 
  geom_point(aes(x = delay, y = iv),
             size = 2) + 
  geom_line(aes(x = delay, y = pred)) +
  geom_text(aes(x = 250, y = 800, label = R2_exp),
            size = 5, parse = TRUE) +
  # windowsの場合はbase_family引数は削除
  theme_gray(base_size = 16, base_family = "HiraKakuPro-W3") + 
  xlab("遅延時間") + ylab("IV/IV_predicted") + 
  facet_wrap(~id)

Qexp

dat %>%
  filter(meanR2 %in% c(max(meanR2), min(meanR2))) %>%
  arrange(meanR2) %>%
  select(id, d0_iv, d2_iv, d30_iv, d180_iv, d365_iv, R2_qexp) %>%
  mutate(id = paste("ID: ", id, sep = "") %>%
           factor() %>%
           fct_inorder()) %>%
  gather(key = delay, value = iv, -id, -R2_qexp) %>%
  mutate(pred = dat %>%
           filter(meanR2 %in% c(max(meanR2), min(meanR2))) %>%
           arrange(meanR2) %>%
           select(id, R2_qexp,
                  d0_iv_qexp_pred, d2_iv_qexp_pred, d30_iv_qexp_pred,
                  d180_iv_qexp_pred, d365_iv_qexp_pred) %>%
           gather(key = delay, value = pred, -id, -R2_qexp) %>%
           pull(pred),
         delay = recode(delay,
                        "d0_iv" = 0,
                        "d2_iv" = 2,
                        "d30_iv" = 30,
                        "d180_iv" = 180,
                        "d365_iv" = 365),
         # expression表現にするために値を作成
         R2_qexp = paste("R^2", "==", round(R2_qexp, 3), sep = "")) %>%
  ggplot() + 
  geom_point(aes(x = delay, y = iv),
             size = 2) + 
  geom_line(aes(x = delay, y = pred)) +
  geom_text(aes(x = 200, y = 800, label = R2_qexp),
            size = 5, parse = TRUE) +
  # windowsの場合はbase_family引数は削除
  theme_gray(base_size = 16, base_family = "HiraKakuPro-W3") + 
  xlab("遅延時間") + ylab("IV/IV_predicted") + 
  facet_wrap(~id)

AUCを用いた相関分析

AUCを計算する自作R関数を読み込む

source( "../R/compute_auc.R" )

AUCを計算

AUCを算出し元データに結合する。

dat <- dat %>%
  split(.$id) %>%
  map_dfr(calcu_auc)

相関分析

cor.test(dat$auc, dat$shaps)

dat %>%
  ggplot(aes(x = auc, y = shaps)) +
  geom_point(colour = "dimgray", alpha = 0.7) +
  geom_smooth(method = "lm",
              colour = "dimgray", fill = "dimgray", alpha = 0.2) +
  # windowsの場合はbase_family引数は削除
  theme_gray(base_size = 12, base_family = "HiraKakuPro-W3") +
  xlab("Area Under the Curve") + ylab("SHAPS")

ykunisato/somatori_kunisato_2019_replication_study documentation built on Nov. 5, 2019, 1:20 p.m.