2012/07/31 現在で，まだ用意されていないようだ
“セキュリティ”環境設定でインストールが許可されているのは、Mac App Store と確認済みの開発元からのアプリケーションのみです。
といわれる。システム環境設定を弄ってもよいが，そこに以下の示唆が表示されているので，それに従うのが無難。
開発元が不明なアプリケーションを開きたい場合は、Control キーを押したままアプリケーションアイコンをクリックして“開く”を選択することにより、個別に許可することもできます。
X11 is not included with Mountain Lion, but X11 server and client libraries for OS X Mountain Lion are available from the XQuartz project: http://xquartz.macosforge.org. You should use XQuartz version 2.7.2 or later.
「ダメ出し」シリーズをお送りしておりますが，ほとんどの人は「なんで，そんなこというの？」，「どこに問題があるっていうの？いいんじゃない？」ということだと思います。
そんな事じゃないんです。悪いんです。悪いと思わないことが悪いんです。
とだけ，言っておきましょう。
分からない人は，今後も，ダメなプログラムを書き続けるだけのことでしょう。
嘆かわしい(^_^;)
「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 が間違ってんじゃね?」という意見は，即時却下。
「割り付けがうまく行く確率」にて
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
「上昇トレンドの検定 ヨンクヒール・タプストラ検定」にて
SAGx というライブラリ（パッケージ）は普通にはないので，代替法を探していて，
http://aoki2.si.gunma-u.ac.jp/R/Jonckheere.html
にあるとか，cor.test で method="kendall" とするのと同じであるとかの収穫があった。
で，ブログ著者は，「この本の例をやってみる。」...「本ではものすごい小さいp値になるらしいけど1になる。保留。」としているが，
JT.test(y[1, ], rep(1, 7), alternative="increasing")
というのは，書式を間違えていると思う。だから，変な結果になっただけだろう。実際の所，以下のように確かに P 値はとっても小さな値になる。
> y <- rbind(c(1.06, 1.03, 1.10, 1.66, 1.72, 1.82, 1.91),
+ c(0.98, 1.02, 1.08, 1.58, 1.68, 1.78, 1.89),
+ c(1.09, 1.10, 1.11, NA, 1.71, 1.83, 1.92))
> y <- c(y)
> g <- rep(1:7, each=3)
> Jonckheere(y, g)
ヨンキー検定
data: y by g
J = 164.5000, E(J) = 85.5000, V(J) = 231.5250, Z-value = 5.1919,
p-value = 1.041e-07
alternative hypothesis: control <= treat-1 <= ... <= treat-n
> cor.test(y, g, method="kendall", alternative="greater")
Kendall's rank correlation tau
data: y and g
z = 5.1919, p-value = 1.041e-07
alternative hypothesis: true tau is greater than 0
sample estimates:
tau
0.8788771
cor.test(x, g, method="kendall") だということ...
http://aoki2.si.gunma-u.ac.jp/R/Jonckheere.html
の例題
> x <- c(
+ 153, 153, 152, 156, 158, 151, 151, 150, 148, 157, # 第 1 群のデータ
+ 158, 152, 152, 152, 151, 151, 157, 147, 155, 146, # 第 2 群のデータ
+ 153, 146, 138, 152, 140, 146, 156, 142, 147, 153, # 第 3 群のデータ
+ 137, 139, 141, 141, 143, 133, 147, 144, 151, 156 # 第 4 群のデータ
+ )
> g <- rep(1:4, each=10)
> Jonckheere(x, g, alternative="decreasing")
ヨンキー検定
data: x by g
J = 446.5000, E(J) = 300.0000, V(J) = 1705.0776, Z-value = 3.5479,
p-value = 0.0001942
alternative hypothesis: control >= treat-1 >= ... >= treat-n
は
> cor.test(x, g, method="kendall", alternative="less")
Kendall's rank correlation tau
data: x and g
z = -3.5479, p-value = 0.0001942
alternative hypothesis: true tau is less than 0
sample estimates:
tau
-0.4391269
と同じだということ。知らなかった。
「R で Excel っぽい色を出す」にて
excel.like.color <- function(n) {
n <- as.integer(n)
:
residue <- (n - 1) %% 6 + 1
quotient <- floor((n - 1) / 6)
という件があるが，整数定数，整数演算，演算順序を考えると以下の test2 で示すようにするのがよいだろう。
test <- function(n) {
n <- as.integer(n)
residue <- (n - 1) %% 6 + 1
quotient <- floor((n - 1) / 6)
return(c(i, residue, quotient))
}
test2 <- function(n) {
n <- as.integer(n)
quotient <- (n - 1L) %/% 6L
residue <- n - quotient * 6L
return(c(n, residue, quotient))
}
> for (i in 1:15) {
+ cat(test(i), " ", test2(i), "\n")
+ }
1 1 0 1 1 0
2 2 0 2 2 0
3 3 0 3 3 0
4 4 0 4 4 0
5 5 0 5 5 0
6 6 0 6 6 0
7 1 1 7 1 1
8 2 1 8 2 1
9 3 1 9 3 1
10 4 1 10 4 1
11 5 1 11 5 1
12 6 1 12 6 1
13 1 2 13 1 2
14 2 2 14 2 2
15 3 2 15 3 2
「幾何平均」にて
max.g <- 0
min.g <- 0
for(i in 1:n.m){
max.g <- max.g + max(D[[i]])
min.g <- min.g + min(D[[i]])
}
max.g
min.g
当たり前に，
max.g <- sum(sapply(D, max))
min.g <- sum(sapply(D, min))
と書くようにしたいものだ。
「幾何平均」で
AllelesId<-NULL
AllelesId[[1]]<-c("7",9,10,11,12,13,14,15,16,17,18)
AllelesId[[2]]<-c("27",28,28.2,29,30,30.2,30.3,31,31.2,32,32.2,33,33.1,33.2,34,34.2)
AllelesId[[3]]<-c("7",6,9,9.1,10,10.1,10.3,11,12,13,14,15)
AllelesId[[4]]<-c("7",8,9,10,11,12,13,14,15,16)
AllelesId[[5]]<-c("12",13,14,15,16,17,18,19,21)
AllelesId[[6]]<-c("5",6,7,8,9,9.3,10)
AllelesId[[7]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[8]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[9]]<-c("16",17,18,19,20,21,22,23,24,25,26,27,28)
AllelesId[[10]]<-c("9.2",10.2,11,11.2,12,12.2,13,13.2,14,14.2,15,15.2,16,16.2,17.2,18)
AllelesId[[11]]<-c("13",14,15,16,17,18,19,20,21,22)
AllelesId[[12]]<-c("7",8,9,10,11,12,13,14)
AllelesId[[13]]<-c("10",11,12,13,14,15,16,17,17.1,18,19,20,21,22,23,24,25,26,27)
AllelesId[[14]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[15]]<-c("17",18,19,20,21,22,22.2,23,23.2,24,24.2,25,25.2,26,27,28)
などとしているけど，そもそも，15 個の要素を持つ list を作りたいなら，vector("list", 15) とすべきだし，1000 個もあったらどうするつもり？ AllelesId[[1]]～AllelesId[[1000]] まで書くの？
プログラムというのは，「繰り返し」をそのまま書かないで済むようにすべき。だれも，1から1000までの和をとるときに 1+2+3+…+1000 なんて書かないでしょ？それと同じ。
以下のように書けばよいだけ。
AllelesId2 <- list(c("7",9,10,11,12,13,14,15,16,17,18),
c("27",28,28.2,29,30,30.2,30.3,31,31.2,32,32.2,33,33.1,33.2,34,34.2),
c("7",6,9,9.1,10,10.1,10.3,11,12,13,14,15),
c("7",8,9,10,11,12,13,14,15,16),
c("12",13,14,15,16,17,18,19,21),
c("5",6,7,8,9,9.3,10),
c("7",8,9,10,11,12,13,14,15),
c("7",8,9,10,11,12,13,14,15),
c("16",17,18,19,20,21,22,23,24,25,26,27,28),
c("9.2",10.2,11,11.2,12,12.2,13,13.2,14,14.2,15,15.2,16,16.2,17.2,18),
c("13",14,15,16,17,18,19,20,21,22),
c("7",8,9,10,11,12,13,14),
c("10",11,12,13,14,15,16,17,17.1,18,19,20,21,22,23,24,25,26,27),
c("7",8,9,10,11,12,13,14,15),
c("17",18,19,20,21,22,22.2,23,23.2,24,24.2,25,25.2,26,27,28))
なぜ，第1要素だけ文字列にしているのかな？全部文字列にするために全てに " " を付けるのをきらったからでしょう。
これで同じものができたかどうかは，identical 関数で調べれば分かる。
> identical(AllelesId, AllelesId2)
[1] TRUE
「完全に同じだ」ということがわかる。
ついでだからだけど，なぜ，以下のように書かないのかはわからない。単にテストデータだからか？だったらこそ，簡潔に書きたいと思うのでは？
AllelesId3 <- list(c("7",9,10,11,12,13,14,15,16,17,18),
c("27",28,28.2,29,30,30.2,30.3,31,31.2,32,32.2,33,33.1,33.2,34,34.2),
c("7",6,9,9.1,10,10.1,10.3,11,12,13,14,15),
c("7",8:16),
c("12",13:19,21),
c("5",6:9,9.3,10),
c("7",8:15),
c("7",8:15),
c("16",17:28),
c("9.2",10.2,11,11.2,12,12.2,13,13.2,14,14.2,15,15.2,16,16.2,17.2,18),
c("13",14:22),
c("7",8:14),
c("10",11:17,17.1,18:27),
c("7",8:15),
c("17",18:22,22.2,23,23.2,24,24.2,25,25.2,26:28))
当然だが，これも，最初の書き方でできるものと完全に同じだ。
> identical(AllelesId, AllelesId3)
[1] TRUE
ちなみに，【幾何平均をもって、「代表値」とすることの意義はあるのか、あるとしたらどのあたりにあるのか…はまだ不明】ということの真意はよく分からないけど，対数正規分布の代表値は幾何平均だよね。
PVアクセスランキング にほんブログ村
