- dput( )で変数名を書き出す
- 表示したい変数名を指定する
- カテゴリー変数として扱いたいものを指定する
- 表データを作成する
- 連続変数の詳細結果を表示する
- カテゴリー変数の詳細結果を表示する
- 体裁を整えて結果を表示する
- おわりに
- 参考資料
ここではsurvivalパッケージに入っているMayo Clinicの原発性胆汁性肝硬変(Primary Biliary Cirrhosis, PBC)のデータを使う。
> library(survival) > str(pbc) 'data.frame': 418 obs. of 20 variables: $ id : int 1 2 3 4 5 6 7 8 9 10 ... $ time : int 400 4500 1012 1925 1504 2503 1832 2466 2400 51 ... $ status : int 2 0 2 2 1 2 0 2 2 2 ... $ trt : int 1 1 1 1 2 2 2 2 1 2 ... $ age : num 58.8 56.4 70.1 54.7 38.1 ... $ sex : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ... $ ascites : int 1 0 0 0 0 0 0 0 0 1 ... $ hepato : int 1 1 0 1 1 1 1 0 0 0 ... $ spiders : int 1 1 0 1 1 0 0 0 1 1 ... $ edema : num 1 0 0.5 0.5 0 0 0 0 0 1 ... $ bili : num 14.5 1.1 1.4 1.8 3.4 0.8 1 0.3 3.2 12.6 ... $ chol : int 261 302 176 244 279 248 322 280 562 200 ... $ albumin : num 2.6 4.14 3.48 2.54 3.53 3.98 4.09 4 3.08 2.74 ... $ copper : int 156 54 210 64 143 50 52 52 79 140 ... $ alk.phos: num 1718 7395 516 6122 671 ... $ ast : num 137.9 113.5 96.1 60.6 113.2 ... $ trig : int 172 88 55 92 72 63 213 189 88 143 ... $ platelet: int 190 221 151 183 136 NA 204 373 251 302 ... $ protime : num 12.2 10.6 12 10.3 10.9 11 9.7 11 11 11.5 ... $ stage : int 4 3 4 4 3 3 3 3 2 4 ...
dput( )で変数名を書き出す
変数名を書き写すのが面倒なのでdput()を使う。
> dput(names(pbc)) c("id", "time", "status", "trt", "age", "sex", "ascites", "hepato", "spiders", "edema", "bili", "chol", "albumin", "copper", "alk.phos", "ast", "trig", "platelet", "protime", "stage")
表示された結果を次にコピペで使う。
表示したい変数名を指定する
all_vars <- c("time", "status", "age", "sex", "ascites", "hepato", "spiders", "edema", "bili", "chol", "albumin", "copper", "alk.phos", "ast", "trig", "platelet", "protime", "stage")
患者IDは不要なのでidを除外した。治療別に集計するのでtrtも除外。
カテゴリー変数として扱いたいものを指定する
- 文字列は自動的にカテゴリー変数扱いされる。
cat_vars <- c("status", "sex", "ascites", "hepato", "spiders", "edema", "stage")
表データを作成する
以下の引数を指定する。
vars：表示する変数factorVars：カテゴリー変数として扱いたい変数strata：比較したい群includeNA：TRUEとすると、カテゴリー変数で欠測も1つのカテゴリーとして扱うtestApprox：カテゴリー変数でサンプル数が多い時に用いる近似を使った検定方法を指定。デフォルトはchisq.test（= カイ2乗検定）argsApprox：上記の検定を行うときの設定。デフォルトはlist(correct = TRUE)となっており、カイ2乗検定で連続値補正を行う。testExact：カテゴリー変数でサンプル数が少ないときに用いる正確な検定方法を指定。デフォルトはfisher.test（= Fisher正確検定）testNormal：正規分布を仮定したときに用いられる検定方法。デフォルトはoneway.test（= 一元配置分散分析）で、これは2群比較ならばt.test（= t検定）に相当argsNormal：上記の検定を行うときの設定。デフォルトはlist(var.equal = TRUE)となっており、等分散を仮定した検定を行う。testNonNormal：正規分布を仮定しないときに用いられる検定方法。デフォルトはkruskal.test（= Kruskal-Wallis検定）で、これは2群比較ならばwilcox.test（= Mann-Whitney U検定）に相当
library(tableone) table1 <- CreateTableOne(vars=all_vars, factorVars=cat_vars, strata="trt", data=pbc, includeNA=TRUE, argsNormal=list(var.equal=FALSE))
連続変数の比較では等分散を仮定しない方がいいかも。
> table1 Stratified by trt 1 2 p test n 158 154 time (mean (SD)) 2015.62 (1094.12) 1996.86 (1155.93) 0.883 status (%) 0.894 0 83 (52.5) 85 (55.2) 1 10 ( 6.3) 9 ( 5.8) 2 65 (41.1) 60 (39.0) age (mean (SD)) 51.42 (11.01) 48.58 (9.96) 0.018 sex = f (%) 137 (86.7) 139 (90.3) 0.421 ascites (%) NaN 0 144 (91.1) 144 (93.5) 1 14 ( 8.9) 10 ( 6.5) NA 0 ( 0.0) 0 ( 0.0) hepato (%) NaN 0 85 (53.8) 67 (43.5) 1 73 (46.2) 87 (56.5) NA 0 ( 0.0) 0 ( 0.0) spiders (%) NaN 0 113 (71.5) 109 (70.8) 1 45 (28.5) 45 (29.2) NA 0 ( 0.0) 0 ( 0.0) edema (%) 0.877 0 132 (83.5) 131 (85.1) 0.5 16 (10.1) 13 ( 8.4) 1 10 ( 6.3) 10 ( 6.5) bili (mean (SD)) 2.87 (3.63) 3.65 (5.28) 0.133 chol (mean (SD)) 365.01 (209.54) 373.88 (252.48) 0.747 albumin (mean (SD)) 3.52 (0.44) 3.52 (0.40) 0.874 copper (mean (SD)) 97.64 (90.59) 97.65 (80.49) 0.999 alk.phos (mean (SD)) 2021.30 (2183.44) 1943.01 (2101.69) 0.747 ast (mean (SD)) 120.21 (54.52) 124.97 (58.93) 0.460 trig (mean (SD)) 124.14 (71.54) 125.25 (58.52) 0.886 platelet (mean (SD)) 258.75 (100.32) 265.20 (90.73) 0.554 protime (mean (SD)) 10.65 (0.85) 10.80 (1.14) 0.199 stage (%) NaN 1 12 ( 7.6) 4 ( 2.6) 2 35 (22.2) 32 (20.8) 3 56 (35.4) 64 (41.6) 4 55 (34.8) 54 (35.1) NA 0 ( 0.0) 0 ( 0.0)
連続変数の詳細結果を表示する
summary(table1$ContTable)
strata別に以下の項目：
n：観察数miss：欠測数p.miss：欠測割合mean,sd,median：平均, 標準偏差, 中央値p25,p75：第1四分位, 第3四分位min,max：最小値, 最大値skew：歪度(skewness)=分布の非対称性を示す指標kurt：尖度(kurtosis)=分布のピークの鋭さの指標
strata間を比較したときの以下の項目：
p-values：P値pNormal：正規分布が仮定できるときの検定結果pNonNormal：正規分布を仮定しないときの検定結果
Standardized mean difference：標準化平均差
カテゴリー変数の詳細結果を表示する
summary(table1$CatTable)
strata別に以下の項目（連続変数と重複している項目は省略）：
level：カテゴリー内の水準freq：度数percent：割合cum.percent：累積割合
strata間を比較したときの以下の項目：
p-values：P値pApprox：サンプル数が大きいときに行う検定pExact：サンプル数が小さいときに行う検定
体裁を整えて結果を表示する
以下の引数を指定する。
cramVars：2値変数で両方のレベルについて表示したい変数を指定。exact：正確検定の結果を表示したいカテゴリー変数を指定nonnormal：正規分布に従わない変数を指定。正規分布を仮定しない検定で求められたP値が用いられるほか、mean(SD)ではなくmedian[IQR]で記述要約される。test：FALSEとすると検定結果（つまりP値）が表示されないdropEqual：TRUEとすると、2値変数で2番目のレベル名が表示されないminMax：TRUEとするとnonnormalに指定した変数でIQRではなく最小-最大が表示される。
q <- print(table1, cramVars="sex", exact=c("ascites","hepato","spiders","stage"), nonnormal=c("time","alk.phos"))
csv形式で書き出すことも可能。
write.csv(q, "table1.csv")
おわりに
dput(names(data))を知って変数名を抽出すると便利だと知りました。- キャットタワーが部屋を占拠しはじめてます。
参考資料
- パッケージ開発者の吉田 和樹先生の解説ページ（英語） cran.r-project.org
