pacman::p_load(tidyverse,readxl,ggplot2,dplyr,tableone,ggpubr,ggrepel,DBI,skimr,dbplyr,pryr,lubridate, survival, survminer,car,svglite,MatchIt, cobalt)

#データ読み込み

user.data<- read.csv(file= "/Users/abekohei/Downloads/1_前処理Try用/script_for_provision/semi_processed_data/user.data.semi_processed.destructed.csv", stringsAsFactors = F)  # 破壊済みデータ
step.data<- read.csv(file= "/Users/abekohei/Downloads/1_前処理Try用/script_for_provision/semi_processed_data/step.data.semi_processed.csv"           , stringsAsFactors = F)
step.monthly.data <- read.csv(file= "/Users/abekohei/Downloads/1_前処理Try用/script_for_provision/semi_processed_data/step.monthly.data.semi_processed.csv"   , stringsAsFactors = F)
exam.data<- read.csv(file= "/Users/abekohei/Downloads/1_前処理Try用/script_for_provision/semi_processed_data/exam.data.semi_processed.destructed.csv", stringsAsFactors = F)  # 破壊済みデータ
cancer.data<- read.csv(file= "/Users/abekohei/Downloads/1_前処理Try用/script_for_provision/semi_processed_data/cancer.data.semi_processed.csv", stringsAsFactors = F)
dialysis.data<- read.csv(file= "/Users/abekohei/Downloads/1_前処理Try用/script_for_provision/semi_processed_data/dialysis.data.semi_processed.csv", stringsAsFactors = F)

head(user.data)          # 条件1に合致するkencomユーザーのデータ。仮個人id単位

head(step.data)          # 条件1に合致するkencomユーザーの歩数データ。仮個人id単位（歩数も個人単位に要約済み）

head(step.monthly.data)  # 条件1に合致するkencomユーザーの月単位の歩数データ。仮個人id x 年月単位（歩数も年月単位に要約済み）

head(exam.data)          # 条件1に合致するkencomユーザーの、2015年4月 ~ 2018年6月の健診データ。仮個人id x 健診受診年月単位。

head(cancer.data)        # がんの診断が出ている人のデータ。レセプト年月は診断が出た年月を示す。仮個人id x レセプト年月単位（条件1に合致していない人も含む）

head(dialysis.data)      # 透析の実施がある人のデータ。レセプト年月は透析を実施した年月を示す。仮個人id x レセプト年月単位（条件1に合致していない人も含む）

exam.data %>% skimr::skim()

Data summary
Name	Piped data
Number of rows	381352
Number of columns	9
_______________________
Column type frequency:
character	2
numeric	7
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
仮個人id	0	1	10	10	0	114787	0
健診受診年月	0	1	10	10	0	48	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
体重_nz	615	1.00	63.99	12.64	31.39	54.72	63.24	71.68	144.68	▃▇▂▁▁
収縮期血圧_nz	731	1.00	119.19	16.16	74.15	107.80	118.29	129.13	213.66	▂▇▃▁▁
拡張期血圧_nz	727	1.00	74.06	11.86	38.01	65.63	73.40	81.76	134.02	▁▇▆▁▁
ldlコレステロール_nz	3322	0.99	122.44	30.97	28.51	100.67	120.56	141.94	283.48	▁▇▅▁▁
hdlコレステロール_nz	3218	0.99	63.55	16.35	23.82	51.64	61.64	73.35	153.16	▃▇▂▁▁
中性脂肪_nz	3218	0.99	105.94	82.28	19.96	59.20	84.43	125.94	1588.09	▇▁▁▁▁
hba1c_nz	16042	0.96	5.49	0.57	3.99	5.18	5.41	5.67	13.23	▇▁▁▁▁

条件２

「資格取得 ~ 資格喪失」が「kencom登録月􏴁15ヶ月前 ~ kencom登録月前月(total 15ヶ月)」を包含している

user.data %>% head()

user.data_2<-
  user.data %>% 
  mutate(kencom登録年月=ymd(kencom登録年月)) %>%
  mutate(資格取得年月=ymd(資格取得年月)) %>% 
  filter(資格取得年月<=kencom登録年月-months(15)&資格喪失年月>=(kencom登録年月-months(1)))

条件３

「レセプトデータ開始 ~レセプトデータ終了」が「kencom登録月􏴁15ヶ月前 ~ kencom登録月前月(total 15ヶ月)」を包含している

user.data_3<-
  user.data_2 %>% 
  mutate(レセプトデータ開始年月=ymd(レセプトデータ開始年月)) %>% 
  mutate(レセプトデータ終了年月=ymd(レセプトデータ終了年月)) %>% 
  filter(レセプトデータ開始年月<=kencom登録年月-months(15)&レセプトデータ終了年月>=(kencom登録年月-months(1)))

条件４

kencom登録月􏴁12ヶ月前 ~ kencom登録月􏴁前月(total 12ヶ月) に健診受診あり ※上記期間中に複数回健診受診がある場合􏴂最古を採用

#カラム名を変更
colnames(exam.data)

## [1] "仮個人id"             "健診受診年月"         "体重_nz"             
## [4] "収縮期血圧_nz"        "拡張期血圧_nz"        "ldlコレステロール_nz"
## [7] "hdlコレステロール_nz" "中性脂肪_nz"          "hba1c_nz"

exam.data<-
  exam.data %>%
  rename(weight_nz=体重_nz) %>%
  rename(systolicBP_nz=収縮期血圧_nz) %>%
  rename(diastolicBP_nz=拡張期血圧_nz) %>%
  rename(ldl_c_nz=ldlコレステロール_nz) %>%
  rename(hdl_c_nz=hdlコレステロール_nz) %>%
  rename(TG_nz=中性脂肪_nz) %>%
  rename(HbA1c_nz=hba1c_nz)

user.data_4<-
  user.data_3 %>% 
  left_join(exam.data,by="仮個人id") %>% 
  mutate(健診受診年月=ymd(健診受診年月)) %>%
  filter(健診受診年月>=kencom登録年月-months(12)&健診受診年月<=kencom登録年月-months(1)) %>% 
  group_by(仮個人id) %>%
  arrange(健診受診年月) %>%
  mutate(番号 = row_number()) %>%
  filter(番号==1) %>%
  select(-番号)  #ベースラインの健診データを抽出
user.data_4 %>% head()

条件5

kencom登録月􏴁12ヶ月前 ~ kencom登録月􏴁前月(total 12ヶ月) に毎月1日以上歩数データあり。

user.data_5<-
  user.data_4 %>% 
  left_join(step.data,by="仮個人id") %>% 
  filter(データ有り月数_pre==12)
user.data_5 %>% head()

# 歩数変化量
user.data_5<-
  user.data_5 %>% 
  mutate(step.diff    = step_mean_post    - step_mean_pre,
         step.diff.wd = step_mean_post.wd - step_mean_pre.wd,
         step.diff.we = step_mean_post.we - step_mean_pre.we)

条件６

「kencom登録月􏴁15ヶ月前 ~ kencom登録月前月(total 15ヶ月)」に「がん􏴁診断と投薬」なしかつ「透析􏴁実施」なし

cancer.data<-cancer.data %>% rename(がん診断年月=レセプト年月)
dialysis.data<-dialysis.data %>% rename(透析実施年月=レセプト年月)

user.data_5_cancer<-
  user.data_5 %>% 
  left_join(cancer.data,by="仮個人id") %>% 
  filter(!is.na(がん診断年月)) %>% 
  filter(がん診断年月>=kencom登録年月-months(15)&
           がん診断年月<=kencom登録年月-months(1))

user.data_5_dialysis<-
  user.data_5 %>% 
  left_join(dialysis.data,by="仮個人id") %>% 
  filter(!is.na(透析実施年月)) %>%
  filter(透析実施年月>=kencom登録年月-months(15)&
           透析実施年月<=kencom登録年月-months(1))

user.data_5_cancer_dialysis<-
  user.data_5_cancer %>% 
  left_join(dialysis.data,by="仮個人id") 　　#該当者誰もいないので除外者なし

exclude_user.data_5<-user.data_5_cancer_dialysis$仮個人id
user.data_6<-
  user.data_5 %>% 
  filter(!仮個人id %in% exclude_user.data_5) 

#user.data_6<-user.data_5 %>% anti_join(user.data_5_cancer_dialysis, by = "仮個人id")でも可能

条件7

kencom登録月 ~ kencom登録月􏴁11ヶ月後(total 12ヶ月)に毎月 1日以上歩数データがある。また、kencom登録前後を通して月あたりの平均データあり日数が 20日以上。また、kencom登録月の12ヶ月前 ~ kencom登録月の前月の平均歩数のSDが99%tileより大きい人を除外する。 ※「歩数データが未来12ヶ月にあり、月平均20日以上あり」に絞り込んだあとに99%tileを算出する。

colnames(user.data_6)

##  [1] "仮個人id"               "性別_sfl"               "年齢_sfl"              
##  [4] "kencom登録年月"         "資格取得年月"           "資格喪失年月"          
##  [7] "レセプトデータ開始年月" "レセプトデータ終了年月" "健診受診年月"          
## [10] "weight_nz"              "systolicBP_nz"          "diastolicBP_nz"        
## [13] "ldl_c_nz"               "hdl_c_nz"               "TG_nz"                 
## [16] "HbA1c_nz"               "データ有り月数_total"   "データ有り月数_pre"    
## [19] "データ有り月数_post"    "月ごとの平均データ数"   "step_mean_pre"         
## [22] "step_mean_post"         "step_sd_pre"            "step_sd_post"          
## [25] "step_mean_pre.wd"       "step_mean_post.wd"      "step_sd_pre.wd"        
## [28] "step_sd_post.wd"        "step_mean_pre.we"       "step_mean_post.we"     
## [31] "step_sd_pre.we"         "step_sd_post.we"        "step.diff"             
## [34] "step.diff.wd"           "step.diff.we"

#提供スクリプトはデータ有り月数_total == 24 だったが、データ有り月数_post == 12 
#user_exam_monthly_stepは条件5でデータ有り月数_pre==12に絞っているので結局結果は同じ。

user.data_6_post<-
  user.data_6 %>% 
  filter(データ有り月数_post==12,月ごとの平均データ数>=20) 

user.data_6 %>%
  filter(データ有り月数_total == 24 & 月ごとの平均データ数 >= 20)

user.data_6_post %>% head()

kencom登録月の12ヶ月前 ~ kencom登録月の前月の平均歩数のSDが99%tileより大きい人を除外する。

user.data_7<-
user.data_6_post %>% 
  filter(step_sd_pre<quantile(user.data_6_post$step_sd_pre,0.99,na.rm=T)) 
user.data_7 %>% head()

#　条件8 kencom登録月􏴁12ヶ月後 ~ kencom登録月􏴁23ヶ月後(total 12ヶ月)に受診した健診 ※上記期間中に複数回健診受診がある場合􏴂最古を採用

exam.data %>% head()

user.data_7 <-
  user.data_7 %>% 
  rename(baseline健診受診年月=健診受診年月)

user.data_8 <-
  user.data_7 %>% 
  left_join(exam.data,by="仮個人id",suffix=c(".b",".f")) %>% 
  mutate(followup健診受診年月=ymd(健診受診年月)) %>%
  filter(followup健診受診年月>=kencom登録年月+months(12)&健診受診年月<=kencom登録年月+months(23)) %>% 
  group_by(仮個人id) %>%
  arrange(followup健診受診年月) %>%
  mutate(番号 = row_number()) %>%
  filter(番号==1) %>%
  select(-番号,-健診受診年月)　　　　　　　#フォローアップの健診データを抽出
user.data_8 %>% head()

ベースライン健診データとフォローアップ健診データを結合

#user_exam<-
  #user_exam_b %>% 
  #left_join(user_exam_f,by="仮個人id",suffix=c(".b",".f"))

user.data_8_detect<-
user.data_8 %>% 
  mutate(weight_nz.diff=weight_nz.f-weight_nz.b) %>%
  mutate(systolicBP_nz.diff=systolicBP_nz.f-systolicBP_nz.b) %>%
  mutate(diastolicBP_nz.diff=diastolicBP_nz.f-diastolicBP_nz.b) %>%
  mutate(ldl_c_nz.diff=ldl_c_nz.f-ldl_c_nz.b) %>%
  mutate(hdl_c_nz.diff=hdl_c_nz.f-hdl_c_nz.b) %>%
  mutate(TG_nz.diff=TG_nz.f-TG_nz.b) %>%
  mutate(HbA1c_nz.diff=HbA1c_nz.f-HbA1c_nz.b) 

user.data_8_detect %>% head()

条件9

baseline、followupで体重・収縮期血圧・拡張期血圧・LDL-c・HDL-c・中性脂肪・HbA1cすべてに欠損がない

user.data_9<-
  user.data_8_detect %>% 
  filter(!is.na(weight_nz.b)&!is.na(weight_nz.f)&
           !is.na(systolicBP_nz.b)&!is.na(systolicBP_nz.f)&
           !is.na(diastolicBP_nz.b)&!is.na(diastolicBP_nz.f)&
           !is.na(ldl_c_nz.b)&!is.na(ldl_c_nz.f)&
           !is.na(hdl_c_nz.b)&!is.na(hdl_c_nz.f)&
           !is.na(TG_nz.b)&!is.na(TG_nz.f)&
           !is.na(HbA1c_nz.b)&!is.na(HbA1c_nz.f))

user.data_9 %>% head()

歩数変化量の5分位点（quintile）で参加者を5つの集団に分け、集団ごとに1~5のラベルをつける

user.data_9 <-
  user.data_9 %>% 
  mutate(step.diff.quintile=case_when(
    step.diff < quantile(user.data_9$step.diff, 0.20) ~ 1,
    quantile(user.data_9$step.diff, 0.20) <= step.diff & 
    step.diff < quantile(user.data_9$step.diff, 0.40) ~ 2,
    quantile(user.data_9$step.diff, 0.40) <= step.diff & 
    step.diff < quantile(user.data_9$step.diff, 0.60) ~ 3,
    quantile(user.data_9$step.diff, 0.60) <= step.diff & 
    step.diff < quantile(user.data_9$step.diff, 0.80) ~ 4,
    quantile(user.data_9$step.diff, 0.80) <= step.diff  ~ 5,TRUE ~ 0))

quintileで区切られる各区分ごとに中央値を算出する。

# quintileで区切られる各区分ごとに中央値を算出する。
# データセット作成
q1 <- user.data_9  %>% filter(step.diff.quintile == 1) %>% select(step.diff)

## Adding missing grouping variables: `仮個人id`

q2 <- user.data_9  %>% filter(step.diff.quintile == 2) %>% select(step.diff)

## Adding missing grouping variables: `仮個人id`

q3 <- user.data_9  %>% filter(step.diff.quintile == 3) %>% select(step.diff)

## Adding missing grouping variables: `仮個人id`

q4 <- user.data_9  %>% filter(step.diff.quintile == 4) %>% select(step.diff)

## Adding missing grouping variables: `仮個人id`

q5 <- user.data_9  %>% filter(step.diff.quintile == 5) %>% select(step.diff)

## Adding missing grouping variables: `仮個人id`

各区分の中央値を格納した連続変数を作成する。

user.data_9 <-
  user.data_9 %>%
  mutate(
    step.diff.median.by.group =
      as.numeric(
        case_when(
          step.diff.quintile == 1 ~ median(q1$step.diff),
          step.diff.quintile == 2 ~ median(q2$step.diff),
          step.diff.quintile == 3 ~ median(q3$step.diff),
          step.diff.quintile == 4 ~ median(q4$step.diff),
          step.diff.quintile == 5 ~ median(q5$step.diff)
        )
      )
  )

条件7合致データと月次歩数データを元に、作図用月次歩数データを作成

user_monthly_data <- 
  user.data_7 %>% select(仮個人id, 性別_sfl, 年齢_sfl)

step.montyly.data.7 <- 
  step.monthly.data %>% 
  inner_join(user_monthly_data, by='仮個人id') 

step.montyly.data.7 %>% skim()

Data summary
Name	Piped data
Number of rows	302448
Number of columns	11
_______________________
Column type frequency:
character	5
numeric	6
________________________
Group variables	None