# ←この記号の後は次の改行までコメントになります。

# 今回の目標はふたつあります。
# 1. パッケージの導入を学ぶ。
# 2. tableone パッケージをインストールして、記述統計量を「一気に」求める。

# 1. 
# Rのコアはウィーン工科大学で開発されて配布されています( https://www.r-project.org/ )
# コアでも十分な機能を持っています。例えば、Lesson 1-11の内容は全てコアの内容のみを
# 用いています。
# さらに、世界中の研究者が様々な機能を追加するアドオンプログラムを開発しています。
# 例えば、生存分析を行うにはsurvivalパッケージを追加します。

# まずリポジトリ(repository)が正しく設定されているかを確認します。
# リポジトリとは、この場合はパッケージの倉庫です。
getOption("repos") # ここでhttps://cran.rstudio.com/ が表示されればOKです。
##     CRAN 
## "@CRAN@"
# このパッケージは、以前亀田の内科医だった吉田和樹先生(現・ハーバード大学公衆
# 衛生大学院)が開発したものです。
# まず、installed.packages()を実行してみて下さい。
# 今までもしパッケージをインストールしていたら、そのリストが表示されます。
# Rstudioの右下のウィンドウからPackagesタグをクリックしても同じようにパッケージ
# の一覧が出ます。こちらのほうが説明も付いていて便利ですね。

# 2. tableone パッケージをインストールして、記述統計量を「一気に」求める。

# tableone パッケージをインストールします。
# dependencies=Tとは、tableone パッケージと依存関係にある(tableone パッケージが
# 利用しているパッケージも併せてインストールするオプションです。通常はT(RUE)にし
# ます。下記はコメント行にしていますが、コメント行の記号(先頭の#)を外して一度
# は実行して下さい。
# install.packages("tableone",dependencies=T)   

# もう一度installed.packages()を実行するか右下ウィンドウのPackageタグをクリック
# すると、tableoneとその依存パッケージが増えているのがわかります。

# パッケージは時々updateされます。定期的に下記のコマンドを実行して下さい。
# update.packages()

# 再びデータセット作成から
# まずはデータ空間をクリアしておいてから、
rm(list=ls())
# 第3回講習会で使ったデータセットを使います。
# ここから77行目まではサンプルデータセットの作成です。
set.seed(1111)
Dataset.0<-
  data.frame(
    Sex=factor(rnorm(n=800,mean=0.5,sd=1)>0,labels=c("Male","Female")),
    Age=50+rpois(n=800,lambda=25)+round(rnorm(n=800,mean=0,sd=5),0),
    Severity=rpois(n=800,lambda=3))
Dataset.0$Comorbidity=
  rpois(n=800,lambda=1)+round(Dataset.0$Age*rnorm(n=800,mean=0.03,sd=0.08)^2,0)
Dataset.0$UnknownConfounder=
  factor(
    with(Dataset.0,
         rnorm(n=800,mean=0)+
           I(Sex=="Male")*rnorm(n=800,mean=0.02)*0.01+
           I(Age-80)*rnorm(n=800,mean=0.002)*0.001+
           I(Severity>3)*rnorm(n=800,mean=0.05)*0.01)>1,
    labels=c("No","Yes"))
Dataset.0$Intervention<-
  factor(
    with(Dataset.0,
         rnorm(n=800,mean=0)+
           I(Sex=="Female")*rnorm(n=800,mean=0.2)+
           I(Age-80)*rnorm(n=800,mean=-0.1)*0.1+
           I(Severity-2)*rnorm(n=800,mean=0.4)+
           I(Comorbidity-1)*rnorm(n=800,mean=-1.2)*0.6+
           I(UnknownConfounder=="Yes")*rnorm(n=800,mean=0.2)*0.2)>0,
    labels=c("No","Yes"))
Dataset.0$Outcome<-
  factor(
    with(Dataset.0,
         rnorm(n=800,mean=0)*10+
           I(Sex=="Male")*rnorm(n=800,mean=0.01)+
           I(Age-75)*rnorm(n=800,mean=0.3)*5+
           I(Severity-1)*rnorm(n=800,mean=0.7)*6+
           I(Comorbidity-1)*rnorm(n=800,mean=0.9)*5+
           I(Intervention=="Yes")*rnorm(n=800,mean=-0.6)*4+
           I(UnknownConfounder=="Yes")*rnorm(n=800,mean=0.3)*0.7)>20,labels=c("No","Yes"))
Dataset.0$UnknownConfounder<-NULL

# まずはいきなり全ての記述統計量を含むTable 1をつくってみましょう。
# 最初にパッケージを読み込みます。library関数を使います。
library(tableone)

# Table1を作るには、CreateTableOne関数を使ってデータセットを指定するだけです。非常に簡単です。

Table.1<-
  CreateTableOne(data=Dataset.0)
print(Table.1)
##                          
##                           Overall      
##   n                         800        
##   Sex = Female (%)          546 (68.2) 
##   Age (mean (sd))         75.02 (7.15) 
##   Severity (mean (sd))     3.03 (1.74) 
##   Comorbidity (mean (sd))  1.48 (1.30) 
##   Intervention = Yes (%)    429 (53.6) 
##   Outcome = Yes (%)         289 (36.1)
# 群わけをします。介入(Intervention)で分けてみましょう。
Table.1<-
  CreateTableOne(data=Dataset.0,strata="Intervention")
print(Table.1)
##                          Stratified by Intervention
##                           No            Yes            p      test
##   n                         371           429                     
##   Sex = Female (%)          256 (69.0)    290 ( 67.6)   0.727     
##   Age (mean (sd))         74.97 (7.32)  75.07 (7.01)    0.837     
##   Severity (mean (sd))     2.77 (1.68)   3.25 (1.75)   <0.001     
##   Comorbidity (mean (sd))  1.87 (1.38)   1.14 (1.13)   <0.001     
##   Intervention = Yes (%)      0 ( 0.0)    429 (100.0)  <0.001     
##   Outcome = Yes (%)         149 (40.2)    140 ( 32.6)   0.033
# 2群にわけたので、検定も追加されました。

# InterventionはもはやTable.1には必要ありません。また、転帰(Outcome)はベースライン
# 変数ではないので除きます。
# CreateTableOneでは、varsに何も指定しないと全変数の記述統計量を出しますが、指定すると
# 指定した変数だけのものになります。
Table.1<-
  CreateTableOne(
    data=Dataset.0,
    vars=c("Sex","Age","Severity","Comorbidity"),
    strata="Intervention")
print(Table.1)
##                          Stratified by Intervention
##                           No            Yes           p      test
##   n                         371           429                    
##   Sex = Female (%)          256 (69.0)    290 (67.6)   0.727     
##   Age (mean (sd))         74.97 (7.32)  75.07 (7.01)   0.837     
##   Severity (mean (sd))     2.77 (1.68)   3.25 (1.75)  <0.001     
##   Comorbidity (mean (sd))  1.87 (1.38)   1.14 (1.13)  <0.001
# 次にprint関数(実際にはprint.TableOne関数)のオプションを使ってみます。
# contDigitsに0を指定して、有効数字を0桁(整数表示)にしてみます。
Table.1<-
  CreateTableOne(
    data=Dataset.0,
    vars=c("Sex","Age","Severity","Comorbidity"),
    strata="Intervention")
print(Table.1,contDigits=0)
##                          Stratified by Intervention
##                           No          Yes         p      test
##   n                       371         429                    
##   Sex = Female (%)        256 (69.0)  290 (67.6)   0.727     
##   Age (mean (sd))          75 (7)      75 (7)      0.837     
##   Severity (mean (sd))      3 (2)       3 (2)     <0.001     
##   Comorbidity (mean (sd))   2 (1)       1 (1)     <0.001
# 重症度と合併症は非正規分布を仮定します。
Table.1<-
  CreateTableOne(
    data=Dataset.0,
    vars=c("Sex","Age","Severity","Comorbidity"),
    strata="Intervention")
print(Table.1,contDigits=0,nonnormal=c("Severity","Comorbidity"))
##                             Stratified by Intervention
##                              No          Yes         p      test   
##   n                          371         429                       
##   Sex = Female (%)           256 (69.0)  290 (67.6)   0.727        
##   Age (mean (sd))             75 (7)      75 (7)      0.837        
##   Severity (median [IQR])      3 [2, 4]    3 [2, 4]  <0.001 nonnorm
##   Comorbidity (median [IQR])   2 [1, 3]    1 [0, 2]  <0.001 nonnorm
# この結果をcsvファイルに保存することもできます。
# 保存してものを表計算ソフトで開いてみて下さい。
write.csv(print(Table.1,contDigits=0,nonnormal=c("Severity","Comorbidity")),"Table1.csv")
##                             Stratified by Intervention
##                              No          Yes         p      test   
##   n                          371         429                       
##   Sex = Female (%)           256 (69.0)  290 (67.6)   0.727        
##   Age (mean (sd))             75 (7)      75 (7)      0.837        
##   Severity (median [IQR])      3 [2, 4]    3 [2, 4]  <0.001 nonnorm
##   Comorbidity (median [IQR])   2 [1, 3]    1 [0, 2]  <0.001 nonnorm
# もっと詳しいことが知りたければ、ヘルプファイルを参照して下さい。
help(CreateTableOne)
help(print.TableOne)