Forest plotをStataで作成してみる

グラフ

このブログでは、統計解析ソフトStataのプログラミングのTipsや便利コマンドを紹介しています.

Facebook groupでは、ちょっとした疑問や気づいたことなどを共有して貰うフォーラムになっています. ブログと合わせて個人の学習に役立てて貰えれば幸いです.

さて、以前の記事で交互作用の一つの表現方法としてForest plotをちらっと説明しました.

今回はそれを描くためのプログラムをご紹介します.メタ解析の結果を表示するための forestplot とは別に、回帰分析などの出力結果を単純にグラフ化する方法をお示しします.

1.回帰分析の結果をinputする

必要なデータを手入力します.input – endで挟んでいきます.

inputのあとは変数名を並べていきますが、numericの場合は変数の種類を指定する必要はありません.しかし、stringの場合にはstr6などのような記号を変数名の前につけます.

そしてstringとわかるように入力する値は””で括ります.

input group_num coef lci uci
	 7 2.55 1.96 3.32
	 6 1.47 1.14 1.90
	 5 1.31 0.93 1.84
	 4 2.58 1.73 3.84
	 3 1.94 1.15 3.25
	 2 4.28 3.53 5.20
	 1 2.73 2.05 3.63
end 

これに群の名称を入れるなどすることもできます.しかし残念ながらこれを自動的にグラフ上に表示することはできません.

input str6 groupnm group_num coef lci uci
	"group1" 7 2.55 1.96 3.32
	"group2" 6 1.47 1.14 1.90
	"group3" 5 1.31 0.93 1.84
	"group4" 4 2.58 1.73 3.84
	"group5" 3 1.94 1.15 3.25
	"group6" 2 4.28 3.53 5.20
	"group7" 1 2.73 2.05 3.63
end 

2.twowayのグラフで表現する

データをインプットした後は、グラフを作るのですが、肝となるのは、scatterとrcapです.scatterは点推定値の値にマーカーを置くことですが、rcapで信頼区間の上下限を結んで行くことができます.

.rcap 下限値の変数 上限値の変数 y軸の値

といれて、horizontalのオプションをつけることで縦に横棒を並べることができます.

twoway (rcap lci uci group_num, horizontal) || (scatter group_num coef)

これが原型になります.ここからの加工は、

  1. グラフの色を統一
  2. legend, y軸タイトルの除去
  3. y軸ラベルの名称変更
  4. x軸を対数軸に
  5. 背景を白く

いろいろつけるときは、delimitですっきりした見た目にしておくといいですね.改行のことはこちらで復習してください.ここでは1と2をやってみます.

#delimit;
 twoway 
   rcap lci uci group_num, horizontal lcolor(black)
   || 
   scatter group_num coef, mcolor(black) 
   legend(off) 
ytitle("")
   ;
#delimit cr 

rcapのオプションでlcolor(black)としているのですが、ytitle()のところを””としてやるとblankになります.legendはlegend(off)とすると何も表示されません.

このあとは縦軸のラベルにそれぞれのグループ名をいれることにします.残念ながら軸ラベルオプションで一つずつ入れていくしかありません.ここでラベルを横向きに表示させるには、angleを0度にしておくことが必要で、ang(0)というオプションが必要になります.軸ラベルの文字の大きさはこの中でlabsize()で指定できます.後々グリッドラインが邪魔になるのでここで消しておきましょう.

#delimit;
 twoway 
   rcap lci uci group_num, horizontal lcolor(black)
 
   ylabel(1 "Group 7" 2 "Group 6" 3 "Group 5" 4 "Group 4" 5 "Group 3" 6 "Group 2" 7 "Group 1", 
   ang(0) labsize(small) nogrid noticks) 
   || 
   scatter group_num coef
, mcolor(black)
   legend(off) 
ytitle("")  
   ;
#delimit cr 

さらにここからX軸を対数軸にして0.5~5までの範囲指定を行います.

#delimit;
 twoway 
   rcap lci uci group_num, horizontal lcolor(black)
 
   ylabel(1 "Group 7" 2 "Group 6" 3 "Group 5" 4 "Group 4" 5 "Group 3" 6 "Group 2" 7 "Group 1", 
   ang(0) labsize(small) nogrid noticks) 
   || 
   scatter group_num coef
, mcolor(black)
   legend(off) 
ytitle("")  
   xscale(range(0.5 5) log)
   xlabel(0.5 1 2 5, format(%9.1f) labsize(small)) xtick(0.5(0.5)5) 
   xline(1,lpattern(dash) lcolor(black))
   ;
#delimit cr 

最後に背景を白くします.graphregionで白くすると言えばよいだけです.

#delimit;
 twoway 
   rcap lci uci group_num, horizontal lcolor(black)
 
   ylabel(1 "Group 7" 2 "Group 6" 3 "Group 5" 4 "Group 4" 5 "Group 3" 6 "Group 2" 7 "Group 1", 
   ang(0) labsize(small) nogrid noticks) 
   || 
   scatter group_num coef
, mcolor(black)
   legend(off) 
ytitle("")  
   xscale(range(0.5 5) log)
   xlabel(0.5 1 2 5, format(%9.1f) labsize(small)) xtick(0.5(0.5)5) 
   xline(1,lpattern(dash) lcolor(black))
   graphregion(color(white) lcolor(white))
   ;
#delimit cr 

自分の場合は最終的に表と組み合わせたりということを考える場合にはPPTにできあがったグラフを貼り付けて、その上から加工しています.

3.”meta” 機能を活用する

これはFacebook groupの方で今村先生から教えていただきました.あるメタ解析の論文の結果を再現したものになります.まずは以下の様にデータを入れていただきます.(発表された論文をもとに数値を再現しています.転記ミス等あるかもしれませんがご容赦ください)

clear
input str24 study N rr ll ul c
"Wakai, 2005" 129	0.63	0.38	1.04	129
"Vatten, 1990" 122	0.7	0.44	1.11	122
"Gago-Dominguez premenopausal, 2003" 221	0.71	0.49	1.02	221
"Gago-Dominguez postmenopausal, 2003" 93	0.89	0.48	1.66	93
"Cho, 2003" 714	0.92	0.73	1.15	714
"Folsom, 2004" 1885	0.92	0.76	1.12	1885
"Holmes postmenopausal, 2003" 2936	1	0.89	1.12	2936
"Toniolo, 1994" 180	1.02	0.61	1.71	180
"Key, 1999" 376	1.05	0.82	1.35	376
"Engeset premenopausal, 2006" 2700	1.1	0.95	1.28	2700
"Engeset postmenopausal, 2006" 786	1.11	0.84	1.46	786
"Holmes premenopausal, 2003" 854	1.17	0.92	1.49	854
"Stripp, 2003" 424	1.47	1.1	1.97	424
"Mills, 1989" 207	1.54	1.14	2.07	207
end 

これをメタ解析用のプログラムに入れ込むだけです.meta setとするのが第一歩なのですが、その前に推定値と標準誤差を入手しておきます.ここではごく簡単に、SEは信頼区間の上限と下限の差を1.96の2倍で割ったものとして設定してしまいます.

gen b = ln(rr)
gen se = (ln(ul)-ln(ll))/(2*1.96)

meta set b se, random(dlaird) studylabel(study) 

#delimit;
  meta forestplot  _id  c _plot _es _esci _weight,  /* 表示するカラム. 研究名、N数、プロット、推定値、重み */
  random(dlaird) eform  sort(b)  /* 推定値を得る */
  crop(0.6 1.6) xtitle("Relative risk") xscale(range(0.6 1.6)) 
  xlabel(0.6 0.8 1.0 1.3 1.6, format(%6.1f))  /* X軸の指定 */ 
  nullrefline nohrule noohomtest noghetstats nogwhomtests nogbhomtests noosigtest nonotes 
                     /* 表示するエレメントのon/off指定 */
  columnopts(_id, title("{bf:Study}"))     /* カラムのタイトル(第1) */
  columnopts(c, supertitle("{bf:N}") title("{bf:cases}") )    /* カラムのタイトル(第2) */
  columnopts(_es _esci, format(%6.2f)  supertitle("{bf:Relative risk}") title("{bf:(95% CI)}")) 
                     /* カラムのタイトル(推定値) */
  columnopts(_weight, format(%6.1f) supertitle("{bf:Weight}") title("{bf:%}")) /* カラムのタイトル(重み)
 */
  omarkeropts(mcolor(black))   markeropts(mcolor(gs6))  insidemarkeropts(mcolor(black)) /* マーカーの指定 */
  cibind(parentheses) ciopts(color(gs6)) /* 信頼区間表示方法 
 */
  title("{bf:Fish consumption and breast cancer risk}") /* グラフタイトル */
  ;
#delimit cr 

これで以下の様な綺麗なグラフが描けると思います.

まとめ

Forest plotの描き方について説明しました.メタ解析の全貌については項を改めていずれ説明したいと思います.(勉強しないといけないな~)

コメント

タイトルとURLをコピーしました