あまりにも，ひどすぎるんです

「ダメ出し」シリーズをお送りしておりますが，ほとんどの人は「なんで，そんなこというの？」，「どこに問題があるっていうの？いいんじゃない？」ということだと思います。

そんな事じゃないんです。悪いんです。悪いと思わないことが悪いんです。

とだけ，言っておきましょう。

分からない人は，今後も，ダメなプログラムを書き続けるだけのことでしょう。

嘆かわしい(^_^;)

ダメ出し：間違えてるよ～

「R によるデータの統計的取り扱い/標準誤差」において

したがって、回帰直線の「切片」の95%信頼区間の上限値・下限値は

resultcoefficients[1] + qt(0.975, length(x)-1) * resultcoefficients[3]resultcoefficients[1] - qt(0.975, length(x)-1) * resultcoefficients[3]

同様に、説明変数 x の回帰係数の95%信頼区間の上限値・下限値は

resultcoefficients[2] + qt(0.975, length(x)-1) * resultcoefficients[4]
resultcoefficients[2] - qt(0.975, length(x)-1) * resultcoefficients[4]

で求めることができます。

うそで～～す。qt の自由度は length(x)-2 でなければなりません。

resultcoefficients は result$coefficients ですし。

confint を使うと，一発で答えが出るヨン。

この標準誤差は何のためにあるのでしょうか。答えは、ズバリ、信頼区間の計算のためです。

これも半分ウソ。それだけのためにある訳じゃあない。

===============================

追記 2012/08/01

コメントで，「length(x)-1 でいいんじゃないの?」ともらったので，そうではないということをはっきりさせよう。

元の記事のプログラムでは，以下のようになる。

> x <- c(61, 72, 84, 95, 97, 98, 100, 113, 126, 130)
> y <- c(83, 82, 99, 96, 115, 108, 95, 111, 114, 135)
> d.fr <- data.frame(x, y)
> RegModel <- lm(y ~ x, data = d.fr)
> summary(RegModel)

Call:
lm(formula = y ~ x, data = d.fr)

Residuals:
    Min      1Q  Median      3Q     Max
-10.362  -5.865   0.100   4.024  11.591

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  40.2680    12.3491   3.261 0.011514 *  
x             0.6509     0.1238   5.259 0.000765 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.104 on 8 degrees of freedom
Multiple R-squared: 0.7756,    Adjusted R-squared: 0.7476
F-statistic: 27.66 on 1 and 8 DF,  p-value: 0.0007653

> result <- summary(RegModel)
> result$coefficients[1]
[1] 40.26801
> result$coefficients[2]
[1] 0.6509425
> result$coefficients[3]
[1] 12.34913
> result$coefficients[4]
[1] 0.1237739

切片の信頼区間が

> result$coefficients[1] + qt(0.975, length(x)-1) * result$coefficients[3]
[1] 68.20369
> result$coefficients[1] - qt(0.975, length(x)-1) * result$coefficients[3]
[1] 12.33233

回帰係数の信頼区間が

> result$coefficients[2] + qt(0.975, length(x)-1) * result$coefficients[4]
[1] 0.9309385
> result$coefficients[2] - qt(0.975, length(x)-1) * result$coefficients[4]
[1] 0.3709466

となるとしているのだが，間違いである。それは，confint を使ってみればわかる。

> confint(RegModel)
                 2.5 %     97.5 %
(Intercept) 11.7908555 68.7451654
x            0.3655195  0.9363656

値が違うのは，qt の自由度が length(x)-1 というのが間違えているから。これを length(x)-2 にしてやってみると，

> result$coefficients[1] + qt(0.975, length(x)-2) * result$coefficients[3]
[1] 68.74517
> result$coefficients[1] - qt(0.975, length(x)-2) * result$coefficients[3]
[1] 11.79086
> result$coefficients[2] + qt(0.975, length(x)-2) * result$coefficients[4]
[1] 0.9363656
> result$coefficients[2] - qt(0.975, length(x)-2) * result$coefficients[4]
[1] 0.3655195

となり一致することがわかる。

「confint が間違ってんじゃね?」という意見は，即時却下。

ダメ出し：R の特徴を生かしたプログラミング

「割り付けがうまく行く確率」にて

k <- 5
p <- runif(k)
p <- p/sum(p)

library(gtools)

cmb <- combinations(k,3)

P <- rep(0,length(cmb[,1]))
for(i in 1:length(cmb[,1])){
    P[i] <- prod(p[cmb[i,]])
}

for を使わないようにというのではなく，コンパクトに書こうということ。
コンパクトになるだけでなく，なにを計算しているかがすぐにわかるというメリットがある。
ついでに，length(cmb[,1]) というのは，nrow(cmb)。これも，関数名から何をしようとしているのかすぐにわかる。
以下のようにすれば，前もって P のメモリ確保をする必要もない（メモリ確保も P <- numeric(nrow(cmb)) とすべし）。1 行だけで終わる。

P2 <- apply(cmb, 1, function(x) prod(p[x]))

同じ結果になることも確かめておこう。

> identical(P, P2)
[1] TRUE