実験計画法 交絡と相関
実験計画法における、因子の交絡関係と相関関係に関して。
実験計画関連の業務ではJMPをよく使う。JMPでは、作成した計画の交絡関係を定量的に把握できる行列が出力される。
その交絡関係の出力値と相関係数の関係が不明だったので、ちょっと考えてみた。
スクリーニングのための実験計画は、実験回数を少数に絞りつつ、沢山の因子の中から応答に影響を与える少数の因子を発見することを目的としている。
スクリーニングのための実験計画では、高次項（交互作用項、2次効果項、またはそれ以上）の存在を無視し、低次の項（主効果、交互作用効果項）のみをモデルに加え、これらと応答の関係を把握する。
実験計画法では、線形モデル（重回帰、ANOVA）をあてはめ因子と応答の関係を把握する。
線形モデル：
各ベクトル、行列のサイズは以下。
Xは計画行列と呼ばれ、各列が変数（主効果、交互作用効果、2次効果）。通常、2次効果まで含め、変数の数がp。
nはデータの行数（実験数）。
βはパラメータベクトル。
p変数の内、モデルに含める変数の数をr（r < p）、それ以外の変数の数をp – rとする。
スクリーニング実験なので、モデルに含める変数は低次の項、それ以外の変数は高次の項で構成される。
モデルに含める変数で構成される行列とそのサイズ。
それ以外の変数の行列とそのサイズ。
と表記する。
モデル：
パラメータベクトルのサイズは以下。
実験データに実際にあてはめるモデルはモデルに含める変数で構成されたもの。
最小二乗法による推定値は以下になる。
この推定値の期待値を計算する。
応答の期待値（E(Y)）は、モデルに含めた変数から構成されるものからだけではなく、モデルに含めなかった変数も含まれる。要するに、真のモデルは①であるが、スクリーニングを目的とした実験なので、あてはめに用いるモデルは②といったイメージである。
上の計算より、
があてはめに用いたモデルパラメータへの推定値にバイアスを与えている項となる。
あてはめに用いたモデルパラメータ推定値が、モデルに含めなかった項から受けるバイアスと考えることができる。
と置く。
サイズは以下。
よって、本当はモデルに含めていない変数が応答へ影響があった場合（β2の要素が0でない場合）、あてはめに用いたモデルパラメータ推定値（β1）は、β2の影響のいくらかの割合をバイアスとして受け取っていることになる。
その影響のいくらかの割合を定量的に表すのが、行列Aとなる。
Aはモデルに含める変数とそれ以外の変数の交絡関係の度合いを定量的に表す行列となる。
ここまでは、以下の文献に取り上げられている。
Optimal Design of Experiments: A Case Study Approach
- 作者: Peter Goos,Bradley Jones
- 出版社/メーカー: Wiley
- 発売日: 2011/08/15
- メディア: ハードカバー
- この商品を含むブログを見る
この行列は、相関係数行列と似ているようで異なる。
重回帰モデルは、変数間の相関が強いデータにあてはめを行うと多重共線性の問題が発生し、パラメータ推定値が不安定になり、応答と因子間の関係把握という実験計画法の目的を達成することができない。
よって、因子間に相関関係がないか、もしくは小さく押さえられているかなどの確認が必要となる。
p個の各変数が標準化（平均0、分散1）されている場合、因子の変数で構成される行列Xの相関係数行列は以下になる。
相関係数行列Rのサイズは以下。
ただし実験計画法では通常、各変数は平均0に中心化されているが、分散は1となっていない。
各変数は -1、0や1の水準を取り、分散は1にならない。
平均は0であるが、分散は1でない場合の相関係数行列は以下。
Dは、p個の変数の分散を対角要素に持つ対角行列である。
と、モデルに含める変数と、それ以外の変数の行列に分割する。
Dも、rまでとr+1からを、2つの対角行列と分割する。
サイズは以下。
右上の行列に注目し、Bと置く。
サイズはAと同じ。
Bは、変数間の相関係数の、モデルに含めた部分とそれ以外の変数の相関係数に対応する部分。
Aは、モデルに含める変数の行列掛け算の逆行列部分がややこしそうであるが、実験計画法による計画では通常、モデルに含める変数間に相関は存在しない（最適計画で無理に作った計画や殆直行計画は別）。
よって、相関がない場合、
は各変数の各要素の2乗和を対角要素にとる対角行列となる。
は、第i変数（n行）の2乗和。
A、B共にサイズはr×(p - r)。これのs行t列の箇所のAとBの値を考える。要するにモデルに含めた変数のs番目とそれ以外の変数のt番目（モデルに含めたものの最初から数えるとr+t番目）の関係ということになる。
Aのs行t列の要素。
Bのs行t列の要素。
は、変数sの2乗和で、s含めすべての変数は平均0に中心化が行われているので、
となる。
となる場合、
とs番目の変数とt番目の変数の分散が等しいとなる。
よって、行列AとBの要素が一致するパターンは、以下。
① ※が0かどうか。0であれば共に0となる
② ※が0でない場合、変数sとtの分散が同じかどうか。同じであれば一致
因子の変数は平均が0に中心化されている、モデルに含める変数間の相関は0としているので、そうでない場合（最適計画を用いるとそういった場合が考えられる）は、もっと関係はややこしくなると思われる。
気になっている分布をいくつか
変数yがカウントデータ（y = 0, 1, 2, 3, ….）のときの分布に関していくつか。
Poisson分布
http://ja.wikipedia.org/wiki/%E3%83%9D%E3%82%A2%E3%82%BD%E3%83%B3%E5%88%86%E5%B8%83
パラメータλのみによって形状が定まる。
平均と分散は等しい
E(y) =λ
Var(y) =λ
λ=1,2,3,4,5と変化させつつ、y = 0,1,2,3,…,10のときの発生率をプロット。
lam <- c(1,2,3,4,5)
y <- 0:10
# 行列に値を格納する。行がλ=1,2,3,4,5、列がy=0,1,2,...,10。
Res <- 1:(length(lam)*length(y))
dim(Res) <- c(length(lam), length(y))
for(i in 1:length(lam)){
for(j in 1:length(y)){
# dpois()はPoisson分布の確率密度を返す関数。
Res[i,j] <- dpois(y[j], lambda=lam[i])
}
}
# 重ね合わせプロット。λ=1のときのプロットにλ=2,3,4,5を重ねるのでif()を使用。
for(i in 1:length(lam)){
if(i == 1){
plot(y, Res[i,], type="o", col=i, ylim=c(0,0.45), main="Poisson分布, λ=1,2,3,4,5",xlab="y", ylab="prob(y)" )
}else{
par(new=TRUE)
plot(y, Res[i,], type="o", col=i, ylim=c(0,0.45), xlab="", ylab="")
}
}
黒；λ= 1 赤：λ= 2 緑：λ= 3 青：λ= 4 水色：λ= 5
Zero-infrated Poisson分布
http://en.wikipedia.org/wiki/Zero-inflated_model
0カウントの発生率を多めに見積もる分布。
Prob(y) = π+ (1-π) f( y ) if y = 0
Prob(y) = (1-π) f( y ) if y = 1, 2, 3, 4,…
f( y )は、Poisson分布などの確率密度分布を示す。
πは、上リンクでは“probability of extra zeros”とされており、y=0の発生率がy>0の場合に比べπだけ高くする、という役割を果たすパラメータ。0から1の間をとり、π= 0とすると、y = 0, 1, 2, 3, 4,… において単純にProb(y) = f( y )となる。
Zero-infrated Poisson分布も上リンクに記載されている。
E(y) =λ(1 - π)
Var(y) =λ(1 - π)(1 + πλ)
と定義される。
Zero-infrated Poisson分布をπ=0.1, 0.2, 0.4. 0.7別にプロットしてみる。
pi1 <- c(0.1, 0.2, 0.4, 0.7)
lam <- c(1,2,3,4,5)
y <- 0:10
# 行列に値を格納する。行がλ=1,2,3,4,5、列がy=0,1,2,...,10。
Res <- 1:(length(lam)*length(y))
dim(Res) <- c(length(lam), length(y))
# さらに、π=0.1,0.2,0.4,0.7別に行列をリストで管理。
for(i in 1:length(pi1)){
if(i == 1){
LRes <- list(Res)
}else{
# append()でリストに要素を追加
LRes <- append(LRes, list(Res))
}
}
# ZI Poisson分布の確率密度を取得。π=0.1,0.2,0.4,0.7のときのそれぞれの値がリストの要素。
for(k in 1:length(pi1)){
for(i in 1:length(lam)){
for(j in 1:length(y)){
# y=0のとき、πだけ発生率が水増しされる。
if(y[j] == 0){
LRes[ [k] ][i,j] <- pi1[k] + (1-pi1[k]) * dpois(y[j], lambda=lam[i])
}else{
LRes[ [k] ][i,j] <- (1-pi1[k]) * dpois(y[j], lambda=lam[i])
}
}
}
}
# π別にプロットを描く。
par(mfrow=c(2,2))
for(k in 1:length(pi1)){
# グラフタイトルを別で作っておく。
title <- paste("ZI Poisson分布 π=", pi1[k] , ", λ=1,2,3,4,5")
for(i in 1:length(lam)){
if(i == 1){
plot(y, LRes[ [k] ][i,], type="o", col=i, ylim=c(0,0.9), main=title,xlab="y", ylab="prob(y)" )
}else{
par(new=TRUE)
plot(y, LRes[ [k] ][i,], type="o", col=i, ylim=c(0,0.9), xlab="", ylab="")
}
}
}
y=0の発生に意図的に重みを置くことになるので、πがある程度大きい場合は、高いy=0の発生から、y=1の発生が一旦落着き、まただんだんと発生が高くなる。
y=0の発生が高い理由が既知の場合に用いる分布であろう。
ガンマPoisson分布（負の二項分布）
過大分散（over-dispersion）が起きているカウントデータに対する分布。
過大分散とは、分散が理論より大きい状態を指す。
Poisson分布は、理論的に平均 = 分散となるが、分散がこの仮定より大きい場合に対応するのがガンマPoisson分布と言われる。
確率密度関数は以下。
E(y) = μ
Var(y) = μ + σ (μ^2)
平均と分散が、σ（過大分散パラメータ）の指定によって独立している。
http://ja.wikipedia.org/wiki/%E8%B2%A0%E3%81%AE%E4%BA%8C%E9%A0%85%E5%88%86%E5%B8%83
負の二項分布の確率関数を、パラメータを置き換え、Combinationを階乗で表現し、階乗をガンマ関数に置き換えると上のガンマPoissonに変形できる。
ガンマPoisson分布をσ=0.01, 1, 2, 4別にプロットしてみる。
mu <- c(1, 2, 3, 4, 5)
sig <- c(0.01, 1, 2, 4)
y <- 0:10
# 行列に値を格納する。行がμ=1,2,3,4,5、列がy=0,1,2,...,10。
Res <- 1:(length(mu)*length(y))
dim(Res) <- c(length(mu), length(y))
# さらに、σ=0.01,1,2,4別に行列をリストで管理。
for(i in 1:length(sig)){
if(i == 1){
LRes <- list(Res)
}else{
LRes <- append(LRes, list(Res))
}
}
# Gamma Poissonを関数（GP）を定義しておく。μ、σ、yが引数で、確率密度を返す。標準関数で準備されているかは不明。
GP <- function(mu1, sig1, y1){
f <- ( gamma(y1+(1/sig1)) / (gamma(y1+1)*gamma(1/sig1)) )*( (mu1*sig1)^y1 / (1+mu1*sig1)^(y1+(1/sig1)) )
return(f)
}
# ガンマPoisson分布の確率密度を取得。σ=0.01,1,2,4のときのそれぞれの値がリストの要素。
for(k in 1:length(sig)){
for(i in 1:length(mu)){
for(j in 1:length(y)){
LRes[ [k] ][i,j] <- GP(mu[i], sig[k], y[j])
}
}
}
# σ別にプロットを描く。
par(mfrow=c(2,2))
for(k in 1:length(sig)){
title <- paste("Γ Poisson分布 σ=", sig[k] , ", μ=1,2,3,4,5")
for(i in 1:length(mu)){
if(i == 1){
plot(y, LRes[ [k] ][i,], type="o", col=i, ylim=c(0,0.8), main=title,xlab="y", ylab="prob(y)" )
}else{
par(new=TRUE)
plot(y, LRes[ [k] ][i,], type="o", col=i, ylim=c(0,0.8), xlab="", ylab="")
}
}
}
σを小さくした場合（0.01）通常のPoisson分布に非常に近くなる。
また、σを大きくすると、Zero-infrated Poisson分布のようにy=0の発生を強調した分布になるようである。
ランダムフォレスト1
前回の投稿「とりあえず決定木モデルを実行のためのメモ」で、決定木モデルを使ってみた。
http://highschoolstudent.hatenablog.com/entry/2013/07/03/193603
これはどんなモデルかと言うならば、目的変数が数値データの場合（回帰木）、その目的変数の平均値に最も大きく差が出るような分岐を説明変数の中から探す（説明変数自体と分岐の水準）。要するにt検定で最も有意となる分岐を見つけて行くのを繰り返す。目的変数がカテゴリカルデータの場合（分類木）、その目的変数の割合に最も大きく差が出るような分岐を説明変数の中から探す。要するにカイ2乗検定で最も有意となる分岐を見つけて行くのを繰り返す。とイメージしておけば問題ない。
では、決定木の応用であるランダムフォレスト。
ランダムフォレストは、ブートストラップサンプリングにより複数のデータセットを作成し、それら各データに決定木をあてはめ、すべての結果を平均（カテゴリカル目的変数の場合は多数決の場合も）して最終結果とする手法（バギングと呼ばれる手法）。ただし、分岐毎に候補変数がランダムに選ばれる。たとえば、説明変数が10あり、決定木を作成する際、分岐の候補となる変数が3つランダムに選ばれ、その中から目的変数をよく区別する一つの変数と水準が決定される。
ランダムフォレストは、単純な決定木より予測精度が高いと言われる。
これは、沢山の決定木の予測値を集めて平均することにより、予測値の分散が小さくなる、と説明される。決定木は、ある程度分岐させるとhigh variance、low biasなモデル（要するにオーバーフィッティング気味なモデル）になるので、沢山の結果の平均をとることによって分散が小さくなり、予測精度が向上する。
（メモ： バギングは平均を取ることにより最終結果のvarianceを減らせる。ただし、biasは平均を取るために使われる各モデルのbiasの期待値と変わらない。対し、ブースティングはbiasを取り除きながら最終的な予測結果を得る手法。）
ただし、平均値の分散は、平均値を構成する統計量（予測値）間の相関にも影響を受けるので相関も小さく抑える工夫が必要になってくる。
Cov(X1, X2) = r
Var(X1) = Var(X2) = v
μ= X1 + X2 / 2
Var(μ) = 2v + 2r / 4 = v + r / 2 ← rを小さく抑える必要がある。
ランダムフォレストでは、分岐のたびに候補となる変数をランダムに選択することにより、似すぎたツリーが生成されるのを防いでいる、ということだと思われる。
次回、ランダムフォレスト、randomForest()、を使ってみる。
とりあえず決定木モデルを実行のためのメモ
Rでの決定木分析（分類木、回帰木）の実行に関して、こうではないかとのメモを記す。
CARTアルゴリズムによる決定木分析を行うパッケージはrpartとmvpartがある。
mvpartはrpartに機能拡張を加えた上位パッケージとのこと。
2重に読み込んだ場合、以下のメッセージが出力される。
library(rpart)
library(mvpart)
次のパッケージを付け加えます: 'mvpart'
The following object(s) are masked from 'package:rpart':
meanvar, na.rpart, path.rpart, plotcp, post, printcp, prune, prune.rpart,
rpart, rpart.control, rpconvert, rsq.rpart, snip.rpart, xpred.rpart
決定木モデルの作成にはrpart()関数を用い、plot()やplotcp()関数などで結果を確認する。mvpartパッケージを使用した場合（上書きでも、単独でも）、rpartパッケージの出力（主に見た目）が異なる箇所がある。
上のメッセージはrpartパッケージで定義されている（表示されている）各関数が、後で読み込んだmvpartパッケージ中にも含まれています、ということだと思う。
それで、後から読み込んだものが適用されるのだと思われる。
が、mvpartのマニュアルとrpartのマニュアルの中のそれらの関数の説明は内容が違うように見えない…。
決定木モデルは、分岐数の多さによってモデルの複雑さが決定される。
回帰木（目的変数が連続値）では；
グループAの平方和とグループBの平方和の和が最小になるように、分岐に使用する変数と分割点を決定する。
分類木（目的変数がカテゴリカル）では；
ジニ係数（Gini Index）、もしくはエントロピー（Entropy, Cross-entropy）を最大化するように、分岐に使用する変数と分割点を決定する。
rpart()関数の引数にparms=list(split="gini")と指定するとジニ係数が用いられ、parms=list(split="information")とするとエントロピーが用いられるようである。
また、Loss行列なども指定できるらしい。
どこまで細かく分割するかはrpart() 関数の引数で指定できる。
CARTは、指定した条件まで木を一旦分割し、その後適当な大きさに剪定し、最終的なモデルを決定する。
データの作成。
y <- rnorm(100)
x1 <- y + 0.5*rnorm(100)
x2 <- y + 0.75*rnorm(100)
x3 <- y + rnorm(100)
Data1 <- data.frame(y,x1,x2,x3)
cor(Data1)
y x1 x2 x3
y 1.0000000 0.8810408 0.8086425 0.6582551
x1 0.8810408 1.0000000 0.6947170 0.5983025
x2 0.8086425 0.6947170 1.0000000 0.5468263
x3 0.6582551 0.5983025 0.5468263 1.0000000
yとの相関がx1>x2>x3となるデータを意図的に作成した。
決定木（回帰木）の実行。
result1 <- rpart( y~., data=Data1, control=rpart.control(minsplit=5, cp=0.01) )
control=rpart.control(minsplit=5, cp=0.01)において、どこまで細かく分岐するかを決定する。
minsplitはマニュアルによると、その枝にあるデータ数が指定した数以下になると、分岐が行われない。
cp(complexity parameter)は、木のサイズとあてはまりの良さをトレードオフする指標で、要するに、木がある大きさのときのモデルの複雑さの指標のようなものと考えられる。小さく指定すると細かく分岐が行われる。詳細はよく分かっていない。
result1
n= 100
node), split, n, deviance, yval
* denotes terminal node
1) root 100 82.74123000 -0.08323216
2) x1< -0.118119 50 20.52873000 -0.76585870
4) x2< -1.121383 18 5.31604000 -1.36086100
8) x1< -1.620781 5 0.99465230 -2.05751900
16) x1< -2.258496 2 0.05018492 -2.55654600 *
17) x1>=-2.258496 3 0.11437500 -1.72483500 *
9) x1>=-1.620781 13 0.96139660 -1.09291600 *
5) x2>=-1.121383 32 5.25564900 -0.43116960
10) x2< 0.3988502 27 3.16579100 -0.53866160
20) x1< -0.5092024 12 1.31139200 -0.73501100 *
21) x1>=-0.5092024 15 1.02165200 -0.38158210 *
11) x2>=0.3988502 5 0.09323544 0.14928700 *
3) x1>=-0.118119 50 15.61460000 0.59939430
6) x1< 1.336195 40 8.11672700 0.42480680
12) x2< 1.20844 34 4.89382400 0.32000810
24) x3< 0.4031585 21 1.99932900 0.16314720 *
25) x3>=0.4031585 13 1.54309500 0.57339870 *
13) x2>=1.20844 6 0.73347440 1.01866600 *
7) x1>=1.336195 10 1.40171400 1.29774500 *
詳細な結果はsummary(result1)で確認できる。
10のグループに分けられた。
木のプロットを書くと理解しやすい。
plot(result1)
text(result1, use.n=TRUE, cex=0.5)
引数cexで、ラベルの大きさを調整できる。
ここで、path.rpart(result1)を実行すると、上プロット上にマウスを持っていくと＋マークになる。ノードの箇所でクリックすると分岐履歴が確認できるので面白い。
また、分岐にどの変数が多く使われているかを見るとやはりx1である、これはyとの相関をx1>x2>x3としたデータであるためである。
変数重要度のような指標を出力する関数があるかは、不明。
モデルの複雑性を確認し、剪定する箇所を決定し、最終的なモデルを作成する。
cp統計量を確認。
printcp(result1)
Regression tree:
rpart(formula = y ~ ., data = Data1, control = rpart.control(minsplit = 5,
cp = 0.01))
Variables actually used in tree construction:
[1] x1 x2 x3
Root node error: 82.741/100 = 0.82741
n= 100
CP nsplit rel error xerror xstd
1 0.563176 0 1.00000 1.01980 0.140235
2 0.120340 1 0.43682 0.46874 0.069506
3 0.073677 2 0.31648 0.40655 0.048833
4 0.040608 3 0.24281 0.37062 0.046885
5 0.030087 4 0.20220 0.36497 0.046381
6 0.024131 5 0.17211 0.34970 0.045379
7 0.016333 6 0.14798 0.30970 0.039702
8 0.010064 7 0.13165 0.28018 0.038337
9 0.010032 8 0.12158 0.27143 0.036761
10 0.010000 9 0.11155 0.27143 0.036761
プロットも行う。
plotcp(result1)
k分割交差検証法による結果である様子。ちなみに、kの数は、control=rpart.control(xval=k)で指定できるそう。
結果の出力はいまいち不明。XerrorはRMSE (Root Mean Square Error：誤差の標準偏差)かな、と思うのだが、xstdがなにかも分からない。rel errorもいまいち不明..。
プロットは、黒い下のラインが結果の出力のrel errorと一致する。上の青いラインが検証用データによるものだと推測される。
縦軸（X-Val Relative Error）が0.3の箇所に引かれている横線がMin+1SEとのことで、最大分岐モデルのxerror+1*xstdで、0.27143 + 0.036761 = 0.308191とのこと。この横線の下側に最も近い検証用データによる結果の分岐数にモデルを決定するのが基準とのこと。詳細はますます不明..。
では、cp=0.010032、分岐数が8にモデルを剪定する。
prune() 関数を用いるとのこと。
result2 <- prune(result1, cp=0.01005)
plot(result2)
text(result2, use.n=TRUE, cex=0.5)
予測は、predict(result2, 予測データ) として行う。
と、不明な点は多いのだが、とりあえずお手本に従った使い方をしてみた。
今回はCARTを使ったが、決定木モデルはCHAIDなど含めいろいろと派生型がある。
個人的には、JMPソフトウェアの「分岐」、「剪定」ボタンで実行できる対話的な方式も使い勝手が良いと思う。
決定木であるが、長所は、モデルが単純な分岐の繰り返しなので解釈が容易、と言われる。
こういった、グループ分けの解釈を決定木に求めるのであればJMPの対話形式は非常に良い。
ただし、あてはまりの良いレベルまで木を分岐させると、たいていの場合、解釈が簡単でないまで木が成長してしまう..。
それに、決定木モデルは不安定とも言われる。あてはめを行うそのデータ（学習データ）にとても依存するモデルで、外れ値などの影響で、ある変数がたまたま上位で分岐されると、その下位の分岐は、その分岐のもとに成り立ってしまうことになる。
この不安定性をバギングと分岐毎の分岐の候補変数の無作為抽出によって克服し、予測精度を高めた決定木モデルの応用がランダムフォレストである。
ランダムフォレストによる予測精度の向上だけではなく、前回（「重回帰と変数選択」）でも触れたが、沢山変数がある場合の重要な変数の抽出に関しても次回以降触れたい。
Rによるデータサイエンス - データ解析の基礎から最新手法まで
- 作者: 金明哲
- 出版社/メーカー: 森北出版
- 発売日: 2007/10/13
- メディア: 単行本（ソフトカバー）
- 購入: 36人 クリック: 694回
- この商品を含むブログ (60件) を見る
重回帰と変数選択
重回帰とその変数選択に関するメモ。
RにはF値による検定でのステップワイズ変数選択法はないのでしょうか？
AICによる変数選択法である、step() 関数などは見当たるが検定による変数選択は見当たらない。Web検索でもヒットしない。
AICによる変数選択と検定による変数選択、どちらが良いとかそういうのは置いておいて、検定で行った時の変数の取り込みに関する基本的性質に関してメモしておきたい。
AICは非常に興味深い指標なので、そのうち詳しく勉強できたらいいと思う。
- 作者: 樺島祥介,北川源四郎,甘利俊一,赤池弘次,下平英寿,土谷隆,室田一雄
- 出版社/メーカー: 共立出版
- 発売日: 2007/07/06
- メディア: 単行本
- 購入: 4人 クリック: 74回
- この商品を含むブログ (12件) を見る
統計の専門家でなければ、AICの原理や思想を理解している人はいないと思う。実務家であれば、分かりやすい検定によって変数を選択している場合も多いと思う。
では、検定を用いた場合の変数選択法（増加法）における変数が選ばれる順序を見ていく。
疑似データを作成。
y <- rnorm(n=50,mean=0,sd=5)
x1 <- y + rnorm(n=50,mean=0,sd=3)
x2 <- x1 + rnorm(n=50,mean=0,sd=3)
x3 <- y + rnorm(n=50,mean=0,sd=10)
標準化しておく。
y <- scale(y)
x1 <- scale(x1)
x2 <- scale(x2)
x3 <- scale(x3)
相関を見る。
X <- data.frame(y,x1,x2,x3)
cor(X)
y x1 x2 x3
y 1.0000000 0.8730817 0.7581724 0.4387364
x1 0.8730817 1.0000000 0.8529398 0.3161180
x2 0.7581724 0.8529398 1.0000000 0.2672547
x3 0.4387364 0.3161180 0.2672547 1.0000000
yとx1の相関が一番高く、次にyと相関が高いのはx2、ただしx1とx2との相関も高い。 x3はyとの相関が一番低い、そしてx1、x2との相関も低くなっている。
x1、x2、x3個別にみると、x1 → x2 → x3の順番でyに対して影響（相関）が大きい。
しかし、重回帰モデルにすべて取り込んだ場合、その影響度の順序は異なる。
結論から言うと、x1 → x3 → x2の順番となる。
full <- lm(y~x1+x2+x3)
summary(full)
Call:
lm(formula = y ~ x1 + x2 + x3)
Residuals:
Min 1Q Median 3Q Max
-1.25633 -0.24685 0.01037 0.28646 1.03247
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.456e-17 6.650e-02 0.000 1.000
x1 7.723e-01 1.307e-01 5.908 3.97e-07 ***
x2 5.107e-02 1.287e-01 0.397 0.693
x3 1.809e-01 7.081e-02 2.555 0.014 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4702 on 46 degrees of freedom
Multiple R-squared: 0.7924, Adjusted R-squared: 0.7789
F-statistic: 58.53 on 3 and 46 DF, p-value: 9.728e-16
x2の係数は小さい。（p値が大きい）
個別に結果をみると、以下のようにx1 → x2 → x3の順番でyに対して影響度が高い。
y_x1 <- lm(y~x1)
summary(y_x1)
Call:
lm(formula = y ~ x1)
Residuals:
Min 1Q Median 3Q Max
-1.30111 -0.31311 0.01479 0.33869 1.08676
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.887e-17 6.967e-02 0.00 1
x1 8.731e-01 7.038e-02 12.41 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4926 on 48 degrees of freedom
Multiple R-squared: 0.7623, Adjusted R-squared: 0.7573
F-statistic: 153.9 on 1 and 48 DF, p-value: 2.2e-16
y_x2 <- lm(y~x2)
summary(y_x2)
Call:
lm(formula = y ~ x2)
Residuals:
Min 1Q Median 3Q Max
-1.8576 -0.4300 0.0155 0.4407 1.3690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.831e-17 9.317e-02 0.000 1
x2 7.582e-01 9.412e-02 8.056 1.81e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.6588 on 48 degrees of freedom
Multiple R-squared: 0.5748, Adjusted R-squared: 0.566
F-statistic: 64.89 on 1 and 48 DF, p-value: 1.814e-10
y_x3 <- lm(y~x3)
summary(y_x3)
Call:
lm(formula = y ~ x3)
Residuals:
Min 1Q Median 3Q Max
-1.5212 -0.5209 -0.1517 0.5198 2.4679
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.116e-17 1.284e-01 0.000 1.00000
x3 4.387e-01 1.297e-01 3.383 0.00144 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.9079 on 48 degrees of freedom
Multiple R-squared: 0.1925, Adjusted R-squared: 0.1757
F-statistic: 11.44 on 1 and 48 DF, p-value: 0.001437
変数選択を行った場合、まずyと一番相関が強いx1が取り込まれる。
x1のyに対する回帰を行い、残差を求める。①
resi.y_x1 <- y_x1$residuals
ようするに、x1がyを説明した部分を取り除いたことになる。
x1のx2、x3に対する回帰を行う、残差を求める。
x1のx2に対する回帰と残差。②
x2_x1 <- lm(x2~x1)
resi.x2_x1 <- x2_x1$residuals
x1のx3に対する回帰と残差。③
x3_x1 <- lm(x3~x1)
resi.x3_x1 <- x3_x1$residuals
ようするに、x1がx2とx3を説明していた部分を取り除いたということだ。
②の①に対する回帰と③の①に対する回帰を行う。
resi.y_x1_resi.x2_x1 <- lm(resi.y_x1~resi.x2_x1)
summary(resi.y_x1_resi.x2_x1)
Call:
lm(formula = resi.y_x1 ~ resi.x2_x1)
Residuals:
Min 1Q Median 3Q Max
-1.33597 -0.31891 0.03114 0.34418 1.02701
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.679e-17 6.957e-02 0.000 1.000
resi.x2_x1 4.949e-02 1.346e-01 0.368 0.715
Residual standard error: 0.4919 on 48 degrees of freedom
Multiple R-squared: 0.002808, Adjusted R-squared: -0.01797
F-statistic: 0.1351 on 1 and 48 DF, p-value: 0.7148
resi.y_x1_resi.x3_x1 <- lm(resi.y_x1~resi.x3_x1)
summary(resi.y_x1_resi.x3_x1)
Call:
lm(formula = resi.y_x1 ~ resi.x3_x1)
Residuals:
Min 1Q Median 3Q Max
-1.22041 -0.23250 0.02731 0.30354 1.03938
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.758e-17 6.521e-02 0.000 1.0000
resi.x3_x1 1.808e-01 6.944e-02 2.604 0.0122 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4611 on 48 degrees of freedom
Multiple R-squared: 0.1238, Adjusted R-squared: 0.1055
F-statistic: 6.78 on 1 and 48 DF, p-value: 0.01223
③の①に対する回帰の方が係数が大きくなる。変数選択によりモデルにはx3が取り込まれることになる。
ちなみに、上で推定した係数は
lm(y~x1+x2)
Call: lm(formula = y ~ x1 + x2) Coefficients: (Intercept) x1 x2 -5.836e-17 8.309e-01 4.949e-02
と
lm(y~x1+x3)
Call: lm(formula = y ~ x1 + x3) Coefficients: (Intercept) x1 x3 -5.509e-17 8.159e-01 1.808e-01
に等しくなる。
今、モデルにはx1とx3が取り込まれている。
y_x1x3 <- lm(y~x1+x3)
summary(y_x1x3)
Call:
lm(formula = y ~ x1 + x3)
Residuals:
Min 1Q Median 3Q Max
-1.22041 -0.23250 0.02731 0.30354 1.03938
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.509e-17 6.590e-02 0.000 1.0000
x1 8.159e-01 7.017e-02 11.628 1.98e-15 ***
x3 1.808e-01 7.017e-02 2.577 0.0132 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.466 on 47 degrees of freedom
Multiple R-squared: 0.7917, Adjusted R-squared: 0.7828
F-statistic: 89.32 on 2 and 47 DF, p-value: < 2.2e-16
x1とx3でyを回帰したときの残差を求める。④
resi.y_x1x3 <- y_x1x3$residuals
x1とx3でx2を回帰し、残差を求める。⑤
x2_x1x3 <- lm(x2~x1+x3)
resi.x2_x1x3 <- x2_x1x3$residuals
上2つの残差で回帰を行う。（目的変数：④、説明変数：⑤）
resi.y_x1x3_resi.x2_x1x3 <- lm(resi.y_x1x3~resi.x2_x1x3)
summary(resi.y_x1x3_resi.x2_x1x3)
Call:
lm(formula = resi.y_x1x3 ~ resi.x2_x1x3)
Residuals:
Min 1Q Median 3Q Max
-1.25633 -0.24685 0.01037 0.28646 1.03247
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.543e-18 6.510e-02 0.000 1.000
resi.x2_x1x3 5.107e-02 1.260e-01 0.405 0.687
Residual standard error: 0.4603 on 48 degrees of freedom
Multiple R-squared: 0.003412, Adjusted R-squared: -0.01735
F-statistic: 0.1643 on 1 and 48 DF, p-value: 0.687
この係数は、x1、x2、x3をモデルに取り込んだときの、x2の係数に等しい。
lm(y~x1+x2+x3)
Call: lm(formula = y ~ x1 + x2 + x3) Coefficients: (Intercept) x1 x2 x3 -5.456e-17 7.723e-01 5.107e-02 1.809e-01
ある応答に影響のある要因をざっと見つけたいといった場合があったとする。
重回帰モデルを検定による変数選択で作成し、取り込まれた変数が応答に影響がありますと言ってしまうのは少々問題がある。取り込まれた変数はそれらが同時に（重回帰モデルによって）応答をよく説明しているのであって、それら各変数が応答に強い影響（もしくは相関）を持つということではないからである。
では、多数の要因から応答へ影響のあるものを拾い上げたい、しかも相関（すなわち直線関係）といった制約なしで、といった場合、「ランダムフォレスト」はどうだろうと考えている。
使ったことはあるが、そこまで分かっていない、今後の課題としてランダムフォレストを勉強する。
ロジット変換に関して
ロジット変換に関して。
割合を目的変数とする厳密な解析では、ロジスティックモデルをあてはめたり、割合の発生率が小さな場合はポアソンモデルをあてはめたりする。目的変数が二項分布やポアソン分布に従うと仮定するからである。
だが、そういった目的変数に対しても、回帰モデルや分散分析モデルで解析を行う場合は多々あると思う。こういったモデルの方が分かりやすいだとか、ソフトウェアの都合上、こういったモデルの出力結果の方が充実している場合等考えられる。
目的変数は二項分布に従うと考えられるが、それにロジット変換を施して回帰/分散分析モデルで解析を行う場合に関してのメモを記す。
本業で工業系のデータを扱うことも多いので、不良品発生率を目的変数として書く。
不良品発生率の他に、化学品製造だと収率を目的変数とするような場合もある。
不良品発生率が真ん中あたりでばらついている場合（例えば50%）のデータだと変換を行わなくても問題がない場合も十分あると思う。
だが、不良品発生が0%付近に分布しているだとか、収率が100%付近に分布しているだとかの偏った分布の場合、そのまま回帰/分散分析モデルをあてはめると、残差の等分散性が全く満たされていなかったり、非線形な関係だったりする。こういった場合、変換を施すと改善される場合がある。
n個のサンプルを調査し、f個の不良品を発見したとする。
fは二項分布に従うと仮定する。
f ～ Bin( n , π)
π：母不良率
fの期待値は、
E( f ) = n π
fの分散は、
V( f ) = n π(1 - π)
となる。
不良品発生率は、p = f / n と定義できる。
よって、pの期待値は
E( p ) = π
Pの分散は、
V( p ) = π(1 - π) / n
と計算できる。
nがある程度大きいときは、中心極限定理によりpは正規分布に従うとする。
p ～ N( π, π(1 - π) / n )
---------- xのロジット変換 ----------
Logit(x) = ln{ x / (1 - x) }
オッズ：x / (1 - x)
対数オッズ：ln{ x / (1 - x) }
-----------------------------------------------
pにロジット変換を施す。
Logit(p) = ln{ p / (1 - p) }
デルタ法により、Logit(p)の期待値、分散は近似的に以下に計算できる。
E( Logit(p) ) = ln{ π/ (1 - π) }
V( Logit(p) ) = 1 / nπ(1 - π)
正規分布に従うとする。
Logit(p) ～ N( ln{ π/ (1 - π) }, 1 / nπ(1 - π) )
---------- デルタ法 ----------
確率変数xがE(x) =μ、V(x) =σ2のとき、微分可能な関数f()に対して以下が成り立つ。
y = f(x)
E(y) = f(μ)、V(y) = f ’(μ)^2 σ2
------------------------------------
デルタ法に関しては、こちらの資料など参考になると思う。
http://www012.upp.so-net.ne.jp/doi/biostat/CT39/delta.pdf
不良品発生率p をロジスティックモデルでモデル化する。
p = 1 / 1 + exp{ -(a+bx) }
簡単にするため、説明変数1つ（x）の単ロジスティックとする。
pにロジット変換を行うと直線になる。
Logit(p) = a + bx
それぞれ、
p ～ N( π, π(1 - π) / n )
Logit(p) ～ N( ln{ π/ (1 - π) }, 1 / nπ(1 - π) )
の分布に従い、
a、bを任意の値とし、不良品発生率pとそのロジット変換を縦に、横をxとして線を描き、その周りに「1.96×標準偏差」をプロットしてみる。
πの推定値はp。
n数は30とする。
n <- 30
切片、傾きをa = 0、b = 0.5 とする。
a <- 0
b <- 0.5
横軸ベクトル。
x <- -100:100 / 10
ロジスティックモデルでpをモデル化。
p <- 1 / ( 1 + exp(-(a+b*x)) )
pの「1.96×標準偏差」の上側下側。
pU <- p + 1.96 * sqrt( p*(1-p)/n )
pL <- p - 1.96 * sqrt( p*(1-p)/n )
plot( x, p, type="l", xlim=c(-10,10), ylim=c(0,1), main="n=30、両側に標準偏差の1.96倍", xlab="x", ylab="p" )
pの「1.96×標準偏差」の上側下側を重ね描き。par( new=TRUE )でplot()関数の重ね描きができるよう。
par( new=TRUE )
plot( x, pU, type="l", col="red", xlim=c(-10,10), ylim=c(0,1), xlab="", ylab="" )
par( new=TRUE )
plot( x, pL, type="l", col="red", xlim=c(-10,10), ylim=c(0,1), xlab="", ylab="" )
グリッド線を追加。
abline( h=(1:10)/10, v=(-5:5)*2, lty=3 )
pが0.5のときのpの分散が最も大きく、pが0や1に近いほどその分散は小さくなる。
pにロジット変換を行ったものにも同様のプロットを行う。
Logit(p)、直線になる。
logit_p <- a+b*x
Logit(p)の「1.96×標準偏差」の上側下側。
logit_pU <- logit_p + 1.96 * sqrt( 1/(n*p*(1-p)) )
logit_pL <- logit_p - 1.96 * sqrt( 1/(n*p*(1-p)) )
pの場合と同様にプロットを行う。
plot( x, logit_p, type="l", xlim=c(-10,10), ylim=c(-5,5), main="n=30、両側に標準偏差の1.96倍", xlab="x", ylab="Logit(p)" )
par( new=TRUE )
plot( x, logit_pU, type="l", col="red", xlim=c(-10,10), ylim=c(-5,5), xlab="", ylab="" )
par( new=TRUE )
plot( x, logit_pL, type="l", col="red", xlim=c(-10,10), ylim=c(-5,5), xlab="", ylab="" )
abline( h=-5:5, v=(-5:5)*2, lty=3 )
このように、中心付近でのLogit(p)の分散は小さく、端に行くほどLogit(p)の分散は大きくなる。
ちなみに、p = 1 / ( 1 + exp{ -Logit(p) } )なので、Logit(p)とpの対応は以下。
LogitP <- 5:-5
P <- 1 / ( 1 + exp( -LogitP ) )
data.frame(LogitP, P)
LogitP P 1 5 0.993307149 2 4 0.982013790 3 3 0.952574127 4 2 0.880797078 5 1 0.731058579 6 0 0.500000000 7 -1 0.268941421 8 -2 0.119202922 9 -3 0.047425873 10 -4 0.017986210 11 -5 0.006692851
最小2乗法でLogit(p) = a + bxのa、bを推定することは理論的には正しくない。最小2乗法による推定は直線の周りに等分散を仮定している。
やったことはないのだが、重み付最小2乗法といったものが用いられると思う。
でも、実務では最小2乗法でやってしまうこともあると思う。データにもよるが、十分傾向は把握できる場合もある。でも、注意が必要。
最小2乗法は直線の周りでデータのばらつきを平等に扱うので、目的変数pの割合が0や1近く偏っているデータだと、そのばらつきを過小評価した推定になってしまう、ということだと思う。要するに、0や1近くのデータがいっぱい固まっている箇所をうまくモデルで捉えることができていないのでは、ということだと思う。
ちなみに、p = f / nのロジット変換（Logit(p) = ln{ p / (1 - p) }）だとpが0や1だと変換できないので、修正を施し、修正p（p*）をp* = ( f+0.5 ) / ( n+1 ) とし、そのロジット変換を経験ロジットと呼ぶそう。
Logit(p*) = ln{ p* / (1 – p*) } = ln{ ( f+0.5 ) / ( n-f+0.5 ) }
ロジット、経験ロジット変換に関して取り上げられている文献は以下。
- 作者: 永田靖
- 出版社/メーカー: 日科技連出版社
- 発売日: 1992/04
- メディア: 単行本
- 購入: 7人 クリック: 10回
- この商品を含むブログ (1件) を見る
- 作者: 永田靖
- 出版社/メーカー: 日科技連出版社
- 発売日: 1996/10
- メディア: 単行本
- 購入: 27人 クリック: 226回
- この商品を含むブログ (8件) を見る
とりあえずニューラルネットを使ってみる
今後ニューラルネット（以下、NN）について調べていくが、とりあえず、nnet()関数でモデルを作成してみる。
NNであるが、複雑な非線形モデルと言ってよい。説明変数と目的変数が線形な関係でないとき、線形モデル（重回帰、主成分回帰、PLSなど）ではうまく近似できないが、NNを用いると近似することができる。
データはオリジナルなものを使用。説明変数2つ、目的変数1つの簡単な例。
head( NNexample )
X1 X2 Y 1 0.0 0 39.87536 2 0.5 0 39.99583 3 1.0 0 40.11665 4 1.5 0 40.23488 5 2.0 0 40.34901 6 2.5 0 40.45954
str( NNexample )
'data.frame': 121 obs. of 3 variables: $ X1: num 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 ... $ X2: num 0 0 0 0 0 0 0 0 0 0 ... $ Y : num 39.9 40 40.1 40.2 40.3 ...
plot( NNexample )
scatterplot3dパッケージのscatterplot3d()関数を用い、3Dプロットしてみる。
library( scatterplot3d )
scatterplot3d( NNexample, angle=280 )
引数angleで眺める角度の調整が行える。
説明変数（X1、X2）と目的変数（Y）の関係は明らかに非線形である。
後の重ね描きで必要なのでオブジェクトを作成しておく。
plotNN <- scatterplot3d( NNexample, angle=280 )
データに線形モデルである重回帰モデルをあてはめてみる。
目的変数がYで、説明変数がX1とX2。
res_lr <- lm( Y~X1+X2, data=NNexample )
あてはめ結果とプロットを重ね描きしてみる。
scatterplot3d( NNexample,angle=280, main="重回帰であてはめ" )
ドットを重ねる。
plotNN$points3d(NNexample$X1, NNexample$X2, predict(res_lr), col="red", pch=20)
回帰モデルに関しては平面も重ねられるようなのでついでに。
plotNN$plane3d(res_lr, col="red")
このように、線形モデルではこのデータの現象を捉えることはできない。今回は1次式だったが、2次式をあてはめても同様である。
NNをあてはめてみる。nnetパッケージのnnet()関数を使用。
library(nnet)
関数は、nnet(式, データ, size=ユニット数, linout=出力層への関数（TRUE：logistic、FALSE：linear）)と指定する。
引数lineoutは、出力層への最後の関数を線形とするかロジスティックとするかの選択のよう。判別の問題を扱うときはロジスティックとするのだと思われる。今回は、目的変数（Y）は連続なので線形とする。
引数sizeはユニット数であるが、これの最適数の決定方法に関して私は知らない。
今回は説明変数が2つなので、ユニット数が幾つでも大して変わらないような気がするし、3Dプロットを使ってあてはまりが正しいか簡単に確認できると思う。
まあ試しで、データを3:1で学習とテスト用に分け、学習データでユニットが1から15までのモデルを作成し、テストデータで評価してみる。
評価の指標は、予測残差の平方和をテストデータ数で割り一行あたりの平均としたもの。小さいほうが良い。
乱数（二項分布）を使い、0と1を約3:1の割合でフラグを立てる。
flag <- rbinom( nrow(NNexample), 1, 0.25)
NNexample2 <- data.frame( NNexample, flag)
フラグが0を学習データとする。
NNexampleTrain <- subset( NNexample2, flag=="0", c(X1,X2,Y) )
head( NNexampleTrain )
X1 X2 Y 1 0.0 0 39.87536 2 0.5 0 39.99583 4 1.5 0 40.23488 5 2.0 0 40.34901 8 3.5 0 40.68306 9 4.0 0 40.80690
nrow( NNexampleTrain )
[1] 92
フラグが1をテストデータとする。
NNexampleTest <- subset( NNexample2, flag=="1", c(X1,X2,Y) )
head( NNexampleTest )
X1 X2 Y
3 1.0 0.0 40.11665
6 2.5 0.0 40.45954
7 3.0 0.0 40.56929
14 1.0 0.5 42.04056
21 4.5 0.5 42.60309
27 2.0 1.0 42.64893
nrow( NNexampleTest )
[1] 29
各ユニット数での評価結果を格納するオブジェクト。
SSEMean <- 1:15
各ユニット数での結果を繰り返し計算。
predict()関数は学習データ自身の予測値だけでなく、テストデータへの予測値も返せるようである。
for(i in 1:15) {
res_temp <- nnet(Y~X1+X2, data=NNexampleTrain, size=i, linout=TRUE)
pred_temp <- predict( res_temp, NNexampleTest )
SSEMean[i] <- sum( (NNexampleTest$Y - pred_temp)^2 ) / nrow(NNexampleTest)
}
結果。
data.frame(SSEMean)
SSEMean
1 1.529610362
2 1.529610923
3 1.529614121
4 0.333588701
5 0.006297643
6 0.002305100
7 0.006254313
8 0.008703751
9 0.605880589
10 1.530035349
11 0.003378762
12 0.009971504
13 1.529620984
14 0.049907408
15 0.016740579
数回、以上のプログラム（データの分割から各ユニット数でのモデル作成）を試してみた。
ユニット数が5以上もあれば、安定して小さな値になるようである。まあ、今回のデータは外れ値もないきれいなデータなので、もちろん一般的な結論でもなんでもない。
では、ユニット数は5としてモデルを作成する。
res_nn <- nnet(Y~X1+X2, data=NNexample, size=5, linout=TRUE)
あてはめ結果とプロットを重ね描きしてみる。
scatterplot3d(NNexample,angle=280, main="ニューラルネットであてはめ")
ドットを重ねる。
plotNN$points3d(NNexample$X1, NNexample$X2, predict(res_nn), col="blue", pch=20)
良くあてはまっている。
残差を確認してみる。
重回帰とニューラルで、横軸を予測値、縦軸を実測値でプロットする。
par(mfrow=c(1,2))
plot(predict(res_lr), NNexample$Y, main="重回帰", xlab="予測", ylab="実測")
abline(a=0, b=1, col="red", lwd=2)
plot(predict(res_nn), NNexample$Y, main="ニューラル", xlab="予測", ylab="実測")
abline(a=0, b=1, col="red", lwd=2)
重回帰は使えないであろうことが分かり、ニューラルネットの有効性が分かる。
回帰モデルのパラメータ推定に用いられる通常の最小二乗法は「SSE を最小化」する。
NNのパラメータ推定ではペナルティ付き最小二乗法が用いられる。「SSE + λ *（パラメータの二乗和）を最小化」する。
SSEは残差平方和。λ（λ> 0）はペナルティ。
λは正の値なので、ペナルティ付き最小二乗法は、パラメータが大きく推定されるのを避ける手法だということが分かる。
λが大きいと、推定されるパラメータの絶対値は全体的に小さくなり、λが小さいと、推定されるパラメータの絶対値は全体的に大きくなる。
推定するパラメータの絶対値が大きすぎるとオーバーフィットの危険性が増す。NNモデルはあらゆる応答を近似することができ、オーバーフィットの危険性を心配しなくてはいけないので、ペナルティ付き最小二乗法はこれの回避のために導入されていることになる。
パラメータの二乗和は、中間層の数やユニットの数を増やすと増えるので、それに掛け合わされているλは、それ自身の値に意味はないのではかと思う。
λの決定は、λの値を変化させながらクロスバリデーションを行い、予測誤差が最小になる値に決定されるのだと思われる。
λに関してもそうだが、nnet()の出力などいろいろと不明なので、後々調べて行く。
