- R言語の基本
- データの読み込み・記述統計量の計算
- グラフによる可視化
- 平均値の検定
- 母比率の検定
- 独立性の検定
- 等分散性の検定
- 一元配置分散分析
- 回帰分析
- ノンパラメトリック検定
- 時系列データの扱い
この記事で解説すること
本シリーズでは、データ分析のためのR言語の使い方について解説しています。
プログラミング初心者を意識した記事となっております。
具体的には次のようなことを学びます。
最後には、学習したことのアウトプットとして演習問題も用意しています！
ぜひ取り組んでみてください！
ファイルの読み込み
read.csv()で読み込む
Rでデータ分析をするには、まずデータを取り込まなければなりません。
データは多くの場合CSVファイルとして提供されます。
CSVファイルとは、「Comma Separated Value」の略で、カンマ（,）で値を区切ったテキストファイル・データのことを指します。
CSVファイルを取り込むには次のコードを書きます。
df1 <- read.csv("./datasets/diabetes.csv")
このコードでは「diabetes.csv」というファイルを、データフレームの形式にして変数df1に代入しています。
diabetes.csvは糖尿病に関するデータセットです。各変数が何を意味するかは下の通りです。
age： 年齢
sex： 性別
bmi： BMI（肥満度）
bp： 平均血圧
s1： TC（血液中の総コレステロール値）
s2： LDL（低比重リポタンパク質、悪玉コレステロール）
s3： HDL（高比重リポタンパク質、善玉コレステロール）
s4： TCH（＝TC÷HDL＝総コレステロール値／善玉コレステロール）
s5： LTG（血液中の中性脂肪値の対数）
s6： GLU（血糖値）
response : ベースラインとなる測定時から1年後の糖尿病の進行度を示す値
＊上記の説明はhttps://atmarkit.itmedia.co.jp/ait/articles/2207/11/news016.htmlから引用しました。
read.csv()のカッコ内には、ファイルがある場所を入力します。これは人によってまちまちなので、あなたの環境にあわせて指定する必要があります（下記参照）。
まずはあなたがどのディレクトリ（フォルダ）で作業しているのかを確認しましょう。
getwd()
"C:/Users/*****/ドキュメント"上のような結果が出力されます。（「*****」にはあなたのユーザー名が入っているはずです。）
＊結果は人によってまちまちなので、必ずしも上記と一致するとは限りません。（上の結果はデフォルト）
この文字列が示しているのは「今あなたが作業しているディレクトリ（フォルダ）」です。今回の場合であれば「今、”ドキュメント”というフォルダの中で作業している」ということになります。
ちなみに、この「作業しているディレクトリ（フォルダ）」のことをワーキングディレクトリ（Working Directory）と言います。
ワーキングディレクトリを変更したい場合（というかデフォルトのままだと後々データの管理が大変になるので変更を推奨）は、次のようにします。
RStudioの場合
上のメニューから
Session > Set Working Directory > Choose Directory
R本体の場合
上のメニューから
ファイル > ディレクトリ変更
次に、データファイルの場所を確認しましょう。
まずは分析したいデータファイルを探し、もしワーキングディレクトリの中に入っていなければ、ワーキングディレクトリへファイルを移動（もしくはコピー&ペースト）しましょう。
これで準備完了です。
例えば「sample.csv」という名前のファイルを読み込むには
read.csv("./sample.csv")
とします。
もし「データファイルはフォルダにまとめて入れておきたい！」という場合は、少し書き加える必要があります。
ワーキングディレクトリの中にdatasetsというフォルダをつくり、そのフォルダ内にsample.csvを置いた場合は
read.csv("./datasets/sample.csv")
とします。
str()でデータの概要を確認する
今、データは「df1」という変数に格納されています。
分析に入る前に、中身を確認しておきたいところですよね。データフレームの概要を確認するには次のようにします。
str(df1)
'data.frame': 442 obs. of 11 variables:
$ age : int 59 48 72 24 50 23 36 66 60 29 ...
$ sex : chr "male" "female" "male" "female" ...
$ bmi : num 32.1 21.6 30.5 25.3 23 22.6 22 26.2 32.1 30 ...
$ bp : int 101 87 93 84 101 89 90 114 83 85 ...
$ s1 : int 157 183 156 198 192 139 160 255 179 180 ...
$ s2 : num 93.2 103.2 93.6 131.4 125.4 ...
$ s3 : int 38 70 41 40 52 61 50 56 42 43 ...
$ s4 : int 4 3 4 5 4 2 3 4 4 4 ...
$ s5 : num 4.86 3.89 4.67 4.89 4.29 ...
$ s6 : int 87 69 85 89 80 68 82 92 94 88 ...
$ response: int 151 75 141 206 135 97 138 63 110 310 ...このように、概要が出力されました。何が書かれているか、簡単に見ていきましょう。
1行目：データ構造とデータのサイズ
'data.frame': 442 obs. of 11 variables:1行目にはデータ構造の種類やデータの大きさに関する情報が書かれています。
「df1」にはデータフレーム型のデータ構造が格納されているので「data.frame」となっています。
「442 obs.」というのは「442 objects」の略だと思いますが、「442行ぶんのデータがあるよ」という意味です。
最後の「11 variables」は列の個数を表します。
2行目～：列名・データ型・先頭数行の値
$ age : int 59 48 72 24 50 23 36 66 60 29 ...2行目以降には各列の概要が書かれています。ここでは2行目の「age」という列を例に説明します。
コロン(：)のすぐ右に「int」とあるのは、その列のデータ型がinteger型、つまり整数型であることを表しています。
Rの主なデータ型は次のようにまとめられます。
| integer(numeric) | 整数型 | 1, -1, 20300, … |
| double(numeric) | 実数型（整数も含む） | 3.1415, -0.001, 125.00012, 1, 1e-10, … |
| complex | 複素数型 | 1i, 1 – 4.5i, 3+0i, complex(re=a,im=b), … |
| character | 文字型 | ”A”, “2014/10/14”, “”, … |
| logical | 論理型 | TRUE, FALSE, T, F, NA, … |
そしてその右側に並んでいる数字たちは、データの一部です。
前処理
データ分析を始める前に前処理を行うことは重要です。というか、分析作業のほとんどは前処理にあると言っても過言ではないかもしれません。
今回のデータの場合はそれほど前処理は必要ないですが、やっておくべきことや、知っておくと便利なことがいくつかあるので紹介します。
適切な型に変換する
$ sex : chr "male" "female" "male" "female" ...データフレームの中のこの列に注目してください。
この列には性別に関するデータが入っており、”male”か”female”のどちらかです。
このように、数値や量では測れない変数のことをカテゴリ変数（質的変数）といいます。
今、「chr」と書かれていることから、この性別データは「文字列型（character）」として扱われていることがわかります。
しかしRでデータ分析をする際には、カテゴリ変数は文字列型ではなく因子型（Factor型）と呼ばれるデータ型で扱った方が便利です。
文字列型から因子型へ変換することで、
「男性のみの場合の平均値を求めたい」
といった操作が簡単にできるようになります。
Rではこの「文字列型→因子型」の変換を次のように行います。
df1$sex <- as.factor(df1$sex)
このように、as.factor()を使えば因子型への変換が行えます。
データフレームから「sex」という列を取り出したい場合、次のようにします。
df1$sex
データフレームの後に「＄」を付けた後に列名を付けます。
すると…
[1] "male" "female" "male" "female" "female" "female" "male" "male" "male" "female" "female" "male"
[13] "female" "male" "female" "male" "female" "male" "female" "female" "female" "male" "female" "male"
......
(中略)
......
[433] "female" "female" "female" "female" "female" "male" "male" "male" "female" "female"ほしい列だけを抜き出してくることができました！
＊ちなみにRを実行するといつも出てくる[1]とか[13]とか[433]は何個目のデータなのかを表しています。例えば[13]とあったら、[13]のすぐ右にある”female”は上から13行目のデータだということがわかります。
str(df1)で、ちゃんと因子型へ変換できているか確認してみましょう。
$ sex : Factor w/ 2 levels "female","male": 2 1 2 1 1 1 2 2 2 1 ...上のようになっていたら成功です！
「w/」って何？？
→ 英語の「with」の省略形だそうです。わざわざ省略しているのは、後に続く情報を少しでも多く表示できるようにしているからなのでしょうね…
「levels」って何？？
→ 日本語に直すと「水準」です。性別なら「female」や「male」のことをそれぞれ水準といいます。都道府県なら「北海道」「青森県」…「沖縄県」のことを水準といいます。
「2 1 2 1 1 2 2 2 1…」って何者？？
→ “female”=1, “male”=2と置き換えたときのデータを表しています。なのでこの数字の羅列は結局「”male” “female” “male” “female” …」と同じことです。どちらが1でどちらが2なのかは、「…levels “female”,”male”: …」で書かれている順番でわかります。
＊このようにカテゴリ変数を数値で置き換えるということは（特に機械学習の分野では）よく行われます。
一部の列を取り出す
「s1」「s2」「s3」「s4」「s5」「s6」の列だけを取り出したデータフレームにしたい…
データフレームから一部の列を取り出すにはいくつか方法があります。
おすすめの方法はこれです。
df1_s <- df1[c("s1","s2","s3","s4","s5","s6")]
df1[]の角カッコの中に列名のベクトルを入れれば複数の列を取り出すことができます。上のコードでは、”s1″から”s6″までの列を取り出して、変数df1_sへ代入しています。
結果はこうなります。
str(df1_s)
'data.frame': 442 obs. of 6 variables:
$ s1: int 157 183 156 198 192 139 160 255 179 180 ...
$ s2: num 93.2 103.2 93.6 131.4 125.4 ...
$ s3: int 38 70 41 40 52 61 50 56 42 43 ...
$ s4: int 4 3 4 5 4 2 3 4 4 4 ...
$ s5: num 4.86 3.89 4.67 4.89 4.29 ...
$ s6: int 87 69 85 89 80 68 82 92 94 88 ...他のやりかたとして、数字で指定する方法もあります。
df1_s <- df1[c(5,6,7,8,9,10)]
ほしい列は、先頭から5,6,7,8,9,10番目の列にあるのでこのようにすることでも列を抜き出すことができます。
上のコードでは
c(5,6,7,8,9,10)
と書いていましたが、このような等差数列はseq()を使うと楽ができます。
seq(5,10, by=1)
これは「初項が5、最後が10となる公差1の等差数列」を生成します。byの値を変えたりすれば3つとばしとかもできます。
列名を変更する
列名を”response”から”y”に変えたい…
データフレームの列名を変えるには次のようにします。
names(df1) <- c("age","sex","bmi","bp","s1","s2","s3","s4","s5","s6","y")
names()でデータフレームの列名を取得し、それに変更後の列名ベクトルを代入しています。
ちょっと面倒ですが、変更しない列名も含め全部の列名を代入しないといけません。つまり
names(df1$response) <- c("y")
のようにはできません。
もし列が100個も200個もあったら、こんなことしていられないよ…
そこで、変更しない列名を書かかないで済む方法を紹介します。
col_names <- names(df1)
col_names[length(col_names)] <- "y"
names(df1) <- col_names
まず、col_namesという変数に列名を代入します（1行目）。
col_namesの中身
[1] "age" "sex" "bmi" "bp" "s1" "s2" "s3" "s4" "s5"
[10] "s6" "response"次に、col_nameの要素である”repsonse”を”y”に書き換えます（2行目）。
2行目は
col_names[11] <- "y"
としてもいいのですが、”response”が何番目にあるのかを数えるのはちょっと大変です。そこで、”response”はベクトルの最後にあることを利用して、「11」とする代わりに「length(col_name)」としています（length()でベクトルの要素の個数を求めてくれます）。
col_namesの中身（書き換え後）
[1] "age" "sex" "bmi" "bp" "s1" "s2" "s3" "s4" "s5" "s6" "y" 最後にこのcol_namesをnames(df1)へ代入してあげれば列名変更が完了します。
記述統計量
カウント
このデータフレームには何個のデータが入っているんだ…？
ということを知りたい場合は、次のようなコードを書きます。
dim(df1) # 行数、列数
nrow(df1) #行数だけ
ncol(df1) #列数だけ
> dim(df1)
[1] 442 11
> nrow(df1) # 行数だけ
[1] 442
> ncol(df1) # 列数だけ
[1] 11dim()を使うと、データフレームの行数と列数がわかります。
今回の場合は
442 11
と出力されます。1番目の値は行数を、2番目の値は列数を表しています。
一方で nrow() や ncol() を使うと、それぞれ行数のみ、列数のみを知ることができます。
＊ちなみに nrow()やncol()を使わなくても次のようにして行数だけ、列数だけを出力することもできます。
dim(df1)[1] # 行数だけ
dim(df1)[2] # 列数だけ
合計
データの合計値を知りたい…
という場合は、sum()を使います。
例えば、列「age」の合計を求めたいなら次のように書きます。
sum(df1$age)
21445では、列「sex」の合計を求めるとどうなるでしょう？
sum(df1$sex)
Error in Summary.factor(c(2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, :
‘sum’ は因子に対しては無意味ですエラーが出てきました。
列「sex」が因子型だからsumができません
と言っています。
このように、数値でない列の合計値を求めようとするとエラーが出てしまうので注意です。
全部の列の合計値を一気に求めたい…
という場合には、apply()を使うと便利です。
たとえば、先ほど「一部の列を取り出す」の節で作成したデータフレーム「df1_s」の全ての列の合計値を一気にまとめたい場合は、次のようにします。
apply(df1_s, 2, sum)
s1 s2 s3 s4 s5 s6
83600.000 51024.100 22006.000 1770.000 2051.504 40337.000カッコ内には3つの引数がありますが、次のような意味です。
- 第1引数 : データ。
- 第2引数：行ごとに集計するか、列ごとに集計するか。行ごと＝1、列ごと＝2
- 第3引数：集計方法。合計を求めるならsumをそのまま入力。
引数とは、ひとことで言うと「カッコ内に入れる情報」のことです。
apply() や sum() など便利な計算をしてくれるものを関数と呼びます。
こういった関数を使うには、このカッコ内には情報を入れなければなりません（入れなくても良い関数もある）。
sum()という関数だったら「何を合計したいのか」という情報が必要になりますよね。そのような情報のことを引数というのです。
なお「どんな引数を指定したらいいのか」というのは関数によってさまざまです。これについては、関数が登場するたびに説明していきますのでムリして覚える必要はありません。
また、引数は1つだけではなく2つ以上存在する場合もあります。この場合には、1つ目の引数を第1引数、2つ目の引数を第2引数と呼んだりすることがあるので、知っておくと良いでしょう。
applyは同じ処理をまとめて行える便利な関数です。この後にも何回か出てきますので、しっかりと押さえておきましょう！
平均値
データの平均値を知りたい…
という場合には、mean()を使います。
全体の平均年齢を調べたい場合は、次のようにします。
mean(df$age)
48.5181また、データフレームの全ての列に対して平均値をそれぞれ求めたい場合はapply()を使いましょう。
apply(df1_s, 2, mean)
s1 s2 s3 s4 s5 s6
189.140271 115.439140 49.787330 4.004525 4.641411 91.260181男性の平均年齢を調べたい…
という場合は、次のようにします。
mean(df1[df1$sex == "male", "age"])
50.94203mean()を使うという点では同じですが、引数がちょっと複雑ですね。何をしているか、詳しく見ていきましょう。
まず、
df1[df1$sex == “male”, “age”]
では何をしているかというと、
df1の列「sex」が”male”になっている行で、そこから列「age」のみを取り出す
ということをしています。
少し複雑ですが、こういった操作はデータ分析で頻繁に行われますので必ずマスターしておきましょう！
実は、ちょっとだけ違う書き方もできます。というか、僕はこちらの書き方をよくします。
mean(df1[df1$sex == "male",]$age)
結局は同じことをしているのですが、違いは「列の抽出を先にやるか後にやるか」ということだけです。
結果は全く同じなのでどちらでも良いのですが
「男性のデータだけを集めたデータフレームを作りたい」
と言った場合には、
df1[df1$sex == "male",]
とすればよいことがわかります。
なので、両方の書き方を知っておくと良いかもしれませんね。
分散・標準偏差
分散や標準偏差を求めたい…
という場合には、var()やsd()を使います。
var(df1$age)
sd(df1$age)
171.8466
13.10903また、データフレームの全ての列に対して分散・標準偏差をそれぞれ求めたい場合はapply()を使いましょう。
apply(df1_s, 2, sd)
s1 s2 s3 s4 s5 s6
34.6080517 30.4130810 12.9348629 1.2830581 0.5223906 11.4963347ここで注意しなければならないことが1点あります。
それはvar()やsd()で計算される値は標本分散・標本標準偏差ではなく、不偏分散・不偏標準偏差だということです。
標本分散：\(\displaystyle s^2=\frac{1}{N}\sum_{i=1}^{N}(x_i-\bar{x})^2 \)
不偏分散：\(\displaystyle \hat{s}^2=\frac{1}{N-1}\sum_{i=1}^{N}(x_i-\bar{x})^2 \) ←sd()で出てくる結果はこっち！
じゃあ標本分散を計算したいときはどうする？
標本分散と不偏分散の間には次のような関係が成り立つことに注意しましょう。
$$\displaystyle s^2 = \frac{N-1}{N}\hat{s}^2$$
よって、次のようなコードで標本分散を計算することができます。
n <- nrow(df1)
v <- var(df1$age)*(n-1)/n
v
171.4578また、標本分散の平方根をとることで標本標準偏差も求めることができます。
sqrt(v)
13.09419中央値
中央値を求めたい…
という場合には、median()を使います。
median(df1$age)
50また、データフレームの全ての列に対して中央値を求めたい場合はapply()を使いましょう。
apply(df1_s, 2, median)
s1 s2 s3 s4 s5 s6
186.00000 113.00000 48.00000 4.00000 4.62005 91.00000最大値・最小値
最大値・最小値を求めたい…
という場合には、max()やmin()を使います。
max(df1$age)
min(df1$age)
79
19また、データフレームの全ての列に対して最大値・最小値を求めたい場合はapply()を使いましょう。
apply(df1_s, 2, max)
apply(df1_s, 2, min)
s1 s2 s3 s4 s5 s6
301.000 242.400 99.000 9.000 6.107 124.000 s1 s2 s3 s4 s5 s6
97.0000 41.6000 22.0000 2.0000 3.2581 58.0000四分位数
四分位数を求めたい…
という場合には、quantile()を使います。
これを使うと、
0％点(=最小値)、25％点、50％点(＝中央値)、75％点、100％点(＝最大値)
を求めることができます。
quantile(df1$age)
0% 25% 50% 75% 100%
19.00 38.25 50.00 59.00 79.00例えば「5％点と95％点」とかって求めることはできるの？
はい、できます。第2引数でパーセンタイル（0以上1以下）を指定すれば求めることができます。
quantile(df1$age, c(0.05,0.95))
5% 95%
25 68最頻値
最頻値を求めたい…
残念ながら最頻値を求める関数はもとから入っていないため、少し面倒です。
最頻値を求めるには、次のようにします。
names(which.max(table(df1$sex)))
"female"ここで使われているtable()は、各値の頻度を求める関数です。
table(df1$sex)
female male
235 207which.max()は、最大値がどこにあるかを探す関数です。つまり、最大値235を探し、その場所を示す「female」を返します。
whch.max(table(df1$sex))
"female"最頻値求めるたびにこれを書くのは面倒だよ…。
と思う方が多いと思います。そんなときは「自作関数」を作ってしまえばよいのです。
つまり、今までのmean()やsum()などのようにmode()を自分で定義します。
mode <- function(x) names(which.max(table(x)))
これを実行したあとは、
mode(df1$sex)
と書くだけで最頻値を求めることができてしまいます！
相関係数
2変数の間の相関係数を求めたい…
という場合には、cor()を使います。
bmi と y の間の相関係数を求めたければ、次のようにします。
cor(df1$bmi, df1$y)
0.5864501また、データフレームの全ての列の組み合わせで相関係数を求めたいなら、cor()の中にデータフレームをそのままぶち込みます。
cor(df1)
Error in cor(df1) : 'x' は数値でなければなりませんError in cor(df1) : 'x' must be numericエラーが出てきました。
これはdf1の中に数字でないデータが混じっていることが原因です。
もちろん、数値でない列を除去してぶち込めばよいのですが、もっと便利な方法があります。
Filter()という関数を利用して、数値型の列のみを取り出したデータフレームを作ることができます。
Filter(is.numeric, df1)
これをcor()の中にぶち込めば…
cor(Filter(is.numeric, df1))
age bmi bp s1 s2 s3
age 1.00000000 0.1850847 0.3336003 0.26006082 0.2192431 -0.07520042
bmi 0.18508467 1.0000000 0.3948794 0.24977742 0.2611699 -0.36672527
bp 0.33360025 0.3948794 1.0000000 0.24180485 0.1848695 -0.17747785
s1 0.26006082 0.2497774 0.2418049 1.00000000 0.8966630 0.05154242
s2 0.21924314 0.2611699 0.1848695 0.89666296 1.0000000 -0.19646890
s3 -0.07520042 -0.3667253 -0.1774779 0.05154242 -0.1964689 1.00000000こんなふうに行列のような形で各列同士の相関係数を計算してくれるのです！
＊統計学ではこの行列には相関行列という名前がついています。
まとめて見る
統計量の計算方法は分かったけど、まとめて計算する方法ってないの…？
実はあるんですよ…。しかもめちゃくちゃ簡単に。
summary(df1)
age sex bmi bp
Min. :19.00 female:235 Min. :18.00 Min. : 62.0
1st Qu.:38.25 male :207 1st Qu.:23.20 1st Qu.: 84.0
Median :50.00 Median :25.70 Median : 93.0
Mean :48.52 Mean :26.38 Mean : 94.6
3rd Qu.:59.00 3rd Qu.:29.27 3rd Qu.:105.0
Max. :79.00 Max. :42.20 Max. :133.0 はい、これだけです。データをながめるだけならこれで十分です。
しかし、場合によっては
「ageの列の平均値だけほしい」
といったことも多々あります。
なので、summary()だけではなく上で紹介した方法もしっかり習得しておきましょう。
演習問題
今回学んだことの復習として、演習問題を用意しました。
wine.csvをダウンロードし、次の問いに答えなさい。
alcohol： アルコール度数
malic_acid： リンゴ酸
ash： 灰分（かいぶん）
alcalinity_of_ash： 灰分のアルカリ度
magnesium： マグネシウム
total_phenols： 全フェノール含量
flavanoids： フラボノイド
nonflavanoid_phenols： 非フラボノイドフェノール
proanthocyanins： プロアントシアニン
color_intensity： 色の濃さ
hue： 色相
od280/od315_of_diluted_wines： 希釈ワイン溶液のOD280／OD315（＝280nmと315nmの吸光度の比）
proline： プロリン
＊上記の説明はhttps://atmarkit.itmedia.co.jp/ait/articles/2208/25/news046.htmlから引用しました。
問1
wine.csvをデータフレーム型として読み込み、データの概要を確認しなさい。
問2
列「od280/od315_of_diluted_wines」の名前を「od280/od315」に変更しなさい。
問3
質的変数をFactor型に変換しなさい。
問4
列「class」について、A,B,Cの頻度をそれぞれ調べなさい。
問5
すべての数値型の列について、平均値をすべて調べなさい。
問6
すべての数値型の列について、標本分散をすべて調べなさい。
問7
全ての列について、基本統計量をまとめて表示させなさい。
問8
相関行列を調べなさい。また「total_phenols」と「flavanoids」の相関係数を調べなさい。
【解答】
# (1) ファイルを読み込む
df2 <- read.csv("./datasets/wine.csv")
# 中身を確認する
str(df2)
# 行数と列数を調べる
dim(df2)
# (2) 列名を変更する
col_names <- names(df2)
col_names[length(col_names)-2] <- "od280/od315"
names(df2) <- col_names
# (3) 質的変数をFactor型に
df2$class <- as.factor(df2$class)
# (4) 各classの個数
table(df2$class)
# (5) 平均値
apply(Filter(is.numeric,df2), 2, mean, na.rm=T)
# 【補足】na.rm=T または na.rm=TRUE とすると欠損値を無視して計算してくれる。そうでないと、計算結果がNA(Not Available)になってしまう。
# (6) 標本分散
n <- nrow(df2)
apply(Filter(is.numeric,df2), 2, var, na.rm=T)*(n-1)/n
# (7) まとめて
summary(df2)
# (8) 相関係数
cor(Filter(is.numeric, df2))
cor(df2$total_phenols, df2$flavanoids)
今回はこれで以上となります。
次回はグラフによるデータの可視化を行います！
記事への意見・感想はコチラ