Regression
Kosugitti
2014年11月21日
今回は 回帰分析 です。 回帰分析は，これまでの「まとめる」を目的にしていた多変量解析モデルと少し違って，「説明する」「予測する」と表現されることがあります。
回帰分析の実例
まず適用例から見てみましょう。いつものirisデータで相関の高い二変数を見てみます。
summary(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## plot(iris$Sepal.Length,iris$Petal.Length)
cor(iris$Sepal.Length,iris$Petal.Length)## [1] 0.8717538幾つかの例外はありますが，基本的に右肩上がりの，相関関係にあることがわかります。 この関係を，関数\(y=f(x)\)で表現したい，ということを考えます。 関数の形はなんでもいいのですが，一番簡単なのは一次関数，\(y=ax+b\)です。
ここで\(y,x\)はデータです。\(a,b\)は未知数です。データから最適な係数\(a,b\)を見つけることが今回の目的です。この関数のあてはめを (線形)回帰分析，と言います。
result.lm <- lm(Petal.Length~Sepal.Length,data=iris)
summary(result.lm)##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.47747 -0.59072 -0.00668 0.60484 2.49512
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
## Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8678 on 148 degrees of freedom
## Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
## F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16plot(iris$Sepal.Length,iris$Petal.Length)
abline(result.lm, col = "red")
これが説明・予測と言われるのは，関数の形がわかると，\(x\)に新しい数値を代入すれば\(y\)が算出できるからです。
また，これもCFAのように「データにモデルを当てはめる」ことになるので，適合度で評価することになります。
説明変数が複数の場合
irisデータには他にも変数があります。他の変数を追加すると，説明の精度があがります。
result.lm2 <- lm(Petal.Length~Sepal.Length + Sepal.Width,data=iris)
summary(result.lm2)##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Sepal.Width, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.25582 -0.46922 -0.05741 0.45530 1.75599
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.52476 0.56344 -4.481 1.48e-05 ***
## Sepal.Length 1.77559 0.06441 27.569 < 2e-16 ***
## Sepal.Width -1.33862 0.12236 -10.940 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6465 on 147 degrees of freedom
## Multiple R-squared: 0.8677, Adjusted R-squared: 0.8659
## F-statistic: 482 on 2 and 147 DF, p-value: < 2.2e-16\(R^2\)の変化を確認しましょう。
標準化されたスコア
複数の変数が説明する場合，どちらがより重要な要因なのでしょうか？ 相対的な重要度を，係数から読み取ることが不適切な場合があります。ローデータの係数は単位に依存するからです。そこで，データの標準化を考えます。
iris.z <- data.frame(scale(iris[,-5]))
result.lm3 <- lm(Petal.Length~Sepal.Length + Sepal.Width,data=iris.z)
summary(result.lm3)##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Sepal.Width, data = iris.z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71139 -0.26580 -0.03252 0.25792 0.99472
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.384e-16 2.990e-02 0.00 1
## Sepal.Length 8.329e-01 3.021e-02 27.57 <2e-16 ***
## Sepal.Width -3.305e-01 3.021e-02 -10.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3662 on 147 degrees of freedom
## Multiple R-squared: 0.8677, Adjusted R-squared: 0.8659
## F-statistic: 482 on 2 and 147 DF, p-value: < 2.2e-16\(R^2\)や相対的な大小関係がわからないこと，Interceptが0であること，などに注目しよう。
変数選択
変数が大量にある場合，どのような説明変数の組み合わせが良いのか，いちいち試すと大変である。これを，適合度基準を元に，自動的に行うのが変数選択であります。
変数を特に選ばずに全て代入することを 強制投入法 ,まず全ての変数を入れて，効果の弱い変数を徐々に減らしていく 変数減少法,一つの変数から徐々に増やしていく 変数増加法, 変数増減法とも言われる Stepwise法 などがありあます。
まず全ての変数を投入したモデルを作り，ステップワイズ法で変数を選択していくことにします。
result.lm.all <- lm(Petal.Length~.,data=iris.z)
result.best <- step(result.lm.all)## Start: AIC=-509.36
## Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width
##
## Df Sum of Sq RSS AIC
## <none> 4.7662 -509.36
## - Sepal.Width 1 2.9037 7.6700 -440.00
## - Sepal.Length 1 5.1029 9.8691 -402.18
## - Petal.Width 1 14.9485 19.7148 -298.39summary(result.best)##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width,
## data = iris.z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.56270 -0.10002 -0.00569 0.10513 0.60561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.766e-16 1.475e-02 0.000 1
## Sepal.Length 3.420e-01 2.736e-02 12.502 <2e-16 ***
## Sepal.Width -1.595e-01 1.691e-02 -9.431 <2e-16 ***
## Petal.Width 6.247e-01 2.919e-02 21.399 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1807 on 146 degrees of freedom
## Multiple R-squared: 0.968, Adjusted R-squared: 0.9674
## F-statistic: 1473 on 3 and 146 DF, p-value: < 2.2e-16今回は結局全ての変数が重要であることがわかりました。
lavaanによる表現
この分析方法は，実はCFAの時に使ったlavaanパッケージでも表現できます。
library(lavaan)## This is lavaan 0.5-17
## lavaan is BETA software! Please report any bugs.model <- '
Petal.Length ~ Sepal.Width + Sepal.Length + Petal.Width
'
result.sem <- sem(model,estimator="GLS",data=iris)
summary(result.sem,standardized=T,fit.measures=T)## lavaan (0.5-17) converged normally after 57 iterations
##
## Number of observations 150
##
## Estimator GLS
## Minimum Function Test Statistic 0.000
## Degrees of freedom 0
## Minimum Function Value 0.0000000000000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 73.289
## Degrees of freedom 3
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent Confidence Interval 0.000 0.000
## P-value RMSEA <= 0.05 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter estimates:
##
## Information Expected
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|) Std.lv Std.all
## Regressions:
## Petal.Length ~
## Sepal.Width -0.646 0.068 -9.528 0.000 -0.646 -0.160
## Sepal.Length 0.729 0.058 12.630 0.000 0.729 0.342
## Petal.Width 1.447 0.067 21.617 0.000 1.447 0.625
##
## Covariances:
## Sepal.Width ~~
## Sepal.Length -0.042 0.030 -1.425 0.154 -0.042 -0.118
## Petal.Width -0.122 0.029 -4.197 0.000 -0.122 -0.366
## Sepal.Length ~~
## Petal.Width 0.516 0.067 7.728 0.000 0.516 0.818
##
## Variances:
## Petal.Length 0.100 0.012 0.100 0.032
## Sepal.Width 0.190 0.022 0.190 1.000
## Sepal.Length 0.686 0.079 0.686 1.000
## Petal.Width 0.581 0.067 0.581 1.000覚えておいて欲しいことは，
- 構造方程式モデリングStructural Equation Modeling，という方法が因子分析や回帰分析を下位モデルとする統合的な表現であること
- 因子分析と回帰分析は，説明する変数が 観測変数 か 潜在変数 かという違いだけであること
- SEMをつかうと適合度を見ながら自由なモデリングができること
の三点です。
回帰分析のまとめ
- 変数間関係に関数を当てはめることで，「予測」「説明」のアナロジー表現ができること
- データにモデルを当てはめる，という考え方はSEMで統合，表現されている
- 回帰分析の関数は様々に考えられ，二次式以上の曲線モデルを当てはめることもできるが(非線形回帰)，もっとも単純なものが単回帰・重回帰分析と呼ばれる線形回帰分析である
本日の課題
サンプルデータをつかって，国語の成績を最もうまく説明する，他の教科のスコアをつかった回帰式を算出し，回答してください。lavaanではなくlm関数を使うこと。