* Rの基本的な機能 [#o35047ac] 大学のサーバの引越しを機に、R の部分について整理する このページは引越し前のサーバで、引越し後には影響しないので、 サンドバッグとしての利用(内容はよいはず) #contents ** 文字列の表示:cat() [#x230be7c] cat("\n All Data: Courtship Latency > 4 \n") ** names 属性を使ってデータに名前を付ける:names() [#ka0d7e39] dataframe ではなくても可 names(mydf) <- "Crosses with D. melanogaster males " ** カラム名変換:names() [#b6a7c507] names(mydf2) <- c("cross","time") ** タブ区切りデータの入力をする [#c46a79bf] *** 小さなデータなら直接書き込むほうが簡単 [#q01a228d] 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)) *** 複数のデータを書いて、一辺に取り込ませることもできる [#m7692466] 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出力 [#a5300046] サイズを黄金比にしてみた (7*(1+sqrt(5)) pdf("myout.pdf", width=(7*(1+sqrt(5))/2), height=7, pointsize=5) ** dataframe の操作 [#l5d9a5b8] *** dataframe の全てのデータに代入 [#jd9a7704] mydf$Order <-0 *** dataframe に条件に一致するものだけに代入する(置換) [#g2814d20] mydf[mydf<0] <- NA mydf2[mydf2==0] <- NA *** dataframe の factor を文字列として取り出す:as.character() [#l6952722] mychar1 <- as.character(mydf1$factor1[i]) mychar2 <- as.character(mydf2$factor2[j]) *** dataframe から条件に一致するものだけ取り出す:subset() [#q420ef6f] mydf2 <- subset(mydf, Order>0) *** dataframe から条件に一致するものだけ取り出す:subset() [#k8238b2e] subset= で行の抽出条件~ select= で列の抽出条件~ select= で取り出し順序の変更ができる~ Observed > 5 のデータの行を取り出すとき~ as.matrix(): dataframe を matrix に変換 mydata <- as.matrix(subset( mydf, subset=Observed>5, select=c(Courtship, Observed) ) ) *** dataframe にデータを加える(横):cbind() [#q0d4be6e] ここでは演算していますが、別な dataframe を右に付けることもできる~ サイズは一致していないとヘンなことになるので注意が必要 mydf <- cbind(mydf, mydf$N-mydf$Copulated) *** dataframe にデータを加える(下):rbind() [#ycaf986e] カラム名が一致しないと上手くいかない mydf <- rbind( mydf, data.frame(columnname1="Dummy", columnname2=NA)) *** dataframe からカラムを取り出し、結合する:cbind() [#x3ad4222] mydf3 <- cbind(mydf2[,3:4], mydf2[,6:10]) *** 空の因子レベルを取り除く:droplevels() [#r7ac0646] [RWiki 因子Tips大全>http://www.okadajp.org/RWiki/?因子Tips大全] mydf2$EDnum <- droplevels(mydf2$EDnum) *** カラム名を変更(mydf$N-mydf$Copulated は面倒なので) [#b59ae843] 何番目のカラムかを意識しないといけません names(mydf)[3] <- "NotC" *** dataframe を表示させる:show() [#ieb9ceac] for 文の中では表示されないので、show() で対応する~ http://d.hatena.ne.jp/yusap/20130514/1368486988 ヘのコメントを参照 show( mydf ) *** dataframe の行数のカウント:nrow() [#kf4e8dc0] x <- nrow(mydf) ** Rの小技 [#q2c2e62b] *** 数値の切り捨て [#e4c9706c] 万能ではないかもしれない~ さらに、これは、車輪の再発明になっているかもしれない (が、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 になる *** 文字列の一致を検査 [#m2a4c629] 文字列が一致のときは、match が integer の 1 を返す~ 文字列が不一致のときは、match が integer の NA を返す~ NA は FALSE の代りにはならないので、is.na() を利用~ if(!is.na((match(myed1, myed2))){ # 一致したときの処理をココに } ** 検定を一気に行うセット [#we5281bc] *** n元配置:分散分析、Tukey HSD 検定~ [#y9918bf8] 分散分析を行ってから、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 *** 分割表の独立性の検定 [#cbd5cf1e] 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 検定 [#pc6cd43f] 事後検定として 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() [#ddb7764b] 二項分布を利用~ 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 の出力:グラフの描画ではない [#z11ad082] plot=FALSE がミソ boxplot(mydf2$time~mydf2$cross, plot=FALSE) *** 相関の計算・検定:cor.test() [#r07dbf82] corresults <- cor.test( mydf3$CourtshipFreq, mydf3$Courtshipsec ) *** 相関を計算:cor() [#fce9b4d3] 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 *** 相関を個々に計算 [#o4d04c79] 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) } } ** 重み付けデータの扱い [#r9e4830a] 重み付けデータは、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 ... ** グラフ [#y37e3a39] *** バーグラフに、95%信頼限界の線分を加える [#ha11290e] 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() [#n7aa2181] factor属性に levels引数を渡す mydf2$cross <- factor(mydf2$cross, levels=c( "Data7", "Data2", "Data5", "Data3", ..., "Data1")) *** ダミーデータ追加 [#nb6f1319] 図を描くとき「空き」を作るためなどに~ 場合によっていは、名前の文字列は空白がいいのかもしれない (でもひとつしか空きは作れないはず) mydf2 <- rbind( mydf2, data.frame(cross="Dummy", time=NA)) *** グラフ領域の操作:par(oma= )) [#m58badc7] par(oma=c(9,3,0,0)) # 余白を作る 下左上右 *** 太字、斜体の対応:substitute(paste()) italic() bold() bolditalic() [#lfcc96a9] substitute と paste のコンビネーション なので、適用できない場面も多々ある title(main= substitute(paste( bold('A. '), bolditalic('D. melanogaster'), bold(' males')) ), cex.main=1.5, cex=1.5, ) *** グラフタイトルに変数を入れる:paste() [#f228c0bb] paste を使って文字列を結合すれば、変数がきちんと展開される italic などが使いたい場合は list を利用する title(main=paste("r", " = ", coefficient, ", ", "P", " = ", pvalue)) *** グラフタイトルに変数を入れる:太字、斜体の対応も必要な場合 substitute() paste() list() [#c1ea9f81] 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 の操作 [#j6823453] leged に title を付けることができる legend("topleft", lty=1:5, title="Female", bty ="n", # 枠線なし # cex = 0.8, # 文字サイズ 0.8倍 legend=c( "DATA1", "DATA2", "DATA3", "DATA4", "DATA5" ) ) *** 複数の折れ線グラフ:matplot() [#m40e4a9d] 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() [#gc5a702c] 文字も入れる 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() [#fb57f1f6] 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() [#a0c47540] まずは plot というのが R のグラフ作成らしい plot(mydf) *** グラフ描画領域内に文字を書く:text() [#cbd02cf6] 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() [#a8983a3a] 描画領域のサイズが変わると位置が変わるので注意~ 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() [#fa9f2b41] lines(X座標のベクトル, Y座標のベクトル)~ lines(点1,点2,点3)の形式ではないので注意 lines(c(3.5,5.4),c(95,95)) *** 散布図を描き、相関を計算、グラフの上に表示 [#de725416] 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)" )