統計解析ソフトウエアR 入門

雑誌名 久留米大学コンピュータジャーナル

巻 28

ページ 2‑14

発行年 2014‑03‑01

URL http://hdl.handle.net/11316/530

- 2 -

1. はじめに はじめに はじめに はじめに

コンピュータと統計解析ソフトウエアの、最近の進歩によって各種の統計処理が格段 に楽になった。しかし、一方で、

SAS

や

SPSS

など信頼性の高い市販のソフトウエアは 個人で購入するには高価である。そういう中、フリーソフトウエアで信頼性の高いソフ トウエアがある。「R」はその１つである。

本稿では、データを処理する際の、統計解析ソフトウエア

R

の基本的な扱い方につ いて説明する。データが入力されている

EXCEL

ファイルを読み込み、基本的なデータ 処理するという状況を想定し、具体的な方法を説明する。

R

はすでにインストール済み だとして話を進める。このソフトウエアは、パソコンの複数の

OS、 Windows、 MacOS、

Linux

などさまざまな

OS

に対応したものがつくられている。ここでは、MacOS版をつ

かうことにしよう。

記載している具体的な項目はつぎのものである。

（１） 統計解析ソフトウエア

R

の起動と終了

（２）

EXCEL

ファイルからのデータの読み込み

（３） 度数分布表の作成

（４） ヒストグラムを描く

（５） ２変数間の関連を散布図でみる

（６） 相関係数を求める

（７） 回帰直線を描く

（８） ２群の平均値の比較の検定

（９） コンピュータシミュレーション

２．統計解析ソフトウエア ２．統計解析ソフトウエア ２．統計解析ソフトウエア

２．統計解析ソフトウエア R の操作 の操作 の操作 の操作

†久留米大学 医学部 看護学科 教授

†The School of Nursing, KURUME University.

Professor.

犬塚裕樹^†

Hiroki Inutsuka

^†

統計解析ソフトウエア 統計解析ソフトウエア 統計解析ソフトウエア

統計解析ソフトウエア R 入門 入門 入門 入門

- 3 - 2.1統計解析ソフトウエア統計解析ソフトウエア統計解析ソフトウエア統計解析ソフトウエアRの起動の起動の起動の起動

R

のショートカットキーをクリックする。すると、図のような画面があらわれる。

起動後画面には「

>

」のプロンプトが表示される。これは、ユーザの入力待ちの状態を あらわしている。

終了する場合には

> q()

と入力する。あるいは、メニューの「

R

」をクリックし、サブメニューの「終了」項目 をクリックする。

R

を起動し、「

>

」のプロンプトがでている状態から、四則演算の計算ができるところ がおもしろい。ためしに、足し算をおこなってみよう。

>15+5

とキーをうち、enterキーをおして入力すると、

図１

R

の初期画面

- 4 -

[1] 20

と結果が表示される。かけ算、累乗はそれぞれ「*」と「^」を利用する。「５×６」と

「

2

の３乗」を計算してみよう。

> 5*6 [1] 30

> 2^3 [1] 8

ここで、下線はユーザの入力部分をあらわすことにする。このように簡単に扱えること がわかる。

2.2 EXCELファイルからファイルからファイルからファイルからのデータののデータののデータの読み込みのデータの読み込み読み込み読み込み

最近では、データは、たいてい、テキストファイルや

EXCEL

ファイルのシートに入 力されているだろう。

R

は、テキストファイ

ルや

EXCEL

の

csv

形式のファイルから、書

き込まれているデータを取り込むことがで きる。ファイルを指定するためには、ファイ ルを作業ディレクトリにおいておく。標準で

は「

/User/home/

ユーザ名」が作業ディレク

トリーに設定されている。しかし、この場所 はユーザが簡単に変更することができる。こ

こでは

/User/hinutsu

として扱うことにしよ

う。

EXCEL

ファイル、「練習.csv」という名称

のファイルの内容が図１にしめされている。

項目は「id」「a1」「a2」の３つある。これを

「練習データ」という名前の「データフレー ム」に読み込む。この「データフレーム」と は

R

におけるデータの保存形式のことであ る。これは、数値や文字など異なるデータタ イプのものをまとめて扱うことができる。

「練習.csv」ファイルからデータを読み込むた めに、つぎのように入力する。

図２ 練習.CSV ファイルの内

- 5 -

>

練習データ

<- read.csv(“

練習

.csv”)

これで図１の内容のデータをもつ「練習データ」という名称のデータフレームが作成 された。

このデータフレームの内容を表示しよう。そのためには

>

練習データ

と入力するだけでよい。すると、図２と同じものが表示される。

このデータのうち、項目「a1」のデータを表示しよう。

>

練習データ[,2]

と入力すると

2

列目の「

a1

」の項目のデータが表示される。

[1] 2.0 7.0 3.0 2.5 4.0 4.0 5.0 4.0

、、、

あるいは、項目名「a1」をつかって

>

練習データ[“a1”]

と入力しても同じ結果となる。項目「

a1

」の内容が表示された。

1 2.0 2 7.0 3 3.0 4 2.5 5 4.0

ただし、前者が横に表示されるのに対して、後者は縦に表示される。

2.3 度数分布表の作成度数分布表の作成度数分布表の作成度数分布表の作成

- 6 -

データが読み込まれたので、データの属性をみるために、まず、各変数の度数分布をみ よう。

項目「a1」の度数分布表を作成する。つぎのように入力する。

> table(

練習データ

[,2])

2 2.5 3 4 4.5 5 6 7 7.5 8 2 2 3 3 2 2 1 3 1 2

１行目にデータ内の数値が、２行目に対応する度数が示されている。これは列番号をつ かった方法である。データの列の指定の方法には、いくつかあり、つぎのように項目名 を入力しても同じ結果が得られる。

> table(練習データ[“a1”])

また、項目名をつかった方法で、「$」をつかっても同じ結果となる。

> table(

練習データ

$a1)

2.4 ヒストグラムを描くヒストグラムを描くヒストグラムを描くヒストグラムを描く

つぎに、ヒストグラムを描いてみよう。「

hist

」関数をつかって、つぎのように入力する。

> hist(練習データ["a1"])

以下にエラー

hist.default(

練習データ

["a1"]) :

'x'

は数値でなければなりません

入力をすると、その後にエラーメッセージがでてきた。「hist」関数では、パラメータ指 定がこのような書式ではエラーとなることがわかる。つぎのように入力すると今度はう まくいく。

> hist(

練習データ

$a1)

- 7 -

図３ ヒストグラム

2.5 2変数間の関連を散布図でみる変数間の関連を散布図でみる変数間の関連を散布図でみる変数間の関連を散布図でみる

項目「

a1

」と「

a2

」の間の相関の度合いを調べよう。まず、散布図を描こう。「

plot

」 関数をつかう。横軸を「a1」、縦軸を「a2」とする散布図を描くためにつぎのように入 力する。

>plot(練習データ$a1,練習データ$a2)

こうして左図が表示される。

図４ 散布図

- 8 - 2.6 相関係数を求める相関係数を求める相関係数を求める相関係数を求める

「

a1

」と「

a2

」の間の

Peason

相関係数を求めよう。そのために「

cor

」関数をつかう。

つぎのように入力する。

> cor(練習データ$a1,練習データ$a2)

[1] 0.7393481

これで、相関係数がおよそ

0.74

であり、かなり正の相関があることがしめされた。

2.7 回帰直線を描回帰直線を描回帰直線を描回帰直線を描くくくく

この散布図の上に重ねて回帰直線を描こう。これは

EXCEL

のようにメニューの１つ のボタンをクリックするだけで回帰直線が描けるわけではない。一気にやれない。まず、

回帰分析をおこない、その結果を利用して直線を描くという手順でおこなう。

R

には、

回帰分析で単回帰分析をおこなう「lm」関数がある。従属変数を「練習データ$a2」と し、独立変数を「練習データ

$a1

」とする単回帰分析をおこなう。そのために、つぎの ように入力する。

> lm(

練習データ

$a2~

練習データ

$a1)

すると、つぎのように結果が表示される。

Call:

lm(formula =

練習データ

$a2 ~

練習データ

$a1)

Coefficients:

(Intercept)

練習データ$a1

1.831 0.716

１

回帰直線の式は

- 9 -

「練習データ$a2」=0.716×「練習データ$a1」+ 1.831

であることを示している。回帰直線の傾きが

0.716

で縦軸の切片が

1.831

である。

この結果を利用して直線のグラフを描くのである。一連の入力はつぎのようにする。

>

単回帰結果 <- lm(練習データ$a2~練習データ$a1)

> plot(練習データ$a1,練習データ$a2)

> abline(

単回帰結果

)

単回帰分析の結果を「単回帰結果」という名称のデータフレームに入れた。つぎに散 布図を描き、そのグラフに回帰直線を重ねて描くという手順である。「abline」関数は、

重ねてグラフを描く機能をもっている。

グラフはつぎのように作成される。

図５ 直線回帰を付けた散布図 2.8 ２群における平均値の２群における平均値の２群における平均値の２群における平均値の比較の比較の比較の検定比較の検定検定検定

項目「a1」と「a2」の２群の母集団における平均値が等しいかどうかを検定しよう。

手元にあるデータの「a1」と「a2」の平均値は「mean」関数をつかって、つぎのように 求めることができる。

- 10 -

> mean(練習データ$a1) [1] 4.738095

> mean(

練習データ

$a2) [1] 5.22381

>

手元のサンプルの２群の平均値は明らかに異なるが、これが母集団における母平均値 が異なるということを意味しない。母集団からサンプルされた少ない数からなる標本の ２群の平均値が異なっていたということである。このサンプルのばらつきの程度から、

この母集団の平均値が明らかに異なると言えるかどうかを調べる。２群の平均値が等し いかどうかの検定として、

t

検定が使える。ただし、この検定が正しく使えるためには、

「２群の母分散が等しい」という前提条件がある。

したがって、まず「

a1

」と「

a2

」の２群における等分散の検定をおこなおう。「

var.test

」 関数によってその検定ができる。下に入力とその結果を示している。

> var.test(

練習データ

$a1,

練習データ

$a2)

var.test(

練習データ

$a1,

練習データ

$a2)

F test to compare two variances

data:

練習データ$a1 and 練習データ$a2

F = 1.0661, num df = 20, denom df = 20, p-value = 0.8876 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval:

0.4325804 2.6273559 sample estimates:

ratio of variances 1.066087

上の結果は、帰無仮説「２群の分散は等しい」の検定結果において、

p

値が

0.8876

と なったことで、有意水準５％ではこの帰無仮説は棄却されなかったことを示している。

したがって、２群の分散は明らかに異なるとは言えなかった。この結果を踏まえて

t

検 定をおこなう。等分散が証明されたわけではないが、否定されなければいちおう等分散 であるとみなすことにしよう、と考えるのが慣例となっている。

- 11 -

> t.test(練習データ 練習データ)

t.test(練習データ$a1,練習データ$a2,var.equal=TRUE)

Two Sample t-test

data:

練習データ

$a1 and

練習データ

$a2 t = -0.7905, df = 40, p-value = 0.4339

alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval:

-1.7274745 0.7560459 sample estimates:

mean of x mean of y 4.738095 5.223810

2.9 コンピュータシミュレーションコンピュータシミュレーションコンピュータシミュレーションコンピュータシミュレーション

（１）

（１）

（１）

（１） 乱数の発生

疑似乱数をつかって、さまざまな目的をもったコンピュータシミュレーションをおこ なうことができる。このシミュレーションはモンテカルロ・シミュレーションとよばれ る。

R

で疑似乱数を発生することができるために、モンテカルロ・シミュレーションも 可能である。

一様乱数を発生させるためには、「

runif

」関数を利用する。これは、開区間

(0 1)

での連続の一様乱数を発生する。いま、５個の乱数を発生させよう。つぎのように入力 する。すると、すぐに結果があらわれる。

> runif(n=5)

[1] 0.003309959 0.928781783 0.888635009 0.378752298 0.157205250

>

乱数を２０個発生させ、これらのヒストグラムを描き、乱数の値の発生度数がどれほ ど一様な分布になっているかを確認しよう。つぎのように入力する。

- 12 -

>

乱数

1 <- runif(n=20)

> hist(乱数 1)

すると、図５の表示がなされる。

図５ ヒストグラム 図６ ヒストグラム

つぎに、乱数を３００個に増やして、ヒストグラムを描いてみよう。

>

乱数

2 <- runif(n=300)

> hist(乱数 2)

こんどは図６が表示され、図５に比べ、かなり一様な分布に近づいてきたことがわかる。

乱数の出現範囲を変更するためには、「runif」関数にオプションを追加記述する。つぎ のように入力すると、開区間

(0 10)

の範囲で乱数が発生する。

> runif(n=5,min=0,max=10)

１０個ずつの乱数を２回発生させ、(x₁

,x

2

,・・・,x

₁₀

)と(y

₁

,y

2

,・・・,y

₁₀

)

の２つのグル ープをつくり、２グループの対応する値の組、(x_i,

y

i

), i=1,2,,,10

が変数（X,Y）の実現値 とみなし、

X

と

Y

の間で相関関係があるかどうかを調べてみよう。

- 13 - 以下のように入力する。

>

乱数グループ

1<-runif(10)

>

乱数グループ

2 <-runif(10)

> cor(乱数グループ 1,乱数グループ 2)

[1] -0.2463142

相関係数は、 -0.2463142となった。

つぎにサンプル数を増やし、各グループ５０個の乱数をつくり、

X

と

Y

の間の相関 係数を求めよう。

>

乱数グループ

1<-runif(50)

>

乱数グループ

2 <-runif(50)

> cor(

乱数グループ

1,

乱数グループ

2) [1] 0.03221177

上の結果から、乱数を１０個から５０個に増やすと、相関係数は、およそ、

-0.25

か らおよそ、0.03に変化し、０に近づき、一様分布をもつ

X

と

Y

の変数間の相関関係は ほとんどなくなることがわかった。乱数には一様乱数以外に、正規分布、二項分布、β 分布、指数分布、幾何分布、対数正規分布など、さまざまな確率分布にしたがう乱数を 関数によって発生させることができる。

（２） 事象の生起のシミュレーション

一様乱数をつかって、確率

0.9

で生じる事象を考え、実際にどのような頻度で事象が 生じるかをシミュレーションしよう。そのために、０と１の間の一様乱数

x

を発生させ、

x

が

であれば、注目する事象が生じたとみなし、離散変数

m

が

m=1 であるとする。x

が

であれば、注目した事象が起こらなかったとみなして、m=0 とする。これを

R

ではつ ぎのように入力する。

0 < x ≤ 0.9

0.9 < x ≤ 1

- 14 -

> x <- runif(n=10)

> x

[1] 0.29708019 0.63871612 0.69683403 0.02948449 0.37808815 0.01236702 0.97239158 0.16094847 0.18708040 0.59387856

> m <- cut(x,c(0,0.9,1))

> m

[1] (0,0.9] (0,0.9] (0,0.9] (0,0.9] (0,0.9] (0,0.9] (0.9,1] (0,0.9] (0,0.9] (0,0.9]

Levels: (0,0.9) (0.9,1)

> levels(m) <- c(1,0)

> m

[1] 1 1 1 1 1 1 0 1 1 1 Levels: 1 0

この結果より、７回目にのみ注目している事象は生じなくて、１０回のうち、９回事 象が生じたことがわかった。

参考文献 参考文献 参考文献 参考文献

[1] 山田剛史、杉沢武俊、村井潤一郎、Rによるやさしい統計学 オーム社/開発局