重回帰のパラメーター推定理論とSAS/IMLによる行列実装
統計学の「基本の「き」」である重回帰を復習しながらSAS/IMLで行列で実装していこうと思います。
理論
重回帰を式で表すと、
となります。
これを行列で表すと
要素を書き下すと
となります。
最小二乗法によるの推定
最小二乗法では、誤差の二乗和を最小化します。つまり、番目の要素に対して、
を最小化します。
ベクトル表記ではを最小化することになるので
に関して最小化するための一階の条件は
より、求める推定量は
となります。
SAS/IMLによる実装
データ
Rにデフォルトで入っている attitude データを使います。
RのデータをSASで使う方法は、以前書いた記事をご覧ください。
y-mattu.hatenablog.com
attitude データは、「管理者に対する態度」というデータセットで、以下のものが int 型で含まれています。
- 全般的格付け
- 雇用者の苦情の処理
- 特別な特権の付与
- 学習の機会
- 昇進機会
- 批判的すぎる
- 進歩
今回は、「全般的格付け」を目的変数に、それ以外を説明変数にしてみます。
SASコードは以下の通りです。
推定したは以下のようになりました。
これが正しいのか、reg
プロシージャの推定で確かめます。
全然有意になってないww
まあでも値はあってました。
今回は単なる勉強のためのコードなのでこれで良しとします。
決定係数とかの行列実装は気が向いたらやります。
それにしてもiml
とreg
の行数の差ときたら...
SASのプロシージャにしろRのパッケージにしろ便利ですね。
おまけ:Rによる実装
行列で
data(attitude) X <- as.matrix(attitude[,-1]) Y <- as.vector(attitude[,1]) X <- cbind(rep(1, nrow(X)), X) Beta <- solve(t(X)%*%X) %*% t(X) %*% Y print(Beta)
結果
[,1] 10.78707639 complaints 0.61318761 privileges -0.07305014 learning 0.32033212 raises 0.08173213 critical 0.03838145 advance -0.21705668
lm
関数でやると
> summary(lm(rating~.,data=attitude)) Call: lm(formula = rating ~ ., data = attitude) Residuals: Min 1Q Median 3Q Max -10.9418 -4.3555 0.3158 5.5425 11.5990 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 10.78708 11.58926 0.931 0.361634 complaints 0.61319 0.16098 3.809 0.000903 *** privileges -0.07305 0.13572 -0.538 0.595594 learning 0.32033 0.16852 1.901 0.069925 . raises 0.08173 0.22148 0.369 0.715480 critical 0.03838 0.14700 0.261 0.796334 advance -0.21706 0.17821 -1.218 0.235577 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 7.068 on 23 degrees of freedom Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628 F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05
当たり前ですが同じ結果です。