R 用の indent 変換例 1

# 同一染色体上にある多型の数Ns
Ns<-4
# Ns個のマーカーのそれぞれのアレルの数をNaS
NaS<-c(3,4,5,2)
# アレル頻度をPaS
PaS<-list()
for(i in 1:Ns){
  PaS[[i]]<-c(rdirichlet(1,rep(1,NaS[i])))
}
NaS
PaS
 
hap1<-hap2<-rep(0,Ns)
for(i in 1:Ns){
  hap1[i]<-sample(1:NaS[i],1,prob=PaS[[i]])
  hap2[i]<-sample(1:NaS[i],1,prob=PaS[[i]])
}
hap1
hap2

これを以下のように変換

# 同一染色体上にある多型の数Ns
Ns <- 4
# Ns個のマーカーのそれぞれのアレルの数をNaS
NaS <- c(3, 4, 5, 2)
# アレル頻度をPaS
PaS <- list()
for (i in 1:Ns) {
    PaS[[i]] <- c(rdirichlet(1, rep(1, NaS[i])))
}
NaS
PaS
 
hap1 <- hap2 <- rep(0, Ns)
for (i in 1:Ns) {
    hap1[i] <- sample(1:NaS[i], 1, prob = PaS[[i]])
    hap2[i] <- sample(1:NaS[i], 1, prob = PaS[[i]])
}
hap1
hap2

あまりにも。。。

> legend(x=par()$usr[2],y=par()$usr[4],legend=colnames(dat),pch=22,lty=0,xjust=-.1
+ ,pt.bg=cm.colors(3),col=1:3,pt.cex=2,x.intersp=2,y.intersp=1.5,title="xpd=TRUE")

プログラムの2行目以降（継続行）が "," で始まるというスタイルは，見たことがない。
-0.1 を -.1 と書くのもいやだなあ。
適切なスペースがないのもひどいなあ。
par の値を取り出すのは，普通は par("width") とか par("user")[2] とかの方がわかりやすい。
小さなことがいろいろ積み重なると，読みにくい（読みたくない）プログラムになる。

Sewave-14

figure チャンクを使って図を挿入するのは記述が面倒くさい。
そこで，以下のような関数を作っておき，results=tex を使って作図した結果を挿入する。

xfigure <- function(src, fn, caption="", label=fn, size=320, width=6, height=6,
                 mgp=c(1.8, 0.6, 0), mar=c(3, 3, 1, 1))
{
    pdf(sprintf("%s.pdf", fn), width=width, height=height)
    par(mgp=mgp, mar=mar)
    eval(parse(text=src))
    dev.off()
    header <- "\\begin{figure}[htbp]\\begin{center}"
    trailer <- sprintf("\\end{center}\\caption{%s}\\label{%s}\\end{figure}", caption, label)
    cat(sprintf("%s\\includegraphics[bb=0 0 %sin %sin, clip, width=%sbp]{%s.pdf}%s",
        header, width, height, size, fn, trailer))
}

利用は，以下のように。第一引数は作図プログラムを文字列で与える。

<<echo=false, results=tex>>=
xfigure("
hist(rnorm(1000))
abline(h=1:4*50, col=2, lty=2)
",
"testtest", "テストです", height=4)
@

結果は以下のようになる。いいね。

ダメ出し：if 文

if 文　というタイトルのついたページ

例を挙げているのだから，揚げ足取ってもしようがないのだけど，パズルとして

以下のようなものは，

> set.seed(123)
> x <- sample(9,12,replace=TRUE)
> x <- matrix(x,nrow=2,dimnames=list(1:2,letters[1:6]))
> for(i in 1:ncol(x)) if(diff(x[,i]) > 0) x[,i] <- NA
> x
a b c d e f
1 NA NA 9 NA 5 9
2 NA NA 1 NA 5 5

以下のように書ける

> x[,x[2,] > x[1,]] <- NA
> x
   a  b c  d e f
1 NA NA 9 NA 5 9
2 NA NA 1 NA 5 5

ダメ出し：散布図に記号を使うときは大きさの指定に注意すべし

変わった散布図 symbol 関数，そして引数 pch について少し　にて

> x <- cbind("a"=1:5,"b"=5:1,"c"=6:10,"d"=10:6,"e"=sample(10,5)/10)
> symbols(x[,1:2],circles=x[,1]/10,inches=FALSE,main="circles")

散布図において，半径（直径）が度数を表すように描いてはならない。
面積が度数を表すように描くべし。
つまり，上のだと，circles=sqrt(x[,1])/10 と指定すべし。

x <- cbind("a"=1:5, "b"=5:1, "c"=6:10, "d"=10:6, "e"=sample(10, 5)/10)
symbols(x[, 1:2], circles=x[, 1]/10, inches=FALSE, main="circles ×")
symbols(x[, 1:2], circles=sqrt(x[, 1])/10, inches=FALSE, main="circles ○")

ダメ出し：ダイナマイトプロット

棒グラフのいろいろ・その2 にて

ここに示されているようなグラフを，ダイナマイトプロットと呼ぶ。
描いてはいけないグラフの一つ。

なお，グラフを jpeg デバイスで描くのもどうかと思う。
せめて png で描こう。

また，プログラムで，カンマの後に空白がないと読みづらい。

ダメ出し：truehist ??

ヒストグラムと箱髭図　にて

> truehist関数は総面積が1に基準化されるので、縦軸が異なります。

というのはうそ。freq=FALSE にすれば同じこと。

一番違うのは，hist では right=TRUE がデフォルトになっていて，一般常識と異なる点（事実，リンク先のページに表示されているヒストグラムが若干違うのは，このせい）。

right=FALSE にすれば truehist と同じになる。

layout(matrix(1:2, 1))
hist(iris$Sepal.Length, main="hist", prob=TRUE, right=FALSE)
library(MASS)
truehist(iris$Sepal.Length, main="truehist")

ちなみに，「truehist関数は色がつきます」だけど，ショッキングピンクならぬ蛍光色がかった水色は，目に痛い。

結局，truehist ってどこが true なの?という感じ。

ダメ出し：行列の回転

image関数でイメージする(R Advent Calenderの番がきた)　にて

> (x <- matrix(1:12,ncol=3))
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

> (y <- t(x[nrow(x):1,ncol(x):1])[ncol(x):1,])
     [,1] [,2] [,3] [,4]
[1,]    4    3    2    1
[2,]    8    7    6    5
[3,]   12   11   10    9

というのだが，以下のようにすればよい。

(y <- t(x[nrow(x):1,]))

ついでに，左90度回転

> t(x[,ncol(x):1])
     [,1] [,2] [,3] [,4]
[1,]    9   10   11   12
[2,]    5    6    7    8
[3,]    1    2    3    4

さらについでに，180度回転

> x[nrow(x):1, ncol(x):1]
     [,1] [,2] [,3]
[1,]   12    8    4
[2,]   11    7    3
[3,]   10    6    2
[4,]    9    5    1

Sweave-13

Sweave-12 で書いた ae フォントを使うかどうかの件

やはり，困るという人もいるようで，マニュアルにちゃんと記載がある。

‘Sweave.sty’ also supports the [noae] option, which suppresses the use of the ae package, the use of which may interfere with certain encoding and typeface selections. If you have problems in the rendering of certain character sets, try this option.

マニュアル，ちゃんと読め ＞ オレ

マクネマー検定

「マクネマー検定についてのページ」で epibasix パッケージに mcNemar 関数があるのを知った。
しかし，この関数の実装は余り良くない。

(1) P 値が 0.001 より小さい場合，
McNemar's Chi^2 Statistic (corrected for continuity) = 12.345 which has a p-value of: 0
などと，表示され吃驚する。digit=3 になっているため 0.001 未満は 0 と表示されてしまうのだ。which has a p-value of less than 0.001 とか，せめて 0.000 とでも表示すればいかが？あるいは，デフォルトを digits=7 程度にしておくとか。

(2) "the number of discordant pairs exceeds 30" ということで，30 未満のときには force=TRUE などとやらねば，
以下にエラー mcNemar(data1, digits = 7) :
  The number of discordant pairs, must be greater than 30 to ensure validity.
などと，惨めな結果を思い知らされる。「なんにも計算しない」のではなく，警告メッセージと共に結果を表示すれば良いではないか。

さらにいえば，30 未満のときには二項検定を行ってその結果を表示すべき。さらにさらに，別に 30 未満に限らず，常に二項検定を行えばよいだけ。

> (x <- matrix(c(21, 3, 26, 72), 2))
     [,1] [,2]
[1,]   21   26
[2,]    3   72

> mcNemar(x)
 以下にエラー mcNemar(x) :
  The number of discordant pairs, must be greater than 30 to ensure validity.

> mcNemar(x, force=TRUE)

Matched Pairs Analysis: McNemar's Chi^2 Statistic and Odds Ratio
 
McNemar's Chi^2 Statistic (corrected for continuity) = 16.69 which has a p-value of: 0
 
McNemar's Odds Ratio (b/c): 8.667
95% Confidence Limits for the OR are: [3.319, -41.595]

> mcNemar(x, force=TRUE, digits=5)

Matched Pairs Analysis: McNemar's Chi^2 Statistic and Odds Ratio
 
McNemar's Chi^2 Statistic (corrected for continuity) = 16.68966 which has a p-value of: 4e-05
 
McNemar's Odds Ratio (b/c): 8.66667
95% Confidence Limits for the OR are: [3.31909, -41.59499]

> mcnemar.test(x)

    McNemar's Chi-squared test with continuity correction

data:  x
McNemar's chi-squared = 16.6897, df = 1, p-value = 4.402e-05

> binom.test(3, 29)

    Exact binomial test

data:  3 and 29
number of successes = 3, number of trials = 29, p-value = 1.524e-05
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.02186374 0.27351520
sample estimates:
probability of success
             0.1034483

Sweave-12

Sweave で作った pdf と，LaTeX だけで作った pdf を比べて，前者と後者のフォントの大きさ（当然ながら，版面の大きさとか行当たりの文字数など）が異なることに気づいた。
LaTeX のソースに \usepackage{Sweave} を加えるだけで，変になることがわかった。

消去法で突き止めた原因は，Sweave.sty にある
\newboolean{Sweave@ae}
\setboolean{Sweave@ae}{true}
が元凶のようだ。このようになっていると
\RequirePackage[T1]{fontenc} というfontを使うようになるようだ。そこで，
\setboolean{Sweave@ae}{false}
とすると，おかしなことはなくなった。

もっとも，そのほうがいいか，元のままがいいかは，個人の好みなんだろうなあ。

ダメ出し：apply関数は遅かった・・・

apply関数は遅かった・・・

確かに　apply は遅い。その代わりに colMeans, rowSums  を使って定義通りに書くと，ヘタに行列演算するより速い。
colMeans, rowSums を使って書いたものと，結果の格納領域を事前に確保して行列演算をするのは，ほとんど同じくらいの速度である。
それでもまあ，実行時間は1秒ミマンなのだから，こんな細かいことを言っても，毒にも薬にもならぬ。

> # apply
> Rprof()
> varA <- apply(d, 2, var)
> Rprof(NULL)
> print(a <- summaryRprof()$sampling.time)
[1] 3.56

> # colMeans, rowSums
> Rprof()
> means <- colMeans(d)
> varA2 <- rowSums((t(d)-means)^2)/(n-1)
> Rprof(NULL)
> print(a2 <- summaryRprof()$sampling.time)
[1] 0.1
> all.equal(varA, varA2)　　　# 結果が一致していることを確認
[1] TRUE

行列演算を使うと速いけど，cbind(  ) を使って結果をつなげてゆくなんてことをすると，遅くなってしまう。

> # matrix operation (2)
> Rprof()
> varC <- NULL
> for(i in 1:1000) {
+ varC <- cbind(　　　# cbind でつなぐのは悪手
+   varC,
+   (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)]^2 -
+   (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)])^2 / n)
+   / (n - 1)
+ )}
> varC <- as.numeric(varC)
> Rprof(NULL)
> print(c <- summaryRprof()$sampling.time)
[1] 0.96
>
> # matrix operation (2-2)
> Rprof()
> varC2 <- matrix(0, 100, 1000)　　# 前もって必要なだけ領域を確保しておく
> for(i in 1:1000) {
+   varC2[,i] <- (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)]^2 -
+     (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)])^2 / n)/ (n - 1)
+ }
> varC2 <- as.vector(varC2)
> Rprof(NULL)
> print(c2 <- summaryRprof()$sampling.time)
[1] 0.08
> all.equal(varC, varC2)　　　# 結果が一致していることを確認
[1] TRUE

ダメ出し：例題をやってみる その2

今回は，R プログラミングではなく，統計学の方でダメ出し
例題をやってみる　その2　で，
> ここでは、抗生物質を使うとプラセボより感染率は小さくなるはずだ、という期待が込められているので片側検定を行なっていると思われる。
のに
> （ここではカイ二乗検定を行なっている。）
とは，何を考えているのだろうか?そのようなときには，prop.test を使う。
そして，alternarive="less" を指定する。そうすると，
>ここでp=0.02と返ってきているが、論文では0.01と言っている。
などと，寝言を言わなくても，ちゃんと p= 0.01043 になる。

---

 

> prop.test(c(32, 53), c(488, 484), alternative="less")

    2-sample test for equality of proportions with continuity
    correction

data:  c(32, 53) out of c(488, 484)
X-squared = 5.3389, df = 1, p-value = 0.01043
alternative hypothesis: less
95 percent confidence interval:
 -1.00000000 -0.01212703
sample estimates:
    prop 1     prop 2
0.06557377 0.10950413

ダメ出し：例題をやってみる その1の続き ブートストラップ法

相変わらず，「空白などという無駄なものは置くものか」という，読む気のしないグッチャリプログラム

元のプログラムは 98 秒かかる（作者の環境では 3 分ぐらいかかるそうだ）が，

        pvalue <- sum(apply(boot, 1, mean) > 37)/l

を

        pvalue <- sum(rowMeans(boot) > 37)/l

に変えてやるだけで，6 秒ほどで計算が終わる。まあ，約 15 倍速くなると言うわけだ。
二重の for も，どうにかしたいと思うが，ヘタに手を入れると，かえって遅くなるのでそのまま。

以下はそのプログラム。

system.time({
set.seed(888)
tm <- c(37, 36, 37.1, 37.1, 36.2, 37.3, 36.8, 37, 36.3, 36.9, 36.7, 36.8) # 本の値
L <- (1:10)*1000 # リサンプリングする回数
trials <- 100 # 1リサンプリングあたりのシミュレーション回数
data <- matrix(0, trials, length(L))
for(l in L) {
    for(i in 1:trials) {
        boot <- matrix(sample(tm, size=l*length(tm), replace=TRUE), nc=length(tm))
        pvalue <- sum(rowMeans(boot) > 37)/l
        data[i, which(L == l)] <- pvalue
    }
}
matplot(t(data), type="p", pch=1, cex=0.2, col=1, axes=FALSE)
axis(1, 1:length(L), L)
axis(2)
lines(1:length(L), apply(data, 2, mean), col=4, lwd=3)
})
   ユーザ   システム       経過 
    6.249      0.163      6.400

ついでに，11/28　MIKUセミナーも
同じく apply(..., 1, mean) を rowMeans にするだけで倍速になる

> system.time({
+ set.seed(888)
+ L <- 100000 # リサンプリングする回数
+ tm <- c(37, 36, 37.1, 37.1, 36.2, 37.3, 36.8, 37, 36.3, 36.9, 36.7, 36.8) # 本の値
+ data <- matrix(sample(tm, size=L*length(tm), replace=TRUE), nc=length(tm)) # リサンプリングします
+ datamean <- rowMeans(data) # 平均値
+ sum(datamean > 37)/L # が、帰無仮説37度より高いか
+ hist(datamean, nclass=1000) # ヒストグラム
+
+ res <- cumsum(rowMeans(data) > 37)/(1:L)
+ plot(res, cex=0.1) # リサンプリングごとの、mu>37の確率。
+ })
   ユーザ   システム       経過 
     2.401      0.045      2.448