- 追加された行はこの色です。
- 削除された行はこの色です。
#contents
* Rの基本的な機能 ＆ tips [#ae2236af]
* Rの基本的な機能 ＆ tips [#ad5b75b8]
** R 本体（Windows版）のアップデート [#ba54b669]
** R 本体（Windows版）のアップデート [#c084b1fa]
RIGHT:22 October 2017
[[Updating R from R (on Windows) – using the {installr} package>https://www.r-statistics.com/2013/03/updating-r-from-r-on-windows-using-the-installr-package/]]~
[[Windows用Rのバージョンアップ>http://d.hatena.ne.jp/ryamada/20140522/1400554341]]
install.packages("installr")
library(installr)
updateR()
** Object の調査：str() [#k7529e80]
** Object の調査：str() [#ka949440]
RIGHT:21 March 2018
str() 関数はいろいろできるらしいのだけれど……
str(obj)
** 数の扱い [#f27734fa]
*** 二進数 digitsBase() [#ob4beb00]
** 数の扱い [#vc09be0b]
*** 二進数 digitsBase() [#f7aee0a3]
RIGHT:23November2013
10進数を二進数のベクトルに変換する~
パッケージ sfsmisc のインストールが必要
> library(sfsmisc)
> mydigits <- digitsBase(1101, base=2)
> mydigits
Class 'basedInt'(base = 2) [1:1]
[,1]
[1,] 1
[2,] 0
[3,] 0
[4,] 0
[5,] 1
[6,] 0
[7,] 0
[8,] 1
[9,] 1
[10,] 0
[11,] 1
> mydigits <- digitsBase(1101, base=2,13) ← ベクトルの長さを指定したいとき
二進数に限らず利用できる
*** 重複しないランダムな数 sample() [#zca2c6e9]
*** 重複しないランダムな数 sample() [#cf4f3bc4]
RIGHT:23November2013
> sample(0:20,3)
[1] 0 11 6
0から20のうちランダムに3つの数字を重複せずに取り出す
*** 数値の切り捨て [#cfd42ace]
*** 数値の切り捨て [#cd55ae45]
RIGHT:21 March 2018
万能ではないかもしれない~
さらに、これは、車輪の再発明になっているかもしれない
（が、round down が見つけられなかったのよ）
sprintf を使うと、round されてしまう~
統計量は round ではなく切り捨てがよいことも多いので
（四捨五入したら有意になってしまうかもしれない((そもそも、そんな微妙なデータでものを言うのは危ないので避けるべき。サンプルサイズを増やすか、データの取り方を考え直すかなどすべき))）、
文字列の操作とすることにして round down を実現する~
文字数を数えて、最終的な文字列の長さを決める~
元々の小数点以下桁数がわからないものを、小数点以下 dp桁に切り捨てる~
マイナス記号が付くと文字列が長くなるので……~
関数化しておくとよいのはわかるけれど、
後々のメンテナンスを考え（ブラックボックス化すると面倒なので）
関数化はしないでおく((コピー＆ペーストのミスの可能性は高まるので注意は必要))~
切り捨て桁（dp）がデータより小さい場合でも 0 fill はしません
mynum <- -0.052745
dp <- 3 # decimal places 小数点以下
# dp <- 9
mynumchar <- as.character(mynum)
if(mynum < 0) {
minusvalue <- 1
mynum <- round(-1 * mynum)
} else {
minusvalue <- 0
mynum <- round(mynum)
}
if( (mydigit <- log10(mynum)+1) < 0 ){ mydigit <- 1 }
mynum <- as.double(substr(mynumchar, 1, (minusvalue+mydigit+dp+1)))
# -0.052 になる
# dp <- 9 のときは元の -0.052745 になる
** 文字列 [#w228e223]
*** 文字列の一致を検査 [#qaf29533]
** 文字列 [#fee6b1df]
*** 文字列の一致を検査 [#x5aa2db1]
RIGHT:21 March 2018
文字列が一致のときは、match が integer の 1 を返す~
文字列が不一致のときは、match が integer の NA を返す~
NA は FALSE の代りにはならないので、is.na() を利用~
if(!is.na((match(myed1, myed2))){
# 一致したときの処理をココに
}
*** 文字列の表示：cat() [#qfeeaa1b]
*** 文字列の表示：cat() [#db490404]
RIGHT:21 March 2018
cat("\n
All Data: Courtship Latency > 4
\n")
** 出力 [#c625fa07]
*** ファイル出力：sink() [#x739b474]
*** 文字を置換 [#t5ab8af3]
RIGHT:23 March 2018
mystring は、list でも、dataframeの要素でもよい~
gsub は、同じ文字列に繰り返し適用~
正規表現が使えるらしい~
mystring <- sub("Re", "", mystring)
mystring <- gsub(".", "", mystring) # 全ての文字を消す
** 出力 [#ob797b9a]
*** ファイル出力：sink() [#qaf96617]
RIGHT:21 March 2018
sink("out.txt")
sink()
*** pdf出力：pdf() [#e3c6a480]
*** pdf出力：pdf() [#g680d436]
RIGHT:21 March 2018
サイズを黄金比 \(\left(1:\frac{1+\sqrt{5}}{2}\right)\) にしてみた
pdf("myout.pdf", width=(7*(1+sqrt(5))/2), height=7, pointsize=5)
...
dev.off()
* データの扱い [#a6aaa25e]
* データの扱い [#s74f5dff]
** データの入力 [#vc54b3d0]
** データの入力 [#m5533fed]
*** テキストファイルから取り込む read.table() [#u5b82831]
*** テキストファイルから取り込む read.table() [#ia783219]
RIGHT:21 March 2018
データフレームとして取り込まれる~
header=TRUE １行目がカラム名として取り込まれる~
fill=TRUE 空のデータがあるとき、空欄が取り込まれる（事故が起こらないか注意が必要）
mydf <- read.table("data.txt", header=TRUE)
mydf <- read.table("data.txt", header=TRUE, fill=TRUE)
*** 小さなデータなら直接書き込むほうが簡単 [#y9c58b4a]
*** 小さなデータなら直接書き込むほうが簡単 [#i19818e5]
RIGHT:21 March 2018
データの部分をスプレッドシートからコピー＆ペーストできる利点あり~
textConnection() を使って偽装することで実現できる
myinputmel <- ("
Female[tab]Copulated[tab]N
'Oregon-R'[tab]96[tab]100
'Canton-S'[tab]86[tab]100
'Hikone-R'[tab]92[tab]100
")
#
#names 属性を使ってデータに名前を付ける
names(myinputmel) <- "Crosses with D. melanogaster males "
#
# dataframe に読み込ませる
# textConnection() を使って偽装する
mydf <- read.table(textConnection(myinputmel))
*** 複数のデータを書いて、一辺に取り込ませることもできる [#fa99d589]
*** 複数のデータを書いて、一辺に取り込ませることもできる [#w182e744]
RIGHT:21 March 2018
myinput1 から myinput5 までは上のように作っておく~
for 文にすれば、関数を作らなくても複数処理できる~
for 文の中では表示されないものが出てくるので、show での対応が必要
myinputall <- c(
myinput1,
myinput2,
myinput3,
myinput4,
myinput5
)
#
mydataset <- length(myinputall)
#
for(i in 1:mydataset) {
mydf <- read.table(textConnection(myinputall[i]),
header=TRUE,
row.names=1)
#
# ココにそれぞれの処理
cat("\n------------------------------------------------------------\n",
names(myinputall[i]),
"\n\n")
#
}
** 重み付けデータの扱い [#dd50fd42]
** 重み付けデータの扱い [#q5f44d2f]
RIGHT:21 March 2018
重み付けデータはRでは扱いづらい（データの取り方がこうなることがある）~
survival analysis は別だけど~
そこで、個々のデータの形に集計し直す
プログラム引用：[[R numericとfactorの相互変換>https://siguniang.wordpress.com/2010/09/30/rnumericとfactorの相互変換/]]
weightdata.txt~
ある時間に対してイベントが起きた個体数を記録したデータ
Time Data1 Data2 Data3 ... Data11
1 6 4 5 7
2 1 1 0 1
...
各個体がそのイベントが起きたデータを持つという形に変換する
mydf <- read.table("weightdata.txt", header=TRUE)
#
# matrix に変換
mymat <- data.matrix(mydf)
#
# 頻度データを個々のデータに変換する
mydf2 <-
cbind(
#
# 数値データの数だけ、名前（交配）を並べる
data.frame(c(
rep(colnames(mymat)[2], sum(mymat[,c(2)])),
rep(colnames(mymat)[3], sum(mymat[,c(3)])),
rep(colnames(mymat)[4], sum(mymat[,c(4)])),
rep(colnames(mymat)[5], sum(mymat[,c(5)])),
rep(colnames(mymat)[6], sum(mymat[,c(6)])),
rep(colnames(mymat)[7], sum(mymat[,c(7)])),
rep(colnames(mymat)[8], sum(mymat[,c(8)])),
rep(colnames(mymat)[9], sum(mymat[,c(9)])),
rep(colnames(mymat)[10], sum(mymat[,c(10)])),
rep(colnames(mymat)[11], sum(mymat[,c(11)]))
)),
#
# 数値データを作成
# 頻度の回数だけ反復する
#
# 単純に dataframe にすると
# なぜか数値データもファクター扱いになってしまう
# そこで、数値に変換してから data.frame に変換する
# https://siguniang.wordpress.com/2010/09/30/rnumericとfactorの相互変換/
# as.numeric(as.character(size))
#
data.frame(as.numeric(
c(
rep(mymat[,1], mymat[,c(2)]),
rep(mymat[,1], mymat[,c(3)]),
rep(mymat[,1], mymat[,c(4)]),
rep(mymat[,1], mymat[,c(5)]),
rep(mymat[,1], mymat[,c(6)]),
rep(mymat[,1], mymat[,c(7)]),
rep(mymat[,1], mymat[,c(8)]),
rep(mymat[,1], mymat[,c(9)]),
rep(mymat[,1], mymat[,c(10)]),
rep(mymat[,1], mymat[,c(11)])
)
) )
)
> mydf2
［長い名前］ ［長い名前］
1 Data1 1
2 Data1 1
...
6 Data1 1
7 Data1 2
...
カラム名が変換操作をした履歴を示す「長い変な名前」になってしまう~
そこで、変換しておかないとニンゲンが理解できなくなる
> names(mydf2) <- c("cross","time")
> mydf2
Cross time
1 Data1 1
2 Data1 1
...
6 Data1 1
7 Data1 2
...
** データの操作 [#h6255b1e]
*** 行列の特定の数を置換 [#af542f12]
** データの操作 [#o646cf61]
*** 行列の特定の数を置換 [#h4596056]
RIGHT:23November2013
x[(x==0)]<-- -1
xという名前の行列の要素のうち、0のものを-1に置換
*** names 属性を使ってデータに名前を付ける：names() [#jaec7a1e]
*** names 属性を使ってデータに名前を付ける：names() [#l109b7c9]
RIGHT:21 March 2018
dataframe ではなくても可
names(mydf) <- "Crosses with D. melanogaster males "
*** カラム名変換：names() [#g67940f6]
*** カラム名変換：names() [#eac812e6]
RIGHT:21 March 2018
names(mydf2) <- c("cross","time")
** dataframe の操作 [#t517bfc0]
** dataframe の操作 [#qea0345d]
*** dataframe の全てのデータに代入 [#ubf54b73]
*** dataframe の全てのデータに代入 [#g6e668e9]
RIGHT:21 March 2018
mydf$Order <-0
*** dataframe に条件に一致するものだけに代入する（置換） [#gd235ea2]
*** dataframe に条件に一致するものだけに代入する（置換） [#uf0fc801]
RIGHT:21 March 2018
mydf[mydf<0] <- NA
mydf2[mydf2==0] <- NA
*** dataframe の factor を文字列として取り出す：as.character() [#i12c5d61]
RIGHT:21 March 2018
*** dataframe に条件に一致するものだけに代入する（置換） [#jccc2115]
RIGHT:23 March 2018
これは for文の中の一部分のコードで i や j を使っての操作
mychar1 <- as.character(mydf1$factor1[i])
mychar2 <- as.character(mydf2$factor2[j])
factor（MyFactor）が、Mylist のいずれかと一致する場合、Cat1 に 0 を代入~
for loop を使わずに１行でできてしまう~
すごい~
文字列が一致のときは、match が integer の 1 を返す~
文字列が不一致のときは、match が integer の NA を返す~
NA は FALSE の代りにはならないので、is.na() を利用~
mydf$Cat1[is.na(match(mydf$MyFactor, Mylist))==FALSE] <- 0
*** dataframe から条件に一致するものだけ取り出す：subset() [#hc841358]
*** dataframe から条件に一致するものだけ取り出す：subset() [#e87ff109]
RIGHT:21 March 2018
mydf2 <- subset(mydf, Order>0)
*** dataframe から条件に一致するものだけ取り出す：subset() [#h5fa1ee0]
*** dataframe から条件に一致するものだけ取り出す：subset() [#pf5b105c]
RIGHT:21 March 2018
subset= で行の抽出条件 ← Observed > 5 のデータの行を取り出すとき~
select= で列の抽出条件、取り出し順序の変更もできる~
as.matrix(): dataframe を matrix に変換
mydata <- as.matrix(subset(
mydf,
subset=Observed>5,
select=c(Courtship, Observed)
)
)
*** dataframe にデータを加える（横）：cbind() [#ieed20c7]
*** dataframe の factor を文字列として取り出す：as.character() [#pb680a41]
RIGHT:21 March 2018
これは for文の中の一部分のコードで i や j を使っての操作
mychar1 <- as.character(mydf1$factor1[i])
mychar2 <- as.character(mydf2$factor2[j])
*** dataframe から条件に一致する factor を文字列リストとして取り出す：as.character() [#m51f7b30]
RIGHT:23 March 2018
Cat1 が 0 である factor（Myfactorw）の取り出し
mycharlist <- as.character(mydf$Myfactor[mydf$Cat1==0])
*** dataframe にデータを加える（横）：cbind() [#q11ff28a]
RIGHT:21 March 2018
ここでは演算していますが、別な dataframe を右に付けることもできる~
サイズは一致していないとヘンなことになるので注意が必要
mydf <- cbind(mydf, mydf$N-mydf$Copulated)
*** dataframe にデータを加える（下）：rbind() [#q716dd96]
*** dataframe にデータを加える（下）：rbind() [#ofbc42c7]
RIGHT:21 March 2018
カラム名が一致しないと上手くいかない
mydf <- rbind( mydf, data.frame(columnname1="Dummy", columnname2=NA))
*** dataframe からカラムを取り出し、結合する：cbind() [#f07ccb81]
*** dataframe からカラムを取り出し、結合する：cbind() [#w66c7ff6]
RIGHT:21 March 2018
mydf3 <- cbind(mydf2[,3:4], mydf2[,6:10])
*** 空の因子レベルを取り除く：droplevels() [#r564c650]
*** 空の因子レベルを取り除く：droplevels() [#cdec30af]
RIGHT:21 March 2018
操作した dataframe には、
元の dataframe の因子の情報（因子レベル）が含まれている~
そこで、空の因子レベルを取り除かないと正しく検定できなくなってしまう
[[RWiki 因子Tips大全>http://www.okadajp.org/RWiki/?因子Tips大全]]
mydf2$EDnum <- droplevels(mydf2$EDnum)
*** カラム名を変更：names()[] [#gc6584ba]
*** カラム名を変更：names()[] [#s5b7cddc]
RIGHT:21 March 2018
何番目のカラムかを意識しないといけません
names(mydf)[3] <- "NotC"
*** dataframe を表示させる：show() [#s1ce247c]
*** dataframe を表示させる：show() [#v963d494]
RIGHT:21 March 2018
for 文の中では表示されないので、show() で対応する~
http://d.hatena.ne.jp/yusap/20130514/1368486988 ヘのコメントを参照
show(
mydf
)
*** dataframe の行数のカウント：nrow() [#j655fad5]
*** dataframe の行数のカウント：nrow() [#o148ee1c]
RIGHT:21 March 2018
x <- nrow(mydf)
***データの扱いについてのメモ [#cf8334a5]
***データの扱いについてのメモ [#ge8640fc]
RIGHT:11 July 2011
mytabledf <- data.frame(summary(myfit)$table)
とすると、myfit の出力テーブルがデータフレームになる。
str(mytabledf)で何が含まれているか分かる（records, n.max, n.start, events, median, X0.95LCL, X0.95UCL です）。
そこで、例えば、
mytabledf$median
で、（Kaplan-Meier estimate の）median だけが出てくる。
* 検定を一気に行うなど [#d5d7b443]
* 検定を一気に行うなど [#ye430ac7]
** データのまとめ [#d7fc63d0]
*** データのまとめ [#t8b03581]
** データのまとめ [#q7bafd08]
*** データのまとめ [#v98e858b]
RIGHT:09 May 2004
各$Ageの$Durationについて、
データのまとめ、平均、標準偏差の出力。~
> sink("out.txt")
> mydata <- read.table("data.txt", header=TRUE, fill=TRUE)
> mydata
> summary(mydata)
> mean(mydata, na.rm=TRUE)
> sd(mydata, na.rm=TRUE)
データの形式
Age2 Age4 Age6 Age8 Age10
9 6 7 4 8
.
.
.
- ''Tomaru, M.'', Oguma, Y. & Watada, M. 2006. Courtship in '''Drosophila quadrilineata''' with a unique male behavioral element, abdomen bending. Journal of Ethology. 24: 133--139. [[doi: 10.1007/s10164-005-0172-4:http://dx.doi.org/10.1007/s10164-005-0172-4]]
*** boxplot に使われる summary data の出力：グラフの描画ではない [#a26b62b4]
*** boxplot に使われる summary data の出力：グラフの描画ではない [#w24a0fa9]
RIGHT:21 March 2018
plot=FALSE がミソ
boxplot(mydf2$time~mydf2$cross, plot=FALSE)
** ノンパラメトリック検定 [#e6cae293]
** ノンパラメトリック検定 [#k2f8466e]
*** Kruskal-Wallis 検定 [#f76b3d51]
*** Kruskal-Wallis 検定 [#h1534618]
RIGHT:21 March 2018
事後検定として Dunn 検定
kruskal.test(Courtshipsec ~ EDnum, data = mydf2)
#
library(FSA)
mydt <- dunnTest(Courtshipsec ~ EDnum, data = mydf2, method="bh")
#
cat("\nDunn test\n")
# なぜか、この表示が sink()するとファイルに記録されなかったので書く
cat("Dunn (1964) Kruskal-Wallis multiple comparison\n")
cat(" p-values adjusted with the Benjamini-Hochberg method.\n")
#
mydf
#
library(rcompanion)
#
cat("\nDunn test: Compact letter display\n")
#
mylist <- cldList(comparison = mydt$Comparison,
p.value = mydt$P.adj,
threshold = 0.05)
#
# Produce summary statistics
Summarize(Courtship ~ EDnum,
data = mydf2,
digits=3)
*** Kruskal-Wallis test と multiple comaprison [#nd06fa1a]
*** Kruskal-Wallis test と multiple comaprison [#y1c1259d]
RIGHT:2017年10月11日
Kruskal-Wallis test（ノンパラメトリック分散分析）と
それに引き続く multiple comaprison（多重比較）を行う。
multiple comaprison は、
Zar (2010)((Zar, J. H. (2010) Biostatistical Analysis, 5th ed, Pearson Prentice-Hall, Upper Saddle River, NJ))
に
「ほぼ」((Zar には、P値の補正に Benjamini-Hochberg method を使うべきとは書いていない。実質的には使っているのだと思うけれども。))
従って、Dunn's method を利用。~
R では、
Benjamini-Hochberg method によって補正された P値も表示してくれる。~
Kruskal-Wallis test はデフォルトインストールだけで使える。~
multiple comaprison は
[[An R Companion for the Handbook of Biological Statistics, Kruskal-Wallis Test>https://rcompanion.org/rcompanion/d_06.html]]
に従って、FSA パッケージを利用。
データ（mydata.txt）の形式
numdata [tab] category
1.6059 [tab] Cat1
... [tab] ...
R code
# install.packages("FSA")
# https://rcompanion.org/rcompanion/d_06.html
#
x <- read.table ("mydata.txt", header=TRUE)
kruskal.test(numdata ~ category, data = x)
#
library(FSA)
dunnTest(numdata ~ category, data = x, method="bh")
分布を確認するためにボックスプロットしてみるなら
boxplot(numdata ~ category, data = x)
分布を確認するために蜂群＋ボックスプロットしてみるなら
library(beeswarm)
boxplot(x$numdata ~ x$category)
beeswarm(x$numdata ~ x$category, corral="wrap", add=T)
（beeswarm が、「x$」方式でないと動かないようなので、boxplotもそうしてみた）
*** ２標本検定 [#z5168c2f]
*** ２標本検定 [#a8d3dfbb]
RIGHT:2017年8月7日
- Wilcoxon-Mann-Whitney test または、Exact Wilcoxon rank sum test、U-test
- Brunner-Munzel test
- t検定も続けて、分散の比較もしておく
myx <- c(48, 34, 72, 86, 42)
myy <- c(35, 46 ,78, 29, 31)
#
# Exact Wilcoxon Mann Whitney test
library(exactRankTests) # wilcox.exact
wilcox.exact(myx,myy,exact=TRUE)
#
library(coin)
wilcox_test(
c(myx,myy) ~ factor(c(rep("myx",length(myx)), rep("myy",length(myy)))),
distribution="exact")
#
# Brunner-Munzel test
library(lawstat) # brunner.munzel.test
brunner.munzel.test(myx, myy)
#
# Brunner-Munzel test
library(nparcomp) # npar.t.test
writeLines("\n\nBrunner-Munzel test (npar.t.test)\n")
mydata <- data.frame(myxy=c(myx, myy),
cond=c(rep("myx",length(myx)),rep("myy",length(myy))))
myresults <- npar.t.test(myxy~cond, data = mydata,
method = "t.app", alternative = "two.sided")
summary(myresults)
#
# permuted Brunner-Munzel test
writeLines("\n\nBrunner-Munzel test (studentized permutation test)\n")
myresults <- npar.t.test(myxy~cond, data = mydata,
method = "permu", alternative = "two.sided")
summary(myresults)
#
# t test
t.test(myx, myy, var.equal=TRUE)
# t test with unequal variance, Welch t test
t.test(myx, myy, var.equal=FALSE)
# F test for variance
var.test(myx, myy)
ところが
> library(exactRankTests) # wilcox.exact
Package 'exactRankTests' is no longer under development.
Please consider using package 'coin' instead.
となる
library coin を使うのは
[[奥村さんの解説>https://oku.edu.mie-u.ac.jp/~okumura/stat/wmw.html]]
に従ってはいるが、出力が~
exactRankTests は、W値~
coin を使った場合は、Z値~
となる。
古くにU検定を教わった身には、W値のほうが、正規化されたZ値よりも気分がいい。
それに現時点では開発停止と書いてあるけれども、package `exactRankTests' は
2017-03-01版がある。~
参照：https://cran.r-project.org/web/packages/exactRankTests/index.html
U検定：
exactRankTests の wilcox.exact の結果と
library coin を使った結果とは一致するはず~
//
Brunner-Munzel test 検定：
library lawstat の brunner.munzel.test の結果と
library nparcomp の npar.t.test を使った結果とは一致するはず~
//
t検定：
t.test で、「var.equal=TRUE」を指定しなくてもデフォルトでこれ~
//
ウェルチの t検定：
t.test で、「var.equal=FALSE」を指定
*** Brunner-Munzel検定 [#md492c58]
[[奥村さんの解説>https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html]]
を見ていただくとして。
論文
- Brunner, E., & Munzel, U. (2000). The Nonparametric Behrens-Fisher Problem : Asymptotic Theory and a Small-Sample Approximation. Biometrical Journal, 42, 17–25. Doi: 10.1002/(SICI)1521-4036(200001)42:1<17::AID-BIMJ17>3.0.CO;2-U
- 名取真人. (2014). マン・ホイットニーのU検定と不等分散時における代表値の検定法. 霊長類研究, 30(1), 173–185. Doi: [[10.2354/psj.30.006>https://doi.org/10.2354/psj.30.006]]
*** Kruskal-Wallis test [#s1544387]
*** Kruskal-Wallis test [#q5c193fa]
RIGHT:09 May 2004
3つ以上の$Ageの$Durationについて検定。Kruskal-Wallis検定の後に、ペアワイズでWilcoxon Rank Sum Test（＝U検定）をする。おまけで最後に分散分析。
> sink("out.txt")
> mydata2 <- read.table("data2.txt", header=TRUE)
> mydata2
> summary(mydata2)
> kruskal.test(mydata2$Duration, mydata2$Age)
> pairwise.wilcox.test(mydata2$Duration, mydata2$Age)
> summary(aov(mydata2$Duration~mydata2$Age))
データの形式
Age Duration
Age2 6
.
.
.
- ''Tomaru, M.'', Oguma, Y. & Watada, M. 2006. Courtship in '''Drosophila quadrilineata''' with a unique male behavioral element, abdomen bending. Journal of Ethology. 24: 133--139. [[doi: 10.1007/s10164-005-0172-4:http://dx.doi.org/10.1007/s10164-005-0172-4]]
カテゴリ毎にサマリを出力（追記 22 February 2015）
> by(mydata2, mydata2$Age, summary)
** パラメトリック検定 [#bef85c88]
** パラメトリック検定 [#mc4be4dc]
*** n元配置：分散分析、Tukey HSD 検定 [#ua24f771]
*** n元配置：分散分析、Tukey HSD 検定 [#i0d3c98e]
RIGHT:21 March 2018
分散分析を行ってから、Tukey HSD 検定~
最近は、Tukey HSD 検定だけでもよいみたいですね~
プログラム引用： https://rcompanion.org/rcompanion/d_05.html
library(FSA)
#
cat("\n\nSummary statistics\n")
#
# Produce summary statistics
Summarize(Courtshipsec ~ EDnum,
data = mydf2,
digits=3)
#
cat("\n\nANOVA\n")
#Fit the linear model and conduct ANOVA
model = lm(Courtshipsec ~ EDnum,
data = mydf2)
#
anova(model)
summary(model)
#
cat("\n\nTukey HSD\n")
#
#Tukey comparisons in agricolae package
library(agricolae)
(HSD.test(model, "EDnum")) # outer parentheses print result
** Contingency Table [#t9c03134]
*** Fisher test [#cd55fede]
*** Fisher test [#c6c3c946]
RIGHT:03 Oct 2005
> mydata <- matrix(c(4,6,3,4), nr=2)
> fisher.test(mydata)
*** Fisher test for 2 by n [#e2dfa30b]
*** Fisher test for 2 by n [#d9020889]
RIGHT:17 September 2010
> mydata <- matrix(c(15,19,10,15,3,8,3,3,9,2,2,0,1,3,2,4,9,7,5,7,4,10,4,4,2,2,4,2,2,1), nr=15)
> mydata
> fisher.test(mydata, workspace=20000000)
workspaceを大きく取らないとデータによってはエラーになる場合がある
*** 分割表の独立性の検定 [#eb3bb1bd]
*** 分割表の独立性の検定 [#le0202f9]
RIGHT:21 March 2018
dataframe は 2x2 ではなくて、nx2 のデータが格納されているとき~
まず全体（nx2）の分割表の独立性検定を行い、
続いて、pariwise comparisons を行う~
Fisher Exact Test、カイ二乗検定、G検定で
カイ二乗検定では期待値を表示させて、
検定を行ってよい条件かどうかを見てわかるようにする~
さらに、独立性の検定の期待値が適切かどうかを表示させる((Zar, J. H. (2010) Biostatistical Analysis, 5th ed, Pearson Prentice-Hall, Upper Saddle River, NJ))~
- 期待値は１以上であるべき
- 期待値が５未満のセルの数は20%以下
Fisher Exact Test：サンプルサイズが大きいときは
workspace を大きく取っても計算ができなくなるので注意が必要
プログラム引用：
An R Companion for the Handbook of Biological Statistics~
Salvatore S. Mangiafico~
http://rcompanion.org/rcompanion/b_01.html~
http://rcompanion.org/rcompanion/b_05.html~
http://rcompanion.org/rcompanion/b_06.html
fisher.test(mydata, workspace=20000000)
#
cat("\nPost hoc pairwise comparison (Fisher exact test) adjusted by Holm's method\n")
#
FUN = function(i,j){
fisher.test(matrix(c(mydata[i,1], mydata[i,2],
mydata[j,1], mydata[j,2]),
nrow=2,),
workspace=20000000)$ p.value
}
#
pairwise.table(FUN,
rownames(mydata),
p.adjust.method="holm")
#
#
# カイ二乗検定
mychitest <- chisq.test(mydata)
#
# 期待値を明示する
expecteddf <- data.frame(mychitest$expected)
#
# 期待値が 0 になるようなデータをカウントする
lessthan1 <- expecteddf[expecteddf <1]
if(length(lessthan1)<1){
lessthan1 <- 0
}
# 期待値が 5 より小さなものをカウントする
lessthan5 <- expecteddf[expecteddf <5]
if(length(lessthan5)<1){
lessthan5inpercent <- 0
} else {
lessthan5inpercent <- 100*length(expecteddf[expecteddf <5])/length(expecteddf[expecteddf > 0])
}
#
# 独立性の検定の期待値が適切かどうかを表示させる
cat(
"Check expected values \n",
"If exp < 1.0, omit such data\n",
"Number of exp less than 1.0: ", lessthan1, "\n",
"If exp < 5.0 are more than 20%, omit some data\n",
"exp less than 5.0: ", lessthan5,
" -> ", lessthan5inpercent, "%\n"
)
#
mychitest
#
# G検定もする
GTest(mydata,correct="williams")
#
# 2x2以上のときの事後検定
library(DescTools)
cat("\n\nPost hoc pairwise comparison (chi-square test wite Yates correction) adjusted by Holm's method\n")
#
FUN = function(i,j){
chisq.test(matrix(c(mydata[i,1], mydata[i,2],
mydata[j,1], mydata[j,2]),
nrow=2,
byrow=TRUE),
correct=TRUE)$ p.value
}
#
pairwise.table(FUN,
rownames(mydata),
p.adjust.method="holm")
#
#
cat("\n\nPost hoc pairwise comparison (G test with William's correction) adjusted by Holm's method\n")
#
FUN = function(i,j){
GTest(matrix(c(mydata[i,1], mydata[i,2],
mydata[j,1], mydata[j,2]),
nrow=2,
byrow=TRUE),
correct="williams")$ p.value # "none" "williams" "yates"
}
#
pairwise.table(FUN,
rownames(mydata),
p.adjust.method="holm")
** Goodness of Fit [#fa70f5d3]
** Goodness of Fit [#x3539842]
*** カイ二乗検定 [#h2b72f37]
*** カイ二乗検定 [#a50cb45e]
RIGHT:14 January 2016
適合性の検定
メンデルの法則などのようなものの検定
Zar (1981) Example 5.1 のデータ~
null hypothesis (H&subsc{0};) は 3:1 （＝メンデル分離するかどうかを検定）
参考にしたページ~
An R Companion for the Handbook of Biological Statistics
by Salvatore S. Mangiafico~
[[G-test of Goodness-of-Fit>http://rcompanion.org/rcompanion/b_04.html]]
> obs <- c(84,16)
> hyp <- c(3/4,1/4)
> chisq.test(obs, p=hyp)
手計算と同じことをする場合
> exp = sum(obs)*hyp
> dif <- obs - exp
> dif2 <- dif*dif
> dif2div <- dif2 / exp
> chi <- sum(dif2div)
> pchisq(chi,df=1,lower.tail=FALSE)
しかし ……、自由度が 1 のときは Yates correction を用いるのが望ましい。
ところが、chisq.test は contingency table では、correction=T
オプションを付ければ Yates correction を行うが、
適合性の検定では行わない。~
そこで、手計算と同じようにせざるを得ない。~
差分の絶対値から0.5を引きそれを二乗するように式を変えればよい（dif2corrected）。
> obs <- c(84,16)
> hyp <- c(3/4,1/4)
> exp = sum(obs)*hyp
> dif <- obs - exp
> dif2corrected <- (abs(dif)-0.5)*(abs(dif)-0.5)
> dif2correcteddiv <- dif2corrected / exp
> chic <- sum(dif2correcteddiv)
> pchisq(chic,df=1,lower.tail=FALSE)
*** G検定を行いたいなら [#c8691e2a]
*** G検定を行いたいなら [#m9aeaaa7]
RIGHT:14 January 2016
> G = 2 * sum(obs * log(obs / exp))
> a <- 2
> n <- sum(obs)
> df <- a-1
> q <- 1 +((a*a-1)/(6*n*df))
> Gadj <- G/q
ここでは、Gadj は Williams' correction~
a はカテゴリ数~
log が自然対数、常用対数は log10~
** 頻度データの95%信頼限界の計算：binom.test() [#cd8b820a]
** 頻度データの95%信頼限界の計算：binom.test() [#qb1cf680]
RIGHT:21 March 2018
二項分布を利用~
binom.test() の結果のうち~
$statistic: データ~
$parameter: データのN~
$conf.int[1]: 95%信頼限界下限~
$conf.int[2]: 95%信頼限界上限~
if(myobs > 0){
mybntest <- binom.test(mycourt, myobs)
cat(
100*mybntest$statistic/mybntest$parameter, # rate in percent
"\t",
100*mybntest$conf.int[1],
"\t",
100*mybntest$conf.int[2],
"\t"
)
** survival analysis [#d5696480]
** survival analysis [#g7b217dd]
RIGHT:11 April 2004
ショウジョウバエの配偶行動を見るときに、観察時間内に交尾しないものがある。そこで、交尾したものを「dead」、交尾しなかったものを「survive」とみなしてsurvival analysisを行う。
Kaplan-Meier estimatorを計算
> library(survival)
> mydata <- read.table("data.txt", header=TRUE)
> mydata
> sink("out.txt")
> myfit <- survfit( Surv(Time, Copl)~Age, data=mydata)
> myfit
> summary(myfit)
プロット
> plot(myfit)
Log Rank test
> survdiff( Surv(Time, Copl)~Age, data=mydata)
データの形式~
Copl=1は交尾。Copl=0は交尾せず（ここでは観察時間が120）。
Age Time Copl
2 41 1
.
.
.
2 120 0
.
.
.
- ''Tomaru, M.'', Oguma, Y. & Watada, M. 2006. Courtship in '''Drosophila quadrilineata''' with a unique male behavioral element, abdomen bending. Journal of Ethology. 24: 133--139. [[doi: 10.1007/s10164-005-0172-4:http://dx.doi.org/10.1007/s10164-005-0172-4]]
** 相関の計算・検定：cor.test() [#z79f9b5f]
** 相関の計算・検定：cor.test() [#w8a17c9e]
RIGHT:21 March 2018
corresults <- cor.test(
mydf3$CourtshipFreq,
mydf3$Courtshipsec
)
** 相関を計算：cor() [#vb5aa407]
** 相関を計算：cor() [#da574f03]
RIGHT:21 March 2018
dataframe に３つ以上のデータが含まれている場合
# mycor <- cor(mydf3) # use指定をしないと NA データがあると NA となる
# mycor <- cor(mydf3, use="complete.obs") # これは全体の covariance を使う
mycor <- cor(mydf3, use="pairwise.complete.obs") # これは pairwise
mycod <- mycor^2 # coeffient of determination
** 相関を個々に計算 [#o193ba29]
** 相関を個々に計算 [#v8cc8aed]
RIGHT:21 March 2018
dataframe に３つ以上のデータが含まれている場合
for(i in 1:7) {
for(j in 1:7) {
mycortest <- cor.test(mydf3[, i], mydf3[, j], method="pearson")
cat("data i: ", names(mydf3)[i], "\ndata j: ", names(mydf3)[j])
print(mycortest)
}
}
#
for(i in 1:7) {
for(j in 1:7) {
mycortest <- cor.test(mydf3[, i], mydf3[, j], method="kendall")
cat("data i: ", names(mydf3)[i], "\ndata j: ", names(mydf3)[j])
print(mycortest)
}
}
#
for(i in 1:7) {
for(j in 1:7) {
mycortest <- cor.test(mydf3[, i], mydf3[, j], method="spearman")
cat("data i: ", names(mydf3)[i], "\ndata j: ", names(mydf3)[j])
print(mycortest)
}
}
** Comparing Regressions [#md8494d1]
*** 異なる条件下で得られた複数の回帰直線を比較する [#q77492a7]
*** 異なる条件下で得られた複数の回帰直線を比較する [#y898adae]
RIGHT:01 February 2008
-傾きに差はある？
-切片に差はある？
長くなったのでこちら -> [[複数の回帰直線の比較>R_ComparingRegressions]]
* グラフ [#wb89e66d]
* グラフ [#m7d7fc47]
** バーグラフに、95%信頼限界の線分を加える [#le4431ee]
** バーグラフに、95%信頼限界の線分を加える [#c3d9f906]
RIGHT:21 March 2018
mybar <- barplot(mydf$CourtshipFreq, names=mydf$EDonly,
main="Courtship", cex.main=2.5, cex.lab=1.6, las=2,
main="Courtship", cex.main=2.5, cex.lab=1.6, cex.names=0.8,
las=2,
ylab="Frequency (%)", ylim=c(0,100), axis.lty=1, xaxt="n"
)
#
#グラフに囲みを付ける 21 March 2018
#グラフに囲みを付ける
box(lty=1)
# 信頼限界（上）を描く
arrows(mybar, mydf$CourtshipFreq, mybar, mydf$CourtCIH,angle=90,length=0,lwd=1)
# 信頼限界（下）を描く
arrows(mybar, mydf$CourtshipFreq, mybar, mydf$CourtCIL,angle=90,length=0,lwd=1)
#length = 0.1 などとするとヒゲも付けられる
** 表示順指定：factor() [#q73e5842]
** ボックスプロット [#o9ce125d]
boxplot(mydf$data ~ mydf$cat, las=2, outline=F,
main="Copulation Duration", cex.main=2.5, cex.lab=1.6,
ylab="Copulation Duration (s)", xaxt="n")
# X軸ラベルは boxplot 関数で書かず別途書き込み
axis(side=1, at=c(1:mynumxaxis), label=mydf2$name, las=2, cex.axis=0.8)
# label=mydf2$name 別な dataframe の文字列を X軸ラベルに利用
# label=F でラベル文字列を書かない
** 表示順指定：factor() [#l114b2a1]
RIGHT:21 March 2018
factor属性に levels引数を渡す~
c() はリストなので、順序の情報が含まれている
→ リストの順序が利用される仕組みになっている
mydf2$cross <- factor(mydf2$cross, levels=c( "Data7", "Data2", "Data5", "Data3", ..., "Data1"))
** ダミーデータ追加 [#lbca4082]
** ダミーデータ追加 [#qc488c35]
RIGHT:21 March 2018
図を描くとき「空き」を作るためなどに~
場合によっていは、名前の文字列は空白がいいのかもしれない
（でもひとつしか空きは作れないはず）
mydf2 <- rbind( mydf2, data.frame(cross="Dummy", time=NA))
** グラフ領域の操作：par(oma= )) [#q44275f2]
** グラフ領域の操作：par(oma= )) [#z40b6ab4]
RIGHT:21 March 2018
par(oma=c(9,3,0,0)) # 余白を作る 下左上右
** 太字、斜体の対応：substitute(paste()) italic() bold() bolditalic() [#bb310690]
** 太字、斜体の対応：substitute(paste()) italic() bold() bolditalic() [#fb768fcb]
RIGHT:21 March 2018
substitute と paste のコンビネーション
なので、適用できない場面も多々ある~
グラフタイトルはおっけい
title(main=
substitute(paste(
bold('A. '), bolditalic('D. melanogaster'), bold(' males'))
),
cex.main=1.5, cex=1.5, )
** グラフタイトルに変数を入れる：paste() [#c3c6f9e6]
** グラフタイトルに変数を入れる：paste() [#o1125407]
RIGHT:21 March 2018
paste を使って文字列を結合すれば、変数がきちんと展開される~
italic などが使いたい場合は substitute() paste() list() を利用する
title(main=paste("r", " = ", coefficient, ", ", "P", " = ", pvalue))
** グラフタイトルに変数を入れる：太字、斜体の対応も必要な場合 substitute() paste() list() [#v2440a97]
** グラフタイトルに変数を入れる：太字、斜体の対応も必要な場合 substitute() paste() list() [#i622d18e]
RIGHT:21 March 2018
substitute と paste のコンビネーション~
関数 substitute(expr, env) は exprに変数が入っているだけでは展開しない~
そこで、env に適当な関数を入れて展開させる（計算もできる）~
今回は単に数字に置き換えるだけなので、list() で対応した
title(
main=substitute(
paste(italic("r"), " = ", mycoefficient,
", ", italic("P"), " = ", mypvalue),
list(mycoefficient=coefficient, mypvalue=pvalue)
),
cex=0.7
)
** legend の操作 [#f2e6906d]
** legend の操作 [#gf5b0cc8]
RIGHT:21 March 2018
leged に title を付けることができる
legend("topleft", lty=1:5,
title="Female",
bty ="n", # 枠線なし
# cex = 0.8, # 文字サイズ 0.8倍
legend=c(
"DATA1",
"DATA2",
"DATA3",
"DATA4",
"DATA5"
)
)
** 複数の折れ線グラフ：matplot() [#m52f6509]
** 複数の折れ線グラフ：matplot() [#xfecd2e5]
RIGHT:21 March 2018
matplot(
cbind(
mydf$DATA1,
mydf$DATA2,
mydf$DATA3,
mydf$DATA4,
mydf$DATA5
),
cex.main=1.5, cex.lab=1.5, cex.axis=1.2,
type="l",lwd=1, lty=1:5, ylim=range(0,100), col="black", las=1,
bty ="n", # 枠線なし
# bty ="l", # 枠線が左と下
ylab="Cumulated number of copulated pairs",
xlab="Time (s)")
** boxplot に蜂群を重ねる：boxplot() beeswarm() [#s0fec42b]
** boxplot に蜂群を重ねる：boxplot() beeswarm() [#r21365a6]
RIGHT:21 March 2018
par(oma=c(9,3,0,0)) # 余白を作る 下左上右
boxplot(mydf2$time~mydf2$cross, outline=FALSE, las=2, cex.main=1.5,
cex.lab=1.5, cex.axis=1.2, main="C. Copulation latency",
ylab="Copulation latency (min)", ylim=c(0,60),
names=c("DATA1", "DATA2", "DATA2", "DATA3", "", "DATA4",
"DATA5", "DATA6", "DATA7", "DATA8"),
bty ="n
)
# outline=F 外れ値を表示させない ← beeswarm と重ねるときは外れ値はまずい
# las=2 軸の表示文字の向き X下から上へ、Yは横
# names = X軸の表示文字列を指定
# bty ="n" 枠線なし
library(beeswarm)
beeswarm(mydf2$time~mydf2$cross, method="center", pch = 16,
corral="gutter", add = TRUE)
** Bar graph：barplot() [#i8ff1573]
** Bar graph：barplot() [#n165700e]
RIGHT:21 March 2018
barplot(copulationdata, cex.main=1.5, cex.lab=1.5, cex.axis=1.2,
main="D. Copulation frequency (1 h)", ylim=c(0,100),
ylab="Copulation (%)", las=2
)
** 散布図を描く：plot() [#d2a2983d]
** 散布図を描く：plot() [#j468a1a8]
RIGHT:21 March 2018
まずは plot というのが R のグラフ作成のキモらしい~
複数のデータセットになっていると、ペアワイズの複数のグラフがプロットされる
plot(mydf)
** グラフ描画領域内に文字を書く：text() [#t7e750f4]
** グラフ描画領域内に文字を書く：text() [#kb070bf3]
RIGHT:21 March 2018
adj=c(0.5,0) # テキストをxy座標中央揃え、下揃え~
x 座標については、0 左揃え、0.5 中央揃え、1 右揃え らしい~
y 座標については、0 下揃え、0.5 中央揃え、1 上揃え らしい~
y 座標はデフォルトで上揃えなので、高さの異なるアルファベットが並ぶと、
ガタガタにならないよう下揃えにする
text(x=1, y=61, labels="a", cex=1.1)
text(x=2, y=61, labels="a", cex=1.1)
text(x=3, y=61, labels="a", cex=1.1)
text(x=4, y=61, labels="a", cex=1.1)
#
text(x=6, y=59.6, labels="ab", cex=1.1, adj=c(0.5,0))
text(x=7, y=59.6, labels="ab", cex=1.1, adj=c(0.5,0))
text(x=8, y=59.6, labels="a", cex=1.1, adj=c(0.5,0))
text(x=9, y=59.6, labels="b", cex=1.1, adj=c(0.5,0))
text(x=10, y=59.6, labels="b", cex=1.1, adj=c(0.5,0))
** 軸ラベルの下に文字を書く：mtext() [#y0fb286d]
** 軸ラベルの下に文字を書く：mtext() [#j08de178]
RIGHT:21 March 2018
描画領域のサイズが変わると位置が変わるので注意~
substitute(paste()) で、斜体、太字もできる
mtext(text="Female", at=-1, padj=37, cex=1.2)
mtext(text="Male", at=-1, padj=43, cex=1.2)
mtext(text=substitute(paste(italic('D. melanogaster'))), at=3, padj=34, cex=1.2)
mtext(text=substitute(paste(italic('D. sechellia'))), at=8, padj=42, cex=1.2)
** グラフ内に線分を描く：lines() [#gd3c7012]
** グラフ内に線分を描く：lines() [#o236bd94]
RIGHT:21 March 2018
lines(X座標のベクトル, Y座標のベクトル)~
lines(点１,点２,点３)の形式ではないので注意
lines(c(3.5,5.4),c(95,95))
** 『「蜂群＋平均値」に折れ線グラフを重ねる』をふたつ重ねる [#w6d1e8c6]
** 『「蜂群＋平均値」に折れ線グラフを重ねる』をふたつ重ねる [#p19c0bb8]
RIGHT:2017年8月7日
- 横軸はカテゴリ
- 一番左のカテゴリは少し性質が違うので、左からひとつめと左からふたつめをつなぐ折れ線は他と別の破線で
- 『「蜂群＋平均値」に折れ線グラフを重ねる』をずらしてふたつ重ねる
- 平均値はグラフ作成の途中で計算させる
- 「ずらし」を沢山使うので、カテゴリを factor と score を使って振り分け
&ref(beeswarmmeanplot.png);
# DATA HERE
# cat0
control0 <- c(3986, 4246, 5102, 4920, 6102)
test0 <- c(3055, 2203, 2948, 3108, 5261)
#
# cat1
control1 <- c(1848, 434, 2072, 1786, 2642)
test1 <- c(1220, 1185, 1702, 1212, 1515)
#
# cat2
control2 <- c(2981, 2217, 2880, 2962, 3764)
test2 <- c(1138, 672, 3669, 472, 512)
#
# cat3
control3 <- c(1453, 2846, 2652, 1447, 2261)
test3 <- c(1520, 2082, 1006, 1566, 1224)
#
# End of DATA HERE
###################################################
#
# Draw graphs
#
library(beeswarm)
#
## 軸の余裕分 myxlimadd
myxlimadd <- 1
#
###################################################
#
postscript("out.eps", horizontal=FALSE, height=6, width=5, pointsize=15)
#
###################################################
#
mytitle <- "My Title"
#
# 縦軸の範囲を決めておく：繰り返し使えるように
yminmax <- c(0,7000)
#
# 縦軸の目盛刻み（ラベルにもなる）
ytick <- c(0,1000,2000,3000,4000,5000,6000,7000)
#
# 横軸のカテゴリ ＋ beewarm用隙間（5）
myxlim <- 11 + myxlimadd
#
## 横軸の目盛刻み
xtick <- c(1.5,4.5,7.5,10.5)
#
# 左：control
#
# データを描く横軸の目盛上の位置：左
myx <- c(1,4,7,10)
# data
data0 <- control0
data1 <- control1
data2 <- control2
data3 <- control3
#
# 「1-3」データのみ：「0」データは NA とする
meanright <- c(NA,mean(data1), mean(data2), mean(data3))
#
# 「1-3」データ左端と「0」データ。ほかのデータは NA とする
meanleft <- c(mean(data0),mean(data1),NA,NA)
#
#
# 「1-3」データ折れ線
#
# cex=3, col="gray" サイズ10倍、色グレー
plot(meanright~myx,type="b",pch="-", cex=3, col="gray", ylim=yminmax,
xlim=c(1,myxlim),xlab="",ylab="",axes=FALSE, main=mytitle)
#
# 「1-3」データ左端と「0」データ折れ線
#
par(new=T)
#
# cex=3, col="gray" サイズ10倍、色グレー
plot(meanleft~myx,type="b",pch="-", cex=3, col="gray", ylim=yminmax,
xlim=c(1,myxlim),xlab="",ylab="", lty=2, axes=FALSE)
#
# beeswarm で同じ点を打つために繰り返し数を入れておく
condition <- factor(c(
rep(1, length(data0)), # 「0」
# dummy
2,3,
# dummy
rep(4, length(data1)), # 「1」
# dummy
5,6,
# dummy
rep(7, length(data2)), # 「2」
# dummy
8,9,
# dummy
rep(10, length(data3)), # 「3」
11 # dummy
))
#
score <- c(
data0,
NA,NA,
data1,
NA,NA,
data2,
NA,NA,
data3,
NA
)
#
# 蜂群プロット
# pch=5：ダイヤモンド◇
beeswarm( score~condition, pch=5, xlim=c(1,myxlim),ylim=yminmax,
axes=FALSE, add=TRUE, correl="wrap")
#
###################################################
#右にずらして重ねる
#
par(new=T)
#
# data
data0 <- test0
data1 <- test1
data2 <- test2
data3 <- test3
#
# データを描く横軸の目盛上の位置：右
myx <- c(2,5,8,11)
#
# 「1-3」データのみ：「0」データは NA とする
meanright <- c(NA,mean(data1), mean(data2), mean(data3))
#
# 「1-3」データ左端と「0」データ。ほかのデータは NA とする
meanleft <- c(mean(data0),mean(data1),NA,NA)
#
# 「1-3」データ折れ線
# cex=3, col="gray" サイズ10倍、色グレー
plot(meanright~myx,type="b",pch="-", cex=3, col="gray",
ylim=yminmax,xlim=c(1,myxlim),xlab="",ylab="",axes=F)
#
# 「1-3」データ左端と「0」データ折れ線
#
par(new=T)
#
# cex=3, col="gray" サイズ10倍、色グレー
plot(meanleft~myx,type="b",pch="-", cex=3, col="gray",
ylim=yminmax,xlim=c(1,myxlim),xlab="",ylab="", lty=2, axes=F)
#
# beeswarm で同じ点を打つために繰り返し数を入れておく
condition <- factor(c(
1, # dummy
rep(2, length(data0)), # 「0」
# dummy
3,4,
# dummy
rep(5, length(data1)), # 「1」
# dummy
6,7,
# dummy
rep(8, length(data2)), # 「2」
# dummy
9,10,
# dummy
rep(11, length(data3)) # 「3」
))
#
score <- c(
NA,
data0,
NA,NA,
data1,
NA,NA,
data2,
NA,NA,
data3
)
#
#pch=1：白抜き丸○
beeswarm( score~condition, pch=1, xlim=c(1,myxlim),ylim=yminmax,
axes=FALSE, add=TRUE, correl="wrap")
#
# 縦軸：side=2
# las = 1：X軸、Y軸ともに数値ラベルを水平に書く
# line=2：縦軸を左にずらす
axis(2,ytick,las=1,line=1)
#
# 横軸：side=1
axis(1,xtick,labels=FALSE)
###################################################
#
# postscript (eps) 出力終了
dev.off()
*ESS [#s80ba893]
Emacs Speaks Statistics：Emacs で R を使う~
Emacs には後からインストールをする必要あり~
[[install-memo5.1]]を参照
** ESS起動 [#f3a197ca]
** ESS起動 [#mfcab239]
RIGHT:03 Oct 2005
M-x R
コマンド入力後
C-c C-j
C-m （改行）
で、１行ずつコマンドを実行する。
*コマンドメモ [#n6486b66]
*コマンドメモ [#w67b1ba7]
RIGHT:30 March 2004
> q() ← quit
> source("commands.R") ← コマンドを書き込んだソースを読み実行
> sink("out.txt") ← ファイルに出力
> sink() ← 出力をコンソールに戻す
> objects() ← 現在あるオブジェクトを表示
> mydata <- read.table("data.txt", header=TRUE) ← ファイルからデータを読み込む
> library(ctest) ← ライブラリを使う
> t.test(A, B) ← t検定（unpaired）。※ A、Bにはデータを入れておく
> var.test(A, B) ← F検定
> t.test(A, B, var.equal=TRUE) ← t検定、等分散
> wilcox.test(A, B) ← U検定（＝Wilcoxon test）
> ks.test(A, B) ← Kolmogorov-Smirnov test
> getwd() ← working directoryの表示
> setwd("c:/cygwin/home/tomaru/stats") ← working directoryの変更
追記 12 July 2013
> length(X) ← データの数
> mean(X) ← 平均
> sd(X)/(length(X)^0.5) ← 標準誤差
追記 22 February 2015~
> help.start() ← htmlヘルプを起動
**コマンドラインで利用する（Windows） [#ycaab92d]
**コマンドラインで利用する（Windows） [#u550b7ec]
RIGHT:05September2013
cygwin shell または、コマンドプロンプトで
rscript command.R
rscript command.R & ← バックグラウンドで実行する場合
※ pathを通しておかなければなりません（MS-DOS！のようですが）。
参考：[[R の入出力画面を経由せずに実行：他のプログラムからの呼び出し:http://takenaka-akio.org/doc/r_auto/chapter_14.html]]
*Rのリンク [#fd0daed0]
- [[RjpWiki>http://www.okada.jp.org/RWiki/]]
|Today:&counter(today); |Yesterday:&counter(yesterday); |Total:&counter(); since 04 April 2006|
![[PukiWiki]](image/tomarufront.png "[PukiWiki]")
