第2回mimic演習
Takahiro Mihara
（前回までの資料も一部記載していますが、コードは最低限にしています。復習は前回の資料を参照してください。）
package
pacman::p_load(tidyverse, lubridate, tableone, skimr,
duckdb, arrow, RSQLite)データベース
FドライブのDatabaseフォルダ内のmimic.dbにアクセスする場合は以下のコード（各自の状況に合わせて適宜修正）。 （mimic.dbが無い場合は以下の操作で空のmimic.dbが作成される）
con_mimic <- dbConnect(duckdb::duckdb(), "F:/Database/mimic.db")mimic.dbの中に入っているデータリスト
dbListTables(con_mimic)## [1] "admissions" "apsiii"
## [3] "callout" "caregivers"
## [5] "chartevents" "comorbidities_AHRQ"
## [7] "comorbidities_Quan" "cptevents"
## [9] "d_cpt" "d_icd_diagnoses"
## [11] "d_icd_procedures" "d_items"
## [13] "d_labitems" "datetimeevents"
## [15] "db_HtWt" "db_LVEF"
## [17] "db_LVEF_assess" "db_PEO"
## [19] "db_comorb" "db_ethnicity"
## [21] "db_patients" "db_rrt"
## [23] "db_sofa" "db_urine"
## [25] "db_venti" "db_vital"
## [27] "df_lab_list" "diagnoses_icd"
## [29] "drgcodes" "height_firstday"
## [31] "height_weight" "icustay_details"
## [33] "icustays" "inputevents_cv"
## [35] "inputevents_mv" "labevents"
## [37] "microbiologyevents" "noteevents"
## [39] "outputevents" "patients"
## [41] "pivot_RRT" "pivot_height"
## [43] "pivot_labo" "pivot_urine_output"
## [45] "pivot_vitals" "prescriptions"
## [47] "procedureevents_mv" "procedures_icd"
## [49] "services" "sofa"
## [51] "test" "transfers"
## [53] "urine_output" "urine_output_firstday"
## [55] "ventilation_classification" "ventilation_duration"
## [57] "weight_duration" "weight_firstday"# dbDisconnect(con_mimic)使用データ
作成したデータ
db_patients <- tbl(con_mimic, "db_patients")<第1回講義・課題> 解答例
Exposureについて：H2RAR使用したかどうかを示す列をdb_patientsに追加する。
ヒント：どのデータを利用するのか考えてください。str_detect()など覚えた関数は最大限利用しましょう。
薬剤使用群と非使用群がそれぞれ何名ずつかカウントし、論文の記載と近いことを確認しましょう。Outcomeについて：入院日を起点として30日死亡したかどうかを示す列を上記にさらに追加する。
ヒント：NA処理が必要になるかも。
（多分簡単にできると思います。）
2. 暴露情報；H2RAs
“prescriptions”
db_presc <- tbl(con_mimic,"prescriptions")
db_presc |> colnames()## [1] "ROW_ID" "SUBJECT_ID" "HADM_ID"
## [4] "ICUSTAY_ID" "STARTDATE" "ENDDATE"
## [7] "DRUG_TYPE" "DRUG" "DRUG_NAME_POE"
## [10] "DRUG_NAME_GENERIC" "FORMULARY_DRUG_CD" "GSN"
## [13] "NDC" "PROD_STRENGTH" "DOSE_VAL_RX"
## [16] "DOSE_UNIT_RX" "FORM_VAL_DISP" "FORM_UNIT_DISP"
## [19] "ROUTE"db_presc |> head(10)db_presc |> count()col_presc <- c("SUBJECT_ID", "HADM_ID", "ICUSTAY_ID", "DRUG", "DRUG_NAME_POE", "DRUG_NAME_GENERIC", "FORMULARY_DRUG_CD")
db_presc_drug <- db_presc |> select(col_presc)
db_presc_drug |> head(10)一旦薬剤名のみのユニークなデータにして（8500行程度）collect()する
薬剤名全体を見渡してみる（理解する）
View（）した後に使える全体検索機能が便利
drugs <-
db_presc_drug |>
distinct(DRUG, DRUG_NAME_POE, DRUG_NAME_GENERIC, FORMULARY_DRUG_CD)
drugs |> count()drugs |> head(10)drugs |> arrange(DRUG)df_drugs <- drugs |> collect()
df_drugs |> head(20)# df_drugs |> View()H2RAa: Famotidine, Ranitidine, Nizatidine, Cimetidine の4剤
- Loratadine, Loratidine, CLARITIN, ClaritinはH1拮抗薬なので除外
- Azacitidine は抗悪性腫瘍薬
商品名もチェックする(Viewの画面で薬剤名を入れると商品名もある程度見つかる)
“Zantac 75”, “Zantac Maximum Strength”
FORMURALYには追加情報なさそうなので、今回は使用しない。
H2RASは薬剤名をネット上から探す。全てにtidineがついているのでまずはそこから探す。
cimetidine
ranitidine
famotidine
nizatidine
roxatidine
lafutidine
lavoltidine
niperotidine
col_h2 <-
df_drugs |>
filter(str_detect(DRUG, regex("tidine", ignore_case = T))|
str_detect(DRUG_NAME_POE, regex("tidine", ignore_case = T))|
str_detect(DRUG_NAME_GENERIC, regex("tidine", ignore_case = T))
) |>
select(DRUG) |> pull() |> unique()
col_h2 |> dput()## c("Famotidine (PO)", "NEO*PO*Ranitidine", "Famoti", "Ranitidine",
## "Cimetidine", "Ranitidine Hcl", "famotidine", "Ranitidine (Liquid)",
## "Famotidine (IV)", "Famotidine", "Azacitidine", "Cimetidine HCl",
## "Nizatidine", "Ranitidine HCl", "Loratadine", "Claritin", "CLARITIN",
## "Azacitidine (Subcut)", "Famo")除外と追記をコードで行って薬剤名リストを作るのが理想だが、上記のように大枠をdput()でベクトル形式の出力をして、結果をコピペしたものに削除＆追記しても良い。
その場合は以下のように薬剤名リストに名前を付けて格納しておく。
vec_h2ras <- c("Cimetidine", "Cimetidine HCl", "Famo", "Famoti", "Famotidine",
"Famotidine (IV)", "Famotidine (PO)", "NEO*PO*Ranitidine",
"Nizatidine", "Ranitidine", "Ranitidine (Liquid)", "Ranitidine HCl",
"Ranitidine Hcl", "famotidine",
"Zantac 75", "Zantac Maximum Strength")心不全患者データ（db_patients）にH2RAs使用者かどうかフラグを立てる
その入院中でなくても使用歴のある患者は使用者とする為、患者IDでフラグを立てる
# db_patients |> summarize(HADM_count = n_distinct(HADM_ID),
# SUBJECT_count = n_distinct(SUBJECT_ID),
# count = n())
H2RAs_user <-
db_presc_drug |>
filter(DRUG %in% vec_h2ras) |>
select(SUBJECT_ID) |>
pull()
db_patients_h2 <-
db_patients |>
mutate(h2ras = ifelse(SUBJECT_ID %in% H2RAs_user, 1, 0),
.after = ICUSTAY_ID)
db_patients_h2 |> head(10)論文通りに分類できたっぽい ＜db_patients_h2で完成＞
db_patients_h2 |> count(h2ras)3. アウトカム；院内死亡、30日死亡
入室日と死亡日から死亡経過日数を算出し、30日死亡を算出する。
・HOSPITAL_EXPIRE_FLAG:
患者が入院期間内に死亡したかどうかを示す。
・EXPIRE_FLAG:
患者が死亡したかどうか、すなわちDODがNULLかどうかを示す。病院内での死亡（DOD_HOSP）と、The
social security master death
indexと患者を照合して特定した死亡（DOD_SSN）の両方が含まれる。
“icustays”にはINTIME,OUTTIMEというICU入室日と退室日があるので、ここからICU死亡は算出できる。どのタイミングでICU入室するか等臨床的意義が不明な点もあるので今回のShadowingには含めない。
hospital los は退院日と入院日の差で算出可能
ICU
losもICUSTAYSに入室日と退室日更にはLOSも入っているので算出は容易。
＜db_PEOとする（P,E,Oの絞り込み＆変数追加完了）＞
db_PEO <-
db_patients_h2 |>
mutate(death_elaps = day(DOD - ADMITTIME),
mortality30 = ifelse(death_elaps <= 30, 1, 0),
hospital_death = HOSPITAL_EXPIRE_FLAG) |>
mutate(INSURANCE = ifelse(INSURANCE == "Medicaid", "Medicare", INSURANCE)) |>
replace_na(list(mortality30 = 0))
db_PEO <- db_PEO |> compute()
db_PEO |> head(10)db_PEO |> count()# dbListTables(con_mimic)continuous <- c("age")
factors <- c("ADMISSION_TYPE", "INSURANCE", "GENDER", "mortality30")
tableone <-
db_PEO |>
select(c(continuous, factors, "h2ras"))|>
collect() |>
CreateTableOne(vars = c(continuous, factors), strata = "h2ras", factorVars = factors)
print(tableone, smd = TRUE, missing = TRUE) ## Stratified by h2ras
## 0 1 p test SMD Missing
## n 5949 4441
## age (mean (SD)) 74.35 (14.30) 71.25 (13.97) <0.001 0.219 0.0
## ADMISSION_TYPE (%) <0.001 0.354 0.0
## ELECTIVE 412 ( 6.9) 795 (17.9)
## EMERGENCY 5267 (88.5) 3546 (79.8)
## URGENT 270 ( 4.5) 100 ( 2.3)
## INSURANCE (%) <0.001 0.122 0.0
## Government 80 ( 1.3) 58 ( 1.3)
## Medicare 4830 (81.2) 3400 (76.6)
## Private 1013 (17.0) 969 (21.8)
## Self Pay 26 ( 0.4) 14 ( 0.3)
## GENDER = M (%) 3074 (51.7) 2446 (55.1) 0.001 0.068 0.0
## mortality30 = 1 (%) 1301 (21.9) 546 (12.3) <0.001 0.256 0.0if(dbExistsTable(con_mimic, "db_PEO")) {
dbRemoveTable(con_mimic, "db_PEO")
}
dbWriteTable(con_mimic,
"db_PEO",
db_PEO %>% collect(),
overwrite = TRUE) 第2回講義
４．背景因子・交絡因子（共変量）
- BigQuery
論文独自のPatients、Exposure、Outcomeと違って共変量候補はどの論文でも共通の項目が多い（身長、体重など）
多くの抽出コードは公式が作成したものが公開されている。
concept概念と呼ばれるファイル群。BigQueryに登録すると全てダウンロードできるのでお勧め（むしろ使うべき）。
https://github.com/MIT-LCP/mimic-code/tree/main/mimic-iii/concepts
mimic iii 開発チームの公開しているコードはSQL。これらのコードをコピペして使うか、またはmimiciii_derivedに既に作成されたファイルがあるのでダウンロードする。
ダウンロード方法は授業スライド参照。
BigQueryでPhysioNetデータにログインすると、上記のテーブルは全て作成済みの状態でおかれている。mimiciii_derived。
論文再現に必要なBigQueryからダウンロードすべきファイルは18ファイルある（もっとあるかも）。
ダウンロードしたファイルを保存し、mimic.dbに取り込む。
以下のパスに保存した場合のコードを示す。
F:\BigDatas\mimic-iii-clinical-database-1.4\BigQuery\
ダウンロードしたときにつけた名前に応じて適宜コードを修正する事！
bg_list = list.files(path = "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery",
pattern = ".csv$",
full.names = T)
bg_list## [1] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/apsiii.csv"
## [2] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/comorbidities_AHRQ.csv"
## [3] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/comorbidities_Quan.csv"
## [4] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/height_firstday.csv"
## [5] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/height_weight.csv"
## [6] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/ICU_STAY_v1.3.csv"
## [7] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/icustay_details.csv"
## [8] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/pivot_height.csv"
## [9] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/pivot_labo.csv"
## [10] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/pivot_RRT.csv"
## [11] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/pivot_urine_output.csv"
## [12] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/pivot_vitals.csv"
## [13] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/sofa.csv"
## [14] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/urine_output.csv"
## [15] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/urine_output_firstday.csv"
## [16] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/ventilation_classification.csv"
## [17] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/ventilation_duration.csv"
## [18] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/weight_duration.csv"
## [19] "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery/weight_firstday.csv"↑ICU_STAY_v1.3.csvは不要。
bg_names = list.files(path = "F:\\BigDatas\\mimic-iii-clinical-database-1.4\\BigQuery",
pattern = ".csv$",
full.names = F)
bg_names <- bg_names |> str_remove(".csv")
bg_names## [1] "apsiii" "comorbidities_AHRQ"
## [3] "comorbidities_Quan" "height_firstday"
## [5] "height_weight" "ICU_STAY_v1.3"
## [7] "icustay_details" "pivot_height"
## [9] "pivot_labo" "pivot_RRT"
## [11] "pivot_urine_output" "pivot_vitals"
## [13] "sofa" "urine_output"
## [15] "urine_output_firstday" "ventilation_classification"
## [17] "ventilation_duration" "weight_duration"
## [19] "weight_firstday"今回も100MBを超えるファイルがいくつかある（最大430MB）。このサイズであればメモリにあげてしまっても良いように思えるが、念のため.db化して扱う方針とする。
後から追加する場合はループを使わず一つ一つ確認しながら入れた方が無難。
# for (i in seq_along(bg_list)) {
#
# callback <- function(df, pos) {
# dbWriteTable(con_mimic, bg_names[i], df, append = TRUE)
# }
#
# read_csv_chunked(bg_list[i],
# DataFrameCallback$new(callback),
# chunk_size = 100000)
# }dbListTables(con_mimic)## [1] "admissions" "apsiii"
## [3] "callout" "caregivers"
## [5] "chartevents" "comorbidities_AHRQ"
## [7] "comorbidities_Quan" "cptevents"
## [9] "d_cpt" "d_icd_diagnoses"
## [11] "d_icd_procedures" "d_items"
## [13] "d_labitems" "datetimeevents"
## [15] "db_HtWt" "db_LVEF"
## [17] "db_LVEF_assess" "db_PEO"
## [19] "db_comorb" "db_ethnicity"
## [21] "db_patients" "db_rrt"
## [23] "db_sofa" "db_urine"
## [25] "db_venti" "db_vital"
## [27] "dbplyr_MCdiXrvfqi" "df_lab_list"
## [29] "diagnoses_icd" "drgcodes"
## [31] "height_firstday" "height_weight"
## [33] "icustay_details" "icustays"
## [35] "inputevents_cv" "inputevents_mv"
## [37] "labevents" "microbiologyevents"
## [39] "noteevents" "outputevents"
## [41] "patients" "pivot_RRT"
## [43] "pivot_height" "pivot_labo"
## [45] "pivot_urine_output" "pivot_vitals"
## [47] "prescriptions" "procedureevents_mv"
## [49] "procedures_icd" "services"
## [51] "sofa" "test"
## [53] "transfers" "urine_output"
## [55] "urine_output_firstday" "ventilation_classification"
## [57] "ventilation_duration" "weight_duration"
## [59] "weight_firstday"- height, weight, BMI
使用するデータ
heightweight
heightfirstday
pivoted_height(pivot_height)
weightfirstday
weightduration
db_heightweight <- tbl(con_mimic, "height_weight")
db_heightFirst <- tbl(con_mimic, "height_firstday")
db_weightFirst <- tbl(con_mimic, "weight_firstday")
db_weightDuration <- tbl(con_mimic, "weight_duration")
db_pivot_height <- tbl(con_mimic, "pivot_height") #subject_id
db_heightweight |> count()db_heightFirst |> count()db_weightFirst |> count()db_weightDuration |> count()db_pivot_height |> count()db_weightDurationとdb_pivot_heightは行数が桁違い（複数測定が縦に入っているに違いない）
db_weightDurationはICU入室期間内で何回も測定されたデータが全て入っている。
ここでは最初のデータを採用する方針とする。
db_weightDuration1 <-
db_weightDuration |>
# arrange(icustay_id, starttime) |>
distinct(icustay_id, .keep_all = T) |>
rename(weight_duration = weight) |>
select(icustay_id, weight_duration) 各データをマージしていく。
（db_pivot_heightはSUBJECT_IDがキーになるので後で処理する）
db_hw_1 <-
db_PEO |>
left_join(db_heightweight, by = c("ICUSTAY_ID" = "icustay_id"), copy = T)
db_hw_1 <- db_hw_1 |> compute()
db_hw_2 <-
db_hw_1 |>
left_join(db_heightFirst, by = c("ICUSTAY_ID" = "icustay_id"), copy = T)
db_hw_2 <- db_hw_2 |> compute()
db_hw_3 <-
db_hw_2 |>
left_join(db_weightFirst, by = c("ICUSTAY_ID" = "icustay_id"), copy = T)
db_hw_3 <- db_hw_3 |> compute()
db_hw_4 <-
db_hw_3 |>
left_join(db_weightDuration1, by = c("ICUSTAY_ID" = "icustay_id"), copy = T)
db_hw_4 <- db_hw_4 |> compute()
db_hw_4 |> count()db_hw_4 |> head(10)身長、体重を1つの値にする。
どの値を優先するのかを決めてコードを組み立てる（論文未記載）
（以下では、height系全てNAでなければ、height_firstが採用される）
db_PEO_HtWt <-
db_hw_4 |>
mutate(
height = case_when(
!is.na(height_first) ~ height_first,
!is.na(height) ~ height,
!is.na(height_echo) ~ height_echo,
!is.na(height_chart) ~ height_chart
),
weight = case_when(
!is.na(weight_first) ~ weight_first,
!is.na(weight_admit) ~weight_admit,
!is.na(weight) ~ weight,
!is.na(weight_duration) ~ weight_duration,
!is.na(weight_echoinhosp) ~ weight_echoinhosp,
!is.na(weight_daily) ~ weight_daily
)
) |>
select(!c("weight_first", "weight_min", "weight_max", "height_first", "height_min",
"height_max", "height_chart", "height_echo", "weight_admit",
"weight_daily", "weight_echoinhosp", "weight_echoprehosp", "weight_duration"))
db_PEO_HtWt |> count()db_pivot_heightは患者ごとのデータ。複数測定が全て縦に入っているので、最初の一つに絞る。
db_pivot_height1 <-
db_pivot_height |>
distinct(subject_id, .keep_all = T) |>
rename(height_pivot = height) |>
select(subject_id, height_pivot)
db_pivot_height1マージする
db_PEO_HtWt1 <-
db_PEO_HtWt |>
left_join(db_pivot_height1, by = c("SUBJECT_ID" = "subject_id"), copy = T) |>
mutate(height = ifelse(!is.na(height), height, height_pivot)) |>
select(!c(height_pivot))
db_PEO_HtWt1 <- db_PEO_HtWt1 |> compute()連続値は必ずヒストグラムを確認する
cols_continuous <- c("height", "weight")
db_PEO_HtWt1 |>
select(cols_continuous) |>
pivot_longer(cols = cols_continuous, names_to = "name", values_to = "value") |>
ggplot()+
geom_histogram(aes(x = value), color = "black")+
facet_wrap(~ name, scales = "free", ncol = 2) +
theme_bw()
極端な値は入力ミスとみなしてNAに変換する
身長250cm以上、80cm以下
体重10kg以下
最後にBMIを計算する
db_PEO_HtWt2 <-
db_PEO_HtWt1 |>
mutate(height = ifelse(height < 250 & height > 80, height, NA),
weight = ifelse(weight > 10, weight, NA)) |>
mutate(BMI = weight/(height/100)^2)
db_PEO_HtWt2 <- db_PEO_HtWt2 |> compute()2群間でのBMI平均値・SDの算出。
論文未記載のため独自に判断した部分が多かったため、まったく一致することはないが、ある程度近い値になっている。
vars <- c("BMI", "height", "weight")
tableone <-
db_PEO_HtWt2 |>
select(c(vars, "h2ras"))|>
collect() |>
CreateTableOne(vars = vars, strata = "h2ras", factorVars = factors)
print(tableone, smd = TRUE, missing = TRUE) |> write.csv("test.csv")## Stratified by h2ras
## 0 1 p test SMD Missing
## n 5949 4441
## BMI (mean (SD)) 28.63 (8.03) 28.99 (8.11) 0.036 0.044 13.9
## height (mean (SD)) 167.67 (10.97) 168.20 (10.78) 0.020 0.049 12.6
## weight (mean (SD)) 80.16 (24.65) 81.98 (24.62) <0.001 0.074 3.0# summary(tableone)連続値は必ず全てヒストグラムで確認する事。
cols_continuous <- c("age","height", "weight", "BMI")
db_PEO_HtWt2 |>
select(cols_continuous) |>
pivot_longer(cols = cols_continuous, names_to = "name", values_to = "value") |>
ggplot()+
geom_histogram(aes(x = value), color = "black")+
facet_wrap(~ name, scales = "free", ncol = 2) +
theme_bw()
最後のマージ用ファイル
背景因子・交絡因子はたくさんあるので、それぞれHADM_IDとだけ紐づけたファイルを作っておき、最後にすべてをマージする。
db_HtWt <-
db_PEO_HtWt2 |>
select("HADM_ID", "weight", "height", "BMI")データベース上に書き込む
if(dbExistsTable(con_mimic, "db_HtWt")) {
dbRemoveTable(con_mimic, "db_HtWt")
}
dbWriteTable(con_mimic,
"db_HtWt",
db_HtWt %>% collect(),
overwrite = TRUE) # dbDisconnect(con_mimic, shutdown=TRUE)<第2回講義・課題>
BigQueryから必要なファイルをダウンロードし、自分の環境下でmimic.dbに入れる（スライドにある18ファイルは参考までに）
（希望者のみでOK。希望しない人は個別に連絡ください。）
入手したテーブルを利用して、対象患者の以下の項目を抽出する
（論文再現が目的であるが、全く同じにならないことに注意）
（授業と同様、HADM_IDと目的変数のみのファイルを作成する）