DD-IV-3reg1

3　メトリックな変数を別の変数で説明する（重回帰分析）

内容
1)分析の目的と手順
2)単回帰分析
３)重回帰分析
4)ノンメトリックな変数を説明変数として入れた回帰分析
５)回帰診断

1)分析の目的と手順

　例）マーケティングミックス：市場反応（認知率の増加、売上の増加など）を引き起こすことが目的。
　市場の反応とマーケティングミックスとを関連づけた分析が必要。

　基準変数（従属変数、被説明変数）→市場の反応（認知率、ブランドへの好意度、売上　など）
　独立変数（説明変数）　→マーケティングミックス（製品の属性、価格、広告コピー、露出量など）
→　基準変数（被説明変数）を説明変数で説明する。

重回帰分析の手順

理論、仮説の決定
　　↓
変数の決定　　　　　理論や仮説を検証するために妥当な変数を選択。
　　↓
予備分析　　　　　　説明変数を被説明変数との関係、説明変数間の関係、はずれ値の有無などをチェック。
　　↓
重回帰分析の実施　　プログラムを書き推定を行う。
　　↓
モデルの全体的なあてはまりのチェック　　R²やF検定の結果よりモデルの全体的なあてはまりをチェック。
　　↓
モデルの部分的なあてはまりのチェック　　個々の説明変数のt検定の結果を解釈
　　↓
変数の追加・削除　　説明していない変数を削除したり、説明しそうな変数を加える。
　　↓
回帰診断　　　　　　残差分析などにより、重回帰分析の仮定が成立しているか、はずれ値はないかなどをチェック。
　　↓
（必要ならば）モデルの改訂　　変数の組み合わせを変更したり、はずれ値を除外して再度推定。
　　↓
　結論

2)単回帰分析：演習用データの例

　開発者数と関係のある変数は？

　被説明変数：開発者数
　説明変数：？？？

・仮説例　
　ダウンロード数が増えると（ユーザーが増えて）開発者に参加する人も増える。
　バグ報告が増えると、開発者に参加する人も増える。

　まずは目で見て相関係数をみる。

plot(ossdata$nofdl,ossdata$npcntb)

cor(ossdata$nofdl,ossdata$npcntb)

#[1] 0.5172402

図　プロットの結果

plot(ossdata$nmsgbug,ossdata$npcntb)

cor(ossdata$nmsgbug,ossdata$npcntb)

#[1] 0.5172402

　それぞれ比例関係があるようにみえる（その他の変数についても同様にできるが、省略）。
　→データにあてはまるような直線をみつけるのが、(重)回帰分析。

・単回帰分析
　このグラフで、データにもっともよくあてはまるような「直線」をみつける。

　つまり、
　開発者数　=　β₁ダウンロード数　+　切片
　　　　 ↓
　　　y　　　　=　β₁x　　　　　　+　β₀　と定式化して、β₁、β₀推定することになる。

　β₁ ：偏回帰係数　partial regression coefficient　→直線の傾き
　β₀：切片 intercept→x= 0のときの yの値

lm( 非説明変数~説明変数,オブジェクト名）　　回帰分析（を含む線形モデルの推定を）行う。

例）y=β₀+β₁x₁+β₂x₂+β₃x₃+β₄x₄;の重回帰モデルを推定する場合には下のように。

lm (y1~x1+x2+ x3+ x4)
　　　　　　　　切片β0は指定しなくてもよい

　 result<-lm(y~x)
　のように指定するとオブジェクトの中に下記変数が格納される。

fitted 推定された係数とデータを用いて計算した従属変数の値。
　　　　　　y^{^} =β₀^{^}+β₁^{^}x₁+β₂^{^}x₂+β₃^{^}x₃+β₄^{^}x₄ residual 残差=y-y^
coefficient 推定された回帰係数

・単回帰分析の実際

#回帰分析の結果を　オブジェクト　result に入れる
result<-lm(ossdata$npcntb~ossdata$nofdl)

#参考　, のあとにオブジェクト名を指定すれば、モデルの指定式の中ではオブジェクト名は不要。
result<-lm(npcntb~nofdl,ossdata)

#結果を表示させる
summary(result)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 3.870e+00 6.231e-01 6.211 1.46e-08 ***
# 切片の　推定値　　　標準誤差　　t値　　p値　と有意水準
#ossdata$nofdl 1.716e-04 2.944e-05 5.828 8.01e-08 ***
#　　この変数の　推定値　　　標準誤差　　t値　　p値　と有意水準
#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
#
#Residual standard error: 5.358 on 93 degrees of freedom
#残差の自由度=サンプル数-推定したパラメーターの数
#　　　　　　=95-2
#Multiple R-Squared: 0.2675, Adjusted R-squared: 0.2597
# R²と修正R²
#F-statistic: 33.97 on 1 and 93 DF, p-value: 8.01e-08
# モデル全体についてのF検定の結果(帰無仮説すべてのβ=0)
#-----------------

#このように計算された結果オブジェクトには推定（結果）についての情報が格納されている。names関数で、オブジェクトに含まれている変数の名前を見ることができる。

names(result)

#いつもと同じように　オブジェクト名$変数名とすれば、それぞれの変数を利用できる。

#推定結果をあてはめた予測値：推定した回帰係数と説明変数を用いて計算した被説明変数の値
result$fitted.values

#残差=実際の値-fitted.values
result$residuals

#見やすいように必要な部分を別のデータセットに入れる
#　　　　　　　被説明変数　　　,説明変数　　予測値　　　　　　　　残差
ys<-data.frame(ossdata$npcntb,ossdata$nofdl,result$fitted.values,result$residuals)
ys
#residual =実際の値-fitted.valuesとなっていることがわかる。
# ossdata.npcntb ossdata.nofdl result.fitted.values result.residuals
#1 6 31496 9.273904 -3.27390415
#2 13 8907 5.398010 7.60198988
#3 3 126 3.891338 -0.89133761
#4 1 2524 4.302794 -3.30279424
#5 2 122 3.890651 -1.89065127

#分散分析表も出力してみる。
aov(result)

#Call:
# aov(formula = result)

#Terms:
# nofdl　と　残差Residual
# nofdl Residuals
#Sum of Squares 975.2163 2669.9416
#変数によって説明される部分と、残差の二乗和
#Deg. of Freedom 1 93
#自由度

#Residual standard error: 5.358083
#Estimated effects may be unbalanced

#実際の被説明変数と、説明変数をプロットしてみる
plot(ossdata$nofdl,ossdata$npcntb)

#推定された偏回帰係数を用いて直線を描く
#abline(切片,傾き）
#のように指定すると y=切片+傾き*xという直線を現在の散布図に加える。

result$coefficient
#ここに推定された係数が格納されている。
#(Intercept) nofdl
#3.8697181174 0.0001715832
#それぞれ次のようにして指定。
abline(result$coefficient[1],result$coefficient[2])

#実際の被説明変数と予測値をプロットしてみる
plot(ossdata$npcntb,result$fitted.values)

○出力される情報についての

・「決定係数 coefficient of determination　R²」
　モデルのあてはまりの良さを表す指標¹ 。

　　　

：データの（平均からの）分散
　　　

：回帰モデルによって説明されるデータの分散
　　　

：回帰モデルと実績値との乖離→回帰モデルによって説明されないデータの分散

　とすると、次式が成り立つ。

　R²：データの分散のうち、回帰モデルによって説明されている分散の割合

　R²は 0～1の間の値をとる（1に近いほどモデルがよくあてはまっている）。

・自由度修正済み決定係数 adjusted coefficient of determination
→変数を増やすと、それが有意に被説明変数を有意に説明していなくても、一貫してR²は向上。
　決定係数は説明変数の数を増やすと見かけ上向上してしまう。

　統計学では、parsimonyの原理（節約の原理=おなじ程度の当てはまりならば、より少ない変数、簡単なモデルの方がよい）を採用しているので、使った変数の数(自由度　degree of freedom)を考慮して決定係数を修正したもの。
　　n:サンプル数　　k：パラメーターの数とすると次式で与えられる（パラメーターの数=1のときには決定係数と自由度修正済み決定係数は一致する。）。

・F検定
　複数の説明変数がある時に意味がある（一つのときには偏回帰係数のt検定で充分）。

　帰無仮説β1=β2=...β_k-1=0　（どの説明変数も従属変数を説明していない。）　を検定

　k:推定されたパラメーターの数（定数項を含む）
　n：サンプル数
　n-k：残差の自由度
　モデルによって説明される分散：
　モデルによって説明されない分散（残差の分散）：

　とすると、次式のFは自由度 k-1、n-kのF分布に従う。

　この値をＦ（kー1,nーk)と比較。αを有意水準として、

　　　F≧Ｆ_α（kー1,nーk)のときには、帰無仮説β1=β2=...β_k-1=0を棄却する。

　　→少なくともどれかは0ではない。少なくとも１つの変数が従属変数を有意に説明している。

・t検定：推定されたパラメーター=β_dという帰無仮説を検定。
　帰無仮説：=β_d　次の統計量 tは、自由度 n-k-1のt分布に従う（n：サンプル数。　k：パラメータの数）。

　　：パラメータβの標準誤差。
　　：パラメータβの推定値。

　　　　|t|>t_α/2(n-k-1)のとき仮説は棄却される。
　　　　|t|<t_α/2(n-k-1)のとき仮説は棄却されない。
→回帰分析においてもっとも重要なのは、説明変数が従属変数を有意に説明しているか否か（推定された偏回帰係数が0か否か→β₀=0）：次式で t値を算出² 。

3）重回帰分析
　単回帰分析では、一つの変数によって説明したが、重回帰分析では複数の変数によって説明する。
　つまり、

　開発者数　=　β₁ダウンロード数　+　β₂バグ報告の数　+切片
　　　　 ↓
　　　y　　　　=　β₁x₁　　　　+β₂x₂　　　　　　+　β₀

#2変数の場合は右辺に　+　で二つの変数を追加すればよい
result2<-lm(npcntb~nofdl+nmsgbug,ossdata)
summary(result2)

#-----------------
#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 3.640e+00 4.986e-01 7.300 9.99e-11 ***
#ossdata$nofdl 8.954e-05 2.604e-05 3.439 0.00088 ***
#ossdata$nmsgbug 2.774e-02 3.782e-03 7.333 8.55e-11 ***
#---
#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

#Residual standard error: 4.28 on 92 degrees of freedom
#Multiple R-Squared: 0.5377, Adjusted R-squared: 0.5277
#F-statistic: 53.51 on 2 and 92 DF, p-value: 3.841e-16

aov(result2)

#Call:
# aov(formula = result2)

#Terms:
# nofdl nmsgbug Residuals
#Sum of Squares 975.2163 984.9564 1684.9851
#Deg. of Freedom 1 1 92

#Residual standard error: 4.279609
#Estimated effects may be unbalanced

#1変数の場合と同様にプロットしてみる。ただし説明変数が二つあるので２枚のグラフ。
#実際の被説明変数と、説明変数をプロットしてみる
plot(ossdata$nofdl,ossdata$npcntb)

#推定された偏回帰係数を用いて直線を描く
result2$coefficient
abline(result2$coefficient[1],result2$coefficient[2])

#もう一つの説明変数についても同様
plot(ossdata$nmsgbug,ossdata$npcntb)
abline(result2$coefficient[1],result2$coefficient[3])

#実際の被説明変数と予測値をプロットしてみる
plot(ossdata$npcntb,result2$fitted.values)

○F検定による、変数を加えたことの効果の検定
　変数を追加したことによって、有意にモデルの説明力が上がったか？を下記によって検定できる。
　ただし、後述のAICを用いる方が容易。

　もとになるモデルの　残差平方和を S_{e1 。残差の}自由度を　df_e1
　変数を追加したモデルの残差平方和を S_{e2 。残差の}自由度を　df_e2
　　変数を追加した事による残差の減少　　ΔS_e=S_e1-S_e2
　　変数を追加した事による自由度の減少　Δdf_e=df_e1-df_e2

　とする。ΔS_eが大きければ変数を追加することによって説明力が高くなった（残差が減少した）と結論図蹴ることができる。
　厳密には、次のFを求めて、F分布表のパーセント点　F(ΔS_e,S_e1;α)と比較する。

#上の例では
#変数が一つの場合
# nofdl Residuals
#Sum of Squares 975.2163 2669.9416
#Deg. of Freedom 1 93

#残差平方和は　2669.9416
#残差の自由度は　93

#変数を追加した場合
#Terms:
# nofdl nmsgbug Residuals
#Sum of Squares 975.2163 984.9564 1684.9851
#Deg. of Freedom 1 1 92
#残差平方和は　1684.9851
#残差の自由度は　92

#Fを計算
delta.Se<- 2669.9416-1684.9851
delta.dfe<- 93-92

F<-delta.Se/delta.dfe/(2669.9416/93)
F

#[1] 34.30822

pf(F,delta.dfe,93,lower.tail=FALSE)
#pf(F値,自由度1,自由度2):F分布のパーセント点を求める関数。
#7.051216e-08
#p値
#ここでの帰無仮説
#2つのモデルで残差は減少していない　が棄却される。

#もう一つ変数を入れてみる
result3<-lm(npcntb~nofdl+nmsgbug+nrelease,ossdata)
summary(result3)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 4.0290647 0.5902670 6.826 9.51e-10 ***
#nofdl 0.0001013 0.0000277 3.658 0.000425 ***
#nmsgbug 0.0270758 0.0038111 7.105 2.62e-10 ***
#nrelease -0.0480682 0.0393045 -1.223 0.224500
#---
#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

#Residual standard error: 4.268 on 91 degrees of freedom
#Multiple R-Squared: 0.5452, Adjusted R-squared: 0.5302
#F-statistic: 36.37 on 3 and 91 DF, p-value: 1.538e-15

4)ノンメトリックな変数を説明変数として入れた回帰分析

　基本的にダミー変数（dummy variable: 0もしくは1の値となる変数）を定義することによって、ノンメトリックな変数を説明変数として入れることが可能。

(a)シングルアンサーの変数（一つしか選べない場合）

　例1.1　性別ダミー
　　　　d=0　男性の時
　　　　d=1　女性の時

　例1.2　OSダミー
　 oslinux　　　そのプロジェクトの開発しているソフトがLinux用だと1。そうでないとき0
　osmicros　　そのプロジェクトの開発しているソフトがMicrosoftのOS(MS-DOS, WINxx)用だと1。そうでないとき0
　 ososind 　　そのプロジェクトの開発しているソフトがOSに依存しない場合1。そうでないとき0
　などはそのまま説明変数として利用可能。

　例2.1　職業ダミー
　　各回答者とも「会社員、学生、自由業、無職、その他」のうち、どれか一つを選んでいる場合。
この場合も次のようにダミー変数を定義できる。

　　会社員ダミー dkaisya=1 会社員のとき
　　 dkaisya=0 会社員でないとき
　　学生ダミー dgakusei=1 学生のとき
　　 dgakusei=0 学生でないとき
　　自由業ダミー djiyuu=1 自由業のとき
　　 djiyuu=0 自由業でないとき
　　無職ダミー dmusyoku=1 無職のとき
　　 dmusyoku=0 無職でないとき
　　その他ダミー dsonota=1 その他のとき
　　 dsonota=0 その他でないとき

　ただし、回帰分析にもちいるときは、これら5つのダミー変数のうち、４つのみを用いる。
　選択肢をどれか一つしか選べないのだから、
　例えば、会社員ダミー=0、学生ダミー=0、自由業ダミー=0、無職ダミー=0のときには、「その他」を選んだはず。ということがいえる。

　例2.2　 devstage 　　開発段階を示す変数
　　　　　1:計画段階　　2:α版以前　3:α版(試作版)のソフト公開　
　　　　　4:β版(正式版の一歩手前)のソフト公開　5:安定版のソフト公開　6:成熟
　6つの値をとるので、5つのダミー変数を定義する。

　注）どの変数を分析に用いるか？は、分析の目的による。
　　　通常は、「その他」のようなものがある場合には、それ以外の変数を用いる。

　このようにシングルアンサー変数の場合には、その変数の水準の数-1個のダミー変数を定義して回帰分析の説明変数に用いる。
　　　「水準」　変数がとることができる値が何通りあるか？
　　　　　　　　　変数が　0、１、２の値をとるときは　３水準
　　　　　　　　　　　　　0、１の値をとるときは２水準。

○ダミー変数の入れ方

　回帰分析のパラメーターは「傾き」と「切片」。それぞれにダミー変数を入れる、入れないで次の４通りのモデルがあり得る。

	切片	傾き
1	ダミー変数を入れない	ダミー変数を入れない	グループ間でパラメーターに違いがない。 y　　　　=　β*x　 + β₀
2	ダミー変数を入れる	ダミー変数を入れない	傾きは同じだが、切片が異なる。 y　　　　=　β₁x　+ β₀ + β_{d0 }dum
3	ダミー変数を入れない	ダミー変数を入れる	切片は同じだが、傾きが異なる。 y　　　　=　(β₁+β_{1d }dum)x　+ β₀
4	ダミー変数を入れる	ダミー変数を入れる	切片、傾きとも異なる。 y　　　　=　(β₁+β_{1d }dum)x　+ β₀ + β_{d0 *}dum

　切片

　ベースとなるモデル

モデル1)切片、傾きともダミー変数に関わらず同じ

　　　y　　　　=　β*x　 + β₀

モデル2）傾きは同じだが、切片が異なる。
　　　y　　　　=　β₁*x　+ β₀ + β_{d0 *}dum
　dum=0のとき
　　　y_dum=1　　　　=　β₁*x　+ β₀ + β_d0
　dum=0のとき
　　　y_dum=0　　　　=　β₁*x　+ β₀

モデル3）切片は同じだが、傾きが異なる。
　　　y　　　　=　(β₁+β_{1d *}dum)*x　+ β₀
　dum=0のとき
　　　y_dum=1　　　　=　(β+β_d)*x　+ β₀
　dum=0のとき
　　　y_dum=0　　　　=　β*x　+ β₀

　モデル4）切片、傾きとも、ダミー変数によって異なる（ダミー変数=1と=0のグループによって異なる）
　切片、傾きに、ダミー変数を入れたモデル

　　　y　　　　=　(β₁+β_{1d *}dum)*x　+ β₀ + β_{d0 *}dum
ただし、
dum:ダミー変数　　0もしくは1
x：メトリックな変数

　dum=0のとき
　　　y_dum=1　　　　=　(β+β_d)*x　+ β₀ + β_d0
　dum=0のとき
　　　y_dum=0　　　　=　β*x　+ β₀
　のようになる。
　つまり、傾きも、切片も異なると考えている。
　β_dダミー変数=0のグループを基準としたときの傾きの変化分β_d0ダミー変数=0のグループを基準としたときの切片の変化分

・上のどれを用いるか？

　(1)仮説などにもとづいて、決定する場合
　　例　ダウンロード数が0の場合は、OSがlinuxであっても、なくても、バグの報告は0
　　　　この場合、切片は共通の0となるはず。
　　　　　よって、すくなくとも切片についてはダミー変数を入れない。

　(2)あてはまりの良さから決定する場合
　　上の1-4のモデル、それぞれについて推定。あてはまりの良さを比較。最もあてはまりがよいモデルを選ぶ。

○実例
　以下では
　　ダウンロード回数nofdlを
　　　　リリース回数　nrelease
　　　ソフトウエアのトピックがゲームか否かダミー　　tpcgames　
　によって説明。４つのモデルを推定して、あてはまりの良さを比較する。

モデル　１　ダミー変数を用いずにメトリックな変数のみで説明。

result1<-lm(nofdl~nrelease,ossdata)
summary(result1)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 4914.1 2397.8 2.049 0.04324 *
#nrelease 504.2 154.3 3.267 0.00152 **
#---
#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

#Residual standard error: 17870 on 93 degrees of freedom
#Multiple R-Squared: 0.103, Adjusted R-squared: 0.09333
#F-statistic: 10.68 on 1 and 93 DF, p-value: 0.001521

#AICの算出（後述）

AIC(result1)
#[1] 2133.893

モデル0　ダミー変数のみで説明してみる。

result0<-lm(nofdl~tpcgames,ossdata)
summary(result0)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 7745 2010 3.854 0.000214 ***
#tpcgames 15040 5235 2.873 0.005040 **
#---
#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

#Residual standard error: 18090 on 93 degrees of freedom
#Multiple R-Squared: 0.08151, Adjusted R-squared: 0.07163
#F-statistic: 8.253 on 1 and 93 DF, p-value: 0.00504
#0という仮説は棄却できるので説明変数として用いても意味がある。

AIC(result0)
#[1] 2136.139

・モデル2（傾きは同じだが、切片が異なる）　
　次のように、ダミー変数を説明変数として入れるだけ。

result2<-lm(nofdl~nrelease+tpcgames,ossdata)
summary(result2)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 3058.4 2402.9 1.273 0.20631
#nrelease 481.2 149.0 3.229 0.00172 **
#tpcgames 14158.7 4995.9 2.834 0.00565 **
#---
#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

# nofdl = 481.2 nrelease + 14158.7* tpcgames +3058.4
#が推定された直線。つまり、次の二つの直線を表す。
#tpcgames=1のとき
# nofdl = 481.2 nrelease + 14158.7 +3058.4
#tpcgames=0のとき
# nofdl = 481.2 nrelease +3058.4

#Residual standard error: 17230 on 92 degrees of freedom
#Multiple R-Squared: 0.175, Adjusted R-squared: 0.1571
#F-statistic: 9.758 on 2 and 92 DF, p-value: 0.0001435

#生のデータと推定された二つの直線を描く。
plot(ossdata$nrelease,ossdata$nofdl)
abline(result2$coef[1],result2$coef[2])
abline(result2$coef[1]+result2$coef[3],result2$coef[2])

AIC(result2)
#2127.941

・モデル3　切片は同じだが、傾きが異なる。
　　　y　　　　=　(β₁+β_{1d *}dum)*x　+ β₀
　説明変数に　メトリックな変数*ダミー変数として生成した新しい変数も入れる。

　

#二つの変数の積をとった新しい変数 relgame を作成
ossdata$relgame<-ossdata$nrelease*ossdata$tpcgames

#nreleaseと作成した変数 relgame によって説明
result3<-lm(nofdl~nrelease+relgame,ossdata)
summary(result3)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 4299.0 2338.6 1.838 0.06925 .
#nrelease 413.9 153.7 2.693 0.00842 **
#relgame 890.9 342.3 2.603 0.01077 *
---
#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

#Residual standard error: 17340 on 92 degrees of freedom
#Multiple R-Squared: 0.1645, Adjusted R-squared: 0.1463
#F-statistic: 9.057 on 2 and 92 DF, p-value: 0.0002566

#上のようにわざわざ新変数を生成しなくても、変数1:変数2 のように、二つの変数をコロン(:)で結んだものを説明変数として指定すれば、上と同じ結果が得られる。
result3<-lm(nofdl~nrelease+nrelease:tpcgames,ossdata)
summary(result3)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 4299.0 2338.6 1.838 0.06925 .
#nrelease 413.9 153.7 2.693 0.00842 **
#nrelease:tpcgames 890.9 342.3 2.603 0.01077 *
#---
#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

#Residual standard error: 17340 on 92 degrees of freedom
#Multiple R-Squared: 0.1645, Adjusted R-squared: 0.1463
#F-statistic: 9.057 on 2 and 92 DF, p-value: 0.0002566

AIC(result3)
#2129.142

・モデル4(傾き、切片とも異なる）について
　　　y　　　　=　(β₁+β_{1d *}dum)*x　+ β₀ + β_{d0 *}dum　　　　　　　=　β₁*x+β_{1d *}dum*x　+ β₀ + β_{d0 *}dum　となるので、新しい変数　dum*xを作成。これと、ダミー変数そのものを説明変数に入れればよい。

result4<-lm(nofdl~nrelease+tpcgames+nrelease:tpcgames,ossdata)
summary(result4)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 3346.9 2449.7 1.366 0.17523
#nrelease 451.5 156.1 2.893 0.00477 **
#tpcgames 10072.6 7967.8 1.264 0.20940
#nrelease:tpcgames 357.8 542.4 0.660 0.51116

#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

#Residual standard error: 17290 on 91 degrees of freedom
#Multiple R-Squared: 0.1789, Adjusted R-squared: 0.1519
#F-statistic: 6.61 on 3 and 91 DF, p-value: 0.0004323

# nofdl = 451.5 nrelease +357.8 nrelease*tpcgames + 10072.6* tpcgames +3346.9

#が推定された直線。つまり、次の二つの直線を表す。
#tpcgames=1のとき
# nofdl = (451.5 +357.8 )nrelease + 10072.6 +3346.9
#tpcgames=0のとき
# nofdl = (451.5 )nrelease +3346.9

#生のデータと推定された二つの直線を描く。
plot(ossdata$nrelease,ossdata$nofdl)
abline(result4$coef[1],result4$coef[2])
abline(result4$coef[1]+result4$coef[3],result4$coef[2]+result4$coef[4])

AIC(result4)
#2129.488

推定のまとめ表

	モデル１	モデル0	モデル2	モデル3	モデル4
nrelease	504.2** (3.267)		481.2** (3.229)	413.9** (2.693)	451.5** (2.893)
傾きへのダミーnrelease*tpccntb				890.9* (2.603)	357.8 (0.660)
切片へのダミーtpccntb		15040** (2.873)	14158.7** (2.834)		10072.6 (1.264)
切片	4914.1* （2.049）	7745*** (3.854)	3058.4 (1.273)	4299.0 (1.838)	3346.9 (1.366)
R²	0.103	0.08151	0.175	0.1645	0.1789
修正R²	0.09333	0.07163	0.1571	0.1463	0.1519
残差平方和	29713596953	30424668920	27327764555	27675351736	27197720513
残差自由度	93	93	92	92	91
AIC	2133.893	2136.139	2127.941	2129.142	2129.488

○AIC(Akaike's Information Criterion: 赤池の情報量基準）

　誤差項が正規分布するという回帰モデルの場合、AICが小さいほどモデルがよいといえる。

AIC=n*log(RSS/n)+2*k

RSS 残差平方和
k 推定したパラメーター数
n サンプル数

　よって、複数のモデルがある場合、AICが最小となるものが、最もあてはまりがよいといえる。この場合はモデル2

参考）ノンメトリックな変数の値に応じてダミー変数を作成する。

この場合のダミー変数はノンメトリックな変数のとる水準の数-1個定義する。
　どの水準にダミー変数を割り当てるかは、分析の目的による。注目する水準には、ダミー変数を割り当てなければ、その効果も計算されない。

　選択肢が設けてあっても、実際にはそれを選んだ回答者がいないこともあるので、あらかじめその変数を単純集計しておくとよい。
　

table(ossdata$devstage)
# 2 3 4 5 6
# 8 16 25 42 4
#この変数は本来は1,2,3,4,5,6の値（6水準）をとるが、サンプルには、これが１になるものがない。
#よって5水準-1=4つのダミー変数をとる。
#今回は、2を基準とすることにして、次の４つのダミー変数を定義

#devst3 devstage=3のとき1　そうでないとき0
#devst4 devstage=4のとき1　そうでないとき0
#devst5 devstage=5のとき1　そうでないとき0
#devst6 devstage=6のとき1　そうでないとき0

#matrix(値,行の数,列の数)
#devstageに対応するダミー変数を定義して0を入れておく。

ossdata$devst3<-matrix(0,length(ossdata$devstage),1)
#ossdata$devstage==3のサンプルのみ、ダミー変数を1にする。
ossdata$devst3[ossdata$devstage==3]<-1

ossdata$devst4<-matrix(0,length(ossdata$devstage),1)
#ossdata$devstage==4のサンプルのみ、ダミー変数を1にする。
ossdata$devst4[ossdata$devstage==4]<-1

ossdata$devst5<-matrix(0,length(ossdata$devstage),1)
#ossdata$devstage==5のサンプルのみ、ダミー変数を1にする。
ossdata$devst5[ossdata$devstage==5]<-1

ossdata$devst6<-matrix(0,length(ossdata$devstage),1)
#ossdata$devstage==6のサンプルのみ、ダミー変数を1にする。
ossdata$devst6[ossdata$devstage==6]<-1

#出力して確かめてみる
a<-data.frame(ossdata$devstage,ossdata$devst3,ossdata$devst4,ossdata$devst5,ossdata$devst6)
a

#定義した変数を切片に入れて回帰。

res<-lm(nofdl~nrelease+devst3+devst4+devst5+devst6,ossdata)
summary(res)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -631.0 6172.0 -0.102 0.91879
#nrelease 465.2 158.1 2.943 0.00414 **
#devst3 1626.1 7536.0 0.216 0.82966
#devst4 2799.0 7276.8 0.385 0.70142
#devst5 8630.0 6814.4 1.266 0.20866
#devst6 26367.0 10723.3 2.459 0.01587 *

#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

#Residual standard error: 17400 on 89 degrees of freedom
#Multiple R-Squared: 0.1869, Adjusted R-squared: 0.1412
#F-statistic: 4.091 on 5 and 89 DF, p-value: 0.002172

AIC(res)
#2132.564

#定義した変数を傾きに入れて回帰。

res2<-lm(nofdl~nrelease+devst3:nrelease+devst4:nrelease+devst5:nrelease+devst6:nrelease,ossdata)
summary(res2)

#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 4664.6 2501.4 1.865 0.0655 .
#nrelease -781.5 1632.0 -0.479 0.6332
#nrelease:devst3 810.8 1704.4 0.476 0.6354
#nrelease:devst4 1122.1 1616.2 0.694 0.4893
#nrelease:devst5 1394.3 1606.1 0.868 0.3877
#nrelease:devst6 3594.8 1711.9 2.100 0.0386 *
#---
#Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

#Residual standard error: 16860 on 89 degrees of freedom
#Multiple R-Squared: 0.2359, Adjusted R-squared: 0.193
#F-statistic: 5.497 on 5 and 89 DF, p-value: 0.0001865

AIC(res2)
#2126.650

#AICはこちらの方が小さいが、nreleaseの符号が（有意ではないが）マイナスなので多重共線性（後述）の恐れあり。

¹ R=√R²は、重相関係数と呼ばれる。