* Rの基本的な機能 [#o2dc5063] 大学のサーバの引越しを機に、R の部分について整理する このページは引越し前のサーバで、引越し後には影響しないので、 サンドバッグとしての利用(内容はよいはず) #contents ** 文字列の表示 [#g5f8ec12] cat("\n All Data: Courtship Latency > 4 \n") ** names 属性を使ってデータに名前を付ける [#c9d1ccfa] dataframe ではなくても可 names(mydf) <- "Crosses with D. melanogaster males " ** カラム名変換 [#fbcd2bb9] names(mydf2) <- c("cross","time") ** タブ区切りデータの入力をする [#o796d63e] *** 小さなデータなら直接書き込むほうが簡単 [#gc7fad50] 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)) *** 複数のデータを書いて、一辺に取り込ませることもできる [#xad6a093] 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") # } ** pdf出力 [#rbfb9e14] サイズを黄金比にしてみた (7*(1+sqrt(5)) pdf("myout.pdf", width=(7*(1+sqrt(5))/2), height=7, pointsize=5) ** dataframe の操作 [#h2305982] - dataframe の全てのデータに代入 mydf$Order <-0 *** dataframe に条件に一致するものだけに代入する(置換) [#z55e5a76] mydf[mydf<0] <- NA mydf2[mydf2==0] <- NA *** dataframe の factor を文字列として取り出す [#z0bf4719] mychar1 <- as.character(mydf1$factor1[i]) mychar2 <- as.character(mydf2$factor2[j]) *** dataframe から条件に一致するものだけ取り出す [#x23e2cd5] mydf2 <- subset(mydf, Order>0) *** dataframe から条件に一致するものだけ取り出す [#ba2f14a5] subset= で行の抽出条件~ select= で列の抽出条件~ select= で取り出し順序の変更ができる~ Observed > 5 のデータの行を取り出すとき~ as.matrix(): dataframe を matrix に変換 mydata <- as.matrix(subset(mydf, subset=Observed>5, select=c(Courtship, Observed)) ) *** dataframe にデータを加える(横) [#ob7a7210] ここでは演算していますが、別な dataframe を右に付けることもできる~ サイズは一致していないとヘンなことになるので注意が必要 mydf <- cbind(mydf, mydf$N-mydf$Copulated) *** dataframe にデータを加える(下) [#led44ee7] カラム名が一致しないと上手くいかない mydf <- rbind( mydf, data.frame(columnname1="Dummy", columnname2=NA)) *** dataframe からカラムを取り出し、結合する [#mea51feb] mydf3 <- cbind(mydf2[,3:4], mydf2[,6:10]) *** 空の因子レベルを取り除く [#r4dfa170] http://www.okadajp.org/RWiki/?因子Tips大全 mydf2$EDnum <- droplevels(mydf2$EDnum) *** カラム名を変更(mydf$N-mydf$Copulated は面倒なので) [#y060442c] 何番目のカラムかを意識しないといけません names(mydf)[3] <- "NotC" *** dataframe を表示させる [#ld68f362] for 文の中では表示されないので、show() で対応する~ http://d.hatena.ne.jp/yusap/20130514/1368486988 ヘのコメントを参照 show( mydf ) *** dataframe の行数のカウント [#t441c2ed] x <- nrow(mydf) ** Rの小技 [#qe5c0075] *** 数値の切り捨て [#t15e7cf2] 万能ではないかもしれない~ さらに、これは、車輪の再発明になっているかもしれない(が、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 になる *** 文字列の一致を検査 [#ab517ef7] 文字列が一致のときは、match が integer の 1 を返す~ 文字列が不一致のときは、match が integer の NA を返す~ NA は FALSE の代りにはならないので、is.na() を利用~ if(!is.na((match(myed1, myed2))){ 一致したときの処理をココに } ** 検定を一気に行うセット [#yf4c2a11] *** n元配置 [#i3bfdf24] 分散分析を行ってから、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 *** 分割表の独立性の検定 [#kf073163] Fisher Exact Test:サンプルサイズが大きいときは workspace を大きく取っても計算ができなくなるので注意が必要 プログラム引用: An R Companion for the Handbook of Biological Statistics~ Salvatore S. Mangiafico~ 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") # # library(DescTools) # # 2x2以上のときの事後検定 # 参考にしたのは # 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 検定 [#ya77eb97] 事後検定として 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%信頼限界の計算 [#v754aa69] 二項分布を利用~ 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, "\t",100*mybntest$conf.int[1], "\t",100*mybntest$conf.int[2],"\t") *** boxplot に使われる summary data の出力:グラフの描画ではない [#a781c894] plot=FALSE がミソ boxplot(mydf2$time~mydf2$cross, plot=FALSE) *** 相関の計算・検定 [#e075880a] corresults <- cor.test( mydf3$CourtshipFreq, mydf3$Courtshipsec ) *** 相関を計算 [#q225d161] # 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 *** 相関を個々に計算 [#gc907728] 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) } } ** 重み付けデータの扱い [#sb8c5015] 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 ... ** グラフ [#nae8c2bf] *** バーグラフに、95%信頼限界の線分を加える [#t7f122e8] 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 などとするとヒゲも付けられる *** 表示順指定:cross(factor属性に levels引数を渡す) [#gf62f89b] mydf2$cross <- factor(mydf2$cross, levels=c( "Data7", "Data2", "Data5", "Data3", ..., "Data1")) *** ダミーデータ追加 [#wac8ebcf] 図を描くとき「空き」を作るためなどに~ 場合によっていは、名前の文字列は空白がいいのかもしれない (でもひとつしか空きは作れないはず) mydf2 <- rbind( mydf2, data.frame(cross="Dummy", time=NA)) *** グラフ領域の操作 [#zbcd566a] par(oma=c(9,3,0,0)) # 余白を作る 下左上右 *** 複数の折れ線グラフ [#sf0c5b0e] 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)") *** 太字、斜体の対応 [#j40600c8] substitute と paste のコンビネーション なので、適用できない場面も多々ある title(main=substitute(paste(bold('A. '), bolditalic('D. melanogaster'), bold(' males'))), cex.main=1.5, cex=1.5, ) *** グラフタイトルに変数を入れる [#h8aab618] paste を使って文字列を結合すれば、変数がきちんと展開される italic などが使いたい場合は list を利用する title(main=paste("r", " = ", coefficient, ", ", "P", " = ", pvalue)) *** グラフタイトルに変数を入れる:太字、斜体の対応も必要な場合 [#b56ec8de] 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 の操作 [#x5295283] leged に title を付けることができる legend("topleft", lty=1:5, title="Female", bty ="n", # 枠線なし # cex = 0.8, # 文字サイズ 0.8倍 legend=c( "DATA1", "DATA2", "DATA3", "DATA4", "DATA5" ) ) *** boxplot に蜂群を重ねる [#vf6c1756] 文字も入れる 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 [#w5d68042] 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 ) *** 散布図を描く [#r3411ddc] plot(mydf3) *** グラフ描画領域内に文字を書く [#af64be4e] 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)) *** 軸ラベルの下に文字を書く [#i0d813ca] 描画領域のサイズが変わると位置が変わるので注意~ 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) *** グラフ内に線分を描く [#jd10a948] lines(c(3.5,5.4),c(95,95)) *** 散布図を描き、相関を計算、グラフの上に表示 [#q033a8ce] 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)" )