・・・ 式1最終更新日 ： 2007.8.14
モデル、母数、GLM
p47 「BOX 3.1 樹木データセットに対する一般線形モデル」
web上からダウンロードできる"datasheetsch03"のTreesより、HEIGHTとVOLUMEの列のみをRに取り込んで分析する。今度はlm()ではなく、一般線形モデルのための関数glm()を用いる。教科書でいうところの分散分析表や係数表はlm()のときと同様にして得られる。
> dat <- excel.w(2) #Excelからデータの読み込み
Read 62 items
> HEIGHT <- dat[,1] #データを変数へ代入
> VOLUME <- dat[,2]
> result.tre <- glm(VOLUME ~ HEIGHT) #解析結果を変数へ代入
> anova(result.tre, test="F") #分散分析表の作成
Analysis of Deviance Table
Model: gaussian, link: identity
Response: VOLUME
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev F Pr(>F)
NULL 30 8106.1
HEIGHT 1 2901.2 29 5204.9 16.165 0.0003784 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(result.tre) #係数表の作成
Call:
glm(formula = VOLUME ~ HEIGHT)
Deviance Residuals:
Min 1Q Median 3Q Max
-21.274 -9.894 -2.894 12.067 29.852
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -87.1236 29.2731 -2.976 0.005835 **
HEIGHT 1.5433 0.3839 4.021 0.000378 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 179.4791)
Null deviance: 8106.1 on 30 degrees of freedom
Residual deviance: 5204.9 on 29 degrees of freedom
AIC: 252.8
Number of Fisher Scoring iterations: 2
実をいうと、今回の解析は関数lm()を用いた場合と同じである。ただ、関数glm()では応答変数（目的変数）の分布が正規分布（gaussian）、二項分布（binomial）、ポアソン分布（poisson）、逆正規分布（inverse.gaussian）、ガンマ分布（Gamma）、応答変数の分布がハッキリしないときに用いる擬似尤度モデル（quasi）のいずれかを指定できる点で異なる。なお、正規分布の場合は省略しても構わない。
それから、出力結果の下の方にAICという値があるが、これはモデルの当てはまりの良さを評価するためのもので、この値は小さいほどよいモデルであると判断できる。したがって、今回のモデルはあまり良いものとはいえない。
また、以下のようにしてlm()を用いたときと同じ形式の分散分析表を作成することもできる。
> fm <- aov(result.tre)
> summary(fm)
Df Sum Sq Mean Sq F value Pr(>F)
HEIGHT 1 2901.2 2901.2 16.165 0.0003784 ***
Residuals 29 5204.9 179.5
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p48 「BOX 3.2 肥料データセットに対する一般線形モデル」
> dat <- excel.w(2)
Read 60 items
> FERTIL <- as.factor(dat[,1]) #カテゴリカル型として代入
> YIELD <- c(dat[,2])
> result.fer <- glm(YIELD ~ FERTIL)
> summary(aov(result.fer)) #分散分析表の作成
Df Sum Sq Mean Sq F value Pr(>F)
FERTIL 2 10.8227 5.4114 5.7024 0.008594 **
Residuals 27 25.6221 0.9490
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(result.fer)
Call:
glm(formula = YIELD ~ FERTIL)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.3750 -0.6715 -0.0720 0.1740 2.5130
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.4450 0.3081 17.676 2.28e-16 ***
FERTIL2 -1.4460 0.4357 -3.319 0.00259 **
FERTIL3 -0.9580 0.4357 -2.199 0.03663 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.9489685)
Null deviance: 36.445 on 29 degrees of freedom
Residual deviance: 25.622 on 27 degrees of freedom
AIC: 88.404
Number of Fisher Scoring iterations: 2
さて、ここで教科書の解説を少し補足しておこう。p49の上半分に表3.3および表3.4についての解説がされているが、これは要するにどういうことなのかをRで試してみよう。
まず今回得られた回帰式は次のようになった。
・・・ 式1
この式は肥料1（FERTIL1）を基準にどの程度、他の肥料の効果があるかどうかを示しているのである。だから、今回の場合は肥料1を0として、肥料2は-1.4460だけ肥料1よりも効果が低いということになる。同様に肥料3も-0.9580だけ肥料1に比べて効果が低いということがいえる。
ただし、教科書（Minitab）では基準を肥料3として計算しているので最終的に得られた結果がRとは異なる。すなわち、
・・・ 式2
という式が得られているが、この両者は式の係数は異なっていても、結果として同じことを示しているのである。このことを詳しく説明する前に、表3.4についての解説について補足しておこう。
既に分かりきったことではあるが、変数FERTILはカテゴリカル型のデータであり、ソフトウェアの中ではダミー変数として扱っている（計算している）。そのため、FERTIL1もFERTIL2もFERTIL3も0か1の値しかとらず、どれか1つが1である場合、その他の変数は全て0になる。
例えば、式1において、FERTIL2=1の場合、他の変数であるFERTIL3は0となる。
YIELD = 5.4450 - 1.4460 * 1 - 0.9580 * 0 = 3.999
今度はFERTIL3=1の場合を考えてみると、他の変数であるFERTIL2は0であるから、次のようになる。
YIELD = 5.4450 - 1.4460 * 0 - 0.9580 * 1 = 4.487
さらに、FERTIL1=1のときはどうなるかというと、他の変数FERTIL2, 3は0となるから、次のようになる。
YIELD = 5.4450 - 1.4460 * 0 - 0.9580 * 0 = 5.4450
これらの数値は、実は各群の平均値である。
| FERTIL1 = 1のとき | YIELD = 5.4450 （これは肥料1の平均値） |
| FERTIL2 = 1のとき | YIELD = 3.999 （これは肥料2の平均値） |
| FERTIL3 = 1のとき | YIELD = 4.487 （これは肥料3の平均値） |
このような観点で式1と式2を比べてみると、両者が同じものであることが理解できよう。
#式1の場合 > f <- function(FERTIL2, FERTIL3) 5.4450 - 1.4460*FERTIL2 - 0.9580*FERTIL3 > f(0,0) #FERTIL1 = 1, FERTIL2 = 0, FERTIL3 = 0 [1] 5.445 > f(1,0) #FERTIL1 = 0, FERTIL2 = 1, FERTIL3 = 0 [1] 3.999 > f(0,1) #FERTIL1 = 0, FERTIL2 = 0, FERTIL3 = 1 [1] 4.487 #式2の場合 > f2 <- function(FERTIL1, FERTIL2) 4.4870 + 0.9579*FERTIL1 - 0.4880*FERTIL2 > f2(1,0) #FERTIL1 = 1, FER2 = 0, FER3 = 0 [1] 5.4449 > f2(0,1) #FERTIL1 = 0, FER2 = 1, FER3 = 0 [1] 3.999 > f2(0,0) #FERTIL1 = 0, FER2 = 0, FER3 = 1 [1] 4.487
このwebサイト内にもこれに関連した内容のページがあるので参照してもらいたい。
[1] 数量化Ⅰ類はダミー変数を導入した回帰分析と同じであるという解説ページ
[2] 分散分析は回帰分析として分析できるという解説ページ
[3] t検定は回帰分析として分析できるという解説ページ
p50 「BOX 3.3 データセットを生成する」
ここではRによる実行例を示す。手順としては、まず
Y = K3 + K1 * DUM1 + K2 * DUM2 + K4 * error
を関数Y(DUM1, DUM2)として定義する。
Y <- function(DUM1, DUM2) 12.2 + 5.0*DUM1 + -2.5*DUM2 + 1.5*rnorm(1)
つまり、関数Yの引数がDUM1=1, DUM2=0と指定された場合はgroup1のデータが1つ生成される。DUM1=0, DUM2 =1とすればgroup2のデータが1つ生成され、両方とも0とすればgroup3のデータが1つ生成されるというわけである。
次に任意の数のデータを生成させるための関数DS(group1, group2, group3)を定義する。
DS <- function(g1, g2, g3){
group1.dat <- c()
for(i in 1:g1) group1.dat[i] <- Y(1, 0)
group2.dat <- c()
for(i in 1:g2) group2.dat[i] <- Y(0, 1)
group3.dat <- c()
for(i in 1:g3) group3.dat[i] <- Y(0, 0)
dat <- cbind(group1.dat, group2.dat, group3.dat)
dat
}
例えば、以下のように引数を指定した場合、各グループ10個のデータが生成される。
> myData <- DS(10,10,10)
> myData
group1.dat group2.dat group3.dat
[1,] 18.24881 11.216955 14.194820
[2,] 15.10206 9.934522 12.376224
[3,] 15.73394 11.153262 11.102854
[4,] 15.95902 10.090320 13.944442
[5,] 13.10929 9.508876 9.806116
[6,] 15.93356 11.427215 11.782779
[7,] 16.73333 7.446832 13.716463
[8,] 16.27401 13.212243 10.490076
[9,] 17.46408 8.001784 12.735103
[10,] 18.28915 9.942723 11.855761
このようにデータを作成することができたら、あとは今までどおり関数glm()を使用して解析すればよい。これによって得られる分散分析表と係数表の数値が教科書と異なるのは、前述したように、Rでは基準となる変数（ダミー変数）をgroup1としているためである。
ちなみに、これらのデータは擬似乱数によって作成されたものなので、1回1回作成されるデータが異なるのは当然であり、このページの通りにRでデータを作成しても、このページで示しているデータとは異なる。だから教科書と同じ結果も得ることはできない。もっとも、非常に近い値を得ることはできる。
> dat <- c(myData)
> group<- as.factor(rep(c(1,2,3), c(10,10,10)))
> result <- glm(dat ~ group)
> summary(aov(result))
Df Sum Sq Mean Sq F value Pr(>F)
group 2 192.709 96.354 38.926 1.111e-08 ***
Residuals 27 66.833 2.475
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(result)
Call:
glm(formula = dat ~ group)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.1754 -0.6511 -0.1770 1.1404 3.0188
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.2847 0.4975 32.732 < 2e-16 ***
group2 -6.0913 0.7036 -8.657 2.85e-09 ***
group3 -4.0843 0.7036 -5.805 3.54e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 2.475306)
Null deviance: 259.542 on 29 degrees of freedom
Residual deviance: 66.833 on 27 degrees of freedom
AIC: 117.17
Number of Fisher Scoring iterations: 2
BOX 3.2の解析例と同じように教科書（Minitab）とRの結果を比べてみよう。
> kyo <- function(group1, group2) 12.1169 + 4.8073*group1 - 2.6145*group2 > kyo(1,0) [1] 16.9242 > kyo(0,1) [1] 9.5024 > kyo(0,0) [1] 12.1169 > R <- function(group2, group3) 16.2847 - 6.0913*group2 - 4.0843*group3 > R(0,0) [1] 16.2847 > R(1,0) [1] 10.1934 > R(0,1) [1] 12.2004
実際、教科書の式とRによって得られた式により得られる各群の平均値の差の絶対値をみてみると、ほとんど差がないことが分かるだろう。これは図にしてみれば分かりやすい（赤い折れ線が教科書のもの、青い折れ線がRのもの）。
> a <- c(16.9242, 9.5024, 12.1169) #教科書の場合 > b <- c(16.2847, 10.1934, 12.2004) #Rでやった場合 > abs(a-b) [1] 0.6395 0.6910 0.0835 #図を描く > x.jiku <- 1:3 #横軸の値 > plot(x.jiku, a, ylim=c(0,20), xaxt="n", type="b", col="red") > axis(1, at=1:3) > points(x.jiku, b, type="b", col="blue")
 |
p51 「BOX 3.4 推定による母数の復元」
ここでやっていることは、要するに真の値（母数）を自分で決めて、それをもとに擬似データを作り出す。そして、そのデータを解析した結果として得られる推定値が、最初に設定した母数とどれほど一致するかというコンピュータシミュレーションをしてみるということである。
そのためには上で定義したY(DUM1, DUM2)という関数のほかに、DSを少し変えたDS2、DS2.repeatを新たに定義する。
#関数DS2の定義
DS2 <- function(g1, g2, g3){
group1.dat <- c()
for(i in 1:g1) group1.dat[i] <- Y(1, 0)
group2.dat <- c()
for(i in 1:g2) group2.dat[i] <- Y(0, 1)
group3.dat <- c()
for(i in 1:g3) group3.dat[i] <- Y(0, 0)
dat <- c(group1.dat, group2.dat, group3.dat)
group <- as.factor(rep(c(1, 2, 3), c(g1, g2, g3)))
res <- cbind(dat, group)
res
}
#関数DS2.repeatの定義
DS2.repeat <- function(n, vec){
mu <- c()
for(i in 1:n){
dat <- DS2(vec[1], vec[2], vec[3])
Y <- c(dat[,1])
X <- as.factor(dat[,2])
result <- glm(Y ~ X)
intercept <- result$coefficients[1]
x2 <- result$coefficients[2]
x3 <- result$coefficients[3]
func <- function(X2, X3) intercept + x2*X2 + x3*X3
mu[i] <- func(0, 1)
mu
}
}
これらの関数の内容は理解できなくてもよいが、どのような関数であるかは知っておく必要がある。DS2は返される結果が少し違うだけで、基本的にはDSと同じものと考えればよい。DS2.repeatは実験を繰り返す回数をnで指定し、各群のデータ数をvecで指定する。例えば、各群のデータが5で、その解析を10回行う場合はn=10, vec=c(5, 5, 5)と指定する。すると、最終的には母数μ=12.2を推定した値が10個得られるわけである。
教科書と同じことを繰り返すのもつまらないので、ここでは少し工夫して、各群のデータ数が増えるにしたがって、どれほど母数へと推定値が近づくのかを折れ線グラフで表現してみる。
> n5 <- DS2.repeat(10, c(5,5,5)) #各群のデータ数が5の場合 > plot(n5, type="b", ylim=c(5,20)) #n=5のときのデータをプロット > abline(h=12.2) #母数μ=12.2のラインを引く > n10 <- DS2.repeat(10, c(10,10,10)) #各群のデータ数が10 > x.jiku <- 1:10 #x軸の値 > points(x.jiku, n10, type="b", col=2) #n=10のときのデータをプロット > n20 <- DS2.repeat(10, c(20,20,20)) > points(x.jiku, n20, type="b", col=3) > n40 <- DS2.repeat(10, c(40,40,40)) > points(x.jiku, n40, type="b", col=4)
 |
p52 練習問題
母集団のもつ変動は解析にどのような影響をもたらすか
これは単に
> Y <- function(DUM1, DUM2) 12.2 + 5.0*DUM1 + -2.5*DUM2 + 1.5*rnorm(1)
と定義した関数において、K4 = 1.5の部分を2, 4, 6, 8, 16と変更して定義しなおしてから、前回やったのと同じようにすればよいだけである。
例えば、K4 = 2の場合を以下に示す。
#関数Yを定義しなおす > Y <- function(DUM1, DUM2) 12.2 + 5.0*DUM1 + -2.5*DUM2 + 2*rnorm(1) > vec <- c(10, 10, 10) > sig2 <- DS2.repeat(10, vec) #各群のデータ数10、実験回数10 > sig2.sa <- abs(sig2 - 12.2) #母数muとの差 > sig2.sa [1] 0.06974111 0.16662609 0.73873711 0.57698105 0.65905538 1.55222552 [7] 0.47226118 0.71881154 0.64416591 0.47092143
同様にしてK4 = 4, 8, 16の場合を行う。
> Y <- function(DUM1, DUM2) 12.2 + 5.0*DUM1 + -2.5*DUM2 + 4*rnorm(1)
> sig4 <- DS2.repeat(10, vec)
> sig4.sa <- abs(sig4 - 12.2)
> Y <- function(DUM1, DUM2) 12.2 + 5.0*DUM1 + -2.5*DUM2 + 8*rnorm(1)
> sig8 <- DS2.repeat(10, vec)
> sig8.sa <- abs(sig8 - 12.2)
> Y <- function(DUM1, DUM2) 12.2 + 5.0*DUM1 + -2.5*DUM2 + 16*rnorm(1)
> sig16 <- DS2.repeat(10, vec)
> sig16.sa <- abs(sig16 - 12.2)
> sig.dat <- cbind(sig2.sa, sig4.sa, sig8.sa, sig16.sa)
> sig.dat
sig2.sa sig4.sa sig8.sa sig16.sa
[1,] 0.06974111 0.5300125 3.6191772 2.383226
[2,] 0.16662609 0.8747844 0.8173634 12.858148
[3,] 0.73873711 0.1803338 1.0435363 6.100796
[4,] 0.57698105 0.4956233 1.1771887 8.146036
[5,] 0.65905538 0.1991413 2.6654236 4.267054
[6,] 1.55222552 3.0735626 1.7485692 1.480965
[7,] 0.47226118 0.1646721 0.4989577 2.609379
[8,] 0.71881154 0.9679702 1.1464112 3.008222
[9,] 0.64416591 0.8652119 1.9769265 7.060064
[10,] 0.47092143 1.5253496 4.7152296 3.055320
> Mean.sig <- apply(sig.dat, 2, mean) #母数との差の平均
> Mean.sig
sig2.sa sig4.sa sig8.sa sig16.sa
0.6069526 0.8876662 1.9408783 5.0969212
#母数との差の平均をプロットする
> plot(Mean.sig, type="b", xaxt="n")
> axis(1, label=c("sig2", "sig4", "sig8", "sig16"), at=c(1,2,3,4))
 |
上の図は、σの値ごとの推定値の平均値と母数との差をプロットしたものである。このようにしてみると、誤差の標準偏差（σ）が大きくなるほど、母数との差が開いていくのだということがよく分かる。この考察は教科書の回答とは異なるものであるが、シミュレーション結果をどのように解釈するかを自分で考えるのも統計学の醍醐味である。
