R lecture lesson2

動画https://web.microsoftstream.com/video/eeddb908-3629-46ef-a114-8f87f451dca9
資料https://goo.gl/KVGoYH

今日の目標

Table1の作り方と記述統計
- Rへのデータの取り扱いになれる
- 関数を理解する
- 記述統計
- 検定を実行できる

参考資料

Cousera　Data science specialization
https://www.coursera.org/specializations/jhu-data-science
R tips　Rの基本文法など　日本語でわかり安い
http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html
R seek　Rに絞った検索エンジン
http://seekr.jp/
Stack overflow プログラミング界のYahoo知恵袋，自分が躓く問題は誰かが経験している
https://ja.stackoverflow.com/
R pubs
https://rpubs.com/
他
-help()ですんなり解決することが多い
-Rで出力されたエラーをGoogleで直接検索

事前準備

header=T: 最初の行名を列名として指定する．Tを選択，省力しても良い
sep=“,”:列の区切りを“,”と指定している．省略しても良い

#Directoryを確認する
getwd()

## [1] "C:/Users/ryohe/R lecture by Shiraishi"

#使用するcsvをdirectryに入れる
RAdata.import<-file("RA-data-example.csv",encoding="CP932")
RAdata.0<-data.frame(read.csv(RAdata.import,header=T,sep=","))
#演習用データ読み込み
Affairs<-read.csv("affairs.csv")

多変量データの解析

複数の変数が含まれるデータセット
臨床疫学研究では検証に用いる変数の数が多い
複数の変数間の関係に興味がある
- アウトカムー治療
- アウトカムー患者背景
- 患者背景を調整した上でもアウトカムー治療

まず最初にやるのはどんなデータか覗いてみる

#dimentionをみる．データフレームの行，列の大きさを表示
dim(RAdata.0)

## [1] 13285   206

#データフレームの列の名前を表示
colnames(RAdata.0, )

#最初の5行を表示
head(RAdata.0)

#様々な変数にヒストグラムを書く
hist(RAdata.0$age)

#変数のテーブルを表示
table(RAdata.0$Steroid)

## 
##            continuous         no  on demand 
##         59       5723       7316        187

table(RAdata.0$tuberculosis,RAdata.0$Steroid)

##      
##            continuous   no on demand
##          0          0    1         0
##   no    59       5717 7312       187
##   yes    0          6    3         0

Table1の作成

患者背景に必要な情報
* 対象症例数
* 各変数の平均値，割合，頻度，中央値，四分位範囲など
* 各変数の欠損
* 治療割付に関連しそうな因子
* イベントに関連しそうな因子

TableOnep Package

患者背景テーブルを作成するためのパッケージ
亀田総合病院，現ハーバード公衆衛生の吉田和樹先生が作成

#TableOneをインストール
install.packages("tableone", dependencies = T)

#パッケージの呼び出し
library(tableone)

パッケージ名は”"でくくる dependencies=T:他の必要なパッケージも読み込む

パッケージを利用する

CreateTableOne()関数を用いる
- data＝データフレーム名
- vars＝小計した変数名を列挙したベクトル
- strata＝層

library(tableone)

## Warning: package 'tableone' was built under R version 3.5.3

#表示する変数リストを作成
varlist<-c("age","sex","tjc","sjc","DAS28","DAS28CRP","tjc28","sjc28","SDAI","CDAI","RF","ACPA")
#ステージを層としてvarlistにある変数を要約する
CreateTableOne(data=RAdata.0, strata = "stage",var=varlist)

##                       Stratified by stage
##                        1               2               3              
##   n                      2634            3363            2258         
##   age (mean (SD))       59.96 (13.87)   63.54 (13.10)   65.58 (12.16) 
##   sex = M (%)             663 (25.2)      796 (23.7)      437 (19.4)  
##   tjc (mean (SD))        1.90 (4.27)     1.97 (4.20)     2.67 (4.78)  
##   sjc (mean (SD))        1.07 (2.23)     1.44 (2.75)     2.31 (3.56)  
##   DAS28 (mean (SD))      2.70 (1.18)     2.92 (1.19)     3.27 (1.22)  
##   DAS28CRP (mean (SD))   2.16 (1.00)     2.32 (1.03)     2.62 (1.07)  
##   tjc28 (mean (SD))      1.42 (3.09)     1.44 (2.93)     1.95 (3.59)  
##   sjc28 (mean (SD))      0.90 (1.92)     1.18 (2.27)     1.86 (2.85)  
##   SDAI (mean (SD))       6.05 (6.92)     6.97 (7.15)     9.32 (7.88)  
##   CDAI (mean (SD))       5.63 (6.51)     6.43 (6.71)     8.72 (7.44)  
##   RF (mean (SD))        97.77 (253.55) 123.04 (303.90) 144.91 (305.55)
##   ACPA (mean (SD))     200.68 (540.33) 208.21 (415.00) 253.83 (566.62)
##                       Stratified by stage
##                        4               p      test
##   n                      3263                     
##   age (mean (SD))       66.35 (11.27)  <0.001     
##   sex = M (%)             341 (10.5)   <0.001     
##   tjc (mean (SD))        2.99 (5.11)   <0.001     
##   sjc (mean (SD))        1.94 (3.03)   <0.001     
##   DAS28 (mean (SD))      3.50 (1.22)   <0.001     
##   DAS28CRP (mean (SD))   2.79 (1.09)   <0.001     
##   tjc28 (mean (SD))      2.13 (3.71)   <0.001     
##   sjc28 (mean (SD))      1.65 (2.60)   <0.001     
##   SDAI (mean (SD))      10.16 (7.92)   <0.001     
##   CDAI (mean (SD))       9.39 (7.45)   <0.001     
##   RF (mean (SD))       141.34 (286.94) <0.001     
##   ACPA (mean (SD))     259.22 (452.24)  0.014

検定は自動選択されている
名義変数はχ二乗検定
連続変数はt検定
オプションでノンパラメトリックを選択肢マンホイットニーU検定に変更できる

演習

以下のデータを環境に取り込み解析を行う
* “affairs.csv”
+ 昨年1年間の不倫回数とその背景データ
+ 既婚者601名，欠損なし
* 各変数の平均，頻度
* 不倫の有無に影響していそうな因子は？

#演習用データ読み込み
Affairs<-read.csv("affairs.csv")

#データを確認
dim(Affairs)

## [1] 601   9

head(Affairs)

##   affairs gender age yearsmarried children religiousness education
## 1       0   male  37        10.00       no             3        18
## 2       0 female  27         4.00       no             4        14
## 3       0 female  32        15.00      yes             1        12
## 4       0   male  57        15.00      yes             5        18
## 5       0   male  22         0.75       no             2        17
## 6       0 female  32         1.50       no             2        17
##   occupation rating
## 1          7      4
## 2          6      4
## 3          1      4
## 4          6      5
## 5          6      3
## 6          5      5

#Summaryしてみる
summary(Affairs)

##     affairs          gender         age         yearsmarried    children 
##  Min.   : 0.000   female:315   Min.   :17.50   Min.   : 0.125   no :171  
##  1st Qu.: 0.000   male  :286   1st Qu.:27.00   1st Qu.: 4.000   yes:430  
##  Median : 0.000                Median :32.00   Median : 7.000            
##  Mean   : 1.456                Mean   :32.49   Mean   : 8.178            
##  3rd Qu.: 0.000                3rd Qu.:37.00   3rd Qu.:15.000            
##  Max.   :12.000                Max.   :57.00   Max.   :15.000            
##  religiousness     education       occupation        rating     
##  Min.   :1.000   Min.   : 9.00   Min.   :1.000   Min.   :1.000  
##  1st Qu.:2.000   1st Qu.:14.00   1st Qu.:3.000   1st Qu.:3.000  
##  Median :3.000   Median :16.00   Median :5.000   Median :4.000  
##  Mean   :3.116   Mean   :16.17   Mean   :4.195   Mean   :3.932  
##  3rd Qu.:4.000   3rd Qu.:18.00   3rd Qu.:6.000   3rd Qu.:5.000  
##  Max.   :5.000   Max.   :20.00   Max.   :7.000   Max.   :5.000

陥りやすいトラブル

連続良は平均値や中央値に意味がある，カテゴリカル変数は頻度に意味がる
変数がいずれであるかはある程度自動判別してくれるがされないときもある
occupationは職業をカテゴリカルにしたものだが，Rが連続変数を勘違いしている

データの列の型を確かめる必要がある

str(Affairs) #strで型を全て書き出す

## 'data.frame':    601 obs. of  9 variables:
##  $ affairs      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ gender       : Factor w/ 2 levels "female","male": 2 1 1 2 2 1 1 2 1 2 ...
##  $ age          : num  37 27 32 57 22 32 22 57 32 22 ...
##  $ yearsmarried : num  10 4 15 15 0.75 1.5 0.75 15 15 1.5 ...
##  $ children     : Factor w/ 2 levels "no","yes": 1 1 2 2 1 1 1 2 2 1 ...
##  $ religiousness: int  3 4 1 5 2 2 2 2 4 4 ...
##  $ education    : int  18 14 12 18 17 17 12 14 16 14 ...
##  $ occupation   : int  7 6 1 6 6 5 1 4 1 4 ...
##  $ rating       : int  4 4 4 5 3 5 3 4 2 5 ...

class(Affairs$occupation) #classは指定して書き出す

## [1] "integer"

EnviromentのDataからでも見れる

データの型の変換

occupationが連続変数として扱われていたため，カテゴリーに変換したい

summary(Affairs$occupation) #結果は連続変数

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   3.000   5.000   4.195   6.000   7.000

Affairs$occupation<-factor(Affairs$occupation) #因子に変換
summary(Affairs$occupation)　#カテゴリカルに記述されるように

##   1   2   3   4   5   6   7 
## 113  13  47  68 204 143  13

#affirsは何回？
table(Affairs$affairs)

## 
##   0   1   2   3   7  12 
## 451  34  17  19  42  38

#histgram
hist(Affairs$affairs)

#不倫をした人をカテゴリカルにしたい
Affairs$affairs #回数が表示されう

##   [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##  [24]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##  [47]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##  [70]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##  [93]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [116]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [139]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [162]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [185]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [208]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [231]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [254]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [277]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [300]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [323]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [346]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [369]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [392]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [415]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [438]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  3  3  7 12  1  1 12  7  2
## [461]  3  1  7 12 12 12  3  7  7  1  1  7  1 12 12  3  7  1  1  1  3  3  1
## [484]  1  7  7  7 12  7  7  1  2 12 12  1 12 12  7  7  1  3 12  7  1  7  1
## [507]  1  1  7  3  7  7  2  7  1  7  2  7  7  3  1  2 12  1  3 12  7  7  1
## [530]  1 12  7  7  1 12  1 12 12  3  3 12 12  2  1  7  1 12 12  7 12 12  7
## [553] 12  2 12 12  7  2 12  7  2  7  3 12 12  2  1 12  2  7  2 12  7  7  2
## [576]  1  3  1 12  1  1  2  7  3  3  2 12 12  3  7  7 12  7 12  3  7  1  7
## [599]  2  2  1

Affairs$affairs>0　#TRUE or FALSEで表示

##   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [89] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [100] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [111] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [122] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [144] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [155] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [166] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [177] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [188] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [199] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [210] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [221] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [232] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [243] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [254] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [276] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [287] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [298] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [309] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [320] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [331] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [342] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [353] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [364] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [375] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [386] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [408] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [419] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [430] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [441] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [452]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [463]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [474]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [485]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [496]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [507]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [518]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [529]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [540]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [551]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [562]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [573]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [584]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [595]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

Affairs$affairs<-Affairs$affairs>0 #カテゴリカルに書き換え
head(Affairs) # TRUE or FALSEになっているか確認

##   affairs gender age yearsmarried children religiousness education
## 1   FALSE   male  37        10.00       no             3        18
## 2   FALSE female  27         4.00       no             4        14
## 3   FALSE female  32        15.00      yes             1        12
## 4   FALSE   male  57        15.00      yes             5        18
## 5   FALSE   male  22         0.75       no             2        17
## 6   FALSE female  32         1.50       no             2        17
##   occupation rating
## 1          7      4
## 2          6      4
## 3          1      4
## 4          6      5
## 5          6      3
## 6          5      5

#affirsで層を作ってテーブルを作成
varlist2<-c("gender","age","education","occupation","children")
CreateTableOne(data = Affairs,  vars = varlist2, strata = "affairs")

##                        Stratified by affairs
##                         FALSE         TRUE          p      test
##   n                       451           150                    
##   gender = male (%)       208 (46.1)     78 (52.0)   0.248     
##   age (mean (SD))       32.18 (9.49)  33.41 (8.61)   0.160     
##   education (mean (SD)) 16.14 (2.38)  16.25 (2.48)   0.637     
##   occupation (%)                                     0.348     
##      1                     90 (20.0)     23 (15.3)             
##      2                     10 ( 2.2)      3 ( 2.0)             
##      3                     32 ( 7.1)     15 (10.0)             
##      4                     47 (10.4)     21 (14.0)             
##      5                    160 (35.5)     44 (29.3)             
##      6                    104 (23.1)     39 (26.0)             
##      7                      8 ( 1.8)      5 ( 3.3)             
##   children = yes (%)      307 (68.1)    123 (82.0)   0.002

# strataは""でくくらないとでない