FC2ブログ

2値データのロジスティック回帰分析

2007-11-22-Thu-06:13

RのデータToothGrowthを用いて2値データの被説明変数を予測するモデルを作ってみる。


len 歯の長さ
supp 摂取法
dose 投与量

まず、ToothGrowthから5行をランダムサンプリングする。
data(ToothGrowth)
samp<-sample(60,5)
ToothGrowth[samp,]
    len supp dose
26 32.5   VC    2
48 21.2   OJ    1
25 26.4   VC    2
30 29.5   VC    2
27 26.7   VC    2

次に、glm関数で二項分布を選択し解析する
attach(ToothGrowth)
Tooth.glm<-glm(supp~len+dose,family=binomial)
サマリーは以下の通り
summary(Tooth.glm)
===========================================
Call:
glm(formula = supp ~ len + dose, family = binomial)
Deviance Residuals:
     Min        1Q    Median        3Q       Max 
-1.59843  -0.96252   0.04771   1.04486   1.84253 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   1.5377     0.7860   1.956  0.05044 .
len          -0.2127     0.0728  -2.921  0.00348 **
dose          2.0886     0.8497   2.458  0.01397 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 83.178  on 59  degrees of freedom
Residual deviance: 72.330  on 57  degrees of freedom
AIC: 78.33
Number of Fisher Scoring iterations: 4
===========================================
実測値のサンプルは、
real<-supp[samp]
予測値は、
pred<-fitted(Tooth.glm)
これらを表示すると、
data.frame(real,pred[samp])
   real pred.samp.
26   VC  0.2319393
48   OJ  0.2926265
25   VC  0.5249773
30   VC  0.3637015
27   VC  0.5090457


予測値を四捨五入して0,1とする
pred1<-round(pred)
> data.frame(real,pred1[samp])
   real pred1.samp.
26   VC           0
48   OJ           0
25   VC           1
30   VC           0
27   VC           1

table(supp,pred1)
    pred1
supp  0  1
  OJ 17 13
  VC  7 23
つまり、OJ30個のうち17が、またVC30個のうち23が正しく予測されていることになる。


detach(ToothGrowth)

スポンサーサイト



ロジスティック回帰分析

2007-11-14-Wed-13:26

携帯電話や、インターネットの普及率などは以下のグラフに示すようなS字曲線を示すはずである。
e<-seq(from=-5,to=5,length=200)
plot(e,exp(e)/(1+exp(e)),type="l")
logistic1.jpg
もしこのようなモデルに線形回帰分析を適用すると、両端の誤差が大きくなってしまう。
以下にテレビの普及率のデータを示す。
year<-c(1966:1984)
adop<-c(0.003,0.016,0.054,0.139,0.263,0.423,0.611,
 0.758,0.859,0.903,0.937,0.954,0.978,0.978,
 0.982,0.985,0.989,0.988,0.992)

このデータにロジスティック回帰分析を適用してみる。
tv<-glm(adop~year,family=binomial)
そして、実測値と予測値のグラフを作ってみる。
plot(year,adop,type="l")
lines(year,fitted(tv),lty=2,col="red",lwd=2)
legend(1975,0.5,c("Actual","Predict"),col=1:2,lty=1:2)
logistic2.jpg
実測値と予測値は近似しているようだ。

念のため、線形回帰分析も行い、ロジスティック回帰分析との予測値の違いを比較する。
tvl<-lm(adop~year)
predl<-predict(tvl,type="response")
pred<-predict(tv,type="response")
data.frame(year,adop,pred,predl)

   year  adop       pred     predl
1  1966 0.003 0.02794040 0.1108842
2  1967 0.016 0.05108657 0.1734877
3  1968 0.054 0.09160042 0.2360912
4  1969 0.139 0.15886453 0.2986947
5  1970 0.263 0.26131340 0.3612982
6  1971 0.423 0.39852715 0.4239018
7  1972 0.611 0.55377664 0.4865053
8  1973 0.758 0.69919961 0.5491088
9  1974 0.859 0.81321495 0.6117123
10 1975 0.903 0.89076552 0.6743158
11 1976 0.937 0.93855114 0.7369193
12 1977 0.954 0.96622512 0.7995228
13 1978 0.978 0.98167919 0.8621263
14 1979 0.978 0.99013428 0.9247298
15 1980 0.982 0.99470837 0.9873333
16 1981 0.985 0.99716781 1.0499368
17 1982 0.989 0.99848590 1.1125404
18 1983 0.988 0.99919105 1.1751439
19 1984 0.992 0.99956794 1.2377474


線形回帰分析は、実測値のとの乖離が大きいばかりでなく、普及率が、100%を超えるという現象まで起こっている。ロジスティック関数は、二項分布と深く関係している。

一般化線形モデル

2007-11-12-Mon-13:40

分散分析(ANOVA)、共分散分析(ANCOVA)および線形回帰分析は一般線形モデルであり、残差が正規分布に従う、という仮定に基づいている。

しかし、実際にはデータが常に正規分布に従うということはない。
一般化線形モデル(GLM: Generalized Liner Model)とは、正規分布、二項分布、ポアソン分布、ガンマ分布、逆正規分布といった分布族にデータを対応させ、非線形の現象を線形モデルのように扱うデータ解析方法である。

サンプルデータとしてRに用意されているairqualityを用いる。
 data(airquality)
airq<-airquality[,1:4]

まずは、4変数の対散布図を使ってデータを俯瞰してみる。 
> pairs(airq,panel=panel.smooth,lwd=2)  #panel.smoothで、データの傾向を示す線が描ける。lwdはその線の太さ。
glm2.jpg


オゾン量を被説明変数、日射量、温度、風力を説明変数とした重回帰分析を実施すると
> airq.lm<-lm(Ozone~.,data=airq)
> summary(airq.lm)
=============================================
Call:
lm(formula = Ozone ~ ., data = airq)


Residuals:
    Min      1Q  Median      3Q     Max
-40.485 -14.219  -3.551  10.097  95.619


Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) -64.34208   23.05472  -2.791  0.00623 **
Solar.R       0.05982    0.02319   2.580  0.01124 * 
Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
Temp          1.65209    0.25353   6.516 2.42e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 21.18 on 107 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-Squared: 0.6059,     Adjusted R-squared: 0.5948
F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16
=============================================
残差のQ-Qプロットは、
> qqnorm(resid(airq.lm))
> qqline(resid(airq.lm))
glm1.jpg

残差は、正規分布に従っているとはいえない。

モデルの選択基準として、AIC(Akaike's Information Criterion)がよく使われる。
この値が小さい方が、良いモデルであるとされている。
まずは、重回帰モデルから。
> AIC(airq.lm)
[1] 998.7171

次に、glmの正規分布を用いた場合
> AIC(glm(Ozone~Solar.R+Wind+Temp,data=airq,family=gaussian))
[1] 998.7171
重回帰モデルと変わらない。

glmのガンマ分布では、
> AIC(glm(Ozone~Solar.R+Wind+Temp,data=airq,family=Gamma))
[1] 939.8778
Gamma分布によるモデルの当てはめの方が良いようだ。

線形重回帰分析

2007-11-12-Mon-05:22

説明変数が複数である回帰分析を重回帰分析と言う。
ここでもまずは、線形重回帰分析を行う。
前回の、加硫時間とゴムの延びに、「加硫温度」を加えたデータフレームを作成する。

x1<-c(21.5,21.5,22.0,23.0,23.0,23.5,23.5,24.0,24.0,24.5,25.0,25.0,25.5,25.5,26.5,27.0,27.0,27.5,28.5,29.0)
x2<-c(15,13,10,17,6,12,15,11,13,13,8,6,12,11,8,9,10,8,10,13)
y<-c(60,70,65,38,65,50,34,44,34,34,53,43,25,29,32,25,17,28,17,6)
karyu<-data.frame(x1,x2,y)
3変数に対する相関係数を求めてみる。
round(cor(karyu),4) #相関係数行列を求める。
        x1      x2       y
x1  1.0000 -0.3459 -0.8760
x2 -0.3459  1.0000 -0.0609
y  -0.8760 -0.0609  1.0000
加硫時間とゴムの延びには高い相関がある。
加硫時間と加硫温度には若干の負の相関がある。
加硫温度とゴムの延びには大きな相関が認められない。

対散布図を作成してみると、上の関係が視覚的に分かる。 
> pairs(karyu,pch=21,bg="red",cex=2)
Soukan5.jpg
x1とyの単回帰分析では、
Multiple R-Squared: 0.7675,     Adjusted R-squared: 0.7545
決定係数0.7675だった。


では、これに、x2を加えるとどうなるか。 
lm(被説明変数~.,data=,)の~.は説明変数にすべての変数を使うの意味
> karyu2.lm<-lm(y~.,data=karyu)
> summary(karyu2.lm)
=======================================
Call:
lm(formula = y ~ ., data = karyu)


Residuals:
    Min      1Q  Median      3Q     Max
-6.4771 -4.5475 -0.7448  3.8380  9.2359


Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 266.8911    16.9727   15.72 1.46e-11 ***
x1           -8.1148     0.5899  -13.76 1.21e-10 ***
x2           -2.4353     0.4365   -5.58 3.32e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 5.339 on 17 degrees of freedom
Multiple R-Squared: 0.9179,     Adjusted R-squared: 0.9082
F-statistic: 94.99 on 2 and 17 DF,  p-value: 5.935e-10
========================================
決定係数0.9179と単回帰分析に比べ、格段に当てはめがよくなっている。
そして回帰式は、
y=266.8911-8.1148x1-2.4353x2

回帰診断図を以下に示す。
> par(mfrow=c(2,2),oma=c(1,1,2,1),mar=c(4,4,2,1))
> plot(karyu2.lm)
Soukan4.jpg


残差の最大値が、前回の15から10への減少している。


次に、相互作用モデルを考えてみる。
これは、説明変数間の相関関係(相互作用)をを考慮したモデルである。
Rでは、^2をつけることによって実施される。
> karyu2.lm2<-lm(y~.^2,data=karyu)
> summary(karyu2.lm2)
Call:
lm(formula = y ~ .^2, data = karyu)
============================================
Residuals:
   Min     1Q Median     3Q    Max
-6.681 -4.285 -1.033  3.223  8.960


Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 234.7341    71.2928   3.293  0.00459 **
x1           -6.7924     2.9068  -2.337  0.03278 *
x2            0.4090     6.1318   0.067  0.94765  
x1:x2        -0.1180     0.2537  -0.465  0.64814  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 5.467 on 16 degrees of freedom
Multiple R-Squared: 0.919,      Adjusted R-squared: 0.9038
F-statistic: 60.48 on 3 and 16 DF,  p-value: 5.981e-09
============================================

この場合は、相互作用効果を加えても変化がないようだ。


回帰分析

2007-11-11-Sun-06:39

加硫時間(説明変数:x)によって、生産されるゴムの延び(被説明変数:y)がどのように変化するか、といった分析を回帰分析と言う。
ここでは説明変数がひとつの場合である単回帰分析を行う。なお、回帰分析には、線形回帰分析と、非線形回帰分析があるが、まずは線形回帰分析を取り扱う。


x<-c(21.5, 21.5, 22.0, 23.0, 23.0, 23.5, 23.5, 24.0, 24.0, 24.5, 25.0, 25.0, 25.5, 25.5, 26.5, 27.0, 27.0, 27.5, 28.5, 29.0)
y<-c(60, 70, 65, 38, 65, 50, 34, 44, 34, 34, 53, 43, 25, 29, 32, 25, 17, 28, 17, 6)
Karyu.F<-data.frame(x,y)
線形回帰分析は
lm(被説明変数~説明変数,data)で行う。
Karyu.lm<-lm(y~x,data=Karyu.F) #Karyu.lmに分析結果を格納
summary(Karyu.lm) #summaryで、主な分析結果を出力する
===============================================
Call:
lm(formula = y ~ x, data = Karyu.F)


Residuals:(1)
     Min       1Q   Median       3Q      Max
-13.8681  -6.5611   0.5846   5.8642  15.5965


Coefficients:
            Estimate(2) Std. Error t value Pr(>|t|) (3)  
(Intercept) 211.8125    22.5776   9.382 2.36e-08 ***
x            -6.9764     0.9052  -7.707 4.15e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 8.731 on 18 degrees of freedom
Multiple R-Squared: 0.7675 (4),     Adjusted R-squared: 0.7545
F-statistic:  59.4 on 1 and 18 DF,  p-value: 4.149e-07
===============================================
(1) 残差の最小値、第一四分位数、中央値、第二四分位数、最大値
(2)y=a+bxの回帰係数a=211.8125とb=-6.9764
(3)回帰係数が0である、という仮説検定の統計量。a,bともに仮説は棄却される。しかし、この場合、aに対する検定は無意味である。なぜなら、加硫時間が0の時にはゴムは生産されない。
(4)R2=SR/ST。yの全体の変動を表す平方和STのうち、独立変数xによって説明される部分、つまり回帰による平方和SRの占める割合を寄与率、もしくは決定係数と呼ぶ。
SR=R2/ST (SR=ST-Seなので)
Se=(1-R2)ST
つまり、R2が1に近いほど、残差Seは小さくなり、yの変動の大部分がxによって説明されることになる。

予測値は以下のようにして求められる。
pred<-predict(Karyu.lm)
これを実データ、残差とともにリスティングすると
resi<-residuals(Karyu.lm)
data.frame(Karyu.F,pred,resi)
      x  y     pred        resi
1  21.5 60 61.82080  -1.8207953
2  21.5 70 61.82080   8.1792047
3  22.0 65 58.33262   6.6673831
4  23.0 38 51.35626 -13.3562601
5  23.0 65 51.35626  13.6437399
6  23.5 50 47.86808   2.1319183
7  23.5 34 47.86808 -13.8680817
8  24.0 44 44.37990  -0.3799033
9  24.0 34 44.37990 -10.3799033
10 24.5 34 40.89172  -6.8917249
11 25.0 53 37.40355  15.5964535
12 25.0 43 37.40355   5.5964535
13 25.5 25 33.91537  -8.9153681
14 25.5 29 33.91537  -4.9153681
15 26.5 32 26.93901   5.0609887
16 27.0 25 23.45083   1.5491671
17 27.0 17 23.45083  -6.4508329
18 27.5 28 19.96265   8.0373455
19 28.5 17 12.98630   4.0137023
20 29.0  6  9.49812  -3.4981193


Rでは、残差を視覚的に分析できる。 
> par(mfrow=c(2,2))
> plot(Karyu.lm)
>
soukan.jpg
左上<残差とフィット値のプロット>
残差の全体像を把握できる。データ5,7,11の残差が大きいことが見て取れる。

右上<残差の正規Q-Qプロット>
通常の回帰分析は、残差が標準正規分布に従うと言う仮定の下で行っている。
点が直線上に並んでいるため、正規性には問題がなさそうだ。


左下<残差の平方根プロット>
縦軸は、標準化した残差の絶対値の平方根。
これも残差の変動状況を確認するためのものである。


右下<Cookの距離のプロット>
Cookの距離は、1つのデータを除いた後に求めた回帰式による予測値と、すべてのデータを用いた回帰式による予測値との食い違いを表す。
通常0.5以上だと大きいと言われている。

二元配置分散分析

2007-11-10-Sat-07:02

薬剤の評価をするときに、全体では用量反応があるという結果が出ても、施設によってそれがみられていないということがある。
逆に、全体では用量反応がないように見えても、施設間で用量反応を打ち消しているような場合がある。


このような、交互作用を見るために、二元配置分散分析を行う。














































  0mg 10mg 20mg 40mg 60mg
Site1

24.8
23.9
24.1


25.0
26.6
27.9

47.5
32.5
29.5


39.8
36.7
30.7

38.5
35.7
38.7


Site2  48.8
42.6
48.0
48.5
47.1
45.2

46.3
48.2
51.8


48.5
46.5
46.7
51.3
39.4
49.8
Site3 46.4
37.4
29.4
27.9
29.2
36.7

30.2
31.7
29.2


31.7
27.2
25.5
30.3
29.6
41.7
Site 4 30.0
38.7
29.2

25.0
39.9
29.4


30.3
30.9
28.4
29.9
27.3
28.8
33.6
32.0
44.3


まずは、データフレームを作成する。
a1<-c(24.8, 23.9, 24.1, 25.0, 26.6, 27.9, 47.5, 32.5, 29.5, 39.8, 36.7, 30.7, 38.5, 35.7, 38.7)
a2<-c(48.8, 42.6, 48.0, 48.5, 47.1, 45.2, 46.3, 48.2, 51.8, 48.5, 46.5, 46.7, 51.3, 39.4, 49.8)
a3<-c(46.4, 37.4, 29.4, 27.9, 29.2, 36.7, 30.2, 31.7, 29.2, 31.7, 27.2, 25.5, 30.3, 29.6, 41.7)
a4<-c(30.0, 38.7, 29.2, 25.0, 39.9, 29.4, 30.3, 30.9, 28.4, 29.9, 27.3, 28.8, 33.6, 32.0, 44.3)
LB<-data.frame(
A=factor(c(rep("site1",15),rep("site2",15),rep("site3",15),rep("site4",15))),
B=factor(rep(c(rep("0mg",3),rep("10mg",3),rep("20mg",3),rep("40mg",3),rep("60mg",3)),4))
,y=c(a1,a2,a3,a4))

用量ごとに平均を出すと、
tapply(LB$y,LB$B,mean)
     0mg     10mg     20mg     40mg     60mg
35.27500 34.03333 36.37500 34.94167 38.74167
用量相関があるような、ないような・・・

Box Plotを書く
boxplot(y~B,data=LB,main="Efficacy per dose",sub="Fig 1",xlab="Dose", ylab="Efficacy" )
ANOVA.jpg
これでもなんとも言えない・・・


普通の分散分析を行ってみると
summary(aov(y~A+B,data=LB))
            Df  Sum Sq Mean Sq F value    Pr(>F)   
A            3 2588.46  862.82 29.1356 3.421e-11 ***
B            4  157.09   39.27  1.3261    0.2726   
Residuals   52 1539.93   29.61                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A(施設間)に有意な差が見られ、B(用量間)には有意な差が見られなかった。

交互作用を考慮した分散分析を行ってみる。
summary(aov(y~A*B,data=LB))
            Df  Sum Sq Mean Sq F value    Pr(>F)   
A            3 2588.46  862.82 36.5548 1.535e-11 ***
B            4  157.09   39.27  1.6638   0.17741   
A:B         12  595.79   49.65  2.1035   0.03922 * 
Residuals   40  944.14   23.60                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

すると、A:B間に有意水準0.05で交互作用効果がある、と判断される。

対応のある分散分析 Analysis of variance (ANOVA)

2007-11-08-Thu-12:43

前回使ったシミュレーションデータを対応のあるデータとする。
f1<-trunc(rnorm(10,61,5))
f2<-trunc(rnorm(10,67,5))
f3<-trunc(rnorm(10,69,5))
f4<-trunc(rnorm(10,75,5))


temp2<-data.frame(F=factor(c(rep("f1",10),
rep("f2",10),rep("f3",10),rep("f4",10))),
no=factor(rep(1:10,4)),
y=c(f1,f2,f3,f4))

   F no  y
1  f1  1 56
2  f1  2 63
3  f1  3 68
4  f1  4 60
5  f1  5 60
6  f1  6 65
7  f1  7 57
8  f1  8 67
9  f1  9 55
10 f1 10 56
..................
..................
31 f4  1 86
32 f4  2 70
33 f4  3 74
34 f4  4 76
35 f4  5 72
36 f4  6 73
37 f4  7 65
38 f4  8 75
39 f4  9 70
40 f4 10 79


つまり、10台の機器が4つの条件(因子F)でテストされた結果と考えられる。


summary(aov(y~F+no,temp2))
========================================
            Df Sum Sq Mean Sq F value    Pr(>F)   
F            3 895.70  298.57 11.1221 6.237e-05 ***
no           9 128.60   14.29  0.5323     0.838   
Residuals   27 724.80   26.84                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
========================================
因子Fの統計量Prは0に近く、差があるといえる。
一方、機器の種類noの統計量Prは0.838と大きく、機器の差があるとは考えられない。


つまり、因子Fの平均値には差が認められるが、機器noの間には有意差が認められなかった。

対応のない分散分析 Analysis of variance (ANOVA)

2007-11-06-Tue-13:53

得られた各水準の平均値が、因子の影響により変化していると言えるかどうかに関する分析。
また乱数を発生させて検討してみたい。


f1<-trunc(rnorm(10,61,5))
f2<-trunc(rnorm(10,67,5))
f3<-trunc(rnorm(10,69,5))
f4<-trunc(rnorm(10,75,5))


temp<-data.frame(F=factor(c(rep("f1",10),
rep("f2",10),rep("f3",10),rep("f4",10))),
y=c(f1,f2,f3,f4))
データとしては、
     F   y
1  f1 64
2  f1 52
3  f1 73
..............
..............
..............
38 f4 74
39 f4 80
40 f4 80
といった感じになる。
f1,f2,f3,f4毎のyの平均値に違いがあるか、つまり因子の影響があるかどうかを検討する。


箱ひげ図を描くと
boxplot(y~F,data=temp,main="Temparature",sub="Fig 1",xlab="Factor", ylab="Temp." )

 ANOVA.jpg


通常分散分析は、関数aovを使う
aov(y~F,data=temp)
===========================
Call:
   aov(formula = y ~ F, data = temp)
Terms:
                     F Residuals
Sum of Squares  1154.6    1251.8
Deg. of Freedom      3        36


 


Residual standard error: 5.896798
Estimated effects may be unbalanced
===========================
これで群間の変動SSB(1154.6、自由度3)および、群内の変動SSW(1251.8、自由度36)が得られる。


しかし、検定まで行うには以下のようにするとよい。
summary(aov(y~F,data=temp))
===========================
            Df  Sum Sq Mean Sq F value    Pr(>F)   
F            3 1154.60  384.87  11.068 2.707e-05 ***
Residuals   36 1251.80   34.77                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
===========================


これによって、群間の不偏分散(384.87)と群内の不偏分散(34.77)が得られる。
さらに、F=384.87/34.77=11.068および、このF値に対応するP値(2.707e-05)も返されている。
これによると、P値はほぼ0に等しく、「各群の平均値が等しい」という仮説は棄却される。


 

1標本の母平均の検定と2標本の平均値の差の検定

2007-11-02-Fri-12:41

まずは、サンプルデータを発生させる。
身長の平均170cm、標準偏差5のデータを作成する。
height1<-rnorm(25,170,5)

<1標本の母平均の検定>
仮説H0:height1の平均が165cm
=============================
t.test(height1,mu=165,alternative ="two.sided")


        One Sample t-test


data:  height1
t = 2.991, df = 24, p-value = 0.00634
alternative hypothesis: true mean is not equal to 165
95 percent confidence interval:
 165.9856 170.3738
sample estimates:
mean of x
 168.1797
=============================
見事に仮説は棄却された。


<2標本の平均値の差の検定>
身長の平均165cm、標準偏差5のデータを作成する。
height2<-rnorm(25,165,5)
まずは、2標本の母分散が等しいかどうかの検定を行う。
=============================
var.test(height1,height2)

        F test to compare two variances


data:  height1 and height2
F = 0.8606, num df = 24, denom df = 24, p-value = 0.716
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.3792375 1.9529286
sample estimates:
ratio of variances
          0.860595
=============================
2標本の母分散が等しいという仮説は棄却されない。


そこで、var.equal = TRUEとして、t検定を行う。
t.test(height1,height2,var.equal = TRUE)
=============================
        Two Sample t-test


data:  height1 and height2
t = 2.9498, df = 48, p-value = 0.004903
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.467993 7.753758
sample estimates:
mean of x mean of y
 168.1797  163.5688
=============================
というわけで、両群の身長に差がない、という仮説は棄却された。


ところで、仮説検定におけるp 値は、サンプル数に依存する。サンプル数を増やせば幾らでもp 値を小さくできる。それをやってみよう。
n数が5,50,300,1000の正規分布乱数をそれぞれ平均170cmと169cmで作り、p値がどうなるか見てみる。
n<-c(5,50,300,1000)


p<-1:4
for (i in 1:4){
heighta<-rnorm(n[i],170,5)
heightb<-rnorm(n[i],169,5)
t<-t.test(heighta,heightb,var.equal = TRUE)
p[i]<-t$p.value
}
p
[1] 0.6838470460  0.5130966198  0.0621181320  0.0002909941

n数が増加するとともにp値が小さくなり、n=1000では95%有意になった。

カイ二乗分布

2007-11-01-Thu-05:23


平方和Sを母分散σ2で割った値をχ2としたときの分布のことである。


S=(n-1)s2で、s2の期待値がσ2なので、


S/σ2 =(n-1)s22 ≒ n-1となるはずである。


 Rでシミュレーションする


平均0、分散1の標準正規分布乱数を10個発生させその値を


X<-(1:10)
for (i in 1:10){
a<-rnorm(10,0,1)
X[i]<-sum((a-mean(a))^2)/(1^2)
}
X
 [1]  7.822012  8.902314  5.710823  3.754676  5.954704 13.709953 19.746950
 [8]  9.704040  7.766902  7.407105
 range(X)
[1]  3.754676 19.746950
結構ばらついているようだ。


10000回で試してみる
X<-(1:10000)
for (i in 1:10000){
a<-rnorm(10,0,1)
X[i]<-sum((a-mean(a))^2)/(1^2)
}
mean(X)
[1] 9.00005
var(X)
[1] 18.15874 

やはりχ2の期待値はn-1のようだ。
ちなみに、χ2の分散の期待値は2(n-1)となる。


10000回も試行すると変な値もでるもので、
range(X)
[1]  0.7041244            36.7698620
こんな値まで出ている。
p=5%でのχ2値が
qchisq(0.05, 9)
[1] 3.325113
p=95%でのχ2値が
qchisq(0.95, 9)
[1] 16.91898
なので、それぞれ明らかに有意である。
ちなみに最小値のpは、
> pchisq(min(X), 9)
[1] 0.0001307555
最大値のpは、
> pchisq(max(X), 9)
[1] 0.999971

最後に、χ2分布を自由度2,4,6,8,10とした時の確率密度をグラフにする。
x=seq(0,20,0.1)
curve(dchisq(x,2),from=0,to=20,main="カイ2乗分布",xlab="X^2",ylab="f(X^2)")
curve(dchisq(x,4),add=T)
curve(dchisq(x,6),add=T)
curve(dchisq(x,8),add=T)
curve(dchisq(x,10),add=T)


 X2.jpg


 

HOME