![[PukiWiki]](image/tomarufront.png "[PukiWiki]")
#author("2019-07-21T10:25:02+09:00","","") #author("2019-07-21T10:29:22+09:00","","") このページは分割しました 21 July 2019 - [[Rの基本的な機能 ＆ tips>R_base]] - [[Rのおけるデータの扱い>R_data]] - [[Rで検定を一気に行うなど>R_statistics]] - [[Rでグラフを作成する>R_graphs]] ----------------------------- 以下は分割前のページ ----------------------------- #contents * Rの基本的な機能 ＆ tips [#je115281] ** R 本体（Windows版）のアップデート [#jb05f73b] 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() [#l844a153] RIGHT:21 March 2018 str() 関数はいろいろできるらしいのだけれど…… str(obj) ** 数の扱い [#qcb24dbc] *** 二進数 digitsBase() [#dfd4d46b] 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() [#x92342f2] RIGHT:23November2013 > sample(0:20,3) [1] 0 11 6 0から20のうちランダムに3つの数字を重複せずに取り出す *** 数値の切り捨て [#dd710a0a] 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 になる ** 文字列 [#a747c8b5] *** 文字列の一致を検査：!is.na(match()) [#y9bdb470] RIGHT:21 March 2018 文字列が一致のときは、match が integer の 1 を返す~ 文字列が不一致のときは、match が integer の NA を返す~ NA は FALSE の代りにはならないので、is.na() を利用~ if(!is.na(match(myed1, myed2)){ # 一致したときの処理をココに } *** 文字列の表示：cat() [#i2f43c4a] RIGHT:21 March 2018 cat は concatenate （結合、連結）であって猫ではない（← お約束事として書いておく） cat("\n All Data: Courtship Latency > 4 \n") *** 文字を置換：sub() gsub() [#wb40a0fa] RIGHT:23 March 2018 mystring は、list でも、dataframeの要素でもよい~ gsub は、同じ文字列に繰り返し適用~ 正規表現が使えるらしい~ mystring <- sub("Re", "", mystring) mystring <- gsub("-", "/", mystring) # - を / に変える mystring <- gsub(".", "", mystring) # 全ての文字を消す ** 出力 [#y4a37c9a] *** ファイル出力：sink() [#wfe082b3] RIGHT:21 March 2018 sink("out.txt") # この部分のコマンドの結果がファイル出力される sink() *** pdf出力：pdf() [#lb11e9db] 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() * データの扱い [#wc23f00f] ** データの入力 [#f868ff48] *** テキストファイルから取り込む：read.table() [#g476d708] RIGHT:21 March 2018 データフレームとして取り込まれる~ header=TRUE １行目がカラム名として取り込まれる~ fill=TRUE 空のデータがあるとき、空欄が取り込まれる（事故が起こらないか注意が必要） mydf <- read.table("data.txt", header=TRUE) mydf <- read.table("data.txt", header=TRUE, fill=TRUE) *** 小さなデータなら直接書き込むほうが簡単 [#m9558310] RIGHT:21 March 2018 データの部分をスプレッドシートからコピー＆ペーストできる利点あり~ textConnection() を使って偽装することで実現できる myinput <- (" 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(myinput) <- "Crosses of D. melanogaster" # # dataframe に読み込ませる # textConnection() を使って偽装する mydf <- read.table(textConnection(myinput)) *** 複数のデータを書いて、一辺に取り込ませることもできる [#y68ed8a1] 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") # } ** 重み付けデータの扱い [#n65b28bf] 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 ... ** データの操作 [#ee33cc27] *** 行列の特定の数を置換 [#n1143f02] RIGHT:23November2013 xという名前の行列の要素のうち、0のものを-1に置換 x[(x==0)]<- -1 *** names 属性を使ってデータに名前を付ける：names() [#a9de31bd] RIGHT:21 March 2018 dataframe ではなくても可 names(mydf) <- "Crosses with D. melanogaster males " *** カラム名を変更：names() [#o0f1a07a] RIGHT:21 March 2018 dataframe ではなくても可~ 全部変更（または名付け） names(mydf2) <- c("cross", "time") ** dataframe の操作 [#oace8ba3] *** dataframe の全てのデータに代入 [#mc19cb5b] RIGHT:21 March 2018 mydf$Order <- 0 ここで、全データに同じ値を代入しようとして mydf <- 0 とすると痛い目にあう（笑）~ R は自動で型変換をするので、この場合は、各データに代入されず、 同じ名前で dataframe から数値（integer）になってしまうというオチ > mydf <- 0 > mydf [1] 0 *** dataframe に条件に一致するものだけに代入する（置換） [#a18a7df4] RIGHT:21 March 2018 mydf[mydf<0] <- NA mydf2[mydf2==0] <- NA *** dataframe に条件に一致するものだけに代入する（置換） [#z934b76a] RIGHT:23 March 2018 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 から条件に一致するものだけ取り出す [#u6fd86e5] RIGHT:23 March 2018 mydf の Factor1 と同じものを、Factor2 に持つ mydf2 のデータを選び出す~ mydf$Factor1 の代わりに上のように Mylist にもできるはず~ FALSE の後の コンマ は必須~ 一致の判断は上述~ mydf22 <- mydf2[!is.na(match(mydf2$Factor2, mydf$Factor1)),] # 行全て mydf22 <- mydf2[is.na(match(mydf2$Factor2, mydf$Factor1))==FALSE,] # 「!」でなくても mydf22 <- mydf2$Cat1[!is.na(match(mydf2$Factor2, mydf$Factor1)),] # あるカラムだけの取り出し *** dataframe から条件に一致するものだけ取り出す [#ea82d279] RIGHT:25 March 2018 条件は & や | を使ってふたつ以上の使える mydf$Cat1[mydf$Cat1 > 0] # あるカラムについて正の値のみ取り出し mydf$Cat2[mydf$Cat1 > 0] # あるカラムについて正の値の場合の、別なカラムの取り出し mydf$Cat2[!(mydf$Cat1 > 0)] # あるカラムについて正ではない場合の、別なカラムの取り出し mydf$Cat3[!(mydf$Cat1 > 0) & (mydf$Cat2 > 0)] # Cat1 が0以下 かつ Cat2 が正の場合の、Cat3の取り出し mydf$Cat3[!(mydf$Cat1 > 0) | (mydf$Cat2 > 0)] # Cat1 が0以下 または Cat2 が正の場合の、Cat3の取り出し *** dataframe から条件に一致するものだけ取り出す：subset() [#oa885987] RIGHT:21 March 2018 mydf2 <- subset(mydf, Order>0) *** dataframe から条件に一致するものだけ取り出す：subset() [#n0a84a56] 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 の factor を文字列として取り出す：as.character() [#j716403a] RIGHT:21 March 2018 これは for文の中の一部分のコードで i や j を使っての操作 mychar1 <- as.character(mydf1$factor1[i]) mychar2 <- as.character(mydf2$factor2[j]) *** dataframe から条件に一致する factor を文字列リストとして取り出す：as.character() [#i832b128] RIGHT:23 March 2018 Cat1 が 0 である factor（Myfactor）の取り出し mycharlist <- as.character(mydf$Myfactor[mydf$Cat1==0]) *** dataframe にデータを加える（横）：cbind() [#id9fd35d] RIGHT:21 March 2018 ここでは演算していますが、別な dataframe を右に付けることもできる~ サイズは一致していないとヘンなことになるので注意が必要 mydf <- cbind(mydf, mydf$N-mydf$Copulated) *** dataframe にデータを加える（下）：rbind() [#m57d5548] RIGHT:21 March 2018 カラム名が一致しないと上手くいかない mydf <- rbind( mydf, data.frame(columnname1="Dummy", columnname2=NA)) *** dataframe からカラムを取り出し、結合する：cbind() [#ld49835b] RIGHT:21 March 2018 mydf3 <- cbind(mydf2[,3:4], mydf2[,6:10]) *** カラム名を変更：names() [#l8b39ac3] RIGHT:21 March 2018 全部変更（または名付け） names(mydf2) <- c("cross","time") *** カラム名を変更：names()[] [#i35f7692] RIGHT:21 March 2018 何番目のカラムかを意識しないといけません names(mydf)[3] <- "NotC" *** 空の因子レベルを取り除く：droplevels() [#a8afd7a5] RIGHT:21 March 2018 操作した dataframe には、 元の dataframe の因子の情報（因子レベル）が含まれている~ そこで、空の因子レベルを取り除かないと正しく検定できなくなってしまう~ RWiki 因子Tips大全： http://www.okadajp.org/RWiki/?因子Tips大全 mydf2$Factor1 <- droplevels(mydf2$Factor1) *** dataframe を表示させる：show() [#o4144613] RIGHT:21 March 2018 for 文の中では表示されないので、show() で対応する~ http://d.hatena.ne.jp/yusap/20130514/1368486988 ヘのコメントを参照 show( mydf ) *** dataframe の行数のカウント：nrow() [#x979fc62] RIGHT:21 March 2018 x <- nrow(mydf) ***データの扱いについてのメモ [#o15a1bfb] 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 だけが出てくる。 * 検定を一気に行うなど [#c88c95b7] ** データのまとめ [#n70328dd] *** データのまとめ [#kc637fe9] 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 の出力：グラフの描画ではない [#sda357d0] RIGHT:21 March 2018 plot=FALSE がミソ boxplot(mydf2$time~mydf2$cross, plot=FALSE) ** ノンパラメトリック検定 [#ma6b7623] *** Kruskal-Wallis 検定 [#jf7b8007] RIGHT:21 March 2018 事後検定として Dunn 検定~ プログラム引用：[[An R Companion for the Handbook of Biological Statistics, Kruskal-Wallis Test>https://rcompanion.org/rcompanion/d_06.html]] 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 [#t503cf35] 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もそうしてみた） *** ２標本検定 [#f7575e3c] 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 [#acb5ec04] 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) ** パラメトリック検定 [#fbd6e1ba] *** n元配置：分散分析、Tukey HSD 検定 [#s26a9cb0] 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 [#fb8991fe] RIGHT:03 Oct 2005 > mydata <- matrix(c(4,6,3,4), nr=2) > fisher.test(mydata) *** Fisher test for 2 by n [#jf91e2bc] 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を大きく取らないとデータによってはエラーになる場合がある *** 分割表の独立性の検定 [#nbaf3812] 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 [#t557be03] *** カイ二乗検定 [#jf298bb4] 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検定を行いたいなら [#se0794f7] 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() [#h3f4bbf3] 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 [#e1a975ad] 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]] ** 相関 [#q3fbc9d9] *** 相関の計算・検定：cor.test() [#m1238936] RIGHT:21 March 2018 corresults <- cor.test( mydf3$CourtshipFreq, mydf3$Courtshipsec ) *** 相関を計算：cor() [#s8f0954f] 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 *** 相関を個々に計算 [#i5201e46] 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 [#k65cfd16] *** 異なる条件下で得られた複数の回帰直線を比較する [#xae766be] RIGHT:01 February 2008 -傾きに差はある？ -切片に差はある？ 長くなったのでこちら -> [[複数の回帰直線の比較>R_ComparingRegressions]] * グラフ [#hb078894] ** 散布図を描く：plot() [#dbdabe75] RIGHT:21 March 2018 まずは plot というのが R のグラフ作成のキモらしい~ dataframe が複数のデータセットになっていると、ペアワイズの複数のグラフがプロットされる plot(mydf) # plot(mydf$Data1, mydf$Data2, cex = 0.5, # ドットサイズ 0.5倍 cex.lab = 0.8, # ラベル文字サイズ 0.8倍 cex.axis = 0.7, # 軸数値文字サイズ 0.7倍 xaxp = c(0.3,1,7), # X軸目盛を 0.3 から 1 まで、７分割＝目盛を8つ yaxp = c(0,60*60,6), # Y軸目盛を 0 から 3600 まで、6分割＝目盛を7つ xlim = c(0.3,1), # X軸を 0.3 から 1 まで ylim = c(0,60*60), # Y軸を 0 から 3600 まで xlab="Courtship Frequency", ylab="Courtship Latency (s)" ) ** 棒グラフ [#h194665b] *** 棒グラフ：barplot() [#w9d4b218] 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 ) *** 棒グラフに95%信頼限界の線分を加える：barplot() arrows() [#i312077f] RIGHT:21 March 2018 mybar <- barplot(mydf$CourtshipFreq, names=mydf$Factor1, main="Courtship", cex.main=2.5, cex.lab=1.6, las=2, ylab="Frequency (%)", ylim=c(0,100), axis.lty=1, xaxt="n" ) # xaxt="n"：X軸ラベルは barplot 関数で書かず #グラフに囲みを付ける 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 などとするとヒゲも付けられる # mybar <- barplot(mydf$CopulationFreq, names=mydf$Factor1, main="Copulation", cex.main=2.5, cex.lab=1.6, cex.names=0.8, las=2, ylab="Frequency (%)", ylim=c(0,100), axis.lty=1) # X軸ラベルを barplot 関数で書く # cex.names=0.8：X軸のフォントサイズの指定 0.8倍 # axis() で操作するときは、cex.axis=0.8 とするので違いに注意（わかりづらい……） # box(lty=1) arrows(mybar, mydf$CopulationFreq, mybar, mydf$CoplCIH,angle=90,length=0,lwd=1) arrows(mybar, mydf$CopulationFreq, mybar, mydf$CoplCIL,angle=90,length=0,lwd=1) ** 箱ヒゲ図（ボックスプロット）：boxplot() [#q689e1d4] RIGHT:21 March 2018 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") # las=2：軸ラベルの向き 2はそれぞれX、Y軸に対して垂直に（Xは↑、Yは→） # outline=F：外れ値を表示しない # 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 でラベル文字列を書かない ** boxplot に蜂群を重ねる：boxplot() beeswarm() [#h7d76c00] 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) ** 複数の折れ線グラフ：matplot() [#l6bd89b8] 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)") ** 表示順指定：factor() [#k8d49c6d] RIGHT:21 March 2018 factor属性に levels引数を渡す~ c() はリストなので、順序の情報が含まれている~ → リストの順序が利用される仕組みになっている mydf2$cross <- factor(mydf2$cross, levels=c( "Data7", "Data2", "Data5", "Data3", ..., "Data1")) mydf2$cross <- factor(mydf2$cross, levels=Mylist) # Mylist が別にあるとき ** ダミーデータ追加 [#ce04aff0] RIGHT:21 March 2018 図を描くとき「空き」を作るためなどに~ 場合によっていは、名前の文字列は空白がいいのかもしれない （でもひとつしか空きは作れないはず） mydf2 <- rbind( mydf2, data.frame(cross="Dummy", time=NA)) ** グラフ領域の操作 [#bbb7d501] *** グラフの軸ラベルの向き：las= [#m37ac0b4] RIGHT:25 March 2018 las=0：それぞれX、Y軸に対して平行に（Xは→、Yは↑） las=1：X、Y軸とも全て水平に（Xは→、Yは→） las=2：それぞれX、Y軸に対して垂直に（Xは↑、Yは→） las=3：X、Y軸とも全て垂直に（Xは↑、Yは↑） *** グラフのラベルの文字サイズ：cex.***s= [#w212fd35] RIGHT:25 March 2018 デフォルトの○倍と指定 cex.main=2.5 # グラフタイトル cex.lab=1.6 # 軸ラベル（X軸、Y軸とも） cex.axis=1.2 # 軸目盛（X軸、Y軸とも） cex.names=0.8 # 棒グラフの横軸目盛（項目名）（Y軸目盛は cex.axis で指定） *** 太字、斜体の対応：substitute(paste()) italic() bold() bolditalic() [#c339acb2] 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() [#ye194b75] RIGHT:21 March 2018 paste を使って文字列を結合すれば、変数がきちんと展開される~ italic などが使いたい場合は substitute() paste() list() を利用する title(main=paste("r", " = ", coefficient, ", ", "P", " = ", pvalue)) *** グラフタイトルに変数を入れる：太字、斜体の対応も必要な場合 substitute() paste() list() [#i92108fb] 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 ) *** 軸ラベルの調整： [#y9d42d0a] RIGHT:24 March 2018 参考：http://uncorrelated.hatenablog.com/entry/20130708/1373256551 par(mgp=c(2.5, 1, 0)) # デフォルト par(mgp=c(3, 1, 0)) *** グラフ内に線分を描く：lines() [#j0f4e1df] RIGHT:21 March 2018 lines(X座標のベクトル, Y座標のベクトル)~ lines(点１,点２,点３)の形式ではないので注意 lines(c(3.5,5.4),c(95,95)) *** グラフ描画領域内に文字を書く：text() [#b5009860] 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() [#sd660c69] 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) しかし~ これを使うより、領域外への描画を許可しておき （par(xpd=NA) や par(xpd=T) として）、 text() や lines() を使うほうが使いやすいと思う *** legend の操作 [#b69f1af5] RIGHT:21 March 2018 legend に title を付けることができる~ ほかにも色々できます~ https://stats.biopapyrus.jp/r/graph/legend.html~ https://symfoware.blog.fc2.com/blog-entry-1503.html~ legend("topleft", # 位置：座標指定もできる lty=1:5, # legendの中の線種 title="Female", # コレ bty ="n", # 枠線なし # cex = 0.8, # 文字サイズ 0.8倍 bg = "transparent", # 背景を透明に：色の指定ができる legend=c( "DATA1", "DATA2", "DATA3", "DATA4", "DATA5" ) ) ** 描画領域の操作 [#ff761483] *** 描画前のおまじない [#tbf2ca69] RIGHT:24 March 2018 ファイルへの書き込みのときはなくてもよい plot.new() *** 描画領域の分割 split.screen(c()) [#l7cfd9c7] RIGHT:24 March 2018 c(横, 縦) split.screen(c(3,2)) # ３列２行 *** 描画領域を更に分割 split.screen(c(*, *), screen = *) [#ybd9b064] RIGHT:24 March 2018 分割する画面を指定する split.screen(c(1,2), screen = 1) split.screen(c(2,2), screen = 2) *** 描画面の指定 screen() [#hff7e76f] RIGHT:24 March 2018 横向きに数が振られていく~ 再分割した場合は、更に数字が増える screen(1) # ココに 1 の描画内容を screen(6) # ココに 6 の描画内容を *** 描画画面の余白 par(oma=c(9,3,0,0)) [#u9cc9c99] RIGHT:24 March 2018 split.screen() と screen() の間に指定することで、外側に余白を作れる split.screen(c(3,5)) par(oma=c(9,3,0,0)) # 余白を作る 下左上右 screen(1) *** 描画領域内のマージンの指定：下左上右 par(mar=c()) [#ybc34c60] RIGHT:24 March 2018 参考：http://uncorrelated.hatenablog.com/entry/20130708/1373256551 par(mar=c(4, 4, 1.5, 0.5) + 0.1) # デフォルト par(mar=c(5, 4, 4, 2) + 0.1) *** 分割画面の外側の操作 [#i5146bc2] RIGHT:24 March 2018 split.screen(c(3,5)) par(oma=c(0,0,3,0)) # 総合表題のための余白を作る 下左上右 # par(cex.main = 2) # title(main = "My Big Title", outer=TRUE, line=1) title(main = "My Small Title", outer=TRUE, line=2) # line= は位置の指定、文字が大きくなると重なることがあるので調整が必要 # par(cex.main = 1) # 念のためサイズを戻す *** 領域外への描画・文字の書き込み par(xpd= ) [#s3969991] RIGHT:25 March 2018 par(xpd=NA) # 領域外への描画を許可 split.screen()の外側まで par(xpd=T) # グラフ描画領域外への描画を許可 screen()の内側まで # lines(c(0.8,7.2),c(-6,-6)) text(x=-4, y=-5, labels="Control", cex=1.2, pos=1) *** 分割後に描画が終わったときのおまじない close.screen() [#c84f0cb7] RIGHT:24 March 2018 split.screen(c(3,5)) # 画面分割後の操作をここに close.screen(all.screens=TRUE) ** 『「蜂群＋平均値」に折れ線グラフを重ねる』をふたつ重ねる [#l3f48207] 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 [#uf310d46] Emacs Speaks Statistics：Emacs で R を使う~ Emacs には後からインストールをする必要あり~ [[install-memo5.1]]を参照 ** ESS起動 [#nc982a71] RIGHT:03 Oct 2005 M-x R コマンド入力後 C-c C-j C-m （改行） で、１行ずつコマンドを実行する。 *コマンドメモ [#w5e24b5d] 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） [#v5430009] 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|
