「統計学関連なんでもあり」の過去ログ--- 044

No.14030　Cｐ統計量について　　　【MINTO】　2010/12/25(Sat) 14:44

現在，変数選択の勉強をしているのですが，変数選択した後の回帰モデルに対する評価をするためにCp統計量というものがあることを知り，Rのwle関数のmle.Cp()を使用し求めたのですが，自分で式にあてはめてみた値と同じような値になることができません。
以下に自分の実行結果を載せますので，どうかご教授お願いします。

            X1          X2         X3         X4
1  -2.53697979 -0.01904025 -1.8160855 -1.5152303
2   1.02631023  0.99582316  0.6921329  0.3423070
3   0.69811795 -0.35955234  0.7016624  1.1583401
4  -1.55305673 -0.76289973 -0.9117866 -0.9277767
5  -0.79343523 -0.69380337 -0.8674818 -1.1402236
6   1.11335083  0.73702726 -0.2243935 -0.2374776
7  -1.22119563  2.18021153 -1.6622586 -1.5455416
8  -0.01800789  1.05301868 -0.4062520 -0.2478122
9   0.25360118 -1.34385384  0.5825237  0.4440714
10  0.41199244  1.48261321  0.3269431  1.1693431
y -5.8925300  3.6247510  2.9950570 -3.7111170 -3.8911180  2.1963660 -2.6793920  -0.1071691 -0.2337557  3.8609350
y2 <- lm(y~c[,1]+c[,2]+c[,3]+c[,4])
(Intercept)       c[, 1]       c[, 2]       c[, 3]       c[, 4]  
     0.2022       1.2354       0.8936       0.5281       1.4610  
y3 <- c[,1]*1.2354+c[,2]*1.2354+c[,3]*0.5281+c[,4]*1.4610+0.2022
z <- y-y3
t(z)%*%z
 3.027479
var(z)
0.3225098
3.027479/0.3225098+(8-10)
7.387247

mle.Cpを使用した場合，5.000となりました。

No.14032　Re: Cp統計量について　　　【青木繁伸】　2010/12/27(Mon) 14:49

あなたが参照したテキストまたはホームページを明示した方が良いでしょう。あなたの勘違いで，計算が間違えているとしても，どこをどのように指摘したらよいか戸惑います。

たとえば，
http://www.statistics4u.info/fundstat_eng/cc_varsel_mallowscp.html
を例に挙げると，SSres は p-1 個の独立変数を使ったときの残差平方和，MSres は，すべての独立変数を使ったときの残差平均平方，N は例数，p は「独立変数の数プラス1」つまり，従属変数も含めて，分析に使った「変数の個数」であることに注意。
あなたの例の場合，独立変数は全部で 4 つなので，MSres は summary(lm(y~c))$sigma^2
4 個の独立変数全部を使ったとき（p = 4+1 = 5）の SSres は sum(residuals(lm(y~c))^2)
N = 15

sum(residuals(lm(y~c))^2) / summary(lm(y~c))$sigma^2 - N + 2 * p
= 1.556342 / 0.3112685 - 10 + 2 * 5 = 5 ではないですか？

以下のようなプログラム断片が参考になるかも。
v は，独立変数行列 x のどの列をモデルに使うかを示すベクトルです。
y <- c(22, 36, 24, 22, 27, 29, 26, 23, 31, 24, 23, 27, 31, 25, 23)
x <- cbind(
x1=c(28, 46, 39, 25, 34, 29, 38, 23, 42, 27, 35, 39, 38, 32, 25),
x2=c(146, 169, 160, 156, 161, 168, 154, 153, 160, 152, 155, 154, 157, 162, 142),
x3=c(34, 57, 48, 38, 47, 50, 54, 40, 62, 39, 46, 54, 57, 53, 32))
MarrowsCp <- function(y, x, v)
{
	all <- lm(y ~ x)
	p <- ncol(x)+1
	SSres <- sum(residuals(all)^2)
	MSres <- SSres/(N-p)
	N <- length(y)
	part <- lm(y ~ x[, v])
	p <- length(v)+1
	SSres <- sum(residuals(part)^2)
	return(SSres/MSres-N+2*p)
}
MarrowsCp(y, x, 1:3)
MarrowsCp(y, x, 1:2)
MarrowsCp(y, x, 1)

No.14039　Re: Cp統計量について　　　【MINTO】　2010/12/27(Mon) 20:34

説明不足にも関わらず，返信していただきありがとうございます。
参考資料は回帰分析の実際という本で，その中にある式を用いました。
式はCp = (SSE)p/σ ̂^2+(2p-n) という式で，(SSE)pはp項からなる回帰式の残差平方和，σ ̂^2は残差分散の推定値，pは変数の個数です。青木先生のプログラムを試し，うまくいくことができました。その結果，私の計算したSSEの値がRと離れたため，Cpの値が違ったのだと思いました。現在はリッジ回帰を用いた変数選択のためにこれらを勉強しているため係数値からCpを求めています。そのため，lm関数は使えず，上記のような式で求めていたのですが，説明足らずですみませんでした。SSEの求め方に何か問題があったのでしょうか？
長々とすみません。もしよければ返信の方，よろしくお願いします。

No.14042　Re: Cp統計量について　　　【青木繁伸】　2010/12/27(Mon) 20:46

色々なことを総合して，特定の状況の場合のも対応できるようにと工夫することはあり得るでしょう。しかし，それが理論的に正しいかどうかは，統計数理科学的に考察しなければなりません。
今回の場合，あなたの理論が正しいかどうか，私には分かりません。

No.14043　Re: Cp統計量について　　　【MINTO】　2010/12/27(Mon) 20:52

返信ありがとうございます不明だった点に関しては，もう一度考え直していきます。いろいろとご教授いただき，ありがとうございました。