- 追加された行はこの色です。
- 削除された行はこの色です。
* 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)"
)