理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

二元配置分散分析 (1)

二元配置分散分析

投稿日2017.3.8
更新日2017.8.17

 

2要因ともに対応なし、繰り返し数がすべて等しい場合 

f:id:yoshida931:20170816111011p:plain

                                                                                             

例)

要因1(水準A1、水準A2)、要因2(水準B1、水準B2、水準B3)

f:id:yoshida931:20170816111320p:plain

 

                     (仮想データ使用)

 

Rの関数で分散分析表を作成してみます 

各セルにベクトル名をつけて処理してみます

それぞれのベクトルは繰り返し数5になっています

f:id:yoshida931:20170817104436p:plain

s1<-c( 7,11,8,10,9 )

s2<-c( 9,11,10,9,11 )

s3<-c( 21,22,19,18,22 )

s4<-c( 23,21,26,27,22 )

s5<-c( 12,16,16,17,15 )

s6<-c( 9,8,11,10,9 )

 

全データのベクトル

s<-c( s1,s2,s3,s4,s5,s6 ) 

 

水準名をつけます

A<-c(rep("A1",15),rep("A2",15))
B<-c(rep("B1",5),rep("B2",5),rep("B3",5),rep("B1",5),rep("B2",5),rep("B3",5))

 

イメージ

f:id:yoshida931:20170817104553p:plain


summary ( aov ( s ~ A * B) )

Df Sum Sq Mean Sq F value Pr(>F) 
A 1 67.5 67.5 21.89 9.40e-05 ***
B 2 73.3 36.6 11.88 0.000259 ***
A:B 2 850.2 425.1 137.87 6.94e-14 ***
Residuals 24 74.0 3.1

f:id:yoshida931:20170626185319p:plain

 

この分散分析表を平方和から求めてみます

 

平均

二元配置分散分析交互作用と分散分析表を求めるためには

次の4つの平均が必要になります

平均1(総平均)

f:id:yoshida931:20170816140407p:plain

平均2(群Aの平均)

f:id:yoshida931:20170816113234p:plain

平均3(群Bの平均)

f:id:yoshida931:20170816113610p:plain

平均4(群(A,B)の平均)

f:id:yoshida931:20170816120956p:plain

 

 

下の表を使用して平方和の分解を行います

f:id:yoshida931:20170816182726p:plain

注意) 一元配置の表と異なるので注意してください.( Ai , Bj ) の中にそれぞれn個( x1,x2,...,xn  、繰り返し=n )のデータが入っており、xaij は( Ai , Bj ) の平均値を示しています. 

 

2元配置の構造モデル

実測値 = 全体平均 + Aの主効果 + Bの主効果 + ABの相互作用 + 誤差項

xij = xa + ( xai, - xa ) + ( xa-j - xa ) + ( xaij - xai, - xa,j + xa ) + ( xij - xaij )

xij - xa = ( xai, - xa ) + ( xa-j - xa ) + ( xaij - xai, - xa,j + xa ) + ( xij - xaij )

両辺の二乗和

∑(f)∑(g)∑(n) ( xij - xa )^2 

    = ∑(f) g*n*( xai, - xa ) ^2                        Aの主効果

      + ∑(g) f*n*( xa-j - xa )                               Bの主効果  

         + ∑(f)∑(g)n* ( xaij - xai, - xa,j + xa ) ^2      ABの相互作用

            + ∑(f)∑(g)∑(n)xij - xaij ) ^2                            誤差項

 

イメージ

f:id:yoshida931:20170816212841p:plain

 

自由度

A群間の変動:f - 1

B群間の変動:g - 1

交互作用ABの変動: ( f - 1 )*( g - 1 )

群内変動:a*b*( n - 1 )

 

平方和の分解(平均は最後にaを付けることにします)

全体平均

xa<- mean ( s )

各セルの平均

s1a<-mean (s1)

s2a<-mean (s2)

s3a<-mean (s3)

s4a<-mean (s4)

s5a<-mean (s5)

s6a<-mean (s6)

A1,A2の平均ベクトル

Aa<-c ( mean ( c( s1 , s2 , s3 ) ) , mean ( c( s4 , s5 , s6 ) ) )

B1,B2,B3の平均ベクトル

Ba<-c ( mean ( c( s1 , s4 ) ) , mean ( c( s2 , s5 ) ) , mean ( c( s3 , s6 ) ) )

構造モデルから求めた二乗和を算出します

全体 ∑(f)∑(g)∑(n) ( xij - xa )^2

sum ( s - xa )^2

= 1064.967

Aの主効果  ∑(f) g*n*( xai, - xa ) ^2      

3*5*sum((Aa - xa) ^2)  
= 67.5            

Bの主効果 ∑(g) f*n*( xa-j - xa )

2*5*sum((Ba - xa) ^2)  

= 73.26667                            

ABの相互作用 ∑(f)∑(g)n* ( xaij - xai, - xa,j + xa ) ^2

5* ( s1a -  mean ( c( s1 , s2 , s3 ) ) - mean ( c( s1 , s4 ) ) + xa )^2+

5* ( s2a - mean ( c( s1 , s2 , s3 ) ) -mean ( c( s2 , s5 ) ) +xa )^2+

5* ( s3a - mean ( c( s1 , s2 , s3 ) ) - mean ( c( s3 , s6 ) ) +xa)^2+

5* ( s4a -  mean ( c( s4 , s5 , s6 ) )-mean ( c( s1 , s4 ) ) + xa)^2+

5* ( s5a -  mean ( c( s4 , s5 , s6 ) )-mean ( c( s2 , s5 ) ) +xa )^2+

5* ( s6a -  mean ( c( s4 , s5 , s6 ) )- mean ( c( s3 , s6 ) ) +xa)^2

= 850.2

誤差項∑(f)∑(g)∑(n)xij - xaij ) ^2

sum(( s1 - s1a )^2) +
sum((s2 -s2a)^2) +
sum((s3 - s3a)^2 )+
sum((s4 - s4a)^2) +
sum((s5- s5a)^2) +
sum((s6 - s6a)^2)

= 74

 

 これでR関数で出力した分散分析表が作成できます

分散分析表

f:id:yoshida931:20170817191845p:plain

一元配置分散分析

一元配置分散分析

投稿日2016.11.2

更新日2017.8.15

 

言葉の整理

要因:実験結果に影響を与える要素

因子:研究対象となる要因( 母平均に差をもたらすと考えられる )

水準:要因を分ける条件( 要因に含まれる項目 )

群:=水準

平方和:偏差の二乗の合計

変動:=平方和

平均平方和:平方和を自由度で割った値(分散)

〇元配置分散分析:〇は因子の数

 

以下、このブログでは「要因」、「群」を使用します

 

分散分析とは

 分散分析は、3群以上のそれぞれの群の母平均が等しいという帰無仮説のもとに検定を行う分析方法です.しかしt検定を繰り返してはいけません.例えば3群の比較を5%有意水準でそれぞれt検定を行った場合、帰無仮説は「3通り全ての検定において差がみられない」となり、対立仮説は「少なくとも一組には差がみられる」となります.つまり、有意水準は0.05ではなく1- 0.95^3 = 0.142 で検定したことになります.つまり、第一種の過誤が生じる危険性が14.2%となります.

(http://yoshida931.hatenablog.com/entry/2017/02/24/171054)

 分散分析では、各群の母平均の推定を行うのではなく、各群が及ぼす効果(要因の効果)の分散を推定します.そのために群間の平均の変動各群内の誤差変動から平均平方和 ( 分散 ) を求め、要因の効果の有無をF分布から推測します.

 

よくみかける一元配置のデータ 

対応のない要因 繰り返しn回

f:id:yoshida931:20170815162240p:plain

 例えば実験A、B、Cをそれぞれ10回ずつ行ったとします.実験ABCは水準、10回は繰り返し数となります.実験結果のバラつきは分散として得られ、分散の原因として次の2点が考えられます.

原因1)実験A、B、Cの効果によるもの(要因の効果)

原因2)被験者の問題、測定の誤差、偶然などによるもの(誤差項)

 

データの構造を表すモデル

測定値全平均 + 要因の効果 ( 原因1 ) + 誤差項 ( 原因2 )

誤差項は互いに独立で、正規分布に従うことを仮定しています. 

 

分散分析の帰無仮説

 各水準の母平均が等しく( Aの母平均=Bの母平均=Cの母平均 ) 、ABCの母平均の分散は0となる

 「原因1によるばらつきの大きさ(平均平方和)」と「原因2によるばらつきの大きさ(平均平方和)」の比率から、要因の効果による分散が「0」かどうかを検定します.検定の結果、有意水準が棄却された場合には、要因の効果が認められると推測します.つまり 分散分析は実験A、B、Cの母平均を比較検証するのではなく、その3つの「効果」の分散を推定するものです.

 

例題)

水準数は3、各水準の標本サイズ ( 繰り返し数 ) は5
治療後の満足度調査を行った.20名を無作為に分けて治療A,B,C,D後にVASで調査した.治療A,B,C,Dで満足度に差があったかどうかを有意水準0.05で検定せよ.
治療=要因、治療A,B,C,D=水準(群)

A B C D
10 10 5 9
6 5 4 5
10 5 11 2
9 12 4 3
10 4 6 1

 以後、水準を「」と記載します

 

一覧表 

f:id:yoshida931:20170815205053p:plain

 

各群の平均(各群の平均)

( 群平均 <- colMeans ( x ) )

A B C D
9.0 7.2 6.0 4.0

 f:id:yoshida931:20161123211908p:plain

 

全データをベクトルにします

all<-c( x$A , x$B , x$C , x$D  )

[1] 10 6 10 9 10 10 5 5 12 4 5 4 11 4 6 9 5 2 3 1

 

水準をベクトルに変換します

levels <-c ( rep ( "A" , 5 ) , rep ( "B" , 5 ) , rep ( "C" , 5 ) , rep ( "D" , 5 ) )

 

さらに各水準もベクトルに変換します

( levels2 <- factor (levels ) )

この時点で群名を数値に変換します

[1] A A A A A B B B B B C C C C C D D D D D
Levels: A B C D

str(levels2)   #数値になっているか確認してみます

Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1 2 2 2 2 2 ...

 

Rの関数を使用すれば分散分析表の完成です
summary ( aov ( all ~ levels2 ) )

 

                  Df   Sum Sq   Mean Sq    F value    Pr(>F) 

levels2       3      66.15        22.05        2.579    0.0898 .
Residuals 16    136.80         8.55
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 一元配置の分散分析表

f:id:yoshida931:20170816084944p:plain

 

この一覧表から分散分析表を作成してみます

f:id:yoshida931:20170816162212p:plain

 

データの構造を表すモデルから平方和の分解

測定値 = 全平均 + 要因の効果 ( 原因1 ) + 誤差項 ( 原因2 )

xki = xa + ( xa_i - xa ) + ( xki - xa_i )

xki - xa =  ( xa_i - xa ) + ( xki - xa_i )

( xki - xa )^2 = ( ( xa_i - xa ) + ( xki - xa_i ) )^2

( xki - xa )^2 = ( xa_i - xa ) ^2 + ( xki - xa_i ) ^2 + 2 ( xa_i - xa ) * ( xki - xa_i ) 

両辺の二乗和

∑( xki - xa )^2 =  ∑∑  ( xa_i - xa ) ^2  +  ∑∑( xki - xa_i ) ^2 + 2∑∑( xa_i - xa ) * ( xki - xa_i ) 

第三項

∑ ( xki - xa_i ) = 0 より ∑∑( xa_i - xa ) * ( xki - xa_i )   = 0

したがって

∑( xki - xa )^2 =  ∑∑  ( xa_i - xa ) ^2  +  ∑∑( xki - xa_i ) ^2 

全体の平方和 = 群間平方和 + 誤差平方和(群内平方和)

 

rbind ( all , levels2)  #行ベクトルを生成して行列を作る

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
全データ 10 6 10 9 10 10 5 5 12 4 5 4
水準2 1 1 1 1 1 2 2 2 2 2 3 3
[,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
全データ 11 4 6 9 5 2 3 1
水準2 3 3 3 4 4 4 4 4

一覧表

f:id:yoshida931:20170815204922p:plain

 

群平均行列  (平方和を求めるために必要)
( 群平均行列 <- matrix ( rep ( 群平均 , 5 ) , nrow = 5 , ncol = 4 , byrow = TRUE ) )

[,1] [,2] [,3] [,4]
[1,] 9 7.2 6 4
[2,] 9 7.2 6 4
[3,] 9 7.2 6 4
[4,] 9 7.2 6 4
[5,] 9 7.2 6 4

f:id:yoshida931:20161123211847p:plain

 

全データの平均
mean( all )
[1] 6.55

 

(データテーブル)ー(全データの平均
x - mean( all )

A B C D
1 3.45 3.45 -1.55 2.45
2 -0.55 -1.55 -2.55 -1.55
3 3.45 -1.55 4.45 -4.55
4 2.45 5.45 -2.55 -3.55
5 3.45 -2.55 -0.55 -5.55

 f:id:yoshida931:20161123211818p:plain

 

(群平均行列)ー(全データの平均
群平均行列 - mean(全データ)

[,1] [,2] [,3] [,4]
[1,] 2.45 0.65 -0.55 -2.55
[2,] 2.45 0.65 -0.55 -2.55
[3,] 2.45 0.65 -0.55 -2.55
[4,] 2.45 0.65 -0.55 -2.55
[5,] 2.45 0.65 -0.55 -2.55

f:id:yoshida931:20161123212336p:plain

 

(全データ)ー(群平均)
x - 群平均行列

A B C D
1 1 2.8 -1 5
2 -3 -2.2 -2 1
3 1 -2.2 5 -2
4 0 4.8 -2 -1
5 1 -3.2 0 -3 

f:id:yoshida931:20161123212620p:plain

 

 (全データー全平均)^2  この合計が総合平方和
#各データが全平均に近いと小さく、離れると大きくなる
( x - mean ( all ) )^2

A B C D
[1,] 11.9025 11.9025 2.4025 6.0025
[2,] 0.3025 2.4025 6.5025 2.4025
[3,] 11.9025 2.4025 19.8025 20.7025
[4,] 6.0025 29.7025 6.5025 12.6025
[5,] 11.9025 6.5025 0.3025 30.8025 

 

f:id:yoshida931:20161123212724p:plain 

 

群間^2 この合計が群間平方和
各群の平均が全平均に近いと小さく、離れると大きくなる
( 群平均行列 - mean ( all ) )^2

[,1] [,2] [,3] [,4]
[1,] 6.0025 0.4225 0.3025 6.5025
[2,] 6.0025 0.4225 0.3025 6.5025
[3,] 6.0025 0.4225 0.3025 6.5025
[4,] 6.0025 0.4225 0.3025 6.5025
[5,] 6.0025 0.4225 0.3025 6.5025

f:id:yoshida931:20161123212856p:plain

 

郡内^2 この合計が群内平方和
各データが群平均に近いと小さく、離れると大きくなる
( x - 群平均行列 )^2

A B C D
[1,] 1 7.84 1 25
[2,] 9 4.84 4 1
[3,] 1 4.84 25 4
[4,] 0 23.04 4 1
[5,] 1 10.24 0 9

f:id:yoshida931:20161123213016p:plain

 

平方和

総平方和
= ∑∑( 実測値- 全平均 )^2
sum( (x - mean ( all ) )^2)
[1] 202.95

群間平方和(水準間平方和、A間平方和)

= 繰り返し数 n × ∑( 各群の平均 - 全平均 )^2
= 5*sum( ( colMeans ( x ) - mean( all ) )^2)
= sum( ( 群平均行列 - mean ( all )^2)     #全く同じ意味です

[1] 66.15

誤差平方和群内平方和、残差平方和)

= ∑∑ ( 実測値 - 各群の平均 )^2

= sum( ( x$A - 9.0 )^2 ) + sum ( ( x$B - 7.2 )^2 ) +

            sum ( ( x$C - 6.0 )^2 ) + sum ( ( x$D - 4.0 )^2 )

= sum( (x-群平均行列)^2)
[1] 136.8

 

自由度(ncol列数、nrow行数)
群間の自由度:群の数(水準数) - 1
ncol ( x ) - 1
群内の自由度:(各群のデータ数 - 1)× 群の数
( nrow ( x ) - 1 ) * ncol ( x )
全体の自由度:全データ数 - 1
length ( 全データ ) - 1

 

平均平方
群間の平均平方 ( 原因2による分散 ) 
( sum( (群平均行列 - mean( 全データ ) ) ^2 ) ) / ( ncol ( x ) - 1 )

[1] 22.05

群内の平均平方 ( 原因1による分散 )
sum( (x - 群平均行列 ) ^2 ) / ( ( nrow ( x ) - 1 ) * ncol ( x ) )

[1] 8.55

全体の平均平方( 全体の分散)
sum( ( x - mean ( all ) ) ^2 )  / ( length ( all ) - 1 )

= var ( all )

= 10.68158

 

F検定

Fを統計量として要因効果の有無を検定します.

帰無仮説:要因の効果はない.A B C D の母平均が等しく、母平均の分散=0.

要因による分散が0のときF値は1に近づき、0でないときには1より大きくなります.

F値 = 群間の平均平方 ÷ 群内の平均平方
( F値 = ( sum ( ( 群平均行列 - mean ( 全データ ) ) ^2 ) ) / ( ncol ( x ) - 1 ) ) /
( sum ( ( x - 群平均行列 ) ^2 ) / ( ( nrow ( x ) - 1 ) * ncol ( x ) ) )

 [1] 2.578947

 

p値
pf ( F値 , 3 , 16 , lower.tail = FALSE )  

R関数の解と違うので、直接F値を挿入する
pf ( 2.578947 , 3 , 16 , lower.tail = FALSE )

[1] 0.08978798

 

結果

F>2.579でp=0.0898となり治療間で満足度の差はないと言えます

 

参考)

山田 剛史, 杉澤 武俊, 村井 潤一郎; Rによるやさしい統計学, オーム社 , 2008

小野 滋;読めば必ずわかる分散分析の基礎  http://elsur.jpn.org/resource/anova.pdf

標準偏回帰係数の求め方

標準偏回帰係数の求め方

投稿日2017.8.12

更新日2017.8.13

 

忘れないように載せておきます

例)予測変数が二つの場合の重回帰分析

母回帰方程式 Y = b0 + b1*X1i + b2*X2i + Ei

予測方程式 y = β0 + β1*x1i + β2*x2i +ei

 

標準偏回帰係数  sβ1、sβ2

sβ1 = ( yとx1iの相関係数 - ( yとx2iの相関係数 × x1iとx2iの相関係数 ) )  ÷ ( 1 - x1iとx2iの相関係数^2 ) 

sβ2 = ( yとx2iの相関係数 - ( yとx1iの相関係数 × x1iとx2iの相関係数 ) )  ÷ ( 1 - x1iとx2iの相関係数^2 ) 

 

Rで求める場合(とても簡単)

テーブルx( 項目はy,x1,x2 )  があったとします

後で比較するためのそのまま重回帰分析しおきます

summary ( lm ( y ~ x1 + x2 , x ) )

 

次にテーブルを標準化します

scx<- scale( x )              #これでテーブルの標準化が完了
sx<- data.frame( scx )    # データフレーム形式に戻す

重回帰分析をすれば標準回帰係数が算出できます

summary ( lm ( y ~ x1 + x2 , sx ) )

 

r-de-r様からの助言を忘れないように、解説を加えて追記しておきます

データを標準化なしで、重回帰の結果から標準偏回帰係数を算出します.

最小二乗法の結果から得られる次の解を利用します.

β1= sβ1 * ( yの不偏標準偏差÷x1iの不偏標準偏差 )

β2= sβ2 * ( yの不偏標準偏差÷x2iの不偏標準偏差 )

 

例)Rのデータセットirisを使用します.Sepal.Length を目的変数、他の変数(Sepal.Width、Petal.Length、Petal.Width)を説明変数とした重回帰分析です.

iris2 = iris[, 1:4]     #1列目から4列目を使用します

( ans2 = lm(Sepal.Length ~ ., iris2) )

 

lm(formula = Sepal.Length ~ ., data = iris2)

Coefficients:

 (Intercept)   Sepal.Width  Petal.Length   Petal.Width 

      1.8560        0.6508        0.7091       -0.556

 

この重回帰分析の結果から求める方法をr-de-r様から教えてもらいました.

次の式で全ての変数の不偏標準偏差を求めます

 

( s = sapply(iris2, sd) ) 

 

Sepal.Length  Sepal.Width Petal.Length  Petal.Width

   0.8280661    0.4358663    1.7652982    0.7622377

 

次に説明変数の偏回帰係数を取り出します

( coe2 = coefficients(ans2)[-1] )

 

Sepal.Width Petal.Length  Petal.Width

   0.6508372    0.7091320   -0.5564827

 

それぞれの偏回帰係数にそれぞれの不偏標準偏差を掛けて、目的変数の不偏標準偏差で割ります.これで標準回帰係数が算出されます

 

coe2 * s[2:4] / s[1]

 

Sepal.Width Petal.Length  Petal.Width

   0.3425789    1.5117505   -0.5122442

 

まとめて書けば、r-de-r様から助言いただいた式になります

coefficients(ans2)[-1] * s[2:4] / s[1]

 

データ標準化後の重回帰分析の解と比較してみます.

iris3 = as.data.frame(scale(iris2))  #データの標準化

ans3 = lm(Sepal.Length ~ ., iris3)

coefficients(ans3)[-1]

 

Sepal.Width Petal.Length  Petal.Width

   0.3425789    1.5117505   -0.5122442

 

同じ値が得られました.

 

テキスト片手にRに触れ、

助言に耳を傾けながら勉強したら

統計学は楽しくなります!

トービット回帰直線

トービット回帰直線

投稿日2017.8.4

打ち切りデータの場合、つまり天井効果や床効果が生じているデータの場合には、トービット回帰直線で分析します.

参考書はもちろん
豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012

 

Rのサンプルanscombeを使います.( n=11の4セット )

x1<-anscombe$x1

x4<-anscombe$x4
y1<-anscombe$y1
y4<-anscombe$y4

par(mfrow = c(1,2))
plot(x1,y1,xlab = "図A",ylab = "")
plot(x4,y4,xlab = "図B",ylab = "")

f:id:yoshida931:20170804195919p:plain



dev.off()

図A 通常の散布図
lm1<-lm(y1~ x1,data= anscombe)
plot(x1,y1)
abline(lm1$coefficients[1],lm1$coefficients[2])

f:id:yoshida931:20170804195838p:plain

視覚的にも何となくこの直線でフィットしているように見えます
x軸を予測値にして、残差を標準化して更に見やすくしてみます
ei<- y1-( lm1$coefficients[1]+ lm1$coefficients[2]* x1)
eis<-ei/sqrt( (sum(ei^2) )/(11 - 2) )  #標準化残差
plot(lm1$coefficients[1]+ lm1$coefficients[2]* x1, eis,xlab="図A(x軸:予測値)",ylab=" ", xaxt="n")
abline(h=0, lty=2, col=2)

f:id:yoshida931:20170804200055p:plain
図B 打ち切りデータがある場合
トービット回帰直線( 関数vglmを使用 )
通常の回帰直線(実践)
y4x4 <- lm(y4 ~ x4); summary(y4x4)
plot(x4,y4)
abline(summary(y4x4)$coefficients1,summary(y4x4)$coefficients2)
関数vglmを使用してトービット回帰直線を破線で加えます
y4x4t <- vglm(y4 ~ x4, tobit(Lower=8, Upper=19), trace=TRUE); summary(y4x4t)
abline(summary(y4x4t)@coefficients1,summary(y4x4t)@coefficients3,lty=2)

f:id:yoshida931:20170804200433p:plain

床効果で下に引っ張られていた回帰直線をやや上方修正できました.

残差分析

残差分析

Rのサンプル ChickWeight を使用
y<-ChickWeight$weight[1:10]
x<-ChickWeight$Time[1:10]
summary(lm(y~x))
plot(x,y,xlab = "生後日数",ylab = "体重")
abline(30.327 ,7.030)         #回帰直線の挿入

f:id:yoshida931:20170802161419p:plain

残差 ei をプロットすることで、回帰モデルからのズレを視覚化してみます

標準化残差Sie

f:id:yoshida931:20170802154446p:plain

ei<-y-(30.327+7.03*x)
eis<-ei/sqrt( (sum(ei^2) )/(10-2) )  #標準化残差
plot(eis,xlab="生後日数",ylab="残差", xaxt="n")
name<-c(0:9)*2
axis(side=1,1:10,labels=name)
abline(h=0, lty=2, col=2)

f:id:yoshida931:20170802161440p:plain

こうやって見たら、回帰直線からの「ズレ」がより理解できます

単回帰分析

投稿日2017.2.2
更新日2017.8.1


線形回帰分析
説明変数と目的変数を直線関係で傾向を示す。
説明変数と目的変数との関係を直線でモデル化する回帰分析。

 

非線形回帰分析
非線形関係でモデル化する回帰分析。
説明変数と目的変数を非直線的関係で傾向を示す。

以下のように定義します
観測値 ( xi , yi )
観測値の平均値 ( my , mx )

 

母回帰方程式を Yi = B0 + B1 * Xi  + Ei とします
   (Xiは確率変数ではなく確定した値)
母回帰係数 B0、 B1 
誤差項 Ei 期待値=0、分散=σ^2 (一定)

 

最小二乗法による回帰係数の決定
母回帰方程式の誤差項の二乗和(平方和)を最小にする回帰係数 β0 , β1 をもとめま
Yi = B0 + B1 * Xi  + Ei
誤差項 Ei = Yi -  B0 - B1 * Xi 
ここで ∑Ei ^2  =  S   とします
このSを最小にするβ0 , β1B0 , B1の推定値とします.

残差と誤差についてはこちら
27-3. 予測値と残差 | 統計学の時間 | 統計WEB

 

Ei = Yi -  B0 - B1 * Xi    より
S = ∑( Yi - ( B0 + B1 * Xi ) )^2
Sを最小にするB0B1を求めることで回帰直線を得る
SはパラメータB0B1の二次関数で最小値が存在します
そこでB0B1偏微分し、0に等しいとします

f:id:yoshida931:20170317200530p:plain

b0=B0,b1=B1として読んでください
いつか綺麗に書き直します

 

標本回帰係数 β1 , β0  
偏回帰係数   β1 = ∑ ( Xi - X平均 )*( Yi - Y平均 )  ÷  ∑ ( Xi - X平均 )
y切片   β0 = Y平均 - β1 * X平均

標本回帰方程式、標本回帰直線(予測式)
y' = β0 + β1*xi 

f:id:yoshida931:20170801140902p:plain

 

次のサンプルを使って考えていきます
yi<-c(114,124,143,158,166)
xi<-c(35,45,55,65,75)
         出典)統計学入門 ,東京大学出版会,p42,1991

 

平均を求めます
my <- mean ( yi )
141
mx <- mean ( xi )
55

最小二乗法より偏回帰係数を求めます
β1 <- sum ( ( yi - my ) * ( xi - mx ) ) / sum ( ( xi - mx )^2 )
1.38

Y切片
β0<- my - β1  * 55

65.1

回帰モデル
y' = 65.1 + 1.38 xi

 

平方和の分解
観測値の平方和(全変動)をSyとします
Sy = ∑ ( yi - my ) ^2
   = ∑ { ( yi - y' ) + ( y' - my ) } ^2
   = ∑ ( yi - y') ^2 + ∑ ( y' - my ) ^2 + ∑ ( yi - y' ) * ( y' - my )
となります
ここで∑ ( yi - y' ) * ( y' - my ) を考えてみます
予測値と残差の相関は0になります(証明は省略)
∑ ( yi - y' ) * ( y' - my ) = 0 

 

全変動
Sy
=  ∑ ( yi - y' ) ^2 + ∑ ( y' - my ) ^2
=1936

 

回帰により説明されない変動( 残差二乗和、誤差二乗和 
自由度 = 1
Se = ∑ ( yi - y' ) ^2   
=31.6

 

回帰により説明される変動
自由度 = (データ数 - 1) -1 = 5 - 2 = 3
SR∑ ( y' - my ) ^2
=1904.4

 f:id:yoshida931:20170801142946p:plain

 


この回帰方程式に実測値 ( xi , yi ) を当てはめると、
回帰式の値と実測値の間にはズレが生じます
ei = yi - y' = yi - ( β0 + β1*xi  ) 
ei : 回帰残差 (誤差とは異なります)

 

推定値の標準誤差 (Residual standard error)
誤差項 Ei の分散 σ^2 は、回帰方程式のあてはまりの良さ表します.
その値を回帰残差で推定します

                   f:id:yoshida931:20170802082254p:plain
s2 = ( ∑ ei ^ 2  ) ÷ ( n - 2 ) 
s2 <- ( sum ( ei ^ 2 ) ) / ( 5 - 2 )
sqrt ( s2 ) = 3.24551
sのことを推定値の標準誤差といい、この値が小さければ小さいほど回帰式は良く適合していると考えます.

 

単回帰の決定係数
観測値の平方和(全変動)Sy は回帰方程式で説明できる変動と説明できない変動に分けられる.決定係数は回帰方程式で説明できる変動の割合です.ピアソンの相関係数の二乗と同じ値.よいモデルは、残差二乗和 Se = ∑ ( yi - y' ) ^2 が小さく、寄与率が1に近いほど、よいモデルと言える.(重回帰の場合には、必ずしもそうではない)

 

決定係数(寄与率)
R^2 = 1 - ( ∑ ei ^ 2  ) ÷ ( ∑ ( yi - my )^2 )
  = ( ∑ ( y' - my ) ^2 ) / ( ∑ ( yi - my ) ^2 )
  = SR / Sy

      = 回帰方程式で説明できる変動の割合SR ÷ 全変動Sy 

sum( ( 65.1 + 1.38 * xi - my ) ^2 ) / sum ( ( yi - my ) ^2 )
0.9836777

        = (cor(yi,xi))^2      # "pearson" (default),

 

偏回帰係数の検定
帰無仮説:母集団において説明変数が目的変数を全く説明していない.(β1=0

母回帰係数 B0、 B1 の検定を行うためには、回帰係数 β0 , β1 の標本分布について知る必要があります.
β1の標本分布は

f:id:yoshida931:20170801163549p:plain

に従います.(証明省略)

ガウス・マルコフの定理を理解後に証明します

σ^2は未知なので回帰残差を使用してt 統計量を求めます.

f:id:yoshida931:20170814130832p:plain

 

サンプル
yi<-c(114,124,143,158,166)
xi<-c(35,45,55,65,75)
xiが
yiを説明するかどうかを検定してみます.
帰無仮説 H0 : β1=0 、xiはyiと正の相関がありそうなので、対立仮説は H1 : β1>0 となります.
サンプルのt値を求めてみます. 
s= Se / ( n - 2 ) = sum ( (yi - ( 65.1+1.38 * x) )^2) / ( 5 - 2) = 10.53333
xの偏差平方和 = sum( ( xi - mean( xi ) )^2) = 1000
t = 1.38 / sqrt ( 10.53333 / 1000 ) = 13.4461
p値 = (1 - pt ( 13.4461, 3 ) )*2 = 0.000889

有意水準1%で回帰係数 β1 の値が信頼できるので、ほぼ確実に母回帰係数B1は0ではないと言える.したがって「説明変数xが目的変数yに影響を与えており、回帰方程式は有意である」ことが考えられる.結果の解釈が重要になります.この結果から「サンプルデータが直線上にある」という結果にはなりません.

 

分散分析
F分布を利用して求めた回帰式が予測に役立つかどうかを検定します

帰無仮説は H0 「求められた回帰モデル y=65.1+1.38x では、x を説明変数として Y の変動は説明できない」

X^2とY^2 がそれぞれ自由度m1とm2のカイ二乗分布に従う互いに独立な確率変数の場合.F =  ( X^2 / m1 ) / ( Y^2 / m2 )  は自由度( m1 , m2 ) のF分布に従う.

 

 単回帰分析の分散分析表
f:id:yoshida931:20170801173907p:plain
F( 1 , n-2 ) は、自由度 [ n-2 , 1 ] のF分布に従う確率分布です.有意水準5%でp値≦0.05のときには、回帰モデルは目的変数yの変動の説明に有意に役立っていると判定しいます.

 

サンプルのF値を求めてみます
1904.4 / ( 31.6 / ( 5 - 2 ) ) = 180.7975

 

Rの関数を使って確認してみます
yi<-c(114,124,143,158,166)
xi<-c(35,45,55,65,75)
summary(lm(formula=yi~xi) )

Call:
lm(formula = yi ~ xi)
Residuals:
   1    2    3    4    5
 0.6 -3.2  2.0  3.2 -2.6

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  65.1000     5.8284   11.17 0.001538 **
xi            1.3800     0.1026   13.45 0.000889 ***
---
Signif. codes: 
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.246 on 3 degrees of freedom
Multiple R-squared:  0.9837,    Adjusted R-squared:  0.9782
F-statistic: 180.8 on 1 and 3 DF,  p-value: 0.0008894

 

散布図と回帰直線
plot(xi,yi)
abline(lm(yi~xi))

f:id:yoshida931:20170613173429p:plain

 

参考) 豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012

感度、特異度、ROC曲線

感度と特異度

投稿日2017.7.31

f:id:yoshida931:20170731180203p:plain


感 度(陽性反応的中度)= b / ( a + b ) 疾患に罹患していて、検査も陽性  
特異度(陰性反応的中度)= c / ( c + d ) 疾患に罹患してなくて、検査も陰性 

 

偽陽性 a / ( a + b ) 疾患に罹患していて検査が陰性
偽陰性  d / ( c + d ) 疾患に罹患していないのに検査が陽性

感度=1-偽陽性
特異度=1-偽陰性

 

Rを使用してROC曲線を描いてみます
install.packages("ROCR")
library(ROCR)
サンプルデータ
f:id:yoshida931:20170731180401p:plain
    出典)柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010,p12

準備…以下のようにデータを修正します
xは検査点数の個数をベクトルにしたものです
x<-c ( rep( 0,10 ) , rep( 1,2 ) , rep( 2,17 ) , rep( 3,19 ) , rep( 4,12 ) )
yは罹患を1、非罹患を0としてベクトルにしたものです
y<-c ( 1 , rep(0,9) , 0 , 0 , rep(1,6) , rep(0,11) , rep(1,11) , rep(0,8) , rep(1,12) )

 

あとはRの関数を使用します

install.packages("ROCR")
library(ROCR)
pd<- prediction(x,y)
roc <- performance (pd, "tpr", "fpr")
    tpr : True positive rate(感度)、fpr : False positive rate(偽陰性率)
plot( roc )

f:id:yoshida931:20170731181614p:plain

aucの算出

auc.tmp <- performance ( pre , "auc" )
auc <- as.numeric ( auc.tmp@y.values  )

 auc = 0.8327778

 

参考ブログ  k-lab
http://www.nemotos.net/?p=836

参考・・・というか、ほとんどそのまま使用しております.
ありがとうございました.