Rを使用したデータ解析・グラフ描き

データの特徴を大まかにつかむ （記述統計） mean(x$diameter) 平均値μ sd(x$diameter) 不偏標準偏差σy：データのばらつきを表す値 データが正規分布している場合、前データの約68％が区間[μ-σy, μ+σy]の中に入っている。 （出典:http://www.cap.or.jp/~toukei/kandokoro/html/14/14_2migi.htm）

データの特徴を大まかにつかむ （頻度分布・箱ひげ図） hist( x$diameter ) ヒストグラムを描く。 boxplot( x$diameter ) 箱ひげ図（box plot）を描く。 (出典: R-tips）

グラフを仕上げる 箱ひげ図のy軸にタイトルを入れる。 スクリプトエディターからRを実行する。 「ファイル」→「新しいスクリプト」 軸タイトルのない図 →　ダメ スクリプトエディターからRを実行する。 Rの命令を保存し、後から再現・再利用できるようにする。 「ファイル」→「新しいスクリプト」 以下のコマンドを書き、選択し、「…Rコードを実行」 　boxplot( x$diameter, ylab = "幹直径 (cm)" )

データの平均値をもう少し詳しく 300個のデータから計算した平均値は、3個のデータから計算した平均値より信頼がおけそうだ。 →平均値のばらつき（標準偏差）を計算する。 標準誤差（σf）：平均値の標準偏差 sd(x$diameter) / sqrt( length(x$diameter ) )

二つの量の間の関係 樹木の幹の直径と樹高との関係は? 直径が大きいとき樹高も高く、直径が小さいとき樹高も低いという関係がありそうだ！ →グラフに描いてみる 直径が大きいとき樹高も高く、直径が小さいとき樹高も低いという関係がありそうだ！ 両者の関係を簡単な式で要約する。 直径が特定の値であるとき、樹高の期待値は？ 回帰直線を引く。データ点に最も「近い」直線を引く。 回帰分析をして、それを図示する（回帰直線を引く）。

前ページを実行するRコード plot( x$diameter, x$height, xlab = "胸高直径（cm）", ylab = "樹高（m）") xlab, ylabで軸のタイトルを指定する。 result <- lm(height ~ diameter, data = x) lm()は線形モデルをあてはめる関数。ここでは、回帰分析を実行し、結果をresultに収める。 height ~ diameter：heightをdiameterで回帰する。heightは応答変数、diameterは説明変数。 data = x：データフレームを指定する。データフレームを指定しておけば、"diameter"の前に"x$"は必要ない。 summary( result ) resultに格納された結果をコンパクトに表示する。

回帰とは？ (単回帰の場合） y = a + bx 応答変数ｙを説明変数xの一次関数で説明する。 グラフ上で、直線関係をあてはめる。 最小自乗法で一次関数のパラメータ（切片と傾き）を求める。 y = a + bx + cx2 y = a + bx1 + cx2 などを仮定する回帰もある。

最小自乗法 （least square method） Yi:　一次式によるyiの予測値 y の場合 (xi, yi) (xi, Yi) を満たすa, bは x ●はデータ。最小自乗法で回帰直線を決める。■は個々のxに対応したyの予測値。予測値は線上にある。 : 最小自乗法によって決められたa, b : x, yの平均

回帰分析の結果（回帰線）の図示 xForLine <- seq( from = 3, to = 8, length = 10 ) yForLine <- coef(result)[ 1 ] + coef(result)[ 2 ] * xForLine 回帰線を描くためのyの値を作る coef(result)[1]: 回帰線の切片 coef(result)[2]: 回帰線の傾き lines( xForLine, yForLine ) 最後に描いたグラフに線を 描き加える

summaryの中身 Coefficients: Estimate Std. Error t value Pr(>|t|) 直線の切片、直径にかかる係数 左の数字の標準誤差 左の数字=0という帰無仮説を検定した結果：p値 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.01515 1.19657 0.013 0.99021 diameter 1.01515 0.20521 4.947 0.00113 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.054 on 8 degrees of freedom Multiple R-squared: 0.7536, Adjusted R-squared: 0.7228 F-statistic: 24.47 on 1 and 8 DF, p-value: 0.001125 R2＝Rじじょう＝決定係数 Coefficients: 係数 Std. Error: 標準誤差 Intercept: 切片 Estimate: 推定値 モデル全体のp値 このモデルは有意水準1％で有意である。

決定係数（R2）とは？ 決定係数は次の式で与えられる統計量でモデルの当てはまりの良さを表す。 単回帰、重回帰、分散分析、共分散分析など応答変数が正規分布に従うと仮定するモデルで使用する 従属変数の分散の中でモデルによって説明される割合を説明する。 応答変数yの数値の散らばり（分散） ＝モデルによって説明される分散＋モデルによって説明されない分散 モデルによって説明されない分散＝残差平方和 0〜1の値をとり、0で説明力なし、1で誤差のない説明（予測）を意味する Yi:モデルによるyiの予測値

決定係数（R2）とは？ y yの平均値 x 点線の長さの自乗の和＝応答変数yの数値の散らばり

決定係数（R2）とは？ ●のyの値：応答変数の実測値 ■のyの値：予測値（実測値と同じxの値に対応する直線上の点のyの値） (xi, Yi) (xi, yi) ●のyの値：応答変数の実測値 ■のyの値：予測値（実測値と同じxの値に対応する直線上の点のyの値） 点線の長さ：予測値と実測値のずれ（残差） データ全体のずれ：点線の長さの自乗の和（残差平方和） 「真ん中を通るように」：残差平方和を最小にするようにa, bを選ぶ＝最小自乗法 R2=1-（残差平方／データ全体のばらつき）

決定係数（R2）とは？ R2: 大 R2: 小

説明変数のp-値とは 帰無仮説の下で実際にデータから計算された統計量よりも極端な統計量が観測される確率を、p値(p-value)という y = a + bx y: 応答変数 x: 説明変数 a, b: パラメータ、係数 帰無仮説1： a=0 対立仮説1: a≠0 帰無仮説2： b=0 対立仮説2: b≠0 b=0だとするとyとxは無関係（xが増減しても、yの値に影響を与えない）となる。そのような(y, x)の組み合わせを複数（データ）をとってくると、通常上記左のようなグラフとなる。しかし、yとxが無関係であっても、ごくまれに一見y=a+bxの関係がありそうに見える上記右のようなデータが得られることがある。このようなことが起こる確率をp-値という。これが小さいと、この係数の値は帰無仮説のもとでは起こりにくいことを意味する。起こりにくいことがたまたま起こったと考えるよりは、このデータを起こりにくいと判断する根拠（帰無仮説）が間違っていると考える方が合理的である。したがって、p-値が小さいほど、yとxは関係がありそうだと考えることができる。 y x y x

処理の効果（分散分析） 処理（treatment）: 施肥（fertilized）、対照（control） 施肥をすると樹木は大きくなるか？ グラフを描いてみる 大きくなりそうだ！ →効果を検定してみる

分散分析1 例題の分散分析で知りたいこと 樹木の直径が施肥によって変化するという考えを数式で表すと… 施肥をすると樹木の直径はいくつになり、しないといくつになるか？ 施肥をする・しないによって生じる直径の差は偶然生じたとは考えられないほど大きいか（有意差）？ 樹木の直径が施肥によって変化するという考えを数式で表すと… y = A A: 処理の効果を示す質的変数 質的変数とは… 処理がcontrolのときは、A1 処理がfertilizedのときは、A2 の値をとる変数

分散分析2 実際にRが使う式は… y = a + A 処理がcontrolのときは、0 a: 切片 A: 処理の効果を示す質的変数 処理がfertilizedのときは、A2 の値をとる変数 A2は、controlを基準とした時のfertilizedの効果 aは前ページのA1と同じ値になる。 なぜ、切片を考えたか? 質的変数が2つ以上になったとき、こうした方が処理しやすい。 なぜ、fertilizedではなく、controlが0とされているのか？ アルファベット順に最も早いものが0とされる。

グラフ描画・分散分析を実行するRコード result2 <- lm( diameter ~ treatment, data = x ) plot( x$treatment, x$diameter, xlab = "処理", ylab = "胸高直径（cm）" ) result2 <- lm( diameter ~ treatment, data = x ) summary( result2 ) anova( result2 ) lm( diameter ~ treatment, data = x ) という書き方が回帰分析とほぼ同じであることに注意（説明変数名のみ異なる）。説明変数が量的変数か質的変数かの違いがあるが、回帰分析と分散分析はよく似ている。

summaryの中身 Coefficients: Estimate Std. Error t value Pr(>|t|) 切片、controlを基準とした時のfertilizedの効果 左の数字=0という帰無仮説を検定した結果：p値 左の数字の標準誤差 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.8000 0.7071 6.788 0.000139 *** treatmentfertilized 1.6000 1.0000 1.600 0.148266 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.581 on 8 degrees of freedom Multiple R-squared: 0.2424, Adjusted R-squared: 0.1477 F-statistic: 2.56 on 1 and 8 DF, p-value: 0.1483 Coefficients: 係数 Std. Error: 標準誤差 Intercept: 切片 Estimate: 推定値 Rじしょう＝R2＝決定係数

anovaの中身（分散分析表） Analysis of Variance Table Response: diameter 　 　　 Df Sum Sq Mean Sq F value Pr(>F) treatment 1 6.4 6.4 2.56 0.1483 Residuals　 8 20.0 2.5 施肥処理の効果 F検定量 P値 残差 分散分析表：要因（ここでは、施肥処理）の効果を、F検定するために数値を整理した表。 応答変数のデータの散らばり（分散）が、要因で説明される分散と要因によって説明されない分散（残差分散）に分けられ、後者に対する前者の大きさが判定される。判定にはF分布が用いられる。 この表の結果、施肥処理は有意な効果があるとは言い難いと結論づけられる。

多重比較 3群以上の平均値を比較する。 分散分析による結論 2つずつを取り出し、有意差の検定をする→ダメ →「3群の平均値すべてが同じではない。」または、「3群の平均値すべてが同じ。」 「平均値すべてが同じ」が棄却されても、何と何の間に有意な差があるのかわからない 2つずつを取り出し、有意差の検定をする→ダメ 検定を繰り返すことにより、１回のみ検定を行った場合より第一種過誤率が大きくなってしまう。 →多重比較をする！

多重比較の実行 仮想データの作成 d <- c( rnorm( 10, 1, 0.5 ), rnorm( 10, 2, 0.5 ), rnorm( 10, 2.1, 0.5 ) ) t <- c( rep( "control", 10 ), rep( "nitrogen", 10 ), rep( "light", 10 ) ) xx <- data.frame( diameter = d, treatment = t ) 説明 c():列挙された要素からベクトルを作る関数 rnorm( a, b, c ): 平均値b, 標準偏差cの正規分布からランダムにa個数字をサンプルし、ベクトルを作る関数 rep( a, b ): aをb回繰り返したベクトルを作る関数 data.frame(): ベクトルからデータフレームを作る関数。括弧の中に、「データフレーム中の列名　＝　ベクトル名」をならべる。 練習1：平均値0.8、標準偏差3の正規分布から数字を1000個ランダムに取り出し、取り出した数字の頻度分布を描く。

多重比較の実行 グラフの作成 plot( xx$treatment, xx$diameter) 練習2：箱ひげ図に縦軸、横軸のタイトルをつけてください。

多重比較の実行 分散分析 result3 <- lm( diameter ~ treatment, data = xx ) summary( result3 ) anova( result3 ) 質問： controlの処理で、樹木の直径は平均でいくらになりますか？ controlと比べ、nitrogen, lightの平均値はいくら高いですか? 処理の効果は有意でしたか?

多重比較の実行 TukeyHSD( aov( diameter ~ treatment, data = xx ) ) aov(): 分散分析を実行する関数、内部的にlm()を使用している。 TukeyHSD(): Tukey法を使用した多重比較を実行する関数。

TukeyHSD()の結果 全体で95%の信頼水準が確保されている Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = diameter ~ treatment, data = xx) $treatment diff lwr upr p adj light-control 1.01601898 0.4482282 1.5838098 0.0003947 nitrogen-control 0.93956774 0.3717769 1.5073586 0.0009529 nitrogen-light -0.07645124 -0.6442420 0.4913396 0.9405481 処理後との平均間の差 処理後との平均間の差の95%信頼区間の下限 処理後との平均間の差の95%信頼区間の上限 この区間に0を含まなければ、差は95%水準で有意

多重比較の結果をグラフに追加する。 plot( xx$treatment, xx$diameter, xlab = "処理", ylab = "直径（cm）" ) mtext( "a", side = 1, line = -13, at = 1.2 ) mtext( "b", side = 1, line = -22, at = 2.2 ) mtext( "b", side = 1, line = -23, at = 3.2 ) グラフに書き入れる文字 文字列を書き込む余白位置を表す番号(1：下，2：左，3：上，4：右) 垂直方向（文字の方向に対して）の文字位置（行数） 水平方向（文字の方向に対して）の文字位置（x軸の値）

練習の答え 練習1 練習2 z <- rnorm( 1000, 0.8, 3 ) hist( z ) あるいは hist( rnorm( 1000, 0.8, 3 ) ) 練習2 plot( xx$treatment, xx$diameter, xlab = "処理", ylab = "直径（cm）" )