KU-R Session 1

```

本日の内容

• R言語について、R/自由(フリー)ソフトウェアの利点と欠点
• EZRについて(紹介のみ)
• RとRStudioのインストール
• タブ補完とヘルプの見方、付属データの利用
• CRANの設定、追加パッケージの探し方、インストール方法、使用方法
• カレントディレクトリの概念とRのファイルの種類
• RMarkdownによる記録の残る解析方法(reproducible research)
• 外部データファイルの読み込み
• 中間データの保存とデータセットの出力
• 解析に向けてのデータの整理(データフレーム、ベクター、カテゴリ変数、二値化、サブセット)

R言語について、R/自由(フリー)ソフトウェアの利点と欠点

• GNU R: http://www.r-project.org

GNU Rは統計解析に特化したプログラミング言語、および、その実行環境です。開発の開始は1993年ごろですが、ここ数年で急速に注目を集めた印象があります。GNU General Public Licenseというフリーソフトウェアライセンスになっており、無料、かつ、配布および改変の自由が保証されています。R本体の開発はCore Developer Teamと呼ばれている統計学者のグループが行っています。

利点

• 利用が自由。
• 利用者が増加している(Harvard Chan Schoolでは生物統計PhDのデータ分析の授業はRでやっている。ただ疫学系はSAS。)
• プログラミング言語でもあることから拡張性が高い。
• 豊富(6000+)な外部パッケージが存在し、新規統計手法の実装もRで行われていることが多い。
• 情報がインターネットに大量に転がっている http://www.r-bloggers.com など参考。

欠点

• 公式のユーザサポートはない(一応メーリングリストはあるがあまり機能していない)。
• コマンドラインの利用は敷居が高い点もある。
• 処理系の老朽化が目立ち、データ量が増えると速度面で不満がある。

EZRについて

R本体はコマンドでの操作のみしか受け付けません。しかし、Rを初心者にも身近にする試みとしていくつかGUIフロントエンドが存在します。

• R Commander: http://www.rcommander.com
• EZR (医療系向けのR Commander改): http://www.jichi.ac.jp/saitama-sct/SaitamaHP.files/statmed.html
• Deducer (グラフ機能が秀逸だが現在は開発が滞っている印象): http://www.deducer.org/pmwiki/pmwiki.php?n=Main.DeducerManual?from=Main.HomePage
• Rz (現在は開発が滞っている印象): http://m884.jp/RzIntro_Ja.html

EZRは自治医科大学さいたま医療センター血液内科の神田先生が開発をしているため医療系のニーズにあった内容になっています。本もあります。

• http://www.amazon.co.jp/初心者でもすぐにできるフリー統計ソフトEZR-Easy-で誰でも簡単統計解析-神田-善伸/dp/4524261583/ref=sr_1_1?ie=UTF8&qid=1425597448&sr=8-1&keywords=ezr

私が普段使っていないので説明が困難という問題点があるため、以降はコマンドラインでRを使用する方法について説明します。

RとRStudioのインストールと設定

R利用する際のコーディングの利便性を向上するソフトウェアとしてRStudioが存在します。現在はRStudioを介してRを利用するのが一般的なので、(私のようにemacsに強い執着があるなどの場合を除き)RStudioを利用してください。RStudioの概要について、画面を見ながら説明します。

• RStudio: http://www.rstudio.com/products/RStudio/

タブ補完とヘルプの見方、付属データの利用

コマンドの入力というと煩雑な気がしますが、実際はタブキーによる補完を多用するのですべてをタイプする必要はありません。いくつか行ってみましょう。ヘルプ機能が比較的充実していますので試してみます。

## ls関数のヘルプを見る
?ls

## ROCという単語に関連したヘルプを一覧
??ROC

いろいろな付属データがあり、解析の練習をするのに最適です。インストールされているパッケージ内のデータセットをすべて見るには下記のコマンドでできます。

lapply(rownames(installed.packages()), function(pkg) {data(package = pkg)})

## フィッシャーのアヤメ(iris)データセットを読み込む
data(iris)
head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

CRANの設定、追加パッケージの探し方、インストール方法、使用方法

Comprehensive R Archive Network (CRAN)はRの追加パッケージの中央レポジトリです。RStudioの設定で0-Cloudに設定すれば良いです。

Task Viewに分野ごとにまとめられているので特定の目的に沿ったパッケージを探すにはこちらが便利です。あとで使用しますので、XLConnectパッケージをインストールしましょう。Javaがインストールされていないと、そちらのインストールが利用前に必要かもしれません。

• CRAN Task View: http://cran.r-project.org/web/views/

さきほど追加パッケージのインストール方法を説明しましたが、インストールしただけでは利用できるようにならず、読み込みが必要です。

library(XLConnect)

## Loading required package: XLConnectJars
## XLConnect 0.2-11 by Mirai Solutions GmbH [aut],
##   Martin Studer [cre],
##   The Apache Software Foundation [ctb, cph] (Apache POI, Apache Commons
##     Codec),
##   Stephen Colebourne [ctb, cph] (Joda-Time Java library)
## http://www.mirai-solutions.com ,
## http://miraisolutions.wordpress.com

下記のコマンドでインストールされているすべてのパッケージ名が見られます。

rownames(installed.packages())

カレントディレクトリの概念とRのファイルの種類

Linuxなどのシェル(コマンドプロンプト)を利用したことがあるかたは御存じかもしれませんが、Rがいまどこのフォルダを見ているか(カレントディレクトリ)という概念があります。外部のファイルを取り込みするのに理解が必要なので確認します。RStudioではシェルの区画の上の部分に表示されており、メニューで変更できますが、コマンドラインでは下記のように行います。

## カレントディレクトリを表示
getwd()

## [1] "/Users/kazuki/Documents/KyotoUniv/2015KyotoU_R/code"

## カレントディレクトリの中身を表示
dir()

##  [1] "auto"                       "figure"                    
##  [3] "KU-R1.html"                 "KU-R1.R"                   
##  [5] "KU-R1.Rmd"                  "KU-R2.html"                
##  [7] "KU-R2.Rmd"                  "nhefs_book.csv"            
##  [9] "nhefs_book.xls"             "nhefs_codebook.xls"        
## [11] "nhefs_export.csv"           "Presentation1.md"          
## [13] "Presentation1.Rpres"        "SweaveDemo-concordance.tex"
## [15] "SweaveDemo.log"             "SweaveDemo.pdf"            
## [17] "SweaveDemo.Rnw"             "SweaveDemo.synctex.gz"     
## [19] "SweaveDemo.tex"             "temp1.RData"               
## [21] "temp2.RData"                "Untitled.html"             
## [23] "Untitled.R"

## カレントディレクトリを変更(下記は現在のカレントディレクトリに変更しているので実質変更なし)
setwd("./")

Rの作成するファイルには様々な種類があります。

• .R Rスクリプトファイル。通常のコードです。
• .RData 中間データ保存ファイル
• .Rmd RMarkdownファイル。このファイルのようなものです。
• .Rnw Sweaveファイル。RとLaTeXのコードを一緒に書くためのファイル。

New Project…の方ではShinyというフレームワークを使ってWeb appを作成することもできます。

• Ebolaの経時的分布のweb app
• web app https://kaz-yos.shinyapps.io/hackebola/
• コード https://github.com/kaz-yos/hack-ebola/tree/master/hackebola

RMarkdownによる記録の残る解析方法(reproducible research)

この資料自体が、RMarkdownとよばれる、コードと文書を一緒に記載して解析コードと結果結果をいっしょに残すための仕組みを使って作成されています。それほど難しくなく、ちょっとした解析結果をまとめるには大変便利ですので利用してみましょう。ファイルの作成方法を説明します。下記はひとつ前のセクションがどのように書かれているかの解説です。

 ## カレントディレクトリの概念とRのファイルの種類

Linuxなどのシェル(コマンドプロンプト)を利用したことがあるかたは御存じかもしれませんが、Rがいまどこのフォルダを見ているか(カレントディレクトリ)という概念があります。外部のファイルを取り込みするのに理解が必要なので確認します。RStudioではシェルの区画の上の部分に表示されており、メニューで変更できますが、コマンドラインでは下記のように行います。


# ```{r}

## カレントディレクトリを表示
getwd()
## カレントディレクトリの中身を表示
dir()
## カレントディレクトリを変更(下記は現在のカレントディレクトリに変更しているので実質変更なし)
setwd("./")

#  ```

Rの作成するファイルには様々な種類があります。

- .R Rスクリプトファイル。通常のコードです。
- .RData 中間データ保存ファイル
- .Rmd RMarkdownファイル。このファイルのようなものです。
- .Rnw Sweaveファイル。RとLaTeXのコードを一緒に書くためのファイル。


New Project...の方ではShinyというフレームワークを使ってWeb appを作成することもできます。

- Ebolaの経時的分布のweb app
- web app https://kaz-yos.shinyapps.io/hackebola/
- コード https://github.com/kaz-yos/hack-ebola/tree/master/hackebola

外部データファイルの読み込み

Hernan Bookのデータセットを利用しますので、ダウンロードしてこのファイルと同じフォルダに入れてください。

## R内からダウンロードすることも可能です
## download.file(url = "http://cdn1.sph.harvard.edu/wp-content/uploads/sites/1268/2013/01/nhefs_book.xls",
##               destfile = "./nhefs_book.xls")

データセットのcode book xls file に記載があります。

Rは各種の外部データファイルを読みこむことができますが、上手くいかない場合はテキスト形式(.csv)が無難です。読み込んだデタを名前をつけてオブジェクトに格納します。CSVファイル、Stata、SPSS、SASのファイルの読み込みはそれぞれ下記のような形になります。

## *.csv 
dat <- read.csv("./file.csv")

## StataとSPSSはforeignパッケージを使用します。
library(foreign)
## Stata *.dta
dat = read.dta("./file.dta")
## SPSS *.sav (パッケージ名をダブルコロンの前に置けばlibrary(foreign)なしでも可能)
dat <- foreign::read.spss("./file.sav")

## SAS *.sas7bdat
dat <- sas7bdat::read.sas7bdat("./file.sas7bdat")

実際多いと思われるExcelファイルの取り込みはXLConnectパッケージが一番安定している印象です。

## library(XLConnect)か XLConnect:: のどちらかがあれば良いです。
library(XLConnect)
nhefs <- XLConnect::readWorksheetFromFile(file = "./nhefs_book.xls", sheet = 1)

## うまくいかない場合は切り替えて.csvにエクセルで変換して読み込んだ方が良いです。
## nhefs2 <- read.csv(file = "./nhefs_book.csv")

中間データの保存とデータセットの出力

基本的に元データとコードを信頼して、中間データを保存しない、というのが一番安全な管理方法かもしれません。ただ、データ処理に時間がかかる場合はそうも言っていられないので、データを保存するのが良いかもしれません。

中間データの保存はRのデータ形式に保存します。

## 環境に存在するオブジェクトを一覧
ls()

## 全部の中間データを保存 ()
save.image("./temp1.RData")

## 指定の中間データのみを保存
save(list = c("iris"), file = "./temp2.RData")


## 特定のデータセットをエクセルで読める形で保存するには.csvにします
## Rは欠測値にNAを使用しますので出力するときは、na = ""で空欄にするのが良いでしょう。
write.csv(nhefs, file = "./nhefs_export.csv", na = "")

こういった中間データを読み込むときは下記のようにします。

## *.RDataの中間ファイルのデータを読み込みする
load("./temp1.RData")

## *.RDataの中間ファイルのデータを読み込みする
load("./temp2.RData")

解析に向けてのデータの整理(データフレーム、ベクター、カテゴリ変数、二値化、サブセット)

データ形式

他の統計解析環境と比べた場合のユニークな点として、データ形式が非常に豊富ということがあげられます。通常のいゆるデータセット(SASやStataのそれと似たもの)はdata frameと呼ばれる形式です。これはデータの最小単位ではなく、data frameの列を形成しているvectorが最小単位になります。

## nhefsデータセットのデータ形式を確認
class(nhefs)

## [1] "data.frame"

## 一部を閲覧
head(nhefs)

##   seqn qsmk death yrdth sbp dbp sex age race income marital school
## 1  233    0     0    NA 175  96   0  42    1     19       2      7
## 2  235    0     0    NA 123  80   0  36    0     18       2      9
## 3  244    0     0    NA 115  75   1  56    1     15       3     11
## 4  245    0     1    85 148  78   0  68    1     15       3      5
## 5  252    0     0    NA 118  77   0  40    0     18       2     11
## 6  257    0     0    NA 141  83   1  43    1     11       4      9
##         ht  wt71      wt82    wt82_71 birthplace smokeintensity
## 1 174.1875 79.04  68.94604 -10.093960         47             30
## 2 159.3750 58.63  61.23497   2.604970         42             20
## 3 168.5000 56.81  66.22449   9.414486         51             20
## 4 170.1875 59.42  64.41012   4.990117         37              3
## 5 181.8750 87.09  92.07925   4.989251         42             20
## 6 162.1875 99.00 103.41906   4.419060         34             10
##   smkintensity82_71 smokeyrs asthma bronch tb hf hbp pepticulcer colitis
## 1               -10       29      0      0  0  0   1           1       0
## 2               -10       24      0      0  0  0   0           0       0
## 3               -14       26      0      0  0  0   0           0       0
## 4                 4       53      0      0  0  0   1           0       0
## 5                 0       19      0      0  0  0   0           0       0
## 6                10       21      0      0  0  0   0           0       0
##   hepatitis chroniccough hayfever diabetes polio tumor nervousbreak
## 1         0            0        0        1     0     0            0
## 2         0            0        0        0     0     0            0
## 3         0            0        1        0     0     1            0
## 4         0            0        0        0     0     0            0
## 5         0            0        0        0     0     0            0
## 6         0            0        0        0     0     0            0
##   alcoholpy alcoholfreq alcoholtype alcoholhowmuch pica headache otherpain
## 1         1           1           3              7    0        1         0
## 2         1           0           1              4    0        1         0
## 3         1           3           4             NA    0        1         1
## 4         1           2           3              4    0        0         1
## 5         1           2           1              2    0        1         0
## 6         1           3           2              1    0        1         0
##   weakheart allergies nerves lackpep hbpmed boweltrouble wtloss infection
## 1         0         0      0       0      1            0      0         0
## 2         0         0      0       0      0            0      0         1
## 3         0         0      1       0      0            0      0         0
## 4         1         0      0       0      0            0      0         0
## 5         0         0      0       0      0            1      0         0
## 6         0         0      0       0      0            0      0         0
##   active exercise birthcontrol pregnancies cholesterol hightax82  price71
## 1      0        2            2          NA         197         0 2.183594
## 2      0        0            2          NA         301         0 2.346680
## 3      0        2            0           2         157         0 1.569580
## 4      1        2            2          NA         174         0 1.506592
## 5      1        1            2          NA         216         0 2.346680
## 6      1        1            0           1         212         1 2.209961
##    price82     tax71     tax82 price71_82  tax71_82
## 1 1.739990 1.1022949 0.4619751 0.44378662 0.6403809
## 2 1.797363 1.3649902 0.5718994 0.54931641 0.7929688
## 3 1.513428 0.5512695 0.2309875 0.05619812 0.3202515
## 4 1.451904 0.5249023 0.2199707 0.05479431 0.3049927
## 5 1.797363 1.3649902 0.5718994 0.54931641 0.7929688
## 6 2.025879 1.1547852 0.7479248 0.18408203 0.4069824

## qsmk変数のvectorをとりだす
## $はこのdata frameの中でqsmkという名前の付いた部分を取り出すという意味です
head(nhefs$qsmk, 100)

##   [1] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 1 0 1 0 0 0 0 1 0
##  [36] 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0
##  [71] 0 1 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1

下記は、ほぼ同じ意味です。

## このdata frameの中でqsmkという名前の付いた部分を取り出す
nhefs["qsmk"]

## qsmkという列を取り出す
nhefs[, "qsmk"]

Rに特徴的なものとしてはlistと呼ばれるデータ形式があります。これは任意のデータを束ねておける非常に強力なデータ型式です。Rの解析結果は、ほぼ必ずこのlistと呼ばれるデータ型式で返ってきます。

## 適当なデータをlistに格納
lst1 <- list(data1 = head(nhefs$qsmk), data2 = tail(nhefs[,1:7]), number = 3, name = "test list")
## listを表示
lst1

## $data1
## [1] 0 0 0 0 0 0
## 
## $data2
##       seqn qsmk death yrdth sbp dbp sex
## 1741 25014    0     0    NA 115  66   0
## 1742 25016    0     0    NA 124  80   1
## 1743 25024    0     0    NA  NA  NA   1
## 1744 25032    0     0    NA 171  77   0
## 1745 25033    0     0    NA 115  90   0
## 1746 25061    1     0    NA 136  90   0
## 
## $number
## [1] 3
## 
## $name
## [1] "test list"

## listの一部だけ取り出す
lst1[["name"]]

## [1] "test list"

実はdata frameもlistの一種です。同じ長さのvectorが束ねられて表になっているわけです。

## data frameの1:10行目と1:5列目を抜き出し
nhefs[1:10, 1:5]

##    seqn qsmk death yrdth sbp
## 1   233    0     0    NA 175
## 2   235    0     0    NA 123
## 3   244    0     0    NA 115
## 4   245    0     1    85 148
## 5   252    0     0    NA 118
## 6   257    0     0    NA 141
## 7   262    0     0    NA 132
## 8   266    0     0    NA 100
## 9   419    0     1    84 163
## 10  420    0     1    86 184

## デフォルトの表示方法をすると、listとして表示される。
print.default(nhefs[1:10, 1:5])

## $seqn
##  [1] 233 235 244 245 252 257 262 266 419 420
## 
## $qsmk
##  [1] 0 0 0 0 0 0 0 0 0 0
## 
## $death
##  [1] 0 0 0 1 0 0 0 0 1 1
## 
## $yrdth
##  [1] NA NA NA 85 NA NA NA NA 84 86
## 
## $sbp
##  [1] 175 123 115 148 118 141 132 100 163 184
## 
## attr(,"class")
## [1] "data.frame"

データの変換

解析の準備に必要な典型的なデータの変換の例を示します。

## 年齢を二値化して新しい変数としてデータセットに作成
summary(nhefs$age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   25.00   33.00   43.00   43.64   53.00   74.00

## 65歳以上かどうかの0,1変数を作成する。
nhefs$elderly <- as.numeric(nhefs$age >= 65)
table(nhefs$elderly)

## 
##    0    1 
## 1653   93

自動的にカテゴリ変数として扱って欲しい場合はfactorというデータ形式を与えます。名前を付けることもできます。

## カテゴリ変数を数値データからfactorに変換
nhefs$diab_cat <- factor(nhefs$diabetes)
## 元の変数は数値としてサマリーされる
summary(nhefs$diabetes)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.9868  2.0000  2.0000

## 新しい変数はカテゴリとしてサマリーされる。
summary(nhefs$diab_cat)

##   0   1   2 
## 877  15 854

## カテゴリ変数のレベルに名前を付ける
nhefs$diab_cat <- factor(nhefs$diabetes, levels = c(0,1,2), labels = c("Never","Ever","Missing"))
summary(nhefs$diab_cat)

##   Never    Ever Missing 
##     877      15     854

欠測値はNAとしてあらわされます。前述の変数はNAとなっていなかったので変換します。

## nhefs$diab_cat変数で"Missing"と一致する場所をとりだす。
head(nhefs$diab_cat[nhefs$diab_cat == "Missing"])

## [1] Missing Missing Missing Missing Missing Missing
## Levels: Never Ever Missing

## それらをNAで置き換える
nhefs$diab_cat[nhefs$diab_cat == "Missing"] <- NA
## 結果確認
summary(nhefs$diab_cat)

##   Never    Ever Missing    NA's 
##     877      15       0     854

## NAのチェックは等号では行えません。
sum(nhefs$diab_cat == NA)

## [1] NA

## is.na()を用います。
sum(is.na(nhefs$diab_cat))

## [1] 854

データセットの一部のみでデータセットを作成する場合はsubsetを利用します。

## 高齢者の行だけでデータセットを作成
nhefs_elderly <- subset(nhefs, elderly == 1)
## 行数を確認
nrow(nhefs_elderly)

## [1] 93

## 高齢者変数の長さ
length(nhefs_elderly$elderly)

## [1] 93

## 変数の内容を確認
nhefs_elderly$elderly

##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

年齢を10歳ごとに区切ってカテゴリ変数を作成します。

nhefs$ageCat <- cut(nhefs$age,
                    breaks = c(-Inf,30,40,50,60,70,Inf), # カットオフの指定。一番外は無限を使うとよい
                    labels = c("20s","30s","40s","50s","60s","70s"), # カテゴリ名
                    include.lowest = TRUE, # 一番最初の値も含める
                    right = FALSE) # [30, 40) 範囲の左側を閉じる。[以上,未満)にする。
summary(nhefs$ageCat)

## 20s 30s 40s 50s 60s 70s 
## 268 441 452 393 168  24

年齢ときちんと対応しているかチェックします。ソートにはdplyrパッケージのarrange()関数が扱いやすいです。

library(dplyr)
dplyr::arrange(nhefs[c("age","ageCat")], age)