UsingR の DOTplot もいい感じ

barplot よりは DensityPlot

UsingR パッケージにある DensityPlot が，いい感じ

plot.formula には subset パラメータがある

plot(Petal.Width~Petal.Length, data=iris, subset=iris$Species=="setosa")
のようにできる。
points, lines にも formula 用のメソッドがある

関数の設計ポリシー

mat.or.vec という，つまらない関数が定義されている

> mat.or.vec
function (nr, nc)
if (nc == 1L) numeric(nr) else matrix(0, nr, nc)

どこがつまらないかというと，やること自体がつまらないという関数スペックの問題やプログラムの書き方など，気に入らないことばかり。

1. 関数名が馬鹿馬鹿しい。もっと，スマートな名前にすべき。
2. nc が 1 のときに，長さ nr のベクトルを作るという仕様なのだから，引数にデフォルト定義を付けておかねばいつも第 2 引数に 1 を指定しないといけないという仕様が馬鹿馬鹿しい。function(nr, nc=1L) と定義すべき。
3. 行数節約で if - else を 1 行で書く必要はない。
4. この関数では，要素はいつも 0 に固定されている。全部が 1 とか，1 から始まる連続する整数とか，設定する数値が変えられるようにしておけばよいのに。
5. ベクトルと行列しか作れないのは面白くない。配列も作れるようにしておけばよいのに。

ということで，以下のような関数を定義してみる。

mva <- function (nr, nc=1, value=0, dim=NULL) # matrix, vector, array どれでも作れる関数という命名
{
if (is.null(dim)) { # ベクトルか行列を作る
if (nc == 1L) {
if (value == "seq") {
return(1:nr)
} else {
return(rep(value, nr))
}
} else {
if (value == "seq") {
value <- 1:(nr*nc)
}
return(matrix(value, nr, nc))
}
} else { # 配列を作る
if (value == "seq") {
value <- 1:prod(dim)
}
return(array(value, dim=dim))
}
}

関数の仕様を決めるのは，センスが必要

自分の書く関数はセンスがあると，自慢している(^_^;)

どんな風に動くかは，実行してみるとわかる。

> mva(5, 1) # ベクトルを作る
[1] 0 0 0 0 0

> mva(5, value=10) # 第 2 引数は省略すると 1 ，つまり，ベクトルを作る。要素は全て 10 という値
[1] 10 10 10 10 10

> mva(3, 4, value=1) # 要素が全部 1 の行列を作る
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 1 1 1
[3,] 1 1 1 1

> mva(dim=c(3, 4, 2), value=2) # 要素が全部 2 の 3 次元配列を作る
, , 1

[,1] [,2] [,3] [,4]
[1,] 2 2 2 2
[2,] 2 2 2 2
[3,] 2 2 2 2

, , 2

[,1] [,2] [,3] [,4]
[1,] 2 2 2 2
[2,] 2 2 2 2
[3,] 2 2 2 2

> mva(dim=c(3, 4, 2), value="seq") # 要素が 1 から始まる整数値を取る，3 次元配列を作る
, , 1

[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12

, , 2

[,1] [,2] [,3] [,4]
[1,] 13 16 19 22
[2,] 14 17 20 23
[3,] 15 18 21 24

このような関数であっても，前もって作っておこうという気には，あまり，ならない，だろう，なあ。

オブジェクトの名前を取り除く

unname を使う

> x <- c(a=2,b=3)
> x[1]
a
2
> unname(x[1])
[1] 2

もちろん， names(z) <- NULL でもよいけど

集計表の表頭・表側に名前を付けるには

table 関数を使うなら，dnn で指定するのがよい。
with を使うのが，一番よいかも。
deparse.level で 2 を指定するのもお手軽。
または，xtabs 関数による。

> set.seed(1)
> d <- data.frame(a=sample(4, 100, replace=TRUE), b=sample(5, 100, replace=TRUE))
> table(d$a, d$b)

1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
> table(d$a, d$b, deparse.level=0)

1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
> table(d$a, d$b, deparse.level=1)

1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
> table(d$a, d$b, deparse.level=2)
d$b
d$a 1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
>
> with(d, table(a, b))
b
a 1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
> with(d, table(a, b, deparse.level=0))

1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
> with(d, table(a, b, deparse.level=1))
b
a 1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
> with(d, table(a, b, deparse.level=2))
b
a 1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
>
> table(d$a, d$b, dnn=list("Side", "Head"))
Head
Side 1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
> with(d, table(a, b, dnn=list("Side", "Head")))
Head
Side 1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8
> xtabs(~a+b, d)
b
a 1 2 3 4 5
1 4 1 7 6 2
2 4 7 10 7 4
3 5 4 3 5 4
4 3 8 4 4 8

話が見えないというのが，わからないのだろうねぇ

2*2 の被験者内ANOVAの多重比較TukeyHSDがうまくいきません †

kachan? (2009-09-10 (木) 21:33:00)

2*2の被験者内分散分析（被験者は４条件をすべて体験した。欠損値等なし）をRで行い、主効果、交互作用（interaction）ともに有意差が出ました。
2*2なので主効果の方はいいのですが、交互作用の多重比較が問題です。
平均値を見ると、あるひとつの条件だけが低い値になっていて、これが全有意差の元になっているようです。
この「あるひとつの条件だけが低い値であり、他の３条件の間には有意差なし」ということを示したいため、TukeyHSDを行ったのですが、

TukeyHSD(aov.acc)

以下にエラー UseMethod("TukeyHSD") :
"TukeyHSD" に適用可能なメソッドがありません

というエラーが出て、うまくいきません。
何が悪いのでしょうか。
何とぞよろしくお願いします。

あんたのやったことが，なんもわかんね。
だで，答えることができね。
あんた，それでも，中級者か？

ベクトルと行列の掛け算

R では，ベクトルと行列の掛け算も柔軟性があるので，v がベクトル，M が行列の場合でも，

t(v) %*% M %*% v

のような記述は過剰反応。単に，

v %*% M %*% v

でよい 。

> v <- runif(5)
> M <- matrix(runif(25), 5, 5)
> t(v) %*% M %*% v
[,1]
[1,] 2.508485
> v %*% M %*% v
[,1]
[1,] 2.508485

まあ，以下のようにするのは考え物だが

> sum(colSums(v * M) * v)
[1] 2.508485

実行時間はどうだろうかと見てみると

> system.time(replicate(10000, t(v) %*% M %*% v))
ユーザ システム 経過
0.185 0.002 0.187
> system.time(replicate(10000, v %*% M %*% v))
ユーザ システム 経過
0.052 0.001 0.053
> system.time(replicate(10000, sum(colSums(v * M) * v)))
ユーザ システム 経過
0.376 0.004 0.379

v %*% M %*% v　は，t(v) %*% M %*% v の1/3 程度の計算時間であるが，1 万回やって 0.1 秒の差なんだから，屁みたいなもんだ。どっちだってよい。

array を ベクトルにするには

a <- array(1:24, dim=c(2,3,4))
c(a)

as.vector(a) でもよいのだが，c(a) には，負けるなあ。

3次元以上の集計表における独立性の検定

ネットサーフィン中，見つけた情報

例えば，以下のようになった。

> d <- data.frame(a=sample(3, 100, replace=TRUE), b=sample(3, 100, replace=TRUE), c=sample(3, 100, replace=TRUE))
> xtabs(~.,d)
, , c = 1

b
a 1 2 3
1 5 6 3
2 2 2 3
3 5 6 2

, , c = 2

b
a 1 2 3
1 5 6 2
2 2 5 3
3 1 2 4

, , c = 3

b
a 1 2 3
1 6 4 5
2 1 2 6
3 5 5 2

> x <- xtabs(~.,d)
> summary(x)
Call: xtabs(formula = ~., data = d)
Number of cases in table: 100
Number of factors: 3
Test for independence of all factors:
Chisq = 16.6219, df = 20, p-value = 0.67737
Chi-squared approximation may be incorrect

lower.tail=FALSE の効力

pnorm, pchisq, pt, pf では，上側確率が欲しいことのほうが圧倒的に多い。
上側確率を求めるとき，1- pnorm(1.96) などとコーディングするヒトがいるが，これは余り好ましくない。pnorm(1.96, lower.tail=FALSE) とすべし。

例証プログラム

ans <- sapply(1:40, function(i) {
x2 <- 1.1^i
p1 <- pnorm(x2, lower.tail=FALSE)
p2 <- 1-pnorm(x2)
dif <- abs(p1-p2)
err <- dif/p1
return(c(x2, p1, p2, dif, err))
})
ans <-as.data.frame(t(ans))
colnames(ans) <- c("Z値", "望ましい方", "望ましくない方", "絶対誤差", "相対誤差")
ans

結果

検定結果というだけならもちろん両者に差はない。しかし，それ以外の目的で使用する場合には大いに差がある。

Z値 望ましい方 望ましくない方 絶対誤差 相対誤差
1 1.1000000 1.3566606e-01 1.3566606e-01 2.7755576e-17 2.0458747e-16
2 1.2100000 1.1313945e-01 1.1313945e-01 5.5511151e-17 4.9064365e-16
3 1.3310000 9.1594505e-02 9.1594505e-02 4.1633363e-17 4.5453997e-16
4 1.4641000 7.1583314e-02 7.1583314e-02 2.7755576e-17 3.8773807e-16
5 1.6105100 5.3643282e-02 5.3643282e-02 3.4694470e-17 6.4676262e-16
6 1.7715610 3.8233729e-02 3.8233729e-02 6.9388939e-18 1.8148619e-16
7 1.9487171 2.5664609e-02 2.5664609e-02 5.2041704e-17 2.0277614e-15
8 2.1435888 1.6032924e-02 1.6032924e-02 2.0816682e-17 1.2983709e-15
9 2.3579477 9.1881411e-03 9.1881411e-03 3.6429193e-17 3.9648056e-15
10 2.5937425 4.7468786e-03 4.7468786e-03 5.1174343e-17 1.0780630e-14
11 2.8531167 2.1646359e-03 2.1646359e-03 2.1250363e-17 9.8170607e-15
12 3.1384284 8.4928220e-04 8.4928220e-04 2.3527187e-17 2.7702438e-14
13 3.4522712 2.7794433e-04 2.7794433e-04 5.3288537e-17 1.9172377e-13
14 3.7974983 7.3081879e-05 7.3081879e-05 1.4731597e-17 2.0157660e-13
15 4.1772482 1.4752852e-05 1.4752852e-05 4.0352650e-18 2.7352439e-13
16 4.5949730 2.1640267e-06 2.1640267e-06 3.5381837e-17 1.6350000e-11
17 5.0544703 2.1579347e-07 2.1579347e-07 4.2870720e-17 1.9866551e-10
18 5.5599173 1.3495125e-08 1.3495125e-08 1.3068143e-17 9.6836034e-10
19 6.1159090 4.8004038e-10 4.8004034e-10 4.0892510e-17 8.5185562e-08
20 6.7274999 8.6301447e-12 8.6300966e-12 4.8084839e-17 5.5717303e-06
21 7.4002499 6.7964202e-14 6.7945649e-14 1.8552996e-17 2.7298188e-04
22 8.1402749 1.9719065e-16 2.2204460e-16 2.4853958e-17 1.2604025e-01
23 8.9543024 1.7094521e-19 0.0000000e+00 1.7094521e-19 1.0000000e+00
24 9.8497327 3.4363263e-23 0.0000000e+00 3.4363263e-23 1.0000000e+00
25 10.8347059 1.1786076e-27 0.0000000e+00 1.1786076e-27 1.0000000e+00
26 11.9181765 4.7584337e-33 0.0000000e+00 4.7584337e-33 1.0000000e+00
27 13.1099942 1.4431026e-39 0.0000000e+00 1.4431026e-39 1.0000000e+00
28 14.4209936 1.9090303e-47 0.0000000e+00 1.9090303e-47 1.0000000e+00
29 15.8630930 5.7065920e-57 0.0000000e+00 5.7065920e-57 1.0000000e+00
30 17.4494023 1.7392690e-68 0.0000000e+00 1.7392690e-68 1.0000000e+00
31 19.1943425 2.0633595e-82 0.0000000e+00 2.0633595e-82 1.0000000e+00
32 21.1137767 2.9714348e-99 0.0000000e+00 2.9714348e-99 1.0000000e+00
33 23.2251544 1.2683379e-119 0.0000000e+00 1.2683379e-119 1.0000000e+00
34 25.5476699 2.9139943e-144 0.0000000e+00 2.9139943e-144 1.0000000e+00
35 28.1024368 4.5734631e-174 0.0000000e+00 4.5734631e-174 1.0000000e+00
36 30.9126805 4.0342157e-210 0.0000000e+00 4.0342157e-210 1.0000000e+00
37 34.0039486 9.7383892e-254 0.0000000e+00 9.7383892e-254 1.0000000e+00
38 37.4043434 1.6605441e-306 0.0000000e+00 1.6605441e-306 1.0000000e+00
39 41.1447778 0.0000000e+00 0.0000000e+00 0.0000000e+00 NaN
40 45.2592556 0.0000000e+00 0.0000000e+00 0.0000000e+00 NaN

なぜ，こんな簡単な質問が出るのだろうかなぁ

中身を変えて同じ関数を何度も繰り返し実行したい †

tefu? (2009-09-07 (月) 21:17:18)

初歩的な質問で申し訳ありません。

ある関数function(x)のxに、毎回違う数字を入れて何度も繰り返したい場合、どのような方法があるでしょうか？

挿入する数字に規則性はなく、数百回単位で繰り返したいのですが、効率的な方法はありますか？

POSIXlt をデータフレームに追加する方法

石田さん
http://cms.ias.tokushima-u.ac.jp/index.php?R_Primitive

RjpWiki初心者 Q&A に strptime() で作成したリスト（ベクトルではない）を既存の data.frame の要素として追加できないという話題がある．POSIXlt をそのまま代入しようとしても，これはリストであり，その要素数は 9 のリストなのでエラーとなる．

> (t <- strptime(paste("1-1-2001 9", 1:20, sep=":"),
"%m-%d-%Y %H:%M"))
> d1 <- data.frame(a=1:20)
> d1$b <- t
以下にエラー `$<-.data.frame`(`*tmp*`, "b", value = list(sec =
c(0, 0, 0, :
replacement has 9 rows, data has 20

d1 <- data.frame(d1, t) とすればよいようですよ。

if - else

google の R コーディング規則で，いろいろ言われているけど，結論としては

コーディング規則は，統一できるようなものではない（組織としての規則は別）

ということに尽きるだろう。

今回も言われてる，if - else で

a <- 2
if (a == 1) {
4
} else {
6
}

のようにすべし。なぜなら，関数内の記述では

a <- 2
if (a == 1) {
4
}
else {
6
}

も可能なのだが，デバッグのためなどでこの部分を切り取ってコンソールにペーストするとエラーになるというのが根拠になっている。

私は，きれいさの点で，後者を推薦する。

後者はコンソールから入力すると確かにエラーになるけど，そんな場合には，前後を { } でくくってやればよいだけ。

具体的には，まずコンソールに { を入力して，次に，ペーストして，最後に } を入力して，リターンキーを押すだけで問題は解消される。

それくらいの手間は掛けても，バチはあたらないだろう。

{
a <- 2
if (a == 1) {
4
}
else {
6
}
}

インタープリタで動いているのでやむを得ないのだ。

a <- 2+3
+5+6

を入力して，　a が 16 にならないとわめいているのと同じ。

a <- 2+3+
5 + 6

とすれば，問題ないわけだ。要するに，インタープリタの立場に立って，入力するしかないのだ。

大人になりなさい。

なんか，同じ人が釣れているんだろうなあ（大漁祝い）

↑の人、初級者コースと同じ人かどうか知らないけど、厭味だけしか書かないならやめてほしいんですけど。しかもご丁寧に色づけまでして・・・。この問題は別にコードなくても通じるし。で、constrOptim hessianをキーワードにググるとメーリングリストにそのもののトピックがいくつかあるので参考にしてみてください。nlme::fdHessが使えるらしいけど、制約付きの場合はちょっとトリッキーなことしてるらしいですよ。あるいはconstrOptimを書き換えて、optimを呼び出してるところにhessian=Tを入れてあげるとヘッセ行列を返すようになるけど、これが正しいのかどうかはよくわからない。 -- 2009-09-03 (木) 23:28:50