大学のサーバの引越しを機に、R の部分について整理する このページは引越し前のサーバで、引越し後には影響しないので、 サンドバッグとしての利用(内容はよいはず) #contents * Rの基本的な機能 [#e35b0122] ** R 本体(Windows版)のアップデート 22 October 2017 [#r50e5362] [[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() ** 文字列の表示:cat() [#v2250c5c] cat("\n All Data: Courtship Latency > 4 \n") ** names 属性を使ってデータに名前を付ける:names() [#lc9eebef] dataframe ではなくても可 names(mydf) <- "Crosses with D. melanogaster males " ** カラム名変換:names() [#d2e471b5] names(mydf2) <- c("cross","time") ** pdf出力:pdf() [#kf2ecf59] サイズを黄金比 \(\left(1:\frac{1+\sqrt{5}}{2}\right)\) にしてみた pdf("myout.pdf", width=(7*(1+sqrt(5))/2), height=7, pointsize=5) * データの扱い [#m5401e7d] ** タブ区切りデータの入力をする [#z1daf9d3] *** 小さなデータなら直接書き込むほうが簡単 [#s7898ad4] 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)) *** 複数のデータを書いて、一辺に取り込ませることもできる [#r682a4e1] 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") # } ** dataframe の操作 [#n592c952] *** dataframe の全てのデータに代入 [#w3d6ef30] mydf$Order <-0 *** dataframe に条件に一致するものだけに代入する(置換) [#q5af72c8] mydf[mydf<0] <- NA mydf2[mydf2==0] <- NA *** dataframe の factor を文字列として取り出す:as.character() [#i12fbf95] mychar1 <- as.character(mydf1$factor1[i]) mychar2 <- as.character(mydf2$factor2[j]) *** dataframe から条件に一致するものだけ取り出す:subset() [#b0bf3e6d] mydf2 <- subset(mydf, Order>0) *** dataframe から条件に一致するものだけ取り出す:subset() [#aac9bda0] subset= で行の抽出条件~ select= で列の抽出条件~ select= で取り出し順序の変更ができる~ Observed > 5 のデータの行を取り出すとき~ as.matrix(): dataframe を matrix に変換 mydata <- as.matrix(subset( mydf, subset=Observed>5, select=c(Courtship, Observed) ) ) *** dataframe にデータを加える(横):cbind() [#v3e1c8aa] ここでは演算していますが、別な dataframe を右に付けることもできる~ サイズは一致していないとヘンなことになるので注意が必要 mydf <- cbind(mydf, mydf$N-mydf$Copulated) *** dataframe にデータを加える(下):rbind() [#y1c769bf] カラム名が一致しないと上手くいかない mydf <- rbind( mydf, data.frame(columnname1="Dummy", columnname2=NA)) *** dataframe からカラムを取り出し、結合する:cbind() [#d0b73994] mydf3 <- cbind(mydf2[,3:4], mydf2[,6:10]) *** 空の因子レベルを取り除く:droplevels() [#zd44d35d] [[RWiki 因子Tips大全>http://www.okadajp.org/RWiki/?因子Tips大全]] mydf2$EDnum <- droplevels(mydf2$EDnum) *** カラム名を変更:names()[] [#f99bc0bb] 何番目のカラムかを意識しないといけません (mydf$N-mydf$Copulated は面倒なので) names(mydf)[3] <- "NotC" *** dataframe を表示させる:show() [#gc607260] for 文の中では表示されないので、show() で対応する~ http://d.hatena.ne.jp/yusap/20130514/1368486988 ヘのコメントを参照 show( mydf ) *** dataframe の行数のカウント:nrow() [#z2728a94] x <- nrow(mydf) ***データの扱いについてのメモ 11 July 2011 [#u9273631] 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 だけが出てくる。 ** 重み付けデータの扱い [#l614b67f] 重み付けデータは、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 ... ** Rの小技 [#i9005f2d] *** 数値の切り捨て [#qa2ae1bc] 万能ではないかもしれない~ さらに、これは、車輪の再発明になっているかもしれない (が、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 になる *** 文字列の一致を検査 [#q043ed5b] 文字列が一致のときは、match が integer の 1 を返す~ 文字列が不一致のときは、match が integer の NA を返す~ NA は FALSE の代りにはならないので、is.na() を利用~ if(!is.na((match(myed1, myed2))){ # 一致したときの処理をココに } * 検定を一気に行うセット [#h464163b] ** n元配置:分散分析、Tukey HSD 検定~ [#t1d51d2b] 分散分析を行ってから、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 ** 分割表の独立性の検定 [#b0053ebd] 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))~ - 期待値は1以上であるべき - 期待値が5未満のセルの数は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") ** Kruskal-Wallis 検定 [#tbd84eb1] 事後検定として 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) ** 頻度データの95%信頼限界の計算:binom.test() [#p9f2558c] 二項分布を利用~ 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" ) ** boxplot に使われる summary data の出力:グラフの描画ではない [#l94b46ac] plot=FALSE がミソ boxplot(mydf2$time~mydf2$cross, plot=FALSE) ** 相関の計算・検定:cor.test() [#e8d9200f] corresults <- cor.test( mydf3$CourtshipFreq, mydf3$Courtshipsec ) ** 相関を計算:cor() [#lc9be09e] dataframe に3つ以上のデータが含まれている場合 # 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 ** 相関を個々に計算 [#f3502237] dataframe に3つ以上のデータが含まれている場合 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) } } * グラフ [#y052d35c] ** バーグラフに、95%信頼限界の線分を加える [#a4b39829] mybar <- barplot(mydf$CourtshipFreq, names=mydf$EDonly, main="Courtship", cex.main=2.5, cex.lab=1.6, las=2, ylab="Frequency (%)", ylim=c(0,100), axis.lty=1, xaxt="n" ) # #グラフに囲みを付ける 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() [#l58dbdd9] factor属性に levels引数を渡す mydf2$cross <- factor(mydf2$cross, levels=c( "Data7", "Data2", "Data5", "Data3", ..., "Data1")) ** ダミーデータ追加 [#dbe9b408] 図を描くとき「空き」を作るためなどに~ 場合によっていは、名前の文字列は空白がいいのかもしれない (でもひとつしか空きは作れないはず) mydf2 <- rbind( mydf2, data.frame(cross="Dummy", time=NA)) ** グラフ領域の操作:par(oma= )) [#sbab7fd1] par(oma=c(9,3,0,0)) # 余白を作る 下左上右 ** 太字、斜体の対応:substitute(paste()) italic() bold() bolditalic() [#k76d337c] substitute と paste のコンビネーション なので、適用できない場面も多々ある title(main= substitute(paste( bold('A. '), bolditalic('D. melanogaster'), bold(' males')) ), cex.main=1.5, cex=1.5, ) ** グラフタイトルに変数を入れる:paste() [#ebd6fc82] paste を使って文字列を結合すれば、変数がきちんと展開される italic などが使いたい場合は list を利用する title(main=paste("r", " = ", coefficient, ", ", "P", " = ", pvalue)) ** グラフタイトルに変数を入れる:太字、斜体の対応も必要な場合 substitute() paste() list() [#we4acad1] 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 の操作 [#za470856] leged に title を付けることができる legend("topleft", lty=1:5, title="Female", bty ="n", # 枠線なし # cex = 0.8, # 文字サイズ 0.8倍 legend=c( "DATA1", "DATA2", "DATA3", "DATA4", "DATA5" ) ) ** 複数の折れ線グラフ:matplot() [#ed986a4c] 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() [#v1226dba] 文字も入れる 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() [#d3e9e1bf] 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() [#v87729b9] まずは plot というのが R のグラフ作成らしい plot(mydf) ** グラフ描画領域内に文字を書く:text() [#r1e6c65d] adj=c(0.5,0) # テキストをxy座標中央揃え、下揃え~ x 座標については、0 左揃え、0.5 中央揃え、1 右揃え らしい~ y 座標については、0 下揃え、0.5 中央揃え、1 上揃え らしい~ 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() [#s114b6f8] 描画領域のサイズが変わると位置が変わるので注意~ 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() [#td349dce] lines(X座標のベクトル, Y座標のベクトル)~ lines(点1,点2,点3)の形式ではないので注意 lines(c(3.5,5.4),c(95,95)) ** 散布図を描き、相関を計算、グラフの上に表示 [#w9295a19] plot(mydf3$CourtshipFreq, mydf3$Courtshipsec, cex = 0.5, # ドットサイズ 0.5倍 # cex.lab = 0.8, # ラベル文字サイズ 0.8倍 cex.axis = 0.7, # 軸数値文字サイズ ○倍 xaxp = c(0.3,1,7), yaxp = c(0,60*60,6), xlim = c(0.3,1), ylim = c(0,60*60), xlab="", # xlab="Courtship Frequency", ylab="Courtship Latency (s)" ) ** 『「蜂群+平均値」に折れ線グラフを重ねる』をふたつ重ねる 2017年8月7日 [#mc957cf6] - 横軸はカテゴリ - 一番左のカテゴリは少し性質が違うので、左からひとつめと左からふたつめをつなぐ折れ線は他と別の破線で - 『「蜂群+平均値」に折れ線グラフを重ねる』をずらしてふたつ重ねる - 平均値はグラフ作成の途中で計算させる - 「ずらし」を沢山使うので、カテゴリを 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() ** survival analysis (11 April 2004)[#h9e1c007] ショウジョウバエの配偶行動を見るときに、観察時間内に交尾しないものがある。そこで、交尾したものを「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]] ** Comparing Regressions [#md8494d1] *** 異なる条件下で得られた複数の回帰直線を比較する(01 February 2008) [#rda59cbd] -傾きに差はある? -切片に差はある? 長くなったのでこちら -> [[複数の回帰直線の比較>R_ComparingRegressions]] ** Contingency Table [#t9c03134] *** Fisher test (03 Oct 2005) [#te11009a] > mydata <- matrix(c(4,6,3,4), nr=2) > fisher.test(mydata) *** Fisher test for 2 by n (17 September 2010) [#kb93275c] > 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を大きく取らないとデータによってはエラーになる場合がある ** Goodness of Fit (14 January 2016) [#n418c3d5] 適合性の検定 メンデルの法則などのようなものの検定 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]] *** カイ二乗検定 [#o2122b54] > 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検定を行いたいなら [#g69a6ee3] > 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~ **ノンパラメトリック検定 [#t1badf9b] *** Kruskal-Wallis test と multiple comaprison 2017年10月11日 [#c4fb3ebb] 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もそうしてみた) *** 2標本検定 2017年8月7日 [#u79b3da1] - 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 (09 May 2004) [#a267b855] 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) ** パラメトリック検定 [#ta052279] ** データのまとめ (09 May 2004) [#xaf6ba55] 各$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]] *tips [#afeef77e] **二進数 (23November2013) [#hf255012] 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) ← ベクトルの長さを指定したいとき 二進数に限らず利用できる * 重複しないランダムな数 (23November2013) [#n84fa29a] > sample(0:20,3) [1] 0 11 6 0から20のうちランダムに3つの数字を重複せずに取り出す **行列の特定の数を置換 (23November2013) [#tb5633a2] x[(x==0)]<-- -1 xという名前の行列の要素のうち、0のものを-1に置換 *ESS [#s80ba893] ** 03 Oct 2005 [#q15bf9ff] iESS起動 M-x R コマンド入力後 C-c C-j C-m (改行) で、1行ずつコマンドを実行する。 *コマンドメモ (30 March 2004) [#ve4bd0d4] > 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) 05September2013 [#q746f9e7] 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/]]