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

統計学備忘録 since2016

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

統計学的に独立ということ

以前の独立という記事を削除して再投稿

二つの事象A,Bが独立であるということ
一方の事象が起こるかどうかが、他方の事象の起こる確率に影響しないこと
2つの事象AとBが P(A\cap B)= P(A) \times P(B)を満たす場合、事象AとBは統計的に独立である.

このときAとBの余事象も独立である

また

Aの余事象とBの余事象も独立である

条件付き確率
事象A,Bが独立の場合、事象Aと事象Bは互いに影響されないので、
P(A \mid B)= \frac{P(A\cap B)}{P(B)}=P(A)
P(B \mid A)= \frac{P(A\cap B)}{P(A)}=P(B)
という両式が成立します

例題
2個の大小サイコロを振った場合に、P(A)とP(B)は統計学的に独立か?
事象A(サイコロの和が6)、事象B(大きいサイコロの目が4)の場合

<- c(1, 2, 3, 4, 5)<- c(5, 4, 3, 2, 1)
xy <- data.frame(,)

f:id:yoshida931:20171006193628j:plain
P(A \cap B)=\frac{1}{36}
P(A)=\frac{5}{36}
P(B)=\frac{1}{6}
P(A) \times P(B)=\frac{5}{216}
P(A \cap B) \neq P(A) \times P(B)
したがって事象Aと事象Bは独立ではない

事象A(サイコロの和が7)、事象B(大きいサイコロの目が4)の場合

<- c(1, 2, 3, 4, 5, 6)<- c(6, 5, 4, 3, 2, 1)
xy <- data.frame(,)

f:id:yoshida931:20171006194032j:plain
P(A \cap B)=\frac{1}{36}
P(A)=\frac{6}{36}
P(B)=\frac{1}{6}
P(A) \times P(B)=\frac{6}{216}=\frac{1}{36}
P(A \cap B) = P(A) \times P(B)
したがって事象Aと事象Bは独立である

ん~なんか不思議な気分ですが…
「統計的独立」当然のことのようにテキストに記載されていることでありますが、
実に意義深く重要なことのように思えます…
もうちっと掘り下げねば…

柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010,p235

二項分布からの最尤推定

二項分布からの最尤推定
投稿日2017.6.13
更新日2017.10.4 理解してなかったのでもう一度勉強します

二項分布の復習
確率1/2が最も尤もらしい例
右側の数値は同時確率
f:id:yoshida931:20171004180715j:plain:w250
どう見てもP=0.5が最も尤もらしい確率です!

最尤原理

  • 最尤原理は、現在起きている事象の起きる確率が最も大きくなる値を、その母数の推定値として「最も尤もらしい」とする考え方です.

\ 例えば、表の確率P,ウラの確率1-Pとした場合には、表表ウラの同時確率はL(P)=P×P×(1-P)となります.L(P)を尤度関数といいます.尤度関数は尤もらしさを表す関数と表現できます.

尤度関数
> L(P)尤度関数は、Pのいろいろな値における尤もらしさを表す関数です. 尤度関数を最大にするPの推定値が最尤推定値、関数としては最尤推定量になります.

例として表が出る確率をP=?の場合を考えていきます.
「表、表、ウラ」の事象が出た場合の最も尤もらしいPの値を推定してみます.直観的には2/3になるようですが・・・

この事情が起きる確率は
L(P) =P×P× (1-P) =P^2×(1-P)
最も尤もらしい確率Pを探し出すために、P=0.2, 0.5, 0.8 を当てはめて比べてみます.
L(0.2)=0.009L(0.5)=0.125L(0.8)=0.128

最尤原理から考えると、P=0.8が尤もらしい値となります.

正確に求めるためには微分が必要になります.L(P)を微分して0となるPの値がL(P)を最大にします.
L(P) =P×P× (1-P) =P^2×(1-P)微分して=0としたときのPが最尤推定値となります.
 dL(P)/dp = (p^2 – p^3)’ = 2p-3p^2 = 0
最尤推定P=2/3となり直観は正解でした.

ところが…L(P) =P^{12}×(1-P)^{32-12} このように試行回数32回、成功回数12回というよううな場合には非常に面倒になってしいます.そこで対数尤度 ( log likelihood function )が必要になります.

対数尤度
尤度関数を考えるとき、積の形なので数学的に扱うのが不便であるため対数をとり和の形にします.
L(P) =P^x×(1-P)^{n-x}
 logL(P)=xlogP+(n-x)log(1-P)
これを微分した解を0としてPの値を算出します.

 logL(P)の微分 =\frac{x}{P}-\frac{n-x}{1-P}

例題を考えた場合に
logL(P) =log(P^{12}×(1-P)^{32-12})=12*logP -20*log(1-P))
logL(P)の微分 =\frac{12}{P}-\frac{20}{1-P}= 0

 P =\frac{12}{32}=0.375が最も尤もらしいPの値となります.

#pに0.370, 0.375, 0.380を当てはめて確認
> p = 0.370
> p^12*(1-p)^20
[1] 6.386041e-10
> p = 0.375
> p^12*(1-p)^20
[1] 6.396988e-10
> p = 0.380
> p^12*(1-p)^20
[1] 6.386118e-10

参考文献
柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010
東京大学教養学部統計学教室 (編集); 統計学入門 (基礎統計学), 東京大学出版会, 1991

配列 (array) の基本

配列 (array) の基本
投稿日2017.10.4

リスト:異なるデータ構造をまとめて一つのオブジェクトとして保存.解析結果の多くはリスト形式.
配列:次元を自由に設定できるオブジェクト.

RのサンプルHairEyeColorを通して勉強していきます

> HairEyeColor
, , Sex = Male

       Eye
Hair    Brown Blue Hazel Green
  Black    32   11    10     3
  Brown    53   50    25    15
  Red      10   10     7     7
  Blond     3   30     5     8

, , Sex = Female

       Eye
Hair    Brown Blue Hazel Green
  Black    36    9     5     2
  Brown    66   34    29    14
  Red      16    7     7     7
  Blond     4   64     5     8

女性だけ取り出してみます

> HairEyeColor[,,2]
       Eye
Hair    Brown Blue Hazel Green
  Black    36    9     5     2
  Brown    66   34    29    14
  Red      16    7     7     7
  Blond     4   64     5     8

金髪の目の色のデータ

> HairEyeColor[4,,]
       Sex
Eye     Male Female
  Brown    3      4
  Blue    30     64
  Hazel    5      5
  Green    8      8

青い目をした男性の髪の色データ

> HairEyeColor[,2,1]
Black Brown   Red Blond 
   11    50    10    30 

HairEyeColorはどんあデータになっているのか?

> HairEyeColor[1:32]
32 53 10  3 11 50 10 30 10 25  7  5  3 15  7  8 36 66 16  4  9 34  7 64  5 29  7  5  2 14 7  8

それではこのベクトルからHairEyeColorを作ってみましょう

> dat <- c(HairEyeColor[1:32])   #ベクトルdatを作成
> dat
 [1] 32 53 10  3 11 50 10 30 10 25  7  5  3 15  7  8 36 66 16  4  9 34  7 64  5 29  7  5  2 14
[31]  7  8
> class(dat)
[1] "numeric"

arrayを使用して多次元のクロス表を作成します(意外と簡単!)

> dat <- array(dat, dim = c(4,4,2))
> dat
, , 1

     [,1] [,2] [,3] [,4]
[1,]   32   11   10    3
[2,]   53   50   25   15
[3,]   10   10    7    7
[4,]    3   30    5    8

, , 2

     [,1] [,2] [,3] [,4]
[1,]   36    9    5    2
[2,]   66   34   29   14
[3,]   16    7    7    7
[4,]    4   64    5    8

各項目の名前を付けます(日本語で・・・)

> name <- list(Hair=c("黒","茶","赤","金"),Eye=c("茶","青","榛","緑"),Sex=c("男性","女性"))
> dimnames(dat)<-name
> dat
, , Sex = 男性

    Eye
Hair 茶 青 榛 緑
  黒 32 11 10  353 50 25 1510 10  7  73 30  5  8

, , Sex = 女性

    Eye
Hair 茶 青 榛 緑
  黒 36  9  5  266 34 29 1416  7  7  74 64  5  8

HairEyeColorとdatの構成要素を見てみます

> str(HairEyeColor)
 table [1:4, 1:4, 1:2] 32 53 10 3 11 50 10 30 10 25 ...
 - attr(*, "dimnames")=List of 3
  ..$ Hair: chr [1:4] "Black" "Brown" "Red" "Blond"
  ..$ Eye : chr [1:4] "Brown" "Blue" "Hazel" "Green"
  ..$ Sex : chr [1:2] "Male" "Female"

> str(dat)
 num [1:4, 1:4, 1:2] 32 53 10 3 11 50 10 30 10 25 ...
 - attr(*, "dimnames")=List of 3
  ..$ Hair: chr [1:4] "黒" "茶" "赤" "金"
  ..$ Eye : chr [1:4] "茶" "青" "榛" "緑"
  ..$ Sex : chr [1:2] "男性" "女性"

多次元データdatを2次元の表にまとめてみます

> ftable(dat)
         Sex 男性 女性
Hair Eye              
黒   茶        32   3611    910    53    2
茶   茶        53   6650   3425   2915   14
赤   茶        10   1610    77    77    7
金   茶         3    430   645    58    8

石田 基広 ;改訂3版 R言語逆引きハンドブック ,シーアンドアール研究所; 改訂3版,2016

Rのサンプル ~ データフレーム編

Rのサンプル ~ データフレーム編
投稿日2017.10.2
更新日2017.10.4 少し見やすくしました

こういうのが統計学の学習に役立つかも?
サンプルの詳細を知りたいときはここ

間瀬茂 R 基本統計関数マニュアル
https://cran.r-project.org/doc/contrib/manuals-jp/Mase-Rstatman.pdf

head(BOD)

  Time demand
1    1    8.3
2    2   10.3
3    3   19.0
4    4   16.0
5    5   15.6
6    7   19.8

head(ChickWeight)

  weight Time Chick Diet
1     42    0     1    1
2     51    2     1    1
3     59    4     1    1
4     64    6     1    1
5     76    8     1    1
6     93   10     1    1

head(CO2)

  Plant   Type  Treatment conc uptake
1   Qn1 Quebec nonchilled   95   16.0
2   Qn1 Quebec nonchilled  175   30.4
3   Qn1 Quebec nonchilled  250   34.8
4   Qn1 Quebec nonchilled  350   37.2
5   Qn1 Quebec nonchilled  500   35.3
6   Qn1 Quebec nonchilled  675   39.2

head(DNase)

  Run       conc density
1   1 0.04882812   0.017
2   1 0.04882812   0.018
3   1 0.19531250   0.121
4   1 0.19531250   0.124
5   1 0.39062500   0.206
6   1 0.39062500   0.215

head(Formaldehyde)

  carb optden
1  0.1  0.086
2  0.3  0.269
3  0.5  0.446
4  0.6  0.538
5  0.7  0.626
6  0.9  0.782

head(Indometh)

  Subject time conc
1       1 0.25 1.50
2       1 0.50 0.94
3       1 0.75 0.78
4       1 1.00 0.48
5       1 1.25 0.37
6       1 2.00 0.19

head(InsectSprays)

  count spray
1    10     A
2     7     A
3    20     A
4    14     A
5    14     A
6    12     A

head(LifeCycleSavings)

             sr pop15 pop75     dpi ddpi
Australia 11.43 29.35  2.87 2329.68 2.87
Austria   12.07 23.32  4.41 1507.99 3.93
Belgium   13.17 23.80  4.43 2108.47 3.82
Bolivia    5.75 41.89  1.67  189.13 0.22
Brazil    12.88 42.19  0.83  728.47 4.56
Canada     8.79 31.72  2.85 2982.88 2.43

head(Loblolly)

   height age Seed
1    4.51   3  301
15  10.89   5  301
29  28.72  10  301
43  41.74  15  301
57  52.70  20  301
71  60.92  25  301

head(Orange)

  Tree  age circumference
1    1  118            30
2    1  484            58
3    1  664            87
4    1 1004           115
5    1 1231           120
6    1 1372           142

head(OrchardSprays)

  decrease rowpos colpos treatment
1       57      1      1         D
2       95      2      1         E
3        8      3      1         B
4       69      4      1         H
5       92      5      1         G
6       90      6      1         F

head(PlantGrowth)

  weight group
1   4.17  ctrl
2   5.58  ctrl
3   5.18  ctrl
4   6.11  ctrl
5   4.50  ctrl
6   4.61  ctrl

head(Puromycin)

  conc rate   state
1 0.02   76 treated
2 0.02   47 treated
3 0.06   97 treated
4 0.06  107 treated
5 0.11  123 treated
6 0.11  139 treated

head(Theoph)

  Subject   Wt Dose Time  conc
1       1 79.6 4.02 0.00  0.74
2       1 79.6 4.02 0.25  2.84
3       1 79.6 4.02 0.57  6.57
4       1 79.6 4.02 1.12 10.50
5       1 79.6 4.02 2.02  9.66
6       1 79.6 4.02 3.82  8.58

head(ToothGrowth)

   len supp dose
1  4.2   VC  0.5
2 11.5   VC  0.5
3  7.3   VC  0.5
4  5.8   VC  0.5
5  6.4   VC  0.5
6 10.0   VC  0.5

head(USArrests)

           Murder Assault UrbanPop Rape
Alabama      13.2     236       58 21.2
Alaska       10.0     263       48 44.5
Arizona       8.1     294       80 31.0
Arkansas      8.8     190       50 19.5
California    9.0     276       91 40.6
Colorado      7.9     204       78 38.7

head(USJudgeRatings)

               CONT INTG DMNR DILG CFMG DECI PREP FAMI ORAL WRIT PHYS RTEN
AARONSON,L.H.   5.7  7.9  7.7  7.3  7.1  7.4  7.1  7.1  7.1  7.0  8.3  7.8
ALEXANDER,J.M.  6.8  8.9  8.8  8.5  7.8  8.1  8.0  8.0  7.8  7.9  8.5  8.7
ARMENTANO,A.J.  7.2  8.1  7.8  7.8  7.5  7.6  7.5  7.5  7.3  7.4  7.9  7.8
BERDON,R.I.     6.8  8.8  8.5  8.8  8.3  8.5  8.7  8.7  8.4  8.5  8.8  8.7
BRACKEN,J.J.    7.3  6.4  4.3  6.5  6.0  6.2  5.7  5.7  5.1  5.3  5.5  4.8
BURNS,E.B.      6.2  8.8  8.7  8.5  7.9  8.0  8.1  8.0  8.0  8.0  8.6  8.6

head(airquality)

  Ozone Solar.R Wind Temp Month Day
1    41     190  7.4   67     5   1
2    36     118  8.0   72     5   2
3    12     149 12.6   74     5   3
4    18     313 11.5   62     5   4
5    NA      NA 14.3   56     5   5
6    28      NA 14.9   66     5   6

head(anscombe)

  x1 x2 x3 x4   y1   y2    y3   y4
1 10 10 10  8 8.04 9.14  7.46 6.58
2  8  8  8  8 6.95 8.14  6.77 5.76
3 13 13 13  8 7.58 8.74 12.74 7.71
4  9  9  9  8 8.81 8.77  7.11 8.84
5 11 11 11  8 8.33 9.26  7.81 8.47
6 14 14 14  8 9.96 8.10  8.84 7.04

head(attenu)

  event mag station dist accel
1     1 7.0     117   12 0.359
2     2 7.4    1083  148 0.014
3     2 7.4    1095   42 0.196
4     2 7.4     283   85 0.135
5     2 7.4     135  107 0.062
6     2 7.4     475  109 0.054

head(attitude)

  rating complaints privileges learning raises critical advance
1     43         51         30       39     61       92      45
2     63         64         51       54     63       73      47
3     71         70         68       69     76       86      48
4     61         63         45       47     54       84      35
5     81         78         56       66     71       83      47
6     43         55         49       44     54       49      34

head(beaver1)

  day time  temp activ
1 346  840 36.33     0
2 346  850 36.34     0
3 346  900 36.35     0
4 346  910 36.42     0
5 346  920 36.55     0
6 346  930 36.69     0

head(beaver2)

  day time  temp activ
1 307  930 36.58     0
2 307  940 36.73     0
3 307  950 36.93     0
4 307 1000 37.15     0
5 307 1010 37.23     0
6 307 1020 37.24     0

head(cars)

  speed dist
1     4    2
2     4   10
3     7    4
4     7   22
5     8   16
6     9   10

head(chickwts)

  weight      feed
1    179 horsebean
2    160 horsebean
3    136 horsebean
4    227 horsebean
5    217 horsebean
6    168 horsebean

head(esoph)

  agegp     alcgp    tobgp ncases ncontrols
1 25-34 0-39g/day 0-9g/day      0        40
2 25-34 0-39g/day    10-19      0        10
3 25-34 0-39g/day    20-29      0         6
4 25-34 0-39g/day      30+      0         5
5 25-34     40-79 0-9g/day      0        27
6 25-34     40-79    10-19      0         7

head(faithful)

  eruptions waiting
1     3.600      79
2     1.800      54
3     3.333      74
4     2.283      62
5     4.533      85
6     2.883      55

head(freeny)

              y lag.quarterly.revenue price.index income.level market.potential
1962.25 8.79236               8.79636     4.70997      5.82110          12.9699
1962.5  8.79137               8.79236     4.70217      5.82558          12.9733
1962.75 8.81486               8.79137     4.68944      5.83112          12.9774
1963    8.81301               8.81486     4.68558      5.84046          12.9806
1963.25 8.90751               8.81301     4.64019      5.85036          12.9831
1963.5  8.93673               8.90751     4.62553      5.86464          12.9854

head(freeny.x)

     lag quarterly revenue price index income level market potential
[1,]               8.79636     4.70997      5.82110          12.9699
[2,]               8.79236     4.70217      5.82558          12.9733
[3,]               8.79137     4.68944      5.83112          12.9774
[4,]               8.81486     4.68558      5.84046          12.9806
[5,]               8.81301     4.64019      5.85036          12.9831
[6,]               8.90751     4.62553      5.86464          12.9854

head(infert)

  education age parity induced case spontaneous stratum pooled.stratum
1    0-5yrs  26      6       1    1           2       1              3
2    0-5yrs  42      1       1    1           0       2              1
3    0-5yrs  39      6       2    1           0       3              4
4    0-5yrs  34      4       2    1           0       4              2
5   6-11yrs  35      3       1    1           1       5             32
6   6-11yrs  36      4       2    1           1       6             36

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

head(iris3)

[1] 5.1 4.9 4.7 4.6 5.0 5.4

head(longley)

     GNP.deflator     GNP Unemployed Armed.Forces Population Year Employed
1947         83.0 234.289      235.6        159.0    107.608 1947   60.323
1948         88.5 259.426      232.5        145.6    108.632 1948   61.122
1949         88.2 258.054      368.2        161.6    109.773 1949   60.171
1950         89.5 284.599      335.1        165.0    110.929 1950   61.187
1951         96.2 328.975      209.9        309.9    112.075 1951   63.221
1952         98.1 346.999      193.2        359.4    113.270 1952   63.639

head(morley)

    Expt Run Speed
001    1   1   850
002    1   2   740
003    1   3   900
004    1   4  1070
005    1   5   930
006    1   6   850

head(mtcars)

                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

head(pressure)

  temperature pressure
1           0   0.0002
2          20   0.0012
3          40   0.0060
4          60   0.0300
5          80   0.0900
6         100   0.2700

head(quakes)

     lat   long depth mag stations
1 -20.42 181.62   562 4.8       41
2 -20.62 181.03   650 4.2       15
3 -26.00 184.10    42 5.4       43
4 -17.97 181.66   626 4.1       19
5 -20.42 181.96   649 4.0       11
6 -19.68 184.31   195 4.0       12

head(randu)

         x        y        z
1 0.000031 0.000183 0.000824
2 0.044495 0.155732 0.533939
3 0.822440 0.873416 0.838542
4 0.322291 0.648545 0.990648
5 0.393595 0.826873 0.418881
6 0.309097 0.926590 0.777664

head(rock)

  area    peri     shape perm
1 4990 2791.90 0.0903296  6.3
2 7002 3892.60 0.1486220  6.3
3 7558 3930.66 0.1833120  6.3
4 7352 3869.32 0.1170630  6.3
5 7943 3948.54 0.1224170 17.1
6 7979 4010.15 0.1670450 17.1

head(sleep)

  extra group ID
1   0.7     1  1
2  -1.6     1  2
3  -0.2     1  3
4  -1.2     1  4
5  -0.1     1  5
6   3.4     1  6

head(stackloss)

  Air.Flow Water.Temp Acid.Conc. stack.loss
1       80         27         89         42
2       80         27         88         37
3       75         25         90         37
4       62         24         87         28
5       62         22         87         18
6       62         23         87         18

head(swiss)

             Fertility Agriculture Examination Education Catholic Infant.Mortality
Courtelary        80.2        17.0          15        12     9.96             22.2
Delemont          83.1        45.1           6         9    84.84             22.2
Franches-Mnt      92.5        39.7           5         5    93.40             20.2
Moutier           85.8        36.5          12         7    33.77             20.3
Neuveville        76.9        43.5          17        15     5.16             20.6
Porrentruy        76.1        35.3           9         7    90.57             26.6

head(trees)

  Girth Height Volume
1   8.3     70   10.3
2   8.6     65   10.3
3   8.8     63   10.2
4  10.5     72   16.4
5  10.7     81   18.8
6  10.8     83   19.7

head(warpbreaks)

  breaks wool tension
1     26    A       L
2     30    A       L
3     54    A       L
4     25    A       L
5     70    A       L
6     52    A       L

head(women)

  height weight
1     58    115
2     59    117
3     60    120
4     61    123
5     62    126
6     63    129

確率点と累積分布

確率点と累積分布 

投稿日2017.7.22

更新日2017.9.28

 

lower.tail
logical; if TRUE (default), probabilities are P[X ≤ x] otherwise, P[X > x].
                         Rのhelpより

 

TRUE=default

 

正規分布

上側確率が0.025になるZ値
qnorm ( 0.975, lower.tail=TRUE )        

= 1.959964  ( Z=1.959964 より下側の累積確率が0.975 )

qnorm ( 0.975 ,lower.tail=FALSE )  

 =  - 1.959964  ( Z= -1.959964 より上側の累積確率が0.975 )
上側確率が0.025になる累積分

pnorm ( 1.96, lower.tail=TRUE )        

=  0.9750021

下側確率が0.025になる累積分

pnorm ( 1.96, lower.tail=FALSE )   

=  0.0249979

 

t分布

例)自由度10

上側確率が0.025になるt値

qt ( 0.975 ,10 ,lower.tail=TRUE )      

= 2.228139

下側確率が0.025になるt値

qt ( 0.025 ,10 )  

=  qt ( 0.975 ,10 ,lower.tail=FALSE )                      

=  - 2.228139

 

F分布

自由度 ( 6 , 24 ) の確率0.05のF値

f<- qf ( 0.95 , 6 , 24, lower.tail=TRUE )     

f = 2.508189

自由度 ( 6 , 24 ) 、F値がfの場合の累積確率

pf ( f , 6 , 24 , lower.tail=TRUE )                

 = 0.95

その場合のp値

pf ( f , 6 , 24 , lower.tail=FALSE )   

 = 1 - pf ( f , 6 , 24 ) 

 = 0.05

 

二項分布

f:id:yoshida931:20170613163424p:plain

dbinom(x, size, prob)           # 確率密度 (山の高さ)
pbinom(q, size, prob, lower.tail = TRUE)   # 累積分
qbinom(p, size, prob, lower.tail = TRUE)   #確率点

例)コインの裏表

f:id:yoshida931:20170613163444p:plain

dbinom(3,3,0.5)=1/8
dbinom(2,3,0.5)=3/8
dbinom(1,3,0.5)=3/8
dbinom(0,3,0.5)=1/8

Z検定

Z検定
投稿日2017.9.19

記号の定義

母平均μx, 母分散σx, 標本x, サイズm, 標本平均xa
母平均μy, 母分散σy, 標本y, サイズn, 標本平均ya

平均値のZ検定

2標本平均の差の検定 - 統計学備忘録 since2016より
平均値の標準誤差=\sqrt{V(xa)=\frac{σx^2}{m}}\ =\ \frac{サンプルの標準偏差}{\sqrt{サンプルサイズ}}
Z=\frac{(xa-μx)}{\sqrt{\frac{σx^2}{m}}}

問題 xa=38,σx=15, n=100, α=0.05 で次の帰無仮説を検定せよ!
帰無仮説H0:μx = 40 対立仮説H1:μx ≠ 40
(出典:村上 正康,安田 正実; 統計学演習, 培風館,1989,p119)
上記のような問題は次のように書き換えて記憶しておくとよいと思います.(現実的ではないのですが)十分大きなサンプル(n=100)の平均が38、分散が15だった場合に、その母集団の平均は40であることを有意水準0.05で検定せよ.
以下Rで計算してみます

z <- (38 - 40) / sqrt(15^2/100)
- 1.959964 = -1.33
qnorm(0.025,lower.tail = T)
= - 1.959964
z>- 1.959964のため帰無仮説は採択される
pnorm(z,lower.tail = T)
= 0.09121122

f:id:yoshida931:20170919125119j:plain:w500
水色の部分の面積が0.09121122ということになります.したがって、有意水準0.05で母集団の平均は40であると言えます.またサンプル数が少ない場合にはt分布を使用した検定になります.ちなみにこの図の描き方は・・・polygonという関数を使うらしいです.http://cse.naro.affrc.go.jp/takezawa/r-tips/r/51.html

plot(dnorm, -4, 4, xaxt="n")
xvals <- seq(-4, -1.33, length=10)        # -4以上-1.33以下 領域をx軸方向に10個の多角形(台形)に等分割
dvals <- dnorm(xvals)             # 対応するグラフの高さ
polygon(c(xvals,rev(xvals)),
        c(rep(0,10),rev(dvals)),col=5)            # 塗りつぶす
name<-c("-1.96","-1.33","0")
axis(side=1,at=c(-1.96,-1.33,0),labels=name)         #指定した場所にnameを挿入

平均値の差のZ検定

平均値の差の標準誤差
V(xa-xy)=V(xa)+V(xy)   分散の加法性
V(xa)=\frac{σx^2}{m}, V(ya)=\frac{σy^2}{n}

平均値の標準誤差 = {\sqrt{\frac{σx^2}{m}+\frac{σy^2}{m}}}

Z=\frac{(xa-ya)-(μx-μy)}{\sqrt{\frac{σx^2}{m}+\frac{σy^2}{m}}}

問題 xa=22.31, xy=21.54, σx=3.8, σx=3.2, m=50, n=40, α=0.05 で次の帰無仮説を検定せよ.また2群の差の95%信頼区間を求めてよ.
帰無仮説H0:μx = μy 対立仮説H1:μx ≠ μy
(出典:村上 正康,安田 正実; 統計学演習, 培風館,1989,p119)
上記のような問題は次のように書き換えて記憶しておくとよいと思います.
m=50、平均22.31、分散3.8 と n=40、平均21.54、分散3.2 の2群には差がないことを有意水準0.05で検定せよ.

z <- (22.31-21.54)/sqrt(3.8^2/50 + 3.2^2/40)
z = 1.043211
qnorm(0.025,lower.tail = F)
= 1.959964
  #  z<1.959964のため、帰無仮説は採択され、2群には差がないという検定結果になります
2群の差の95%信頼区間は
(22.31-21.54) - 1.959964 * sqrt(3.8^2/50 + 3.2^2/40)
=-0.6766606
(22.31-21.54) + 1.959964 * sqrt(3.8^2/50 + 3.2^2/40) 
= 2.216661
  # -0.6766606 < 2群の差 < 2.216661 より2群には差がないことが推測できます.

比率のZ検定

確率変数Xが二項分布B(n,p)に従う場合、期待値E[X]=np, 分散V[X]=np(1-p)となります.nがある程度大きい場合に、比率pの推定値としてp'=x/n を用います. n=5の場合yesが3回でたらP(YES)=3/5.
p'の期待値はE(p')=p
p'の分散はV(p')=V(x/n)=\frac{V(x)}{n^2}=\frac{np(1-p)}{n^2}=\frac{p(1-p)}{n}

比率の標準誤差={\sqrt{\frac{p(1-p)}{n}}}

z=\frac{p'-p}{\sqrt{\frac{p(1-p)}{n}}}

問題  ある月に1020人を調査したところ内閣支持率は57%でした.内閣支持率の95%信頼区間を推定せよ.

z <- (0.57 - p)/sqrt(0.57*0.43/1020) # ~N(0,1)
pmi <- 0.57-z*sqrt(0.57*0.43/1020) # 下限
pma <- 0.57+z*sqrt(0.57*0.43/1020) # 上限
z <- qnorm(0.025,lower.tail = F) # p値0.025のときのZ値
0.5396178 < p < 0.6003822

比率の差のZ検定

比率の標準誤差={\sqrt{\frac{p(1-p)}{n}}} より比率の分散は\frac{p(1-p)}{n}
例としてYes, Noの割合を考えます.サイズmのA群のyesの比率をp1, サイズnのB群のyesの比率をp2とします.p1-p2の期待値はE(p1-p2), p1-p2の分散はV(p1-p2)=V(p1)+V(p2)・・・分散の加法性より.

比率の差の標準誤差=\sqrt{{\frac{p1*(1-p1)}{m}+\frac{p2*(1-p2)}{n}}}

A群、B群の母集団の比率をP1、P2とした場合

z = \frac{(p1-p2)-(P1-P2)}{\sqrt{\frac{p1*(1-p1)}{m}+\frac{p2*(1-p2)}{n}}}

どんなときに使用するか?まずはクロス表から考えてみます.
問題 男性530人、女性580人にある質問した場合に、yesと答えた割合は男性0.24、女性0.30だった.慌て者は迷わず女性が多いと考えると思います.この回答結果に有意差(α=0.05)はあるのでしょうか?帰無仮説H0:P1=P2

z <- (0.25-0.30)/sqrt(0.25*0.75/530+0.30*0.70/580)
z = -1.868793
qnorm(0.025, lower.tail = T)
= -1.959964
z > -1.959964
-0.1024393 < P1-P2 < 0.002439296
したがって帰無仮説は採択され、男女の割合に差は認められない.

仮説検定のP値や信頼区間は「自分があわて者である」かどうかを教えてくれる.それをどう活かすかというところにこそ、あなたの経験と勘を活かせばいいのである.

西内 啓; 統計学が最強の学問である[実践編]---データ分析のための思想と方法, ダイヤモンド社, 2014

独立性の検定

独立性の検定
投稿日2016.11.4
更新日2017.9.13

f:id:yoshida931:20170913173322j:plain:w500
  
この分割表において独立とは  Ai ∩ Bj の各確率に対して
全ての i j(どの i j でも)に対して  P(Ai ∩ Bj )=P( Ai ) * P( Bj )  となる

AとBが完全に独立な場合の例
  A1 : A2 : A3  がどのBでも  1 : 3 : 6 
  B1 : B2 : B3  がどのAでも  1 : 3 : 5
f:id:yoshida931:20170913173555j:plain:w250

帰無仮説の表現いろいろ
事象Ai とBjは独立である 
A1 , A2 , ... , Ak の起こり方はどのBjに対しても共通である
全ての i jに対して  P(Ai ∩ Bj )=P( Ai ) * P( Bj )

期待度数(理論度数)
Ai の出現確率の推定値は  P( Ai ) = 行の周辺度数 / 総度数 = ni・ / n
Bj の出現確率の推定値は  P( Bj ) = 列の周辺度数 / 総度数 = nj・ / n
したがって セルOij の期待度数Eiは、総度数×セルの確率
Ei  =  n * P( Ai ∩ Bj )  = n * P( Ai ) * P( Bj )  =  n * ( ni・ / n ) * ( nj・ / n )

カイ2乗値
=  ∑∑ ( ( ( 観測度数 - 期待度数 )^2 ) / 期待度数 )

2つの異なる属性A,Pについて度数を集計する.
rn <- c ( "P1" , "P2" , "P3" )
A1 <- c ( 12, 10, 13 )
A2 <- c ( 20, 12, 16 )
A3 <- c ( 20, 15, 19 )
y <-  data.frame (A1,A2,A3)
rownames ( y ) <- paste0 ( rn )
f:id:yoshida931:20170913194511j:plain:w200
周辺度数を算出します

ym <- matrix(c(A1,A2,A3),3,3)
addmargins(ym)

f:id:yoshida931:20170913195030j:plain
次は期待値を求めます

E1A <- (52*35)/137
E1B <- (52*48)/137
E1C <- (52*54)/137
E2A <- (37*35)/137
E2B <- (37*48)/137
E2C <- (37*54)/137
E3A <- (48*35)/137
E3B <- (48*48)/137
E3C <- (48*54)/137

カイ2乗値を求めてみます  自由度 = ( 3-1 ) * ( 3-1 ) = 4

( 12 - E1A )^2 / E1A +
  ( 20 - E1B )^2 / E1B +
  ( 20 - E1C )^2 / E1C +
  ( 10 - E2A )^2 / E2A +
  ( 12 - E2B )^2 / E2B +
  ( 15 - E2C )^2 / E2C +
  ( 13 - E3A )^2 / E3A +
  ( 16 - E3B )^2 / E3B +
  ( 19 - E3C )^2 / E3C
=0.5099473

P値

pchisq(0.5099473, 4, lower.tail = F)      # 累積分布 pchisq ( カイ二乗値 , 自由度 ) 
=0.9725254

帰無仮説は棄却できないため、P1~3におけるAの比率に明確な違いがあるとは言えません.

Rの関数で確認します

chisq.test ( y , correct=F )

f:id:yoshida931:20170913200859j:plain

Yates(イェ-ツ)補正
chisq.test(x)
2行×2列のクロス集計表のデータに対して行われる補正分割表から得られる.分割表から得られるカイ二乗値は跳び跳びの値しかとらないがχ2 分布は連続分布である。 このため,2×2 分割表の場合には連続性の補正をしたほうがよい。検出力は低下するが、より正確な検定が可能になる。

フィッシャーの正確確率検定
fisher.test(x)
分割表において,期待値が 5 以下の桝目が全体の桝目の 20% 以上あるか,期待値が 1 以下の桝目が 1 つでもある場合には,χ2 分布を利用する「独立性の検定」は不適当である。そのような場合に実施する検定。

http://aoki2.si.gunma-u.ac.jp/lecture/Cross/Fisher.html


例)HairEyeColorを使って練習してみましょう!

HairEyeColor
str(HairEyeColor)

table [1:4, 1:4, 1:2] 32 53 10 3 11 50 10 30 10 25 ...
table [行数, 列数, 表数]

y<-addmargins(HairEyeColor) #周辺度数付き
#HairEyeColorの第一表(Sex = Male)にアクセスする
y[,,1]

f:id:yoshida931:20170913184150j:plain
それでは第一表(Sex = Male)のカイ二乗値を求めてみます

addmargins(HairEyeColor) #これで周辺度数を確認します
mh<-HairEyeColor[,,1] #第1表のみ抜き出します
xx<-matrix(c(56,143,34,46),4,4) #行合計マトリックス
yy<-t(matrix(c(98,101,47,33),4,4)) #列合計マトリックス
xy<- xx*yy/279 #各データの期待度数マトリックス
sum((mh-xy)^2/xy) #t値の算出
= 41.28029

行合計マトリックス
f:id:yoshida931:20170913193309j:plain
列合計マトリックス
f:id:yoshida931:20170913193322j:plain
各データの期待度数マトリックス
f:id:yoshida931:20170913193423j:plain
Rで確認してみます

chisq.test ( HairEyeColor[,,1] )

f:id:yoshida931:20170913193701j:plain

残差分析

y <- HairEyeColor[,,1]
ct <- chisq.test ( y )
# 期待値
ct$expected
# 標準化残差
ct$residuals
# 調整済み標準化残差の算出
ct$residuals/sqrt(outer(1-rowSums(y)/sum(y), 1-colSums(y)/sum(y)))