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

統計学備忘録（R言語のメモ）

since2016　ときどきTEXのメモ

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

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

二つの事象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.009、L(0.5)＝0.125、L(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  3
  茶 53 50 25 15
  赤 10 10  7  7
  金  3 30  5  8

, , Sex = 女性

    Eye
Hair 茶 青 榛 緑
  黒 36  9  5  2
  茶 66 34 29 14
  赤 16  7  7  7
  金  4 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   36
     青        11    9
     榛        10    5
     緑         3    2
茶   茶        53   66
     青        50   34
     榛        25   29
     緑        15   14
赤   茶        10   16
     青        10    7
     榛         7    7
     緑         7    7
金   茶         3    4
     青        30   64
     榛         5    5
     緑         8    8

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

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

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　　　#2変数6ケース

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　　#4変数578ケース

head(CO2)　草の耐寒性実験

  Plant   Type  Treatment conc uptake     # 個体識別番号、原産地、処理、CO2濃度、CO21摂取量
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　　#5変数84ケース

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

対応のないｔ検定の検定力分析（事後の分析）

t分布 検定力

検定力(1-β)：帰無仮説が「偽」であるとき（母集団に差があるとき）にサンプルから有意差を得る確率
効果量(effect size)：標準化された平均値差(検定により色々あります...d family, r family)
標本数：研究者が決定しなければならない
α(有意水準)：慣例的に5％、1％で固定されている

検定力分析を実行する場合に、効果量の問題を解決しなければなりません．効果量を表す指標は数多く存在するので、まずどれを使うかが問題となります．今回はt検定の検定力分析に利用されている、次の二つの効果量について勉強していきます．

効果量＝効果の大きさ
実験的操作の効果、変数間の関係性の強さを表す指標(大　中　小)

   Cohoen's d : グループごとの平均値の差を標準化したd family
   r : 変数間の関係の強さを示すr family

効果量の目安
   小  中  大
d 0.2 0.5 0.8
r 0.1 0.3 0.5

まずは記号の定義・・・これ大事
母平均μx, 標本x, サイズm, 標本平均xa, 不偏分散Sx^2
母平均μy, 標本y, サイズn, 標本平均ya, 不偏分散Sy^2

標本効果量(cohensのd効果量)の求め方
d ＝ ｘの平均値とｙの平均値の差の絶対値÷2群の共通の標準偏差
プールした分散
S^2=\frac{\sum(x-xa)^2+\sum(y-ya)^2}{m+n-2}=\frac{(m-1)Sx^2+(n-1)Sy^2}{m+n-2}

例
x <- c(39, 47, 48, 59, 65, 68, 74, 79, 80, 82, 87, 88, 89, 94, 99)
y <- c(25, 27, 32, 37, 38, 44, 45, 52, 54, 60, 65, 67, 68, 69, 71)
pool <- ((length(x) - 1) * var(x) + (length(y) - 1) * var(y)) / (length(x) + length(y) - 2)
t <- (mean(x) - mean(y)) / sqrt((pool / length(x)) + (pool / length(y)))
= 3.6449

d <- (abs(mean(x) - mean(y)))/pool
=1.33 (dは1を超える場合もある)

2標本平均の差の検定 - 統計学備忘録 since2016

Rを使ってt値から算出してみます

t.test(x,t,var.equal = T)
d = t*sqrt((length(x)+length(y))/length(x)*length(y))
  = 3.6449*sqrt(30/(15*15))
  = 1.33

Rを使って検定力を算出します

library(pwr)
pwr.t2n.test(n1=15,n2=15, d =1.33 , sig.level =0.05 , power = NULL,alternative = "greater")

t test power calculation
n1 = 15
n2 = 15
d = 1.33
sig.level = 0.05
power = 0.971769
alternative = greater

power = 0.97より、帰無仮説が偽の場合に有意差が検出される確率は97％です．高い検定力であると言えます．

標本効果量(r)の求め方
r= \sqrt{\frac{t^2}{t^2+自由度}}

t.test(x,y,var.equal = T)
t = 3.6449
r <- sqrt(t^2/(t^2+18))
=0.6516516

効果量は「大」でした

参考）水本篤; 竹内理. 研究論文における効果量の報告のために―基本的概念と注意点―. 2008.