/**** ***** ***** ***** ***** ***** ***** ***** ***** *****
*
* 日本疫学会プレセミナー
* Rで実践！ 美しい Figure ＆ Table を作成しよう
* Stataでの実行に関するフォロー
* by Toshiharu Mitsuhashi @sankyoh
* Drawing Part2 analysis and graph
*
***** ***** ***** ***** ***** ***** ***** ***** ***** ****/

// Part1で作ったdtaファイルの読込みです。
use nhefs, clear

// グラフテーマを設定します。ここで設定すると、後のグラフコマンドでは標準でwhite_tableauになります。
// ただし、今回は、個々のグラフコマンドにもテーマの設定を行なっています。
set scheme white_tableau

// Stataメモリ上の既存グラフを削除
capture graph drop *

/***************************************************************************
*
* Section 2 Data
*
***************************************************************************/

des
sum
list in 1/10

/***************************************************************************
*
* Section 3 Overview & Missing values（欠測値）
*
***************************************************************************/

////////// ////////// ////////// //////////
* Section 3.1 Overview
////////// ////////// ////////// //////////

misstable sum
misstable pattern
misstable tree
misstable nested

////////// ////////// ////////// //////////
* Section 3.2 Missingness at each variable
////////// ////////// ////////// //////////

qui: do missing_lollipop_plot // 定義したProgramの読込み（近日中にado fileにする予定）
missing_lollipop_plot

////////// ////////// ////////// //////////
* Section 3.3 Missingness relationship between variables
////////// ////////// ////////// //////////

///// 3.3.1 全ての変数 /////
* 欠損値プロット
missingplot, scheme(white_tableau) msymbol(oh)

///// 3.3.2 連続変数 × 名義変数 /////
// 生存／死亡におけるalcholohowmachの欠損値の個数
misstable sum alcoholhowmuch if death==0
misstable sum alcoholhowmuch if death==1

///// 3.3.3 連続変数 × 連続変数 /////
* 2つの連続変数間について、一方が欠損の場合にもう一方の変数の分布がどのようになっているかを確認
gen cholesterol2 = cholesterol
recode cholesterol2 (.=0)
capture graph drop *

// 生存者におけるグラフ
twoway (scatter cholesterol2 pack_year_n if death ==0 & cholesterol!=., jitter(2 2)) ///
	(scatter cholesterol2 pack_year_n if death ==0 & cholesterol==., jitter(2 2) ///
	ylabel(0(100)400) xlabel(0(50)150) scheme(white_tableau) legend(off) title("Alive") name(foobar1, replace))

// 死亡者におけるグラフ
twoway (scatter cholesterol2 pack_year_n if death ==1 & cholesterol!=., jitter(2 2)) ///
	(scatter cholesterol2 pack_year_n if death ==1 & cholesterol==., jitter(2 2) ///
	ylabel(0(100)400) xlabel(0(50)150) scheme(white_tableau) legend(off) title("Death") name(foobar2, replace))

// 2つのグラフを結合
graph combine foobar1 foobar2, ///
	cols(2) name(chol_smk, replace) scheme(white_tableau)

/***************************************************************************
*
* Section 4 Continuous variable (連続変数)
*
***************************************************************************/

////////// ////////// ////////// //////////
* Section 4.1 ヒストグラム
////////// ////////// ////////// //////////

twoway histogram pack_year_n, ///
	start(0) width(5) freq ///
	scheme(white_tableau) ///
	xtitle("Pack-Year") ///
	name(hist, replace)

////////// ////////// ////////// //////////
* Section 4.2 密度関数プロット
////////// ////////// ////////// //////////

twoway (kdensity pack_year_n if death==0) ///
	(kdensity pack_year_n if death==1, ///
	scheme(white_tableau) ///
	legend(label(1 "Alive") label(2 "Death")) ///
	xtitle("Pack-Year") ///
	ytitle("Density") ///
	name(dens, replace))

////////// ////////// ////////// //////////
* Section 4.3 箱ひげ図
////////// ////////// ////////// //////////

* 箱ひげ図(overオプションを使うとき)
graph box pack_year_n , ///
	scheme(white_tableau) ///
	ytitle("Pack-Year") ///
	over(death) ///
	name(box_over, replace)

* 箱ひげ図(byオプションを使うとき)
graph box pack_year_n , ///
	scheme(white_tableau) ///
	ytitle("Pack-Year") ///
	by(death) ///
	name(box_by, replace)

* 箱ひげ図にドットプロットを載せる
// ssc install stripplot
stripplot pack_year_n, ///
	scheme(white_tableau) ///
	vertical ///
	ytitle("Pack-Year") ///
	over(death) ///
	box ///
	jitter(30 10) ///
	name(box_dot, replace)

////////// ////////// ////////// //////////
* Section 4.4 ドットプロット
////////// ////////// ////////// //////////

stripplot pack_year_n, scheme(white_tableau) over(death) vertical name(dotplot1, replace)
stripplot pack_year_n, scheme(white_tableau) over(death) vertical msize(0.5) name(dotplot2, replace)
stripplot pack_year_n, scheme(white_tableau) over(death) vertical jitter(5 5) center name(dotplot3, replace)

* Beeswarmっぽいが、Beeswarmではない。
stripplot pack_year_n, scheme(white_tableau) over(death) vertical center stack name(beeswarm, replace)

* stripplotを使わない場合
dotplot pack_year_n , over(death)
dotplot pack_year_n , over(death) msize(0.5)
dotplot pack_year_n , over(death) jitter(5 5) center

////////// ////////// ////////// //////////
* Section 4.5 レインクラウドプロット
////////// ////////// ////////// //////////

// 簡易に行う方法を見つける事が出来ませんでした。
// Statalistに下記スレッドがあります。
// https://www.statalist.org/forums/forum/general-stata-discussion/general/1690233-raincloud-plot-in-stata

/***************************************************************************
*
* Section 5 Nominal variable（名義変数）
*
***************************************************************************/

////////// ////////// ////////// //////////
* Sectin 5.1 棒グラフ
////////// ////////// ////////// //////////

///// 5.1.1 ラベルなし /////
// seqnでなくても良いが、欠損値のない数が良い。
graph bar (count) seqn, scheme(white_tableau) over(marital) ytitle("count") name(bar1, replace)
fre marital // 度数分布表

///// 5.1.2 ラベルあり /////
label define marital 1 "Under17" 2 "Married" 3 "Widowed" 4 "Never" 5 "Divorced" 6 "Separated" 8 "Unknown"
label values marital marital
graph bar (count) seqn, scheme(white_tableau) over(marital) ytitle("count") name(bar2, replace)
fre marital // 度数分布表

////////// ////////// ////////// //////////
* Section 5.2　横並び棒グラフ
////////// ////////// ////////// //////////

graph bar (count) seqn, scheme(white_tableau) over(death) over(marital) ytitle("count") name(bar3, replace)
tab marital death

////////// ////////// ////////// //////////
* Section 5.3 100%積み上げ棒グラフ
////////// ////////// ////////// //////////

label define education 5 "College or more" 4 "College dropout" 3 "HS" 2 "HS dropout" 1 "8th grade or less"
label values education education

bysort education: egen frac_1 = sum(death)
bysort education: gen frac_3 = _N
gen frac   = frac_1/frac_3
gen frac_2 = frac_3 - frac_1

graph hbar (mean) frac_2 frac_1, over(education, sort(frac)) stack percent legend(label(1 "Alive") label(2 "Death")) name(prop_bar, replace)

/***************************************************************************
*
* Section Multivariable（多変量）
*
***************************************************************************/

graph matrix sbp dbp age income, half msymbol(smcircle_hollow) name(multi, replace)

exit