表1 ●分散分析 Analysis of variance (ANOVA) 単因子実験 B1 B2 A1 〇 〇 A2 〇 × 参考文献 鷲尾泰俊著 実験計画法入門 日本規格協会 対象物の特性を調べるために、条件を変えて反応をみることを考える。 どのように条件を変えると少ない条件回数で正しい情報を多く得られるかを 表2 考える方法は実験計画法(design of experiment)と呼ばれる。 要因実験 B1 B2 水準(level)と呼ぶことにする。単因子実験(single factor experiment)は、一度に A1 〇 〇 因子を1つずつ取り上げて最適条件を探す。これに対して要因実験(factorial A2 〇 〇 対象物の変化する特性を実験の特性値、条件を因子(factor)、条件の各場合を experiment)は問題となる要因を同時に取り上げ、因子の水準の全ての組み合わせについて実験を 行い最適条件を探す。各因子の効果が加法的でないとき、交互作用(interaction)が存在するとい う。表1では、A を A1 から A2 に変えた時の効果は B が B1 のときしか知ることができない。とこ ろが表2では1回実験を追加しただけで、全ての場合(4つの場合)の効果を知ることができる。 分散分析では、ある特性についての観測値がいろいろな値に変化しているとき、観測値の変化の 大きさを、それを説明するために取り上げた説明変数(条件)に 水準 A1 A2 ・ Aa よって影響された部分と、それで説明できない部分とに分ける。 1回目 X11 X21 ・ Xa1 Ⅰ.データの構造模型(母集団モデル) 2回目 X12 X22 ・ Xa2 3回目 . ・ ・ ・ ・ ・ ・ ・ X2n ・ Xan xij = m + a i + eij xij :Ai 水準でのj番目の繰り返しデータ m :一般平均(定数) . n回目 X1n a i :Ai 水準の効果、ただし åa i = 0 eij :誤差(変量)で平均値 とする i 0,分散σ2 の正規分布 N 和 (0,σ 2)とする X1s 平均 X2s Xas X1m X2m Xam Ⅱ.1因子分散分析(一元配置の分散分析) a種類の水準の各条件をn回繰り返して得られたデータ X11,.. .Xan が与えられるとき、全デー タのばらつきは個々のデータの全体平均値 Xmm からの偏差平方和としてはかられる。 総平方和(ST)= A 間平方和(SA) + 誤差平方和(Se) a n 2 a a n i =1 i =1 j =1 2 å å ( x ij - x mm) = n å ( x im - x mm) + å å ( x ij - x im) i =1 j =1 変動要因 平方和 自由度 2 平均平方(不偏分散) A間 SA ΦA=a-1 VA=SA /ΦA 誤差 Se Φe =a(n-1) Ve=Se /Φe 総 ST ΦT =an-1 F0 VA/Ve 平均平方 VA, Ve は統計量で、 その期待値(不偏推定量)は [VA] =σ2+nσA2, [Ve] =σ2 ただし σA2= 1 a 2 åa i a - 1 i =1 = 1 a -1 å (x im - xmm )2 検定したい仮説 H0 は、A の水準間に差はない(σA2=0)という仮説である。 期待値は、Ve は平均的にはσ 2 になり、VA はσ2+nσA2 になる。 もし、H0 が正しければ VA も Ve も平均的には同じ値になるはず。 よって比 VA/Ve を計算してみる。 H0 が正しければ VA/Ve は1に近く、H0 が正しくなければ VA/Ve は1より大きくなるであろう。 統計理論から H0 が正しければ VA/Ve は自由度(φA,φe)の F 分布に従うことが分かっている。 よって、検定の有意水準をαとすれば、VA/Ve>F(φ A,φe;α)ならば H0 を棄却する。 VA/Ve≦F(φ A,φe;α)ならば H0 を棄却しない。 31 水準 (例)ある製品のエネルギー効率を測定した。因子として季節を とりあげ、水準として四季の4水準を選んで5回繰り返して測定 した。このデータを分散分析して季節による差があるかを調べよ。 1因子実験(一元配置)を EXCEL の表計算で行う場合 春 夏 秋 冬 1回目 79 80 83 83 2回目 75 86 89 78 3回目 79 81 85 78 4回目 76 84 84 79 5回目 77 81 84 81 1因子実験(一元配置)を EXCEL の分析ツールを使う場合 1.条件を横に、繰り返し数を縦にデータを並べる。 2. [データ]メニューから[データ分析]を選択する(EXCEL2007 の場合) 3.さらに「分散分析:一元配置」を選択する。 32 2.最初にここをクリックしてデータ領域をドラッグして選択する。 危険率の値を入れる。データの先頭がラベルのときは、チェックを入れ最後に OK ボタン押す データ先頭がラベル のときにチェック 危険率(間違える確率)αの値 3.結果の表が別シートとして出力される (計算式から求めた値と一致) FDIST(11.438, 3, 16) から求められる (分散比が 11.4 になる確率値(面積)) A間 誤差 FINV(0.05, 3, 16) 危険率と自由度か ら決まる分散比 F 分布曲線 この線より右側の面積(確率)が 0.05 この面積が 0.000296 3.2388 計算した F 値がこの値より左にくれば、水準 間で差がないとみなす。右にくれば、差があ ると見なす(ただし間違う確率(危険率)が 5%) 33 11.438 実際のデータから 計算した F 値 Ⅲ.2因子要因実験(2元配置) 2要因実験では、要因 A の影響、要因 B の影響を独立に調べることができる他、A と B を同時に 変えることによる相互作用の影響(交互作用)をも調べることができる。 構造模型は、 x ijk = m + a i + b j + (ab )ij + eijk となる。 μ:一般平均,α i:Aの主効果;Σαi=0 ,βi:B の主効果;Σβj=0 a b i j (αβ)ij:AとBの交互作用効果(A×B) ; å (ab ) ij = 0 , å (ab ) ij = 0 A の検定、B の検定、A×B の検定をするために総平方和を3つの成分と誤差成分に分解する。 ST = SA + SB + SAxB + Se a b n 2 すなわち、 SB a b i j 2 2 å å å ( x ijk - x mmm ) = nb å ( x imm - x mmm ) + na å ( x mjm - x mmm ) i j k a b i j 2 a b n i j k 2 + n å å ( x ijm - x imm - x mjm + x mmm) + å åå ( x ijk - x ijm) SAXB Se 最初の添字 i ・・要因 A の水準(レベル)を表す(i=1, 2, ・・, a) 変数の記号 2番目添字 j ・・要因 B の水準(レベル)を表す (j=1, 2, ・・, b) 3番目添字 k ・・同じ条件での繰り返す実験番号 (k=1, 2, ・・, n) xijk x : = ijm 要因 A の水準 i、要因 B の水準 j の k 番目のデータを表す 1 n å n k =1 x ijk : 要因 Ai、Bj での繰り返し実験データの平均 (添字mは平均 mean を意味 する) ximm = x mjm 1 b n åå bn j=1 k =1 xijk 1 a n = åå an i =1 k =1 x ijk : 要因 Ai でのすべての実験データの平均 : 1 a b n = xmmm abn ååå xijk i =1 j =1 k =1 変動要因 平方和 要因 Bj でのすべての実験データの平均 : すべての実験データの平均 自由度 平均平方 分散比 平均平方期待値 SA a-1 VA VA/Ve σ 2+bnσA2 B間 SB b-1 VB VB/Ve σ 2+anσB2 SAXB (a-1)(b-1) VAXB 誤差 Se ab(n-1) 総 ST abn-1 V V A = S A /(a - 1) B = S B /(b - 1) = S AXB /(a - 1)(b - 1) AXB e = S e / ab(n - 1) s A2 = 分析には不要 A間 A×B 間 V V VAXB/Ve Ve σ2+nσAXB2 σ2 1 a 2 åa i a - 1 i =1 a i = ximm - x mmm s B2 = 1 b 2 åb j b - 1 j =1 b j = x mjm - x mmm 2 s AXB = a b 1 (ab ) ij2 åå ( a - 1)(b - 1) i =1 j =1 (ab ) ij = xijm - ximm - x mjm + xmmm 34 二因子の場合の EXCEL を使った分散分析(二元配置)の方法 1.二因子と、繰り返しの3次元データのうち、繰り返しを内部に織り込んだ表を作成する。 2.メニューの「データ」から「分析ツール」を選択し、 さらに、「分散分析:繰り返しのある二元配置」を選択する 3.下記のようにラベル(A1~A4、B1~B4)を含むデータを入力範囲に入れる ラベルも入力範囲に含ませる 危険率 繰り返し数 3.OK ボタンで出力の表が出る。最後の分散分析表のみを使う。 データから計算した 分散比 平方和 危険率がα=0.05 となるときの 分散比の値 データの分散比と なる F 分布の確率値 A 要因 B 要因 3.717×10-6 の意味 【課題 18】2つの要因 A,B の水準を変えて行った実験結果は下図のようになった。 A および B の主効果、A と B の交互作用効果があるかどうか調べよ。 (危険率 5%とする) Xijk A1 A2 A3 B1 B2 B3 B4 28 33 38 36 25 39 43 44 32 32 36 37 28 35 35 34 36 34 32 30 33 33 33 32 a= b= n= 35 3 4 2 [手順] 1.上記原データ表から、繰り返し平均の表(表2) 、Ai、Bj の平均、Ai, Bj の偏差平方和を求め て、SA,SB を求める 2.交互作用の各項の表(表3)を作成して SAXB を求める 3.AiBj の偏差平方和の表(表4)を作成して Se を求める 4.各データ値の総平均からの偏差平方和の表(表5)を作成して、ST を求める 5.A,B の各自由度、各平均平方(不偏分散) 、分散比の順で求めて分散分析表を作成する 6.各 F 値を求めて、分散比と比較して有意差の有無を検定する (下記表は http://www.coins.tsukuba.ac.jp/~fukui/da/4anova.xls 条件Aの水準数:a = 条件Bの水準数:b = 繰り返しデータ数:n = 3 4 2 原データ表 Xijk B1 28 25 32 28 36 33 A1 A2 A3 繰り返しの平均データ表 Xijm B1 A1 A2 A3 Bj の平均 Xm1m Bj の偏差平方 にある) 学籍番号: 記号の意味 Xijk : Xijk はAの水準が i で、Bの水準が j で、繰り返し回の番号がk で表されるデータ X1mm : A1の水準の全てのデータの平均 X2mm : A2の水準の全てのデータの平均 以下同様 Xm1m : B1の水準の全てのデータの平均 Xm2m : B2の水準の全てのデータの平均 以下同様 X12m : A1でB2の水準の繰り返しデータの平均 以下同様 Xmmm : 全てのデータの平均 氏名: B2 33 39 32 35 34 33 B3 38 43 36 35 32 33 B4 36 44 37 34 30 32 B2 B3 B4 Xm2m Xm3m Ai の平均 X1mm = X2mm = X3mm = Xmmm = Xm4m (Xm1m-Xmmm)^2 (Xm2m-Xmmm)^2 (Xm3m-Xmmm)^2 (Xm4m-Xmmm)^2 Ai の偏差平方 (X1mm-Xmmm)^2 = (X2mm-Xmmm)^2 = (X3mm-Xmmm)^2 = 上3つの和:S'A= S A=n*b*S'A= Bj の偏差平方の和 S'B= AXB 交互作用の各項の表 このAi xBj の各項の計算 表 (Xijm-Ximm-Xmjm+Xmmm)^2 B1 B2 A1 A2 SB =n*a*S'B= B3 B4 AxB間偏差平方の和 S' AXB = A3 SAXB =n*S' AXB = Ai・Bj の偏差平方 (Xijk - Xijm)^2 A1 B1 B2 B3 B4 A2 A3 誤差平方和 Se = 原データの全平均からの偏差平方表 (Xijk - Xmmm)^2 B1 A1 B2 B3 B4 A2 A3 全データの偏差平方和 =FINV(危険率, 条件の自由 ST = SA + S B + SA XB + Se = α = 平方和 自由度 平均平方(分散 分散比 条件:A 条件:B 交互作用:AxB 誤差:e 合計: 度, 誤差の自由度) で危険率の確率となるF値 (分散比)が求まる 0.05 確率値 F境界値 =FDIST(分散比, 条件の自由 度, 誤差の自由度)で、分散比 (F値)の値になる確率が求ま る 結論: B の主効果は(ある、または、ない) A×B の交互作用効果は(ある、または、ない) A の主効果は(ある、または、ない) B1 オプション課題 A1 右図の 2 元配置の条件での実験結果のデータ A2 (繰り返し数=2)のデータで、要因 A, 要因 B, 交互効果 A×B の効果を分析せよ。 36 A3 A4 B2 1 0 1 3 4 6 4 2 B3 3 5 8 4 8 11 5 2 4 1 6 5 8 5 1 4 付録 a n a a n i =1 j =1 i =1 i =1 j =1 1. å å ( x ij - x mm)2 = n å ( x im - x mm)2 + å å ( x ij - x im)2 a a n åå ( x i =1 j =1 ij の証明 n - x mm ) 2 = åå ( xij - xim + xim - xmm ) 2 i =1 j =1 = åå ( xij - xim ) 2 + åå ( xim - x mm ) 2 + 2åå ( xij - xim )( xim - xmm ) i a 第 3 項: n åå ( x i =1 j =1 j i j i - xim )( xim - x mm ) = åå xij xim - åå xij x mm - åå xim + åå xim xmm 2 ij i a n åå x i =1 j =1 a x ij im j i i a n a i =1 j =1 i =1 n a n - åå xim = -n å xim i =1 j =1 a i =1 j =1 n j 2 2 i =1 n åå x i となり第3項は 0 となる。 a 2 j 2 i =1 j =1 ここで、上式の各項は、 a j = å xim å xij = n å xim - åå xmm xij = -anxmm 第2項は、 j a x im mm = nxmm å xim = anxmm i =1 a åå ( xim - x mm ) 2 = n å ( xim - xmm ) 2 i =1 j =1 となるのは明らか。 i =1 値 xam xmm x2m x1m x2n xan xa1 x1n x21 x11 データ番号 A1 A2 ・・・ Aa A1 A2 ・・・ Aa 値 結局、分散分析とは 上の図から想定されるように、データ群を、データ模型に沿うとみなした各水準のデータ群と 沿わないと、データ模型から外れたデータ群とに分けて、模型に沿うデータは各水準毎の代表値 の集合のばらつき具合と、データ模型に沿わない部分のばらつき(Se)とがほぼ同じであれば、デ ータ模型は正しいと判断するデータ解析手法。 37 付録 2.<VA> =σ 2+nσA2,<Ve> =σ2 の証明 期待値に関する公式 1.cを定数とするとき、[c] = c ]は期待値を表わす ただし、[ 2.x1,x2,・・,xn が変数で、c1,c2,・・,cn が定数のとき、 [ c1x1+c2x2+・・+cnxn ] = c1[x1] + c2[x2] + ・・ +cn[xn] 3.変数 x1,x2,・・,xn が互いに独立で、いずれも期待値(平均値)が 0、分散がσ 2 のとき、 V (x) = n 1 2 s , [å ( xi - x ) 2 ] =(n-1)σ2 n i =1 ただし、 x は平均、V( )は分散を表わす これらを使って、[VA] =σ2+nσ A2,[Ve] =σ2 を示す まず、 [V A ] = [ = xim = 1 n 1 n xij = å ( m + a i + eij ) = m + a i + eim å n j =1 n j =1 xmm = 1 a 1 a x = ( m + a i + eim ) = m + emm å im a å a i =1 i =1 SA n a n a ]= [å ( xim - xmm ) 2 ] = [å (a i + eim - emm ) 2 ] a - 1 a - 1 i=1 a - 1 i =1 n 2 [å a i ] + 2å a i [(eim - emm )] + [å (eim - emm ) 2 ] a -1 i i i ここで、第2項は [ eij ] = 0 の仮定より 0 となる。 第3項は、a個の独立な変数: e1m , e2 m ,, eam の、それらの平均 emm からの偏差平方和であること から、 [ a å (eim - emm ) 2 ] = (a - 1)V (eim ) = (a - 1) i =1 s2 n a n s2 2 2 {å a i + 0 + (a - 1) } = s 2 + ns A a - 1 i=1 n a n Se 1 [Ve ] = [ ]= [åå ( xij - xim ) 2 ] a(n - 1) a(n - 1) i =1 j =1 [V A ] = = 2 ここで、 s A = 1 a 2 åai a - 1 i=1 a n a n 1 1 [åå (eij - eim ) 2 ] = [å (eij - eim ) 2 ] å a(n - 1) i =1 j =1 a(n - 1) i =1 j =1 n ここで、 を得る。したがって、 å (e j =1 ij - eim ) 2 はn個の独立な変数: ei1 , ei 2 , , ein の、それらの平均 eim からの偏差平方 和であることから、 [ n å (e j =1 したがって、 [Ve ] = ij - eim ) 2 ] = (n - 1)V (eij ) = (n - 1)s 2 a 1 å (n - 1)s 2 = s 2 a(n - 1) i =1 となる。 38 となる。 ●多変量の相関を調べるときの注意 1.相関分析ができないデータ(散布図で調べる) (1)x軸、またはy軸に平行な直線、または、 (a) (b) 特定の直線上にのるデータ、もしくは曲線上の 固定した値しかとらないデータは、確率変数と みなせない (a),(b) (2)明らかに「比例的関係」に無いデータ(c) (3)2つ以上のかたまりがあるデータ(d)(e) (c) (d) (4)異常値があるデータ(f) (5)x軸、またはy軸の計測値がカテゴリー名で あったり、単なる順序尺度値であるデータ(g) (e) 2.相関分析における「要注意」データ(散布図では わかりにくい)(全て例は架空) (1)混在によって、相関がないのにあるように見える 足の大きさと数学力・・単に学年を全部あわせただけ(h) (2)混在によって、相関があるのにないように見える (f) 数学力と英語力に男女差があるのに混ぜると ないように見える(i) (3)切断によって、相関がないのにあるように見える 世界各国の酒類の消費量と GNP の相関(j) (有色人種の国を除外していた) (g) (4)切断によって、相関があるのにないように見える 入社試験成績と入社後の活躍の例(k) (h) (i) (j) (k) 3.データ間の関係は、一般にはわかりにくく、分析で得られた関係がかならずしも 正しいとは限らない。別の分析をすれば、もっともらしい別の関係が出てくることも ありうる。このような過ちを犯す危険性があることを認識しておく必要がある。 39 ●多変量統計解析 1.いろいろな要因によってある項目を推定(説明)したいとき、あるいは 観測されている複数個の項目を代表する総合的指標を求めたいとき 要因 変数 外的 基準 外的基準の有無 有 量的変数 量的変数 質的変数 (単/重)回帰分析 判別分析 数量化Ⅰ類 数量化Ⅱ類 質的変数 無 主成分分析 因子分析 数量化Ⅲ類 2.その他、ものや項目を似たもの同志をまとめるように分類したい 3.項目間の複雑な相関関係を説明する潜在的構造を知りたい -> クラスター分析 -> 潜在構造分析 分析法の特徴 ・(重)回帰分析 ある変数yとそれに影響を及ぼすと考えられる他の変数 x=(x1, x2, .. , xp )に 関するデータに基づいて、新たなxに対して y を予測する1つの方法 ・判別分析 いくつかの変数 x1, x2, …, xp に関して群毎に得られている過去のデータに基づき、 これらの変量の値から、個体がどの群に属するかを判別(予測)する方法 (回帰分析の目的変数yの名義尺度版) ・主成分分析 多くの変量の値を、できるだけ情報の損失なしに少数個の総合的指標で代表させ る方法 ・因子分析 多くの変量の持っている情報を少数個の潜在的因子で説明しようとする方法 潜在構造分析 ・数量化 I 類 質的な要因に関する情報に基づいて、量的に測定された外的基準(目的変数)の 値を説明、予測するための方法 ・数量化 II 類 質的な要因に基づいて、質的な形で与えられた外的基準を予測あるいは判別 する方法 ・数量化 III 類 予測すべき外的基準がないとき、個体の種々のカテゴリーへの反応に基づい て、個体とカテゴリーの両方を数量化し分類を行う方法 ・クラスター分析 異質なものの混ざり合っている対象を、類似度に基づいて似たもの同志を集 めていくつかのクラスター(集落)に分類する方法 ・多次元尺度構成法 類似性を数値で評価した表から、評価した特徴を抽出し、グループ化する 40 ●線形(単)回帰分析 (Regression Analysis) 男女別出生児数及び死亡者数の推移 http://www.stat.go.jp/data/jinsui/2009np/pdf/gaiyou.pdf#pa 人口動態をみるのに、説明変量 ge=1 xを時間(年)とし、人口変化を被説明 変量yで記述する場合の2変量から 数理モデル化を行う。2変量間の関係 は線形関係に近似可能であり、線形 モデルを適用できることがわかる。 〇モデル化の目的 変量xi とyi に最もよく当てはまる 直線の方程式を求める。 〇モデル化の方法 直線と標本点との間の偏差を含んだ何らかの指標により、あてはまりの良い直線を決定する。 〇最小2乗基準 偏差の平方和を最小とする直線を最も当てはまりの良い直線と考える(最小2乗法) 。 回帰分析とは、ある変数yとそれに影響を及ぼすと考えられる他の変数xに関するデータに基づ いて bˆ + bˆ x 0 ® 1 y のように予測する1つの方法である。経済や経営の分野での予測、品質管理 (QC)における最適値の探索など広い分野に応用される。 1.線形回帰モデル(linear regression model) y =b i 0 +b 1 ( i=1, 2, . . . n) ・・・(1) ただし、期待値: E{ei } = 0 ,分散: V {ei } = s x +e i ・・・(2) yˆ i = bˆ0 + bˆ1 xi 回帰式 2 によって説明変数(独立変数)xi から目的変数(従属変数) yi を予測することを考える。このとき予測誤差(残差)ei は式(1)から次のように表される。 ei = yi - yˆ i = yi - ( bˆ0 + bˆ1 xi ) (i=1, 2, . . . ,n) ・・・(3) 予測誤差(残差)ei は正負の値を取りうるので、全体的に小さくするために、その平方和 n n Se = å ei = å{ y 2 i =1 i =1 ˆ ˆ 2 i - ( b 0 + b 1 xi)} ・・・ (4) を最小にするようにβ0, β1 を求める(最小2乗法)。そのために、(4)式をβ0,β1 で偏微分する。 ¶ å ei 2 ¶ bˆ i =1 0 ¶ å ei 2 ¶ bˆ n =å (5)式より、 ¶ bˆ ( y i - bˆ 0 - bˆ 1 x i) 2 n = -2 å ( y - bˆ - bˆ i =1 i 0 1 x)=0 ・・・(5) i 0 n ¶ i =1 ¶ bˆ =å 1 ¶ ( y i - bˆ 0 - bˆ 1 x i) 2 n = -2å x i ( y - bˆ - bˆ i =1 i 0 1 x)=0 ・・・(6) i 1 å y - å bˆ i 0 - bˆ1 å xi = 0 すなわち 41 nbˆ0 + bˆ1 å xi = å yi ・・・(5)' (6)式より、 åx y i i 2 - bˆ0 å xi - bˆi å xi = 0 すなわち bˆ0 åx i 2 + bˆ1 å xi = å xi yi æ å xi 2 - å xi öæ å yi ö ç ÷ç ÷ å xi2 ö÷æç bˆ0 ö÷ = æç å yi ö÷ から、 æç bˆ0 ö÷ = çè - å xi n ÷øçè å xi yi ÷ø ç bˆ ÷ å xi ÷øçè bˆ1 ÷ø çè å xi yi ÷ø n å xi è 1ø 2 å x i å xi æ n çå x i è (5)',(6)'より ç ・・(6)' となる。 各成分は、 bˆ1 = bˆ0 = n å xi y i - å xi å y i n å xi - (å xi ) 2 2 = 1 1 xi × å yi å n n = 1 2 2 å xi - n ( n å xi ) åx y i åx å y - åx åx y nå x - (å x ) 2 i i 2 i å(x i -n× y m å xi - x m å xi y i åx y åx i i 2 i - n × xm y m - n × xm 2 i = 2 i ここで、 S xx = i i åx i 2 i - nxm 2 ただし、xm = - xm ) 2 = å xi - 2 xm å xi + n × xm = å xi - n × xm 2 2 2 2 n 1 n å , ym = n1 å yi n i =1 xi i =1 2 S xy = å ( xi - xm )( yi - y m ) = å xi yi - x m å yi - y m å xi + n × xm y m = å xi yi - n × x m y m とすれば、 また、 bˆ1 = S xy / S xx , bˆ0 = y m - bˆ1 xm が得られる。 S e = å { yi - ( bˆ 0 + bˆ1 xi )}2 = å { yi - y m - bˆ1 ( xi - x m )}2 = S yy - bˆ1 S xy 回帰モデルの誤差 ei の母分散σ2 の推定値は 平方和の分解 sˆ 2 = Ve = S e / f e = S e /( n - 2) 残差平方和 S e 2.寄与率 ・・(8) ・・・(9) 回帰平方和 S R S yy = å ( yi - y m ) 2 = å { yi - ( bˆ0 + bˆ1 xi )}2 + å {( bˆ0 + bˆ1 xi ) - y m }2 = SR + Se 寄与率 ・・・(7) 2 ここで、 S R = bˆ1S xy = S xy / S xx 2 R 2 = S R / S yy = 1 - S e / S yy = S xy /( S xx S yy ) = rxy 自由度調整済寄与率 2 R* = 1 - S e / fe S /( n - 2) = 1- e S yy / fT S yy /(n - 1) 42 2 ・・・(10) ・・・(11) ・・・・(12) 3.EXCEL =(B3-B$14)^2 $ を つ け る とコ ピ ーし て も、$の次のアドレス値は 変わらない。(通常はコピー すると 参照先のアドレ ス は相対的に変化する) =(C3-H3)^2 =average(B3:B12) =count(B3:B12) グラフの作成方法 1.xとyのデータ領域(B3:C12)をマウスをドラッグして選択 2.メニューの「挿入」からグラフの「散布図」の点群表示を選択して作図する 3.回帰式、寄与率(相関係数の2乗)を求めるときには 3.1 グラフの中のデータの1点をマウスでクリックすると「グラフツール」が出る 3.2 「グラフツール」→「レイアウト」→「近似曲線」→「線形近似曲線」を選択する 3.3 「グラフツール」→「レイアウト」→「軸」→「主軸横軸」→「その他のオプション」 で、軸の数値範囲を決めることが可能。縦軸も同様 チェックマークをマウスでクリックして付けて[OK]ボタンを押す EXCEL の分析ツールによる回帰分析の結果(下記の表) ・ 決定係数が 0.93 であるので、 回帰式で y 値の 93%が説明できる 分析ツールが[データ]タブで現れるリ ストに見つからないとき エクセルの左上の Office マーク 概要 回帰統計 重相関 R 0.966413 重決定 R2 0.933955 補正 R2 0.925699 標準誤差 1.47071 観測数 10 から[エクセルのオプション]→[アドイ ン]を選択して分析ツールを選択して、 [設定]ボタンでアドインのリストが出 るので[分析ツール]にチェックマーク を入れて[OK]を押す 分散分析表 回帰 残差 合計 自由度 変動 分散 観測された分散比 有意 F 1 244.6961 244.6961 113.1287324 5.35E-06 8 17.3039 2.162988 9 262 切片 X値1 係数 標準誤差 t 61.85277 1.58778 38.95551 4.523033 0.425249 10.6362 EXCEL 関数で求める場合 回帰係数=linest(y 範囲,x 範囲) y 切片=intercept(y 範囲、x 範囲) P-値 下限 95% 上限 95% 下限 95.0%上限 95.0% 2.07227E-10 58.19134 65.5142 58.19134 65.5142 5.34611E-06 3.542407 5.50366 3.542407 5.50366 43 4.R による単回帰分析 計算例 > x <- c(2.2,4.1,5.5,1.9,3.4,2.6,4.2,3.7,4.9,3.2) > y <- c(71,81,86,72,77,73,80,81,85,74) > result <- lm(y ~ x) > summary(result) 入力 Call: lm(formula = y ~ x) Residuals: Min 1Q -2.3265 -0.7849 lm(y~x)は、 linear model で y を x で 表す意味 Median 3Q Max -0.4219 0.8890 2.4120 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 61.8528 1.5878 38.96 2.07e-10 *** x 4.5230 0.4252 10.64 5.35e-06 *** --Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 出力 y = 4.523 x + 61.8528 が線形回帰モデル Residual standard error: 1.471 on 8 degrees of freedom Multiple R-Squared: 0.934, Adjusted R-squared: 0.9257 F-statistic: 113.1 on 1 and 8 DF, p-value: 5.346e-06 > plot(x,y) > abline(result) データプロットと 回帰直線描画指令 入力 出力グラフ 44 5.推定値の信頼性 ^ 5.1 β1 の変動について S 1 bˆ1 = xy = S xx S xx n å(x i =1 i - x )( yi - y ) = 1 1 {å ( xi - x ) yi - y å ( xi - x )} = S xx S xx å (x i - x ) yi と変形できる。xi は選択した値であり定数、yi は xi に対して抽出された確率変数である。そうす ^ ^ ると、上の式で、β1 は確率変数 yi の一次式となっている。したがって、β1 の偏差平方和 Sβ 1 は、 各項の分散の和になるので、 S b 1 = ^ 5.2 β0 の変動について Se 2 S xx å (x i - x)2 = Se S xx 1 n bˆ0 = å yi - x bˆ1 も同様に、yi の一次式となっているので、 n i =1 1 1 x2 Sb 0 = ( )Se + x 2 Sb 1 = Se ( + ) n n S xx 6.回帰係数の検定と推定 統計量 bˆ1 は平均β1、分散(σ2/Sxx)の正規確率分布に従う確率変数である。 ˆ 標準化を行うと u = b1 - b1 は平均0、分散1の標準正規分布に従う確率変数となる。 s 2 / S xx ˆ 分散σ2 の推定値として(9)式の Ve を代入すると、 t = b1 - b1 は自由度 n-2 のt分布に従う。 Ve / S xx 帰無仮説 対立仮説 H0:β 1=0(xが変化してもyは変化しないので回帰式に意味がない) , bˆ1 H1:β 1≠0(回帰式に意味がある)として検定統計量: t 0 = を求め、 Ve / S xx |t0|≧t(φe,α),(φe=n-2)なら有意水準αで有意(回帰式に意味がある)と判定する。 {t(φ,P)}2=F(1,φ;P)の関係から、|t0|≧t(φe,α)の判定は、以下の分散比の判定と同じになる。 2 F0 = t0 = 2 bˆ1 S = R ³ F (1, f e ; a ) Ve / S xx Ve β 1 の95%信頼区間は右のようになる。 ・・・(13) V bˆ1 ± t (f e ,0.05) e S xx 【演習 19】(回帰分析) 右のデータ(右表)から回帰式、寄与率を求めてみよう。 回帰直線のグラフも描くことでイメージがつきやすい。 Excel と R の両方でやってみて、対応関係を理解しよう。 45 ・・・(14) 説明変数 目的変数 xi yi 42 80 49 85 32 74 27 73 34 77 25 72 22 71 37 81 41 82 55 86 ●重回帰分析 1.説明変数pがp=2個(以上)の場合の解析方法 重回帰モデル: yi = b 0 + b1 x1i + b 2 x2 i + e i ,ただしε i は平均 0、分散σ 2 の正規分布とする。 i 番目の予測値 yˆ i と残差 ei を次で定義する。 yˆ i = bˆ0 + bˆ1x1i + bˆ2 x2i ・・・・(20), ei = yi - yˆ i = yi - ( bˆ0 + bˆ1 x1i + bˆ2 x2 i ) ・・・(21) n 残差平方和 Se = å ei 2 = å { yi - ( bˆ0 + bˆ1 x1i + bˆ2 x2i )}2 ・・(22) i =1 を最小にする bˆ0 , bˆ1 , bˆ2 を求める。 そのために、 S e を bˆ0 , bˆ1 , bˆ2 の各々で偏微分して 0 とおいて連立方程式を作る。 b 0 n + b1 å x1i + b 2 å x2i = å yi b 0 + b1 x1m + b 2 x2 m = y m b1S11 + b 2 S12 = S y1 b 0 å x1i + b1 å x1i + b 2 å x2 i x1i = å x1i yi 2 b 0 å x2 i + b1 å x1i xi 2 + b 2 å x2 i = å x 2i y i 2 b1S12 + b 2 S 22 = S y 2 ここで、 S11 = å ( x1i - x1m )2 , S 22 = å ( x2i - x2 m ) 2 , S12 = å ( x1i - x1m )( x2i - x2 m ) S1 y = å ( x1i - x1m )( yi - ym ) , S 2 y = å ( x2i - x2m )( yi - ym ) , S yy = å ( yi - ym ) 2 これらより次を得る。 bˆ = y - bˆ x - bˆ x 0 m bˆ1 = 1 1m ・・・・(23) ただし、添え字mは平均を表わし、^は推定値を表す。 2 2m S y1 S y2 S12 S 22 S11 S12 S12 S 22 , ˆ b2 = S11 S12 S y1 S y2 S11 S12 S12 S 22 (偏回帰係数) ・・・(24) (以降^は省略) 2.説明変数の重要度の評価 各説明変数が目的変数に及ぼす影響の大きさを調べるためには、その係数である偏回帰係数の大 小では比較できない(単位が異なる)ので、各説明変数を規格化して重回帰分析を適用する。そこで 得られる係数は、標準回帰係数と呼ばれるが、これらは次の式でも求めることができる。 b1¢ = b1 ´ S11 , S b 2¢ = b 2 ´ 22 S yy S yy ここで、 b1¢ , b 2¢ は標準回帰係数、 b1 , b 2 は偏回帰係数 標準回帰係数同士を比較して、各説明変数が目的変数に及ぼす影響力を調べることができる。 式(24)で分母が 0 となる場合は解が存在しない。これは多重共線性(multi co-linearity)が存在する。 多重共線性があれば、 xi の偏回帰係数と、 y と xi との単相関係数の符号が異なる。 重回帰モデルの誤差εの母分散σ 2 の推定値は次式で求まる。 sˆ 2 = Ve = S e / fe = Se /(n - p - 1) ・・・・(25) (pは説明変数の個数) S yy = å ( yi - ym ) 2 = bˆ1S1 y + bˆ2 S 2 y + S e = S R + S e ・・・(26) 回帰式はデータと予測値が良く合うほうが良い。そこで理論値と実測値の重相関係数で評価する。 重相関係数: R = å ( y - y )( yˆ - yˆ ) å ( y - y ) å ( yˆ - yˆ i m i 寄与率(決定係数): i m ・・(27) m 2 i 2 m) R = S R / S yy = ( S yy - S e ) / S yy = 1 - S e / S yy 2 ・・・(28) 説明変数の個数が増えると寄与率は大きくなるので、必ずしも推定精度を表現しなくなる。 そこで、自由度調整済寄与率 2 R* = 1 - S e / fe S /(n - p - 1) =1- e S yy / fT S yy /(n - 1) 46 ・・・(29) を用いる。 表 駅前チェーン店の販売高 (例)駅前チェーン店の1ヶ月あたりの販売高と、 店の間口、取扱品目数、各駅の乗降客数を15点に ついてまとめた表が右にある。 販売高 乗降客数 店の間口 取扱品目数 平均値 標準偏差 195.4 45.3 180.9 56.2 3.9 1.4 163.2 61.8 販売高の変動と、要因の変動との関係はどうなって いるか? 各変量と販売高との相関図を作ると 下のようになっており、各説明変量は販売高と 正の相関関係を持つことがわかる。相関係数は 最大でも 0.76 であるから、複数の変量で販売高 の要因が決まることになるのがわかる。 販売高 乗降客数 店の間口 店番号 (万円/月) (百人/日) (m) 1 131 149 2.3 2 197 188 2.4 3 222 282 2.6 4 170 183 2.9 5 168 221 2.4 6 155 128 4.3 7 205 153 4.8 8 269 197 6.6 9 206 125 6.0 10 126 77 3.6 11 261 210 5.4 12 248 258 5.2 13 250 254 3.8 14 157 114 3.0 15 166 174 3.2 16 ? 175 4.5 取扱品目数 (品数) 99 178 110 95 82 104 190 288 169 156 283 170 228 145 151 200 したがって、4つの変量を同時に分析しなければならない。 R = 0.68 R = 0.63 280 260 260 240 240 240 220 220 220 180 200 販売高 200 180 200 180 160 160 160 140 140 140 120 120 120 100 100 50 150 250 350 R = 0.76 280 260 販売高 販売高 280 100 2 4 乗降客数 6 8 0 100 店の間口 200 300 取扱品目数 重回帰分析により推定式を求めると、次の頁のように、 (販売高)=β0+β 1*(乗降客数)+β2*(間口)+β3*(取扱品数) =7.62+0.51*(乗降客数)+12.43*(間口)+0.288*(取扱品数) となる。各変数の係数は、偏回帰係数であるが、どの変数が一番販売高に影響が強いかを調べる ためには、各変数の変動割合を考慮した標準回帰係数を求めて比較する必要がある。 標準回帰係数が一番大きいのは、β1'(乗降客数)=0.63 であるので、これが一番販売高に影響 を与えることがわかる。 47 400 (x1i-x1m) 2 2 (x2i-x2m) (x3i-x3m) 2 2 (yi-ym) 1015.5 50.9 10228.0 4.6 1610.7 2794.9 776.6 260.3 3121.1 10788.3 848.8 5949.6 5348.5 4471.2 47.2 2.56 2.25 1.69 1 2.25 0.16 0.81 7.29 4.41 0.09 2.25 1.69 0.01 0.81 0.49 4121.6 219.0 2830.2 4651.2 6593.4 3504.6 718.2 15575.0 33.6 51.8 14352.0 46.2 4199.0 331.2 148.8 4147.4 2.6 707.6 645.2 750.8 1632.2 92.2 5417.0 112.4 4816.4 4303.4 2766.8 2981.2 1474.6 864.4 和:S11 47315.7 S22 27.8 S33 57376.4 Syy 30713.6 x3m 163.2 ym 195.4 平均値: x1m 180.9 D= x2m 3.9 S11 S12 S13 S12 S22 S23 S13 S23 S33 = (x1i-x 1m) (x 1i-x1m) (x2i-x2m) (x1i-x1m) (x2i-x2m) (x 3i-x3m) *(x2i-x2m) *(x3i-x 3m) *(x 3i-x3m) *(yi-ym) *(yi-ym) *(yi-ym) 50.99 2045.8 102.7 2052.2 103.04 4134.5 -10.7 105.6 -22.2 11.4 -2.4 23.7 -131.47 -5380.3 69.2 2690.1 -34.58 -1415.12 -2.13 -145.5 68.2 -54.2 25.4 1732.3 -60.2 -3258.8 121.8 -1099.7 41.1 2224.9 -21.15 3129.7 -23.7 2135.8 -16.16 2391.7 -25.08 -746.8 24.1 -267.5 8.64 257.3 43.56 2013.4 337.0 1187.4 198.72 9185.3 -117.32 -324.0 12.2 -592.2 22.26 61.5 31.16 747.8 2.2 7208.3 20.82 499.7 43.7 3490.2 179.7 1911.1 98.4 7858.9 100.3 524.5 8.8 4057.2 68.38 357.7 -7.3 4739.0 -6.5 3993.1 -5.46 3538.1 60.18 1217.0 16.4 2567.7 34.56 698.9 4.81 83.8 8.54 201.9 20.58 358.7 S12 -40.7 47315.7 -40.7 8241.4 -40.7 27.8 898.4 S13 8241.4 S23 898.4 8241.4 898.4 = 57376.4 S1y 26002.8 S2y 583.3 3.46E+10 -40.7 27.8 898.4 8241.4 898.4 /D = 57376.4 S11 S1y S13 β2= S12 S2y S23 S13 S3y S33 /D = 47315.7 -40.7 8241.4 26002.8 583.3 31907.8 8241.4 898.4 /D = 57376.4 12.43 S11 S12 S1y β3= S12 S22 S2y S13 S23 S3y /D = 47315.7 -40.7 8241.4 -40.7 27.8 898.4 26002.8 583.3 /D = 31907.8 0.288 重回帰分析 方法 SR 29710.5 Se 1003.1 0.63 0.37 0.39 寄与率:R2 = 1 - Se/SR = 26002.8 583.3 31907.8 R を用いた 和: S3y 31907.8 標準回帰係数 β1'=β1*√(S11/Syy) = β2'=β2*√(S22/Syy) = β3'=β3*√(S33/Syy) = S1y S12 S13 β1= S2y S22 S23 /D = S3y S23 S33 β0= ym - (β1*x1 m + β2*x2 m + β3*x3 m ) = 予測値 Yi=β0 + β1*x1i + β 2*x2i + β3*x3i (Yi - ym) 2 (yi - Yi) 2 140.75 2986.5 95.09 184.65 115.5 152.50 215.49 403.6 42.38 164.40 960.9 31.34 173.82 465.6 33.89 156.35 1524.9 1.82 200.10 22.1 24.01 273.16 6047.2 17.33 194.69 0.5 127.95 136.62 3455.6 112.70 263.43 4628.3 5.91 252.87 3302.5 23.70 250.13 2995.5 0.02 144.86 2554.5 147.43 179.68 247.2 187.08 0.966 0.510 自由度調整済寄与率: R*2= = 1 - (Se/(n-4))/(Syy/(n-1)) = 0.958 7.62 > x1<- c(149,188,282,183,221,128,153,197,125,77,210,258, 254,114,174) > x2<- c(2.3,2.4,2.6,2.9,2.4,4.3,4.8,6.6,6.0,3.6,5.4,5.2,3.8,3.0,3.2) > x3 <- c(99,178,110,95,82,104,190,288,169,156,283,170,228,145,151) > y <- c(131,197,222,170,168,155,205,269,206,126,261,248,250,157,166) > result <- lm(y ~ x1 + x2 + x3) > summary(result) Call: lm(formula = y ~ x1 + x2 + x3) Residuals: Min 1Q Median 3Q Max -13.678 -5.344 -1.350 6.054 12.349 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.62493 11.10101 0.687 0.506388 x1 0.51007 0.04551 11.207 2.34e-07 *** x2 12.43483 2.64166 4.707 0.000643 *** x3 0.28814 0.05881 4.900 0.000472 *** --Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 9.55 on 11 degrees of freedom Multiple R-Squared: 0.9673, Adjusted R-squared: 0.9584 F-statistic: 108.6 on 3 and 11 DF, p-value: 1.869e-08 48 3.重回帰式の検討 実データ yi と理論値 yˆi との差は残差(または誤差) ei = yi - yˆi である。 残差が小さいほど分析精度が良いので残差平方和 S e = åe i 2 = å ( yi - yˆ i )2 を調べる。 (偏差平方和) S yy =(回帰平方和) S R +(残差平方和)S e S 寄与率(決定係数) R = 1 - e S yy 2 は大きいほうが 残差平方和 S e となるので は小さくなる。 そこで、S R とSe の大きさの比較を行う。そこで、それぞれの不偏分散(平均平方) VR ,Ve の比 V (分散比)FR = R を調べる。 FR がある値より大きい場合に有意として採択する(分散分析)。 Ve 説明変数を選択するときの注意 1.説明変数は、目的変数に効いている(相関係数の高い)変数だけをモデルに含ませる。 2.そのためには、適切な説明変数を選択することが重要で、基本的には、説明変数相互で高い 相関があるもの(0.9 以上)を探し、目的変数との相関係数が低い方の変数を削除する(多重共 線性を避けるため)。 3.または、説明変数の取捨選択のために、変数を1つずつ増やすか減らして自由度調整済寄与 率の大きな変数を選択するように試みるステップが必要となる。 変数増加法についてステップを以下に示す。 1.定数項だけのモデル MODEL0 から出発する。 MODEL0: yi = b 0 + e i に変数 x1 か、変数 x2 のどちらかの変数を取り込んだ単回帰式 MODEL1: yˆ i = bˆ0 + bˆ1 x1i または yˆ i = bˆ0 + bˆ2 x 2 i ・・・(30) において、 yの偏差平方和 Syy(自由度φT=n-1)と、残差平方和 Se(M1)(自由度φ e(M1)=n-2)を使う。 1変数の場合の残差平方和 S e ( M 1) は式(26)より、S e = S yy - S R = S yy - bˆ1- S1 y + bˆ2 S 2 y で 求まる。MODEL0 が正しいとき、 分散比: F0 = ( S yy - S e ( M 1) ) /(fT - fe ( M 1) ) S e ( M 1) / f e ( M 1) = ( S yy - S e ( M 1) )( n - 2) S e ( M 1) は自由度φ T-φe(M1)、 φe(M1)の F 分布に従うので、F0 が大きいとき(例えば2より大きい)場合に式(30)を採択する。 変数 x1 か、変数 x2 のいずれかを使って F0 を計算し、大きい方の変数のモデルを取り入れる。 いずれも大きくない場合は MODEL0 を採択する。 (注)変数 x1 、または変数 x2 の2つの場合でそれぞれの S e ( M 1) の値は異なる。 2.MODEL1 が採択されたとき、次に変数 x1 、変数 x2 の両方を取り入れた MODEL2 を考える。 MODEL2: yˆ = bˆ + bˆ x + bˆ x ・・・(31) i 0 1 1i 2 2i MODEL2 での残差平方和を Se(M2)(自由度をφe(M2)=n-3)として、 F0 = (S e ( M 1) - S e ( M 2) ) /(fe ( M 1) - fe ( M 2 ) ) Se ( M 2) / fe ( M 2 ) = (S e ( M 1) - S e ( M 2) )(n - 3) Se ( M 2) は自由度φe(M1)-φe(M2)、φe(M2)の F 分布に従うので、 F0 値が、前の MODEL1 の場合より大きいときに MODEL2 を採択する。 49 投資効果 投資金額 投資人員 (y) (x1) (x2) 54 72 2 45 63 5 44 53 4 33 57 16 30 51 16 32 38 4 3.単回帰とするときは、いずれの説明変数を使う方が 39 51 11 回帰式の精度が良いか。寄与率から判断しなさい。 60 73 1 54 69 5 45 77 22 【演習 20】(重回帰式) 右表はある投資に関するデータである。 1.これらのデータを元に重回帰式を求めて(投資効果が 目的変数である) 、その時の寄与率、及び、自由度調整済 寄与率を求めてみよう。 2.説明変数の投資金額と投資人員のいずれがより重要と いえるか。 4.重回帰式も含めてどの回帰式を使うのが最も良いかを 自由度調整済寄与率で判断して考察しなさい。 (Excel でやっても、R を用いても良い) (Excel の計算表が http://www.coins.tsukuba.ac.jp/~fukui/da/multiple_regression.xls にある) 50
© Copyright 2024 ExpyDoc