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

No.04947　AICでモデルを検討したいのですが　　【SS5】　2007/12/06(Thu) 17:53

ボーリング孔を利用した地盤調査があります。
1m間隔で地上から各深度までの振動の到達時間を取得し，深度／時間から地盤の速度値を推定するものです。
単純に1m間隔で区分することもできますが，野外データで測定誤差が比較的大きいため通常は数点の平均的な勾配から速度を決定します。
以下の測定データを2つのモデルに速度区分したときに，どちらのモデルが妥当かをAICで見積もれないでしょうか。

深度(m)　時間(msec)　モデル1　モデル2
　 0　　　　 0.00　　　　　　　　　　　
　 1　　　　 2.22　　　　　　　　　　415　　
　 2　　　　 4.94　　　　　　　　　(m/sec)　
　 3　　　　 7.20　　　　　　　　　　　
　 4　　　　 9.65　　 486　　　　＿＿＿＿
　 5　　　　11.52　　 (m/sec)　
　 6　　　　13.10　　　　　　　　　　　
　 7　　　　14.61　　　　　　　　　　　
　 8　　　　16.39　　　　　　　　　　695 　
　 9　　　　17.85　　　　　　　　　(m/sec)　
　10　　　　18.97　　＿＿＿＿　　　　　
　11　　　　19.95　　　　　　　　　＿＿＿＿
　12　　　　20.43　　　　　　　　　　　
　13　　　　20.76　　　　　　　　　　　
　14　　　　21.07　　　　　　　　　　　
　15　　　　21.52　　　　　　　　　　　
　16　　　　21.93　　 3082　　　　　3218
　17　　　　22.21　　(m/sec)　　　　(m/sec)
　18　　　　22.39　　　　　　　　　　　
　19　　　　22.60　　　　　　　　　　　
　20　　　　22.88　　　　　　　　　　　

やってみたことは以下の手順です。

1．測定値と理論到達時間から各速度区分内で直線回帰を求める
2．残差平方和をそれぞれ単純に合計する
3．測定データ数×LN(残差平方和の総和)＋2×速度区分数＝AICとする
4．モデル1とモデル2のAICを比較する

恐らくAICの計算方法が間違っていると思いますが，このようなデータの比較は可能でしょうか。
どなたかご存知でしたらご教示頂きますよう宜しくお願い致します。

No.04951　Re: AICでモデルを検討したいのですが　　【青木繁伸】　2007/12/06(Thu) 20:24

示されたデータで何を何で予測するのかよくわからなかったので，別の例を作ってやってみました。データを折れ線で回帰するとき，2本の折れ線と3本の折れ線どっちが良いかというもの。注釈なしで（間違えているかも知れないけど）

> df <- structure(list(x = c(0.292, 1.97, 3.853, 6.091, 7.973, 9.906, 
11.992, 14.027, 15.909, 17.995, 19.979, 21.81), y = c(21.189, 
16.783, 14.44, 8.44, 5.253, 3.284, 3.378, 5.722, 5.628, 10.503, 
16.127, 16.877), model1 = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 2L), .Label = c("part1", "part2"), class = "factor"), 
model2 = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L), .Label = c("part1", "part2", "part3"), class = "factor")), 
.Names = c("x", "y", "model1", "model2"),
row.names = c(NA, -12L), class = "data.frame")
> df
        x      y model1 model2
1   0.292 21.189  part1  part1
2   1.970 16.783  part1  part1
3   3.853 14.440  part1  part1
4   6.091  8.440  part1  part1
5   7.973  5.253  part1  part2
6   9.906  3.284  part1  part2
7  11.992  3.378  part2  part2
8  14.027  5.722  part2  part2
9  15.909  5.628  part2  part3
10 17.995 10.503  part2  part3
11 19.979 16.127  part2  part3
12 21.810 16.877  part2  part3

> lm.model1 <- lm(y~model1*x, df)
> plot(df$y~df$x, ylim=c(0,25))
> points(lm.model1$fitted.values~df$x, col=2)
> lm.model2 <- lm(y~model2*x, df)
> points(lm.model2$fitted.values~df$x, col=4, pch=19)
> summary(lm.model1)

Call:
lm(formula = y ~ model1 * x, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5094 -0.6978  0.1403  0.7571  1.8694 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    21.1202     1.0330   20.45 3.43e-08 ***
model1part2   -36.9056     3.1118  -11.86 2.34e-06 ***
x              -1.9057     0.1715  -11.11 3.85e-06 ***
model1part2:x   3.4094     0.2414   14.12 6.14e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.401 on 8 degrees of freedom
Multiple R-Squared: 0.9628,	Adjusted R-squared: 0.9489 
F-statistic: 69.04 on 3 and 8 DF,  p-value: 4.637e-06 

> AIC(lm.model1)
[1] 47.28954
> summary(lm.model2)

Call:
lm(formula = y ~ model2 * x, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2162 -0.7945 -0.1353  0.9566  1.7186 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    21.6534     1.1678  18.542 1.59e-06 ***
model2part2   -18.1154     3.5378  -5.120 0.002177 ** 
model2part3   -47.4519     5.9533  -7.971 0.000208 ***
x              -2.1106     0.3124  -6.757 0.000513 ***
model2part2:x   2.1900     0.4317   5.073 0.002282 ** 
model2part3:x   4.1230     0.4376   9.423 8.12e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.349 on 6 degrees of freedom
Multiple R-Squared: 0.9741,	Adjusted R-squared: 0.9526 
F-statistic:  45.2 on 5 and 6 DF,  p-value: 0.0001102 

> AIC(lm.model2)
[1] 46.93006
> anova(lm.model1, lm.model2)
Analysis of Variance Table

Model 1: y ~ model1 * x
Model 2: y ~ model2 * x
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1      8 15.7130                           
2      6 10.9266  2    4.7864 1.3142 0.3363

白丸が実測値。赤丸が2本の折れ線で回帰（model1），青丸が3本の折れ線で回帰（model2）。

No.04957　Re: AICでモデルを検討したいのですが　　【SS5】　2007/12/07(Fri) 10:48

さっそく返信頂きましてありがとうございました。
いただいた内容についてはほとんど理解できてないのですが，ご指摘の方法で解決できそうです。
まずは先生の例で内容をトレースしたいのですが，Rを使えばできるのでしょうか。

No.04960　Re: AICでモデルを検討したいのですが　　【SS5】　2007/12/07(Fri) 15:16

最新のRにアップデートしたら実行できました。
ありがとうございました。