東京国際大学大学院

心理データ解析２（後期）

東京大学大学院農学生命科学研究科　大森宏

４．回帰分析

相関

　標本（サンプル）に対し，２つの変数 x，y が測定されているとする．たとえば，x が身長(m)であり，y が体重(kg)である．大きさ n の標本（サンプル）に対し，２つの変数の組のデータが，

（x₁ ，y₁ ），（x₂ ，y₂ ）， …，（x_n ，y_n ）

であったとする．変数間の関連性の強さを測る量として共分散（Covariance），Cov[x ，y ] がある．これは，変数に対する平均を，

として，

と定義される．

　共分散は測定単位により大きさが変わるので，これをおのおのの変数の標本分散， Var[x]，Var[y]，

で標準化したものがピアソン（Peason）の積率相関係数 r，-1≦r≦1，であり，

と定義される．これは，変数間の線形的関係の強さ，（x が大きいと y も大きく，x が小さいと y も小さい，）を測る指標で，｜r｜＝1 のときは，変数 x，y は完全な直線関係にあり，r ＝0 のときは，線形的な関係がない．r が 1 に近いときは，正の相関関係があるといい， r が -1 に近いときは，負の相関関係があるという．

データ散布図と相関係数
sokan sokan

sokan sokan

相関データダウンロード

# 相関係数の R スクリプト
r <- read.csv("r1.csv") # csv データ読み込み
x <- r[,1]; y <- r[,2] #
plot(x,y, xlim=c(2,8), ylim=c(0,6)) # x，y の散布図
title(main="相関係数 r = 0.23") #
cov(x, y) # x と y の共分散
cor(x, y) # x と y　の相関係数
x <- 10*x #
plot(x,y, xlim=c(2,8), ylim=c(0,6)) # x，y の散布図
cov(x, y) # x と y の共分散
cor(x, y) # x と y　の相関係数

**# 相関係数の R スクリプト**
r <- read.csv("r1.csv")	# csv データ読み込み
x <- r[,1]; y <- r[,2]	#
plot(x,y, xlim=c(2,8), ylim=c(0,6))	# x，y の散布図
title(main="相関係数 r = 0.23")	#
cov(x, y)	# x と y の共分散
cor(x, y)	# x と y　の相関係数
x <- 10*x	#
plot(x,y, xlim=c(2,8), ylim=c(0,6))	# x，y の散布図
cov(x, y)	# x と y の共分散
cor(x, y)	# x と y　の相関係数

入試得点とその後の成績
　入試得点とその後の成績には相関が強いのではないかと考えられるが，実際に測ってみるとそれほどでもないことが多い．このことから，入試を行うことを疑問視する人も出てくる（入試得点と成績に相関がないのだから，入試を行うことに意味がない）．しかしながら入試では，多くの学生が入学できないので，本来入学していたら取れたであろう成績のデータが欠測していると考えるべきである．
　統計学的には以下のように考えればよい．いま，入試得点（x）とその後の成績（y）に r の相関があるとする．そして，x と y は相関 r を持つ 2 変量正規分布に従うとする．しかし現実には，入試得点の低い者（x ＜ x₀）は入学できないので，この 2 変量正規分布を x = x₀ で切断（truncate）した分布のみが観測される．
　いま，入試得点（x）と成績（y）がそれぞれ，平均 50，分散 100（標準偏差 10），相関 0.7 を持つ 2 変量正規分布とし，合格ラインは x = 54 点であるとする．すなわち，x ≧ 54 である 330 名程度が合格する．合格者内での入試得点と成績の相関は，約 0.47 であり，本来あるべき相関 0.7 より低下して，見かけ上相関がそれほど強くないようにみえてしまう．

# 入試得点と成績の R スクリプト
library(MASS) # MASS ライブラリーの読み込み
s <- matrix(c(100, 70, 70, 100), nrow=2) # 分散共分散行列の定義
x <- mvrnorm(1000,c(50,50),s) # 平均 50，分散 s の 2 変量正規分布乱数 1000 個生成
plot(x, xlab="入試得点", ylab="初年度成績") # 分布全体を表示
points(x[x[,1]>54,], col="red") # 入学者のみ赤で表示
abline(v=54) # 合格ライン
title(main="本来の相関と観測される相関") #
cor(x)[1,2] # 本来あるべき相関
cor(x[x[,1]>54,])[1,2] # 実際に観測される相関

**# 入試得点と成績の R スクリプト**
library(MASS)	# MASS ライブラリーの読み込み
s <- matrix(c(100, 70, 70, 100), nrow=2)	# 分散共分散行列の定義
x <- mvrnorm(1000,c(50,50),s)	# 平均 50，分散 s の 2 変量正規分布乱数 1000 個生成
plot(x, xlab="入試得点", ylab="初年度成績")	# 分布全体を表示
points(x[x[,1]>54,], col="red")	# 入学者のみ赤で表示
abline(v=54)	# 合格ライン
title(main="本来の相関と観測される相関")	#
cor(x)[1,2]	# 本来あるべき相関
cor(x[x[,1]>54,])[1,2]	# 実際に観測される相関

直線回帰

モデル

　２つの変数 x ，y に対し，y の値が x の値の動きにつれて線形的に変化すると仮定される，つまり，

y ＝ a ＋ b x

という関係が成り立っていると考えられる場合である．これを y の x に対する直線回帰といい，a ，b を回帰係数という．また，変数 y を従属変数，目的変数といい，変数 x を独立変数，説明変数という．

最小２乗法

　データに最もよくあてはまる直線回帰式を得るには，データ点（x_i ，y_i ），と回帰による推定点，（x_i ，y^_i ）， y^_i ＝ a ＋ b x_i ，の間の距離の２乗和 S が最小になるような回帰係数 a ，b を求める．つまり，

を最小化する a ，b を求める問題に帰着する．これを最小２乗法という．

　これは，S を a ，b で偏微分して 0 とおくことによって得られる．つまり，

の連立方程式を a ，b で解けばよい．これより，

が得られる．

入試得点と初年度成績
　下の表はアメリカのあるビジネススクール（MBA）での女子学生の入試得点と初年度成績である．

入試得点（x） 680 500 600 420 480 630 550 590 610 500 640 570 610

初年度成績（y） 332 265 309 253 276 326 299 310 324 327 334 301 336

入試得点（x）	680	500	600	420	480	630	550	590	610	500	640	570	610
初年度成績（y）	332	265	309	253	276	326	299	310	324	327	334	301	336

　入試得点を説明変数（x），初年度成績を目的変数（従属変数）（y）として回帰式を求めてみる．
　なお一般に，入試ではある得点以下の学生は入学していないので，入試得点と初年度成績の関係は弱くなってしまうが，ここではそのことを無視して話しを進める．

データダウンロード

# 入試得点と初年度成績回帰式の R スクリプト
mba <- read.csv("mbagrade2.csv") # データ読み込み　
mbaf <- mba[mba[,1]=="F",] # 女性データのみ抽出　
mbaf.lm <- lm(初年度成績 ~ 入試得点, data=mbaf) # 回帰：lm(y ~ x, data=zzz), y = ax + b 　
summary(mbaf.lm) # 結果表示　
plot(mbaf[,2:3], xlim=c(400,700), ylim=c(2.5,3.5)) # データ散布図　
abline(mbaf.lm, col="blue") # 回帰直線　
points(mbaf[,2], mbaf.lm$fitted.value, pch=19, col="red") # 回帰推定値　
segments(mbaf[,2],mbaf[,3],mbaf[,2], mbaf.lm$fitted.value) # 回帰残差　
title(main="回帰直線と回帰残差") # 　

legend(locator(1), legend=c("データ","回帰推定値"), pch=c(1,19), col=c("black","red"))

**# 入試得点と初年度成績回帰式の R スクリプト**
mba <- read.csv("mbagrade2.csv")	# データ読み込み
mbaf <- mba[mba[,1]=="F",]	# 女性データのみ抽出
mbaf.lm <- lm(初年度成績 ~ 入試得点, data=mbaf)	# 回帰：lm(y ~ x, data=zzz), y = ax + b
summary(mbaf.lm)	# 結果表示
plot(mbaf[,2:3], xlim=c(400,700), ylim=c(2.5,3.5))	# データ散布図
abline(mbaf.lm, col="blue")	# 回帰直線
points(mbaf[,2], mbaf.lm$fitted.value, pch=19, col="red")	# 回帰推定値
segments(mbaf[,2],mbaf[,3],mbaf[,2], mbaf.lm$fitted.value)	# 回帰残差
title(main="回帰直線と回帰残差")	#
legend(locator(1), legend=c("データ","回帰推定値"), pch=c(1,19), col=c("black","red"))

回帰式の統計モデル

　推定された直線回帰式がどの程度現実のデータに適合しているかを調べるために，回帰式が従う統計モデルを考える．標本の格データ点，（x_i ，y_i ），が，

y_i ＝ a ＋ b x_i ＋ e_i ， e_i ～ N( 0，σ² )

であると仮定する．e_i は誤差（error），あるいは，残差（residual）で，直線回帰式では説明がつかない部分を表し，これが互いに独立に平均 0，分散 σ² の正規分布に従うと仮定する．誤差の大きさが大きいときは，直線回帰式ではデータが説明できないと考える．

残差分散と回帰係数の標準誤差

　回帰で説明がつかない残差平方和 S_e は，

で求められる．これの自由度は n－2 であるので（２つの回帰係数分の自由度を除く），回帰の残差（誤差）分散は，

で求められる．

　一般に，Var(y_i ) ＝ σ² であるとき，その定数倍の分散は，

Var(ay_i ) ＝ a²σ²， Var(Σ_ia_i y_i ) ＝ Σ_ia_i ² σ²

であり，従属変数 y のデータ y_i は，

y_i ～ N( a ＋ b x_i ，σ² )

と分布するので，回帰係数 b の分散は，

となる．この分散の平方根を回帰係数 b の標準誤差という．

回帰係数の標準誤差による t 検定

目的変数 y が説明変数 x との回帰関係にないという帰無仮説，

H₀：b ＝ 0，

を考えてみよう．

回帰係数 b の推定値 b^ の分散は，

と推定できるので，b^ の標準偏差（標準誤差）は， s _b と推定される．これより，回帰係数をその標準誤差で割った t 値が，帰無仮説のもとで，

のように，自由度 n－2 の t 分布に従うことを利用して回帰係数の検定が行える．すなわち，自由度 n－2 の t 分布の 97.5％点を t₀ とすると，

|t | ＞ t₀ → 帰無仮説を有意水準 5 ％で棄却（回帰関係が有意に認められる）

|t | ≦ t₀ → 帰無仮説を棄却しない（回帰関係が認められない）

と定式化できる．

分散分析

平方和分解

　回帰式により，従属変数 y のデータ y_i は，

y_i ＝ y^_i ＋ (y^_i － y_i ) ＝回帰値＋残差

のように分解される．この分解に対応して従属変数データの総平方和 S_T は，

S_T ＝ Σ_i （y _i － y^- ) ² ＝ Σ_i （y^_i － y^- ) ² ＋ Σ_i （y _i － y^_i ) ² ＝ S_R ＋ S_e

総平方和＝回帰平方和＋残差平方和

のように分解される．これを平方和の分解という．この分解に対応して自由度は，

n－1 ＝ 1 ＋ n－2

と分解される．

決定係数（重相関係数の２乗）

　データが直線回帰式でよく説明できるのは，回帰平方和が大きく，残差平方和が小さい場合である．総平方和のうち回帰平方和で説明される割合を決定係数，もしくは重相関係数の２乗といい，

で定義される．なお，重相関係数 R とは，データ y _i と回帰値 y^_i との間の相関係数である．これより，以下の分散分析表ができる．

回帰分析の分散分析表
変動因平方和自由度平均平方 F 値

回帰 S_R 1 S_R F ＝ S_R/s_e²

残差 S_e n－2 s_e² ＝ S_e/n－2 　

全体 S_T n－1 　　

回帰分析の分散分析表
変動因	平方和	自由度	平均平方	F 値
回帰	S_R	1	S_R	F ＝ S_R/s_e²
残差	S_e	n－2	s_e² ＝ S_e/n－2
全体	S_T	n－1

F 検定

従属変数 y が説明変数 x の回帰関係にないという帰無仮説，

H₀：b ＝ 0，

を考える．帰無仮説のもとでは，回帰平均平方 S_R と残差分散 s_e² がともに誤差 σ² の不偏推定量になるので，その比 F 値が，

F ＝ S_R/s_e² ～ F（1，n－2），

という F 分布に従うことを利用して検定ができる．すなわち，分子，分母自由度が 1，n－2 である F 分布 F（1，n－2）の95％点を F₀ とすると，

F ＞ F₀ → 帰無仮説を有意水準 5 ％で棄却（回帰関係が有意に認められる）

F ≦ F₀ → 帰無仮説を棄却しない（回帰関係が認められない）

と定式化できる．

# 回帰分析の詳細の R スクリプト
anova(mbaf.lm) # 分散分析表
n <- nrow(mbaf); n # データ数
st <- var(mbaf[,3])*(n-1); st #　総平方和
sr <- var(mbaf.lm$fitted.values)*(n-1); sr #　回帰平方和
se <- var(mbaf.lm$residuals)*(n-1); se #　残差平方和
s <- sqrt(se/(n-2)); s #　残差標準偏差
sx <- var(mbaf[,2])*(n-1); sx #　x の総平方和
b <- cov(mbaf[,2],mbaf[,3])*(n-1)/sx; b #　回帰係数
sb <- s/sqrt(sx); sb #　回帰係数bの標準誤差
r2 <- sr/st; r2 # 重相関係数の２乗
r <- cor(mbaf[,3], mbaf.lm$fitted.values); r # 目的変数と回帰推定値の相関
fv <- sr/(se/(n-2)); fv # F値
tv <- b/sb; tv # t値

**# 回帰分析の詳細の R スクリプト**
anova(mbaf.lm)	# 分散分析表
n <- nrow(mbaf); n	# データ数
st <- var(mbaf[,3])*(n-1); st	#　総平方和
sr <- var(mbaf.lm$fitted.values)*(n-1); sr	#　回帰平方和
se <- var(mbaf.lm$residuals)*(n-1); se	#　残差平方和
s <- sqrt(se/(n-2)); s	#　残差標準偏差
sx <- var(mbaf[,2])*(n-1); sx	#　x の総平方和
b <- cov(mbaf[,2],mbaf[,3])*(n-1)/sx; b	#　回帰係数
sb <- s/sqrt(sx); sb	#　回帰係数bの標準誤差
r2 <- sr/st; r2	# 重相関係数の２乗
r <- cor(mbaf[,3], mbaf.lm$fitted.values); r	# 目的変数と回帰推定値の相関
fv <- sr/(se/(n-2)); fv	# F値
tv <- b/sb; tv	# t値

回帰式の信頼区間

回帰係数の信頼区間

　回帰係数の標準誤差 s _b を用いて，回帰係数 b の信頼区間がつくれる．すなわち，自由度 n－2 の t 分布の 97.5％点を t₀ とすると，回帰係数 b の 95％信頼区間の幅 d は，d ＝ t₀ s _b となるので， 95％信頼区間は，

b^ － t₀ s _b ＜ b ＜ b^ ＋ t₀ s _b

となる．

回帰直線の信頼区間

　データから推定された回帰直線は，データの平均（x^-，y^-）を通るので，

Y = y^- + b^ （x - x^-）

とおける．これより，Y の分散は，

となる．誤差分散 σ² は，データの残差分散 s_e² で推定されるので，Y の標準誤差は，

となり，これが自由度 n - 2 を持つ．よって，推定回帰式にこの標準誤差の t₀ 倍を加えたものが，回帰式の 95％信頼幅となる．

回帰予測値の信頼幅

　回帰式から得られる予測値 y~ は，回帰式に誤差項が加わって，

y~ = y^- + b^ （x - x^-）+ e

となるので，その分散は，

となる．先ほどと同様に，誤差標準偏差を残差標準偏差で置き換えると，回帰予測値 y^~ の標準誤差は，

となり，回帰式の 95％信頼幅の外側に回帰予測値の 95％信頼幅が描ける．

# 回帰式の信頼幅の R スクリプト
x <- mbaf[,2] # 入試得点
y <- mbaf[,3] # 初年度成績
d <- lm(y ~ x) # 回帰モデル
new <- data.frame(x=seq(400,700,by=2)) # 予測したい範囲の定義
dc <- predict(d, new, interval="confidence", level=0.95) # 回帰推定値（回帰直線）の 95 ％信頼幅
dp <- predict(d, new, interval="prediction", level=0.95) # 回帰予測値の 95 ％信頼幅

matplot(new$x, cbind(dc,dp[,-1]), lty=c(1,2,2,3,3), type="l", col=c("blue","blue","blue","red","red"), xlab="入試得点", ylab="初年度成績")
points(x, y) # データの表示
title(main="回帰式と予測値の 95 ％信頼幅") #

# 信頼幅を詳しく計算
n <- nrow(mbaf); n # データ数
sx <- var(x)*(n-1); sx # x の偏差平方和
se <- summary(d)$sigma; se # 残差標準偏差
b <- d$coefficients[2]; b # 回帰係数
mx <- mean(x) # x の平均
my <- mean(y) # y の平均
t0 <- qt(0.975, df=(n-2)); t0 # 自由度 n - 2 の t 分布 97.5％点
plot(x, y, xlim=c(400,700), ylim=c(2,4)) # データ散布図
abline(d, col="blue") # 回帰直線
sr <- 1/n + (new$x-mx)^2/sx #
y1 <- my + b*(new$x-mx) + t0*se*sqrt(sr) # 回帰直線 95％信頼幅上限
y2 <- my + b*(new$x-mx) - t0*se*sqrt(sr) # 回帰直線 95％信頼幅下限
points(new$x, y1, type="l", lty=2, col="blue") #
points(new$x, y2, type="l", lty=2, col="blue") #
yp1 <- my + b*(new$x-mx) + t0*se*sqrt(sr+1) # 回帰予測値 95％信頼幅上限
yp2 <- my + b*(new$x-mx) - t0*se*sqrt(sr+1) # 回帰予測値 95％信頼幅下限
points(new$x, yp1, type="l", lty=3, col="red") #
points(new$x, yp2, type="l", lty=3, col="red") #
title(main="回帰式と予測値の 95 ％信頼幅") #

**# 回帰式の信頼幅の R スクリプト**
x <- mbaf[,2]	# 入試得点
y <- mbaf[,3]	# 初年度成績
d <- lm(y ~ x)	# 回帰モデル
new <- data.frame(x=seq(400,700,by=2))	# 予測したい範囲の定義
dc <- predict(d, new, interval="confidence", level=0.95)	# 回帰推定値（回帰直線）の 95 ％信頼幅
dp <- predict(d, new, interval="prediction", level=0.95)	# 回帰予測値の 95 ％信頼幅
matplot(new$x, cbind(dc,dp[,-1]), lty=c(1,2,2,3,3), type="l", col=c("blue","blue","blue","red","red"), xlab="入試得点", ylab="初年度成績")
points(x, y)	# データの表示
title(main="回帰式と予測値の 95 ％信頼幅")	#
# 信頼幅を詳しく計算
n <- nrow(mbaf); n	# データ数
sx <- var(x)*(n-1); sx	# x の偏差平方和
se <- summary(d)$sigma; se	# 残差標準偏差
b <- d$coefficients[2]; b	# 回帰係数
mx <- mean(x)	# x の平均
my <- mean(y)	# y の平均
t0 <- qt(0.975, df=(n-2)); t0	# 自由度 n - 2 の t 分布 97.5％点
plot(x, y, xlim=c(400,700), ylim=c(2,4))	# データ散布図
abline(d, col="blue")	# 回帰直線
sr <- 1/n + (new$x-mx)^2/sx	#
y1 <- my + b(new$x-mx) + t0se*sqrt(sr)	# 回帰直線 95％信頼幅上限
y2 <- my + b(new$x-mx) - t0se*sqrt(sr)	# 回帰直線 95％信頼幅下限
points(new$x, y1, type="l", lty=2, col="blue")	#
points(new$x, y2, type="l", lty=2, col="blue")	#
yp1 <- my + b(new$x-mx) + t0se*sqrt(sr+1)	# 回帰予測値 95％信頼幅上限
yp2 <- my + b(new$x-mx) - t0se*sqrt(sr+1)	# 回帰予測値 95％信頼幅下限
points(new$x, yp1, type="l", lty=3, col="red")	#
points(new$x, yp2, type="l", lty=3, col="red")	#
title(main="回帰式と予測値の 95 ％信頼幅")	#

重回帰分析

統計モデル

　説明変数が 2 つ以上になった場合である．いま，p 個の説明変数 x₁， x₂，…，x_p，により目的変数 y が，

y = b₀ + b₁x₁ + b₂x₂ + … + b_px_m + e

と表現できるとする．ここで，b₀b₁，…，b_p は回帰係数，e は，誤差である．
　ここで，n（n ＞ p）個のデータがあり，目的変数が y = (y₁，…，y_n)'，j 番目の説明変数 x_j の値が， x₁ = (x_1j，…，x_nj)' であったとする．回帰係数ベクトルを b = (b₀，…，b_m)'，誤差ベクトルを e = (e₁，…，e_n)' とすると，データが満たす構造モデルは，

となる．ここで，誤差ベクトルは，平均 0，分散共分散行列 σ²I の n 次元正規分布に従うとする．これは，n 個の誤差 e_i が互いに独立に平均 0 分散 σ² の正規分布に従うと考えたとき（普通行う仮定）の多変量表現である．
　最小２乗法により，回帰係数の推定値は，

と求められる．

多変量の線形関数の平均と分散

　p 変量確率変数 x の平均を μ，分散共分散行列を Σ とすると，任意の p×p スカラー行列 A に対し，

E[Ax] = AE[x] = Aμ， Var[Ax] = AVar[x]A' = AΣA'，

が成り立つ．

回帰係数の分散

　重回帰モデルの回帰係数ベクトルの推定量 b^{^} の分散は，

となる．残差分散 σ² の推定値を s_e² とすると，回帰係数の分散は，s_e²(X'X)^-1 で推定される．

ハット行列とてこ比

　さて，回帰推定値ベクトル y^{^} は，

と書ける．ここで，H は，説明変数データベクトル y から回帰推定値ベクトル y^{^} を生成する行列なので，ハット行列（hat matrix）と呼ばれている．ハット行列の対角成分 h_ii が大きいと i 番目のデータは推定値に大きな影響を及ぼすので，これを「てこ比（leverage）」と呼んでいる．てこ比の平均的な値は，(p + 1)/n なので，たとえば，2(p + 1)/n より大きなてこ比をもつようなデータには注意を払う必要がある．すなわち，このデータの存在が回帰係数の推定に大きな影響を与えているからである．

標準化残差

　回帰残差推定量ベクトル e^{^} は，

とハット行列を用いて書ける．これより，その分散は，

となる．なお，I - H がべき等行列であることに注意せよ．従って，個々の回帰残差推定量 e^{^}_i の分布は，

となる．この r_i を標準化残差（standarized residual）と言う．

Cook 距離

　個々のデータが回帰係数の推定値に与える影響を調べるため，i 番目のデータを除いたときの回帰係数推定値を b^{^}_-i として，これが全データを用いたときの回帰係数推定値 b^{^} との違いの大きさで測る．これが Cook 距離で，以下のように定義される．

ここで，y^{^}_-i = Xb^{^}_-i である．また，多少の計算により Cook 距離は標準化残差とてこ比で計算され，最後の式のように表せる．
Cook 距離が 0.5 を超えるとそのデータは影響が「大きい」とされ，1 を超えると影響が「特に大きい」とされる．

回帰診断

　回帰分析を行ったときに，データが回帰モデルによくフィットしているかや，データの中に回帰モデルからはずれたも（異常値）がないか，などを調べた方がより安全である．これを回帰診断（regression diagnostics）という． R などの統計ソフトが普及する以前は，手間がかかるので回帰診断まで行うことはあまりやらなかったが，現在では手軽にできるので，行うのが普通になってきていると思われる．このため，回帰診断の考え方や着目点などを理解する必要がでてきたと言える．　回帰診断は，以下の２点からなる．

誤差分析：回帰残差（誤差）は，独立，等分散，正規性の３つの仮定を置くことが普通であるが，これが満たされているかを調べる．
誤差に何らかのトレンドが見られるときは回帰モデルがデータにフィットしていない可能性があるので，別の回帰モデル（多項式回帰など）を試す必要がある．
極端に誤差の大きな観測値は異常値（out lier）の可能性がある．
感度分析：観測値（データ）が回帰係数に及ぼす影響をみる．少数のデータが回帰係数に大きな影響を与えているときは，このデータの処遇を考える必要がある．

**# 回帰診断の R スクリプト**
mba <- read.csv("mbagrade2.csv") # データ読み込み　 mbaf <- mba[mba[,1]=="F",] # 女性データのみ抽出　 mbaf.lm <- lm(初年度成績 ~ 入試得点, data=mbaf) # summary(mbaf.lm) # 結果表示　 plot(mbaf[,2:3], type="n", xlim=c(400,700), ylim=c(2.5,3.5)) # データ散布図　 text(mbaf[,2:3], rownames(mbaf)) # データ番号表示 abline(mbaf.lm, col="blue") # 回帰直線　 x11() # 新しいグラフウィンドウ plot(mbaf.lm) # 回帰診断はたったこれだけ．（４枚のグラフが出てくる） # # 回帰診断を詳しくみる # てこ比 lev <- hatvalues(mbaf.lm); lev # てこ比 X <- model.matrix(mbaf.lm) # 説明変数行列 H <- X %% solve(t(X) %% X) %% t(X) # ハット行列 diag(H) # てこ比の計算による導出 # 標準化残差 rstd <- rstandard(mbaf.lm); rstd # 標準化残差 se <- summary(mbaf.lm)$sigma # 残差標準偏差 mbaf.lm$residuals/(sesqrt(1-lev)) # 標準化残差の計算による導出 # Cook 距離 x <- mbaf$入試得点 y <- mbaf$初年度成績 cookd <- cooks.distance(mbaf.lm); cookd # Cook 距離 b <- mbaf.lm$coefficients; b # 回帰係数 lm1 <- lm(y[-1] ~ x[-1]) b1 <- lm1$coefficients; b1 # １番目のデータを除いた回帰係数 t(b-b1) %% t(X) %% X %% (b-b1)/se^2/2 # １番目のデータの Cook 距離 rstd^2lev/(1-lev)/2 # Cook 距離の計算による導出 # 回帰推定値 - 残差プロット x0 <- mbaf.lm$fitted.values # 回帰推定値 y0 <- mbaf.lm$residuals # 回帰残差 plot(x0, y0, type="n", xlab="回帰推定値", ylab="残差") # text(x0, y0, rownames(mbaf), cex=0.8) abline(h=0, lty=3) title(main="回帰推定値 - 残差プロット") # 標準化残差 Q - Q プロット a <- qqnorm(rstd, type="n", main="") # 標準化残差の Q - Q プロット（正規分布からのずれ） text(a, rownames(mbaf), cex=0.8) qqline(rstd, col="red") title("標準化残差正規 Q-Q プロット") # S - L プロット yr <- sqrt(abs(rstd)) # 標準化残差絶対値の平方根 plot(x0, yr, type="n", xlab="回帰推定値", ylab="残差絶対値の平方根") # text(x0, yr, rownames(mbaf), cex=0.8) title(main="S - L プロット") # てこ比 - 標準化残差プロット plot(lev, rstd, type="n", xlim=c(0,0.5), xlab="てこ比", ylab="標準化残差") text(lev, rstd, rownames(mbaf), cex=0.8) abline(h=0, lty=3) abline(v=0, lty=3) title(main="てこ比 - 標準化残差プロット") xx <- seq(0, 0.5, by=0.01) yc <- .52(1-xx)/xx # Cook 距離 0.5 となる標準化残差 yc1 <- sqrt(yc) yc2 <- -yc1 points(xx, yc1, type="l", lty=2, col="red") # Cook 距離 0.5 線 points(xx, -yc1, type="l", lty=2, col="red")

# 回帰診断の R スクリプト

mba <- read.csv("mbagrade2.csv")  		# データ読み込み　 
mbaf <- mba[mba[,1]=="F",]  		# 女性データのみ抽出 　 
mbaf.lm <- lm(初年度成績 ~ 入試得点, data=mbaf) # 
summary(mbaf.lm)  				# 結果表示 　 
plot(mbaf[,2:3], type="n", xlim=c(400,700), ylim=c(2.5,3.5))  # データ散布図 　 
text(mbaf[,2:3], rownames(mbaf))		# データ番号表示
abline(mbaf.lm, col="blue")  		# 回帰直線 　 
x11()					# 新しいグラフウィンドウ
plot(mbaf.lm)				# 回帰診断はたったこれだけ．（４枚のグラフが出てくる）
#
# 回帰診断を詳しくみる
# てこ比
lev <- hatvalues(mbaf.lm); lev		# てこ比
X <- model.matrix(mbaf.lm)			# 説明変数行列
H <- X %*% solve(t(X) %*% X) %*% t(X)	# ハット行列
diag(H)					# てこ比の計算による導出
# 標準化残差
rstd <- rstandard(mbaf.lm); rstd		# 標準化残差
se <- summary(mbaf.lm)$sigma		# 残差標準偏差
mbaf.lm$residuals/(se*sqrt(1-lev))		# 標準化残差の計算による導出
# Cook 距離
x <- mbaf$入試得点
y <- mbaf$初年度成績
cookd <- cooks.distance(mbaf.lm); cookd	# Cook 距離
b <- mbaf.lm$coefficients; b 		# 回帰係数
lm1 <- lm(y[-1] ~ x[-1])
b1 <- lm1$coefficients; b1 			# １番目のデータを除いた回帰係数
t(b-b1) %*% t(X) %*% X %*% (b-b1)/se^2/2	# １番目のデータの Cook 距離
rstd^2*lev/(1-lev)/2			# Cook 距離の計算による導出
# 回帰推定値 - 残差プロット
x0 <- mbaf.lm$fitted.values			# 回帰推定値
y0 <- mbaf.lm$residuals			# 回帰残差
plot(x0, y0, type="n", xlab="回帰推定値", ylab="残差")	# 
text(x0, y0, rownames(mbaf), cex=0.8)
abline(h=0, lty=3)
title(main="回帰推定値 - 残差プロット")
# 標準化残差 Q - Q プロット
a <- qqnorm(rstd, type="n", main="")		# 標準化残差の Q - Q プロット（正規分布からのずれ）
text(a, rownames(mbaf), cex=0.8)
qqline(rstd, col="red")
title("標準化残差正規 Q-Q プロット")
# S - L プロット
yr <- sqrt(abs(rstd))			# 標準化残差絶対値の平方根
plot(x0, yr, type="n", xlab="回帰推定値", ylab="残差絶対値の平方根")	# 
text(x0, yr, rownames(mbaf), cex=0.8)
title(main="S - L プロット")
# てこ比 - 標準化残差 プロット
plot(lev, rstd, type="n", xlim=c(0,0.5), xlab="てこ比", ylab="標準化残差")
text(lev, rstd, rownames(mbaf), cex=0.8)
abline(h=0, lty=3)
abline(v=0, lty=3)
title(main="てこ比 - 標準化残差 プロット")
xx <- seq(0, 0.5, by=0.01)
yc <- .5*2*(1-xx)/xx			# Cook 距離 0.5 となる標準化残差
yc1 <- sqrt(yc)
yc2 <- -yc1
points(xx, yc1, type="l", lty=2, col="red")		# Cook 距離 0.5 線
points(xx, -yc1, type="l", lty=2, col="red")

重回帰分析の例

入試得点と初年度成績：性別を加えた重回帰
　前節では，アメリカのあるビジネススクール（MBA）での女子学生の入試得点と初年度成績の回帰分析を行った．男子学生のデータもあるので，これを加え，性別（x₁）と入試得点（x₂）という２つの変数を説明変数として初年度成績（y）の重回帰分析を行ってみる．

# 成績重回帰分析の R スクリプト
mba <- read.csv("mbagrade2.csv") # データ読み込み
mbaf <- mba[mba[,1]=="F",] # 女子データ
mbam <- mba[mba[,1]=="M",] # 男子データ
mba.lm <- lm(初年度成績 ~ 性別 + 入試得点, data=mba) # 重回帰分析
summary(mba.lm) # 結果表示
plot(mba[,2:3], type="n", xlim=c(400,750)) #
points(mbam[,2:3], pch=21) # 男子データ（白丸）
points(mbaf[,2:3], pch=19) # 女子データ（黒丸）
b <- mba.lm$coefficients # 回帰係数
x <- seq(400,750,by=1) #
ym <- b[1] + b[2] + b[3]*x # 男子回帰式
yf <- b[1] + b[3]*x # 女子回帰式
lines(x, ym, type="l", col="blue") #
lines(x, yf, type="l", col="red") #
legend(locator(1), c("男子","女子"), pch=c(21,19)) #
title(main="MBA 入試得点と初年度成績") #

**# 成績重回帰分析の R スクリプト**
mba <- read.csv("mbagrade2.csv")	# データ読み込み
mbaf <- mba[mba[,1]=="F",]	# 女子データ
mbam <- mba[mba[,1]=="M",]	# 男子データ
mba.lm <- lm(初年度成績 ~ 性別 + 入試得点, data=mba)	# 重回帰分析
summary(mba.lm)	# 結果表示
plot(mba[,2:3], type="n", xlim=c(400,750))	#
points(mbam[,2:3], pch=21)	# 男子データ（白丸）
points(mbaf[,2:3], pch=19)	# 女子データ（黒丸）
b <- mba.lm$coefficients	# 回帰係数
x <- seq(400,750,by=1)	#
ym <- b[1] + b[2] + b[3]*x	# 男子回帰式
yf <- b[1] + b[3]*x	# 女子回帰式
lines(x, ym, type="l", col="blue")	#
lines(x, yf, type="l", col="red")	#
legend(locator(1), c("男子","女子"), pch=c(21,19))	#
title(main="MBA 入試得点と初年度成績")	#

入試得点と初年度成績：性別ごとの単回帰
　重回帰分析を行うと，入試得点に対する回帰係数が男子，女子とも同じ値になり，男子と女子の違いは， y 切片の違いのみになってしまう．男子と女子で入試得点に対する回帰係数が異なると考えられる場合には，男子，女子で別々の単回帰を行う．しかし，検定に用いる誤差分散は男女共通である．

# 男女別々の単回帰分析の R スクリプト
mbaf.lm <- lm(初年度成績 ~ 入試得点, data=mbaf) # 女子回帰
summary(mbaf.lm) # 結果表示
mbam.lm <- lm(初年度成績 ~ 入試得点, data=mbam) # 男子回帰
summary(mbam.lm) # 結果表示
plot(mba[,2:3], type="n", xlim=c(400,750)) #
points(mbam[,2:3], pch=21) # 男子データ（白丸）
points(mbaf[,2:3], pch=19) # 女子データ（黒丸）
abline(mbaf.lm, col="red") # 女子回帰式
abline(mbam.lm, col="blue") # 男子回帰式
legend(locator(1), c("男子","女子"), pch=c(21,19)) #
title(main="MBA 入試得点と初年度成績") #
mba2.lm <- lm(初年度成績 ~ 性別 + 性別 : 入試得点, data=mba) # 性別ごとの回帰
summary(mba2.lm) # 結果表示
anova(mba2.lm) # 分散分析表示

**# 男女別々の単回帰分析の R スクリプト**
mbaf.lm <- lm(初年度成績 ~ 入試得点, data=mbaf)	# 女子回帰
summary(mbaf.lm)	# 結果表示
mbam.lm <- lm(初年度成績 ~ 入試得点, data=mbam)	# 男子回帰
summary(mbam.lm)	# 結果表示
plot(mba[,2:3], type="n", xlim=c(400,750))	#
points(mbam[,2:3], pch=21)	# 男子データ（白丸）
points(mbaf[,2:3], pch=19)	# 女子データ（黒丸）
abline(mbaf.lm, col="red")	# 女子回帰式
abline(mbam.lm, col="blue")	# 男子回帰式
legend(locator(1), c("男子","女子"), pch=c(21,19))	#
title(main="MBA 入試得点と初年度成績")	#
mba2.lm <- lm(初年度成績 ~ 性別 + 性別 : 入試得点, data=mba)	# 性別ごとの回帰
summary(mba2.lm)	# 結果表示
anova(mba2.lm)	# 分散分析表示

モデル選択と AIC

　確率モデルのパラメータ推定には，通常，最尤法が用いられる．しかしながら，重回帰分析などで説明変数の個数（パラメータ数）を決めようとすると，一般に，パラメータ数が多いほどデータへのモデルの当てはまり（fitting）が良くなるので，最尤法でパラメータ数を決めるとパラメータ数の多いモデルが「良い」とされてしまう．パラメータ数の多いモデルは，パラメータの値を推定したデータにはよく当てはまるが，同様の状況から得られた別のデータへの当てはまりが悪くなることが知られている．このような現象を解釈しすぎ（over fitting）という．
　これを避けるには，できるだけ単純なモデルを考えるのがよいとされている．これを実現するモデル選択の基準として，

X = (モデルのデータへの当てはまり) + (モデルの複雑さへのペナルティ)

の形式のものがいくつか提案されている．この中で有名な基準の一つが AIC (Akaike Information Criterion) である．AIC は，

AIC = - 2×(モデルの最大対数尤度) + 2×(モデルの自由パラメータ数)

と定義される．モデルの最大対数尤度は，確率モデルの最尤推定値を確率モデルに代入したときの尤度の対数を取ったものであり，モデルのデータへの当てはまりのよさを評価している．モデルの自由パラメータ数は，モデルの複雑さの尺度の一つで，パラメータ数の少ないモデルほど単純でよいものと考えられる．結局，AIC の小さなモデルがよいとされる．

　k 個のパラメータ θ を持つ回帰モデル

y = f(x ; θ) + e

において，残差 e が正規分布に従うモデルでは，n 個のデータから得られた残差分散の最尤推定値を v² とすると，回帰モデルの最大対数尤度は，

l = -(1/2) [n log 2π + n log v² + n ]

となる．これより，回帰モデルでの AIC は，

AIC = (n log 2π + n log v² + n ) + 2k

となり，AIC の小さな回帰モデルがよいとされる．

入試得点と初年度成績：モデル選択
　入試得点と初年度成績の関係を分析するときに，３つのモデルが考えられる．それは，

性別，入試得点重回帰モデル（回帰係数，y 切片 2，残差分散の 4 パラメータ）
残差分散と入試得点に対する傾きの回帰係数は男女共通で，y 切片が性別により異なる．
男女別々回帰，残差分散共通モデル（回帰係数 2，切片 2，残差分散の 5 パラメータ）
残算分散は男女共通であるが，回帰係数，y 切片とも男女で異なる．
男女完全別々モデル（回帰係数 2，切片 2，残差分散 2 の 6 パラメータ）
男子と女子では全く異なり，共通のものがない．

各モデルの分散分析表で，各モデルに意味があるかは答えてくれるが，どのモデルがデータとよく合っているかを判断してくれる数値がない．このようなときに AIC を用いるとモデル間の比較ができる．各モデルの自由パラメータは括弧内に記載されている．
　AIC の結果から男女完全別々モデルがこのデータにはよく適合しているようである．すなわち，女子は入試得点と初年度成績に強い回帰関係が認められる（R ² = 0.698）が，男子は残差分散が大きく，強い回帰関係が無い（R ² = 0.088）ことからもこの結果は支持される．女子は初年度成績は入試得点でかなり予測できるが，男子は初年度成績のばらつきが大きいので，成績を予測することはできない．

# モデル選択の R スクリプト
anova(mba.lm) # 性別，入試得点重回帰の分散分析表
AIC(mba.lm) # 性別，入試得点重回帰モデルの AIC
n <- nrow(mba) # 総データ数
v2 <- anova(mba.lm)[3,2]/n # 残差分散の最尤推定値（n で割る）
aic <- n*log(2*pi) + n*log(v2) + n + 2*4; aic # 自由パラメータ数 4 の AIC
anova(mba2.lm) # 男女別々回帰，残差分散共通の分散分析表
AIC(mba.lm) # 男女別々回帰，残差分散共通モデルの AIC
anova(mbaf.lm) # 男子単回帰分散分析表
anova(mbam.lm) # 女子単回帰分散分析表
AIC(mbaf.lm)+AIC(mbam.lm) # 男女完全別々モデルの AIC

**# モデル選択の R スクリプト**
anova(mba.lm)	# 性別，入試得点重回帰の分散分析表
AIC(mba.lm)	# 性別，入試得点重回帰モデルの AIC
n <- nrow(mba)	# 総データ数
v2 <- anova(mba.lm)[3,2]/n	# 残差分散の最尤推定値（n で割る）
aic <- nlog(2pi) + nlog(v2) + n + 24; aic	# 自由パラメータ数 4 の AIC
anova(mba2.lm)	# 男女別々回帰，残差分散共通の分散分析表
AIC(mba.lm)	# 男女別々回帰，残差分散共通モデルの AIC
anova(mbaf.lm)	# 男子単回帰分散分析表
anova(mbam.lm)	# 女子単回帰分散分析表
AIC(mbaf.lm)+AIC(mbam.lm)	# 男女完全別々モデルの AIC

1 月の日最低気温の月平均値
　1941年～1970年での日本の各都市における 1 月の日最低気温の月平均値（y）が，各都市の緯度（x₁），経度（x₂），標高（x₃）でよく説明できるかを重回帰分析で解析してみる．データは以下の通りである．

データダウンロード

# 1 月最低気温重回帰の R スクリプト
kion <- read.csv("kion.csv") # 気温データ読み込み
pairs(kion[,2:5]) # 変数間散布図一覧
cor(kion[,2:5]) # 変数間相関
kion1.lm <- lm(気温 ~ 緯度 + 経度 + 標高, data=kion) # 3 変数重回帰
summary(kion1.lm) # 結果表示
anova(kion1.lm) # 分散分析表示
n <- nrow(kion) # データ数
x0 <- rep(1,n) #
x <- as.matrix(cbind(x0, kion[,3:5])) # 説明変数行列
se2 <- anova(kion1.lm)[4,3] # 残差分散
v <- se2 * solve(t(x) %*% x) # 回帰係数の分散共分散行列
sqrt(diag(v)) # 回帰係数の標準誤差

**# 1 月最低気温重回帰の R スクリプト**
kion <- read.csv("kion.csv")	# 気温データ読み込み
pairs(kion[,2:5])	# 変数間散布図一覧
cor(kion[,2:5])	# 変数間相関
kion1.lm <- lm(気温 ~ 緯度 + 経度 + 標高, data=kion)	# 3 変数重回帰
summary(kion1.lm)	# 結果表示
anova(kion1.lm)	# 分散分析表示
n <- nrow(kion)	# データ数
x0 <- rep(1,n)	#
x <- as.matrix(cbind(x0, kion[,3:5]))	# 説明変数行列
se2 <- anova(kion1.lm)[4,3]	# 残差分散
v <- se2 * solve(t(x) %*% x)	# 回帰係数の分散共分散行列
sqrt(diag(v))	# 回帰係数の標準誤差

説明変数の選択

すべての組み合わせの AIC

　全変数を使って重回帰分析を行ったが，経度（x₂）の回帰係数の有意確率が小さくないので，経度の情報は気温を説明するのに必要ないかも知れない．これは日本では，緯度が高く標高が高いほど気温が低いと考えられることとも一致している．モデル選択の方法として AIC を利用してみる．
　3 個の説明変数があるので，説明変数の組み合わせは 2³ = 8 通りある．このすべての組合わせに対して AIC の値を計算し，最も小さな値をもつモデルを採用することにする．

# 変数選択の R スクリプト
kion0.lm <- lm(気温 ~ 1, data=kion) # 説明変数無し回帰
anova(kion0.lm) #
kion11.lm <- lm(気温 ~ 緯度, data=kion) # 説明変数：緯度，回帰
anova(kion11.lm) #
kion12.lm <- lm(気温 ~ 経度, data=kion) # 説明変数：経度，回帰
anova(kion12.lm) #
kion13.lm <- lm(気温 ~ 標高, data=kion) # 説明変数：標高，回帰
anova(kion13.lm) #
kion21.lm <- lm(気温 ~ 緯度 + 経度, data=kion) # 説明変数：緯度，経度，回帰
anova(kion21.lm) #
kion22.lm <- lm(気温 ~ 経度 + 標高, data=kion) # 説明変数：経度，標高，回帰
anova(kion22.lm) #
kion23.lm <- lm(気温 ~ 緯度 + 標高, data=kion) # 説明変数：緯度，標高，回帰
anova(kion23.lm) #
kion3.lm <- lm(気温 ~ 緯度 + 経度 + 標高, data=kion) # 説明変数：緯度，経度，標高，回帰
anova(kion3.lm) #
AIC(kion0.lm) # 説明変数無し回帰 AIC
AIC(kion11.lm) # 説明変数：緯度，回帰 AIC
AIC(kion12.lm) # 説明変数：経度，回帰 AIC
AIC(kion13.lm) # 説明変数：標高，回帰 AIC
AIC(kion21.lm) # 説明変数：緯度，経度，回帰 AIC
AIC(kion22.lm) # 説明変数：経度，標高，回帰 AIC
AIC(kion23.lm) # 説明変数：緯度，標高，回帰 AIC
AIC(kion3.lm) # 説明変数：緯度，経度，標高，回帰 AIC

**# 変数選択の R スクリプト**
kion0.lm <- lm(気温 ~ 1, data=kion)	# 説明変数無し回帰
anova(kion0.lm)	#
kion11.lm <- lm(気温 ~ 緯度, data=kion)	# 説明変数：緯度，回帰
anova(kion11.lm)	#
kion12.lm <- lm(気温 ~ 経度, data=kion)	# 説明変数：経度，回帰
anova(kion12.lm)	#
kion13.lm <- lm(気温 ~ 標高, data=kion)	# 説明変数：標高，回帰
anova(kion13.lm)	#
kion21.lm <- lm(気温 ~ 緯度 + 経度, data=kion)	# 説明変数：緯度，経度，回帰
anova(kion21.lm)	#
kion22.lm <- lm(気温 ~ 経度 + 標高, data=kion)	# 説明変数：経度，標高，回帰
anova(kion22.lm)	#
kion23.lm <- lm(気温 ~ 緯度 + 標高, data=kion)	# 説明変数：緯度，標高，回帰
anova(kion23.lm)	#
kion3.lm <- lm(気温 ~ 緯度 + 経度 + 標高, data=kion)	# 説明変数：緯度，経度，標高，回帰
anova(kion3.lm)	#
AIC(kion0.lm)	# 説明変数無し回帰 AIC
AIC(kion11.lm)	# 説明変数：緯度，回帰 AIC
AIC(kion12.lm)	# 説明変数：経度，回帰 AIC
AIC(kion13.lm)	# 説明変数：標高，回帰 AIC
AIC(kion21.lm)	# 説明変数：緯度，経度，回帰 AIC
AIC(kion22.lm)	# 説明変数：経度，標高，回帰 AIC
AIC(kion23.lm)	# 説明変数：緯度，標高，回帰 AIC
AIC(kion3.lm)	# 説明変数：緯度，経度，標高，回帰 AIC

変数増減法の AIC

　説明変数の数が多くなり，すべての組み合わせを調べることが大変な場合には，変数増加法，変数減少法，その組み合わせである変数増減法がある．ここでは，変数増減法もやってみる．

# 変数増減法の R スクリプト
library(MASS) # MASS ライブラリィー読み込み
null <- lm(気温 ~ 1, kion) # 説明変数無し
full <- lm(気温 ~ 緯度 + 経度 + 標高, kion) # 説明変数３つ
result <- stepAIC(null, scope=list(lower=null, upper=full), data=kion) # 変数増減法
summary(result) # 結果表示

**# 変数増減法の R スクリプト**
library(MASS)	# MASS ライブラリィー読み込み
null <- lm(気温 ~ 1, kion)	# 説明変数無し
full <- lm(気温 ~ 緯度 + 経度 + 標高, kion)	# 説明変数３つ
result <- stepAIC(null, scope=list(lower=null, upper=full), data=kion)	# 変数増減法
summary(result)	# 結果表示

多項式回帰

　目的変数 y と説明変数 x との関係が直線関係では説明できないようなときには， x の多項式で回帰することを考える．p 次までの多項式を考えると，

y = b₀ + b₁x + b₂x² + … + b_px^p + e

というモデルとなる．ここで，x₂ = x²，…， x_p = x^p と変数変換を行えば通常の重回帰分析と同じである．

自動車の速度と制動距離
　自動車の速度 x mile/h とブレーキをかけてからの静止するまでの距離 y ft のデータに対し， y を x での多項式回帰を行ってみる．自動車の速度が 0 のときの制動距離が 0 なので，原点を通る多項式回帰が自然であるが，一般の多項式回帰もあてはめてみる．考えたモデルは，

直線回帰： y = b₀ + b₁x
原点を通る直線回帰：y = b₁x
２次回帰： y = b₀ + b₁x + b₂x²
原点を通る２次回帰： y = b₁x + b₂x²
３次回帰： y = b₀ + b₁x + b₂x² + b₃x³
原点を通る３次回帰： y = b₁x + b₂x² + b₃x³

である．よさそうなモデルの選択に AIC を用いたところ，原点を通る２次式回帰のあてはまりがよいとの結果を得た．

# 多項式回帰の R スクリプト
cars # R 付属データ 'cars' の呼び出し
plot(cars, xlim=c(0,30), xlab="速度", ylab="制動距離") # cars の散布図
title(main="制動距離に対する速度の多項式回帰") #
y <- cars$dist # 目的変数（制動距離）
x <- cars$speed # 説明変数（速度）
car1.lm <- lm(y ~ x) # １次回帰
abline(car1.lm, lty=3, col="red") # 赤点線で表示
summary(car1.lm) #
car10.lm <- lm(y ~ 0 + x) # 原点を通る１次回帰
abline(car10.lm, col="red") # 赤実線で表示
summary(car10.lm) #
x2 <- x^2 # x の２乗を新しい変数で定義
car2.lm <- lm(y ~ x + x2) # 重回帰（２次式回帰）
b <- car2.lm$coefficients # 回帰係数
xv <- seq(0, 30, by=0.2) # x の点列の定義
yv <- b[1] + b[2]*xv + b[3]*xv^2 # x の点列に対する２次回帰推定値
lines(xv, yv, lty=3, col="blue") # 青点線で表示
summary(car2.lm) #
car20.lm <- lm(y ~ 0 + x + x2) # 原点を通る２次式回帰
b <- car20.lm$coefficients # 回帰係数
yv <- b[1]*xv + b[2]*xv^2 # x の点列に対する２次回帰推定値
lines(xv, yv, col="blue") # 青実線で表示
summary(car20.lm) #
x3 <- x^3 # x の３乗を新しい変数で定義
car3.lm <- lm(y ~ x + x2 + x3) # ３変数重回帰（３次式回帰）
b <- car3.lm$coefficients # 回帰係数
yv <- b[1] + b[2]*xv + b[3]*xv^2 + b[4}*xv^3 # x の点列に対する３次回帰推定値
lines(xv, yv, lty=3, col="green") # 緑点線で表示
summary(car3.lm) #
car30.lm <- lm(y ~ 0 + x + x2 + x3) #
b <- car30.lm$coefficients # 回帰係数
yv <- b[1]*xv + b[2]*xv^2 + b[3]*xv^3 # x の点列に対する３次回帰推定値
lines(xv, yv, col="green") # 緑実線で表示
summary(car30.lm) #
AIC(car1.lm) # １次回帰の AIC
AIC(car10.lm) # 原点を通る１次回帰の AIC
AIC(car2.lm) # ２次式回帰の AIC
AIC(car20.lm) # 原点を通る２次式回帰の AIC
AIC(car3.lm) # ３次式回帰の AIC
AIC(car30.lm) # 原点を通る３次式回帰の AIC

**# 多項式回帰の R スクリプト**
cars	# R 付属データ 'cars' の呼び出し
plot(cars, xlim=c(0,30), xlab="速度", ylab="制動距離")	# cars の散布図
title(main="制動距離に対する速度の多項式回帰")	#
y <- cars$dist	# 目的変数（制動距離）
x <- cars$speed	# 説明変数（速度）
car1.lm <- lm(y ~ x)	# １次回帰
abline(car1.lm, lty=3, col="red")	# 赤点線で表示
summary(car1.lm)	#
car10.lm <- lm(y ~ 0 + x)	# 原点を通る１次回帰
abline(car10.lm, col="red")	# 赤実線で表示
summary(car10.lm)	#
x2 <- x^2	# x の２乗を新しい変数で定義
car2.lm <- lm(y ~ x + x2)	# 重回帰（２次式回帰）
b <- car2.lm$coefficients	# 回帰係数
xv <- seq(0, 30, by=0.2)	# x の点列の定義
yv <- b[1] + b[2]xv + b[3]xv^2	# x の点列に対する２次回帰推定値
lines(xv, yv, lty=3, col="blue")	# 青点線で表示
summary(car2.lm)	#
car20.lm <- lm(y ~ 0 + x + x2)	# 原点を通る２次式回帰
b <- car20.lm$coefficients	# 回帰係数
yv <- b[1]xv + b[2]xv^2	# x の点列に対する２次回帰推定値
lines(xv, yv, col="blue")	# 青実線で表示
summary(car20.lm)	#
x3 <- x^3	# x の３乗を新しい変数で定義
car3.lm <- lm(y ~ x + x2 + x3)	# ３変数重回帰（３次式回帰）
b <- car3.lm$coefficients	# 回帰係数
yv <- b[1] + b[2]xv + b[3]xv^2 + b[4}*xv^3	# x の点列に対する３次回帰推定値
lines(xv, yv, lty=3, col="green")	# 緑点線で表示
summary(car3.lm)	#
car30.lm <- lm(y ~ 0 + x + x2 + x3)	#
b <- car30.lm$coefficients	# 回帰係数
yv <- b[1]xv + b[2]xv^2 + b[3]*xv^3	# x の点列に対する３次回帰推定値
lines(xv, yv, col="green")	# 緑実線で表示
summary(car30.lm)	#
AIC(car1.lm)	# １次回帰の AIC
AIC(car10.lm)	# 原点を通る１次回帰の AIC
AIC(car2.lm)	# ２次式回帰の AIC
AIC(car20.lm)	# 原点を通る２次式回帰の AIC
AIC(car3.lm)	# ３次式回帰の AIC
AIC(car30.lm)	# 原点を通る３次式回帰の AIC

多変量解析

に続く．

参考文献

心理・教育のための統計法（第 2 版），山内光哉，1998，サイエンス社
工学のためのデータサイエンス入門－フリーな統計環境Rを用いたデータ解析－，間瀬茂ら，2004，数理工学社
実践生物統計学－分子から生態まで－（第 1 章，第 2 章），東京大学生物測定学研究室編（大森宏ら）， 2004，朝倉書店
R で学ぶデータマインニング I －データ解析の視点から－，熊谷悦生・船尾暢男，2007，九天社
R で学ぶデータマインニング II －シミュレーションの視点から－，熊谷悦生・船尾暢男，2007，九天社
生物統計学入門，上村賢治・高野泰・大森宏，2008，オーム社

mba <- read.csv("mbagrade2.csv")	# データ読み込み
mbaf <- mba[mba[,1]=="F",]	# 女性データのみ抽出
mbaf.lm <- lm(初年度成績 ~ 入試得点, data=mbaf)	# 回帰：lm(y ~ x, data=zzz), y = ax + b
summary(mbaf.lm)	# 結果表示
plot(mbaf[,2:3], xlim=c(400,700), ylim=c(2.5,3.5))	# データ散布図
abline(mbaf.lm, col="blue")	# 回帰直線
points(mbaf[,2], mbaf.lm$fitted.value, pch=19, col="red")	# 回帰推定値
segments(mbaf[,2],mbaf[,3],mbaf[,2], mbaf.lm$fitted.value)	# 回帰残差
title(main="回帰直線と回帰残差")	#
legend(locator(1), legend=c("データ","回帰推定値"), pch=c(1,19), col=c("black","red"))