R_statistics の変更点

• 追加された行はこの色です。
• 削除された行はこの色です。
• R_statistics へ行く。
• R_statistics の差分を削除

---

#author("2019-07-21T10:22:36+09:00","","")
#author("2019-07-21T10:30:29+09:00","","")
#contents

* 検定を一気に行うなど [#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&#8211;25. Doi: 10.1002/(SICI)1521-4036(200001)42:1<17::AID-BIMJ17>3.0.CO;2-U
- 名取真人. (2014). マン・ホイットニーのU検定と不等分散時における代表値の検定法. 霊長類研究, 30(1), 173&#8211;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]

*** 等分散の検定：Bartlett 検定など [#x1061cc4]
RIGHT:20 July 2019

> なお，事前検定を行うことが不適切であることはだんだん理解されてきているので，この観点から言えば「等分散検定後に普通の t 検定」というのは好ましくない。分散が等しかろうと等しくなかろうと，最初からズバリ「等分散を仮定しない t 検定」を行うのが正しいやり方である。
https://oku.edu.mie-u.ac.jp/~okumura/stat/ttest.html

ではあるのだけれど、検定をする必要もあるかもしれない。

また、http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc043/12857.html
も参考に。


参考: http://www.cookbook-r.com/Statistical_analysis/Homogeneity_of_variance/

 # Homogeneity of variance
 bartlett.test(Courtshipsec ~ EDnum, data=mydf2)
 library(car)
 leveneTest(Courtshipsec ~ EDnum, data=mydf2)
 fligner.test(Courtshipsec ~ EDnum, data=mydf2)



*** 一元配置の分散分析 oneway.test() [#b2d12591]
RIGHT:20 July 2019

デフォルトで等分散を仮定しない検定を行う~
referee から文句を言われない／標準が示されていない場合は、
等分散の検定を''行わず''、こちらを使うのがよいと思う~

https://oku.edu.mie-u.ac.jp/~okumura/stat/ttest.html

http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc043/12857.html

 oneway.test(Courtshipsec ~ EDnum, data=mydf2)

更に……~
文句を言われないなら、分散分析はせずに Tukey HSD 検定がよさそう((Zar, J. H. (2010) Biostatistical Analysis, 5th ed, Pearson Prentice-Hall, Upper Saddle River, NJ))


*** 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 &20 July 2019

 mydata <- matrix(c(4,6,3,4), nr=2)
 fisher.test(mydata)
 fisher.test(mydata)$p.value

*** Fisher test：mid-P [#xb36b0e9]
RIGHT:20 July 2019

2 &times; 2 分割表の検定においては：~
[[Lydersen et al. (2009) Recommended tests for association in 2×2 tables. Statist. Med. 28:1159--1175>http://www.doi.org/10.1002/sim.3531]]
は現時点では、unconditional test か mid-p を推奨~
[[Ruxton & Neuhauser (2010) Good practice in testing for an association in contingency tables. Behav. Ecol. Sociobiol 64:1505--1513>http://www.doi.org/10.1007/s00265-010-1014-0]]
は、現時点（2010）では、Yates 補正なしカイ自乗検定か mid-p を推奨~
上の２つの論文とも、Barnard exact test はいいかもしれないとも書いている~
実際には、期待値が5より小さいものも多いので、
カイ自乗検定は使えないデータセットもある

 library(exact2x2)
 mydata <- matrix(c(4,6,3,4), nr=2)
 fisher.exact(mydata, midp=TRUE)
 fisher.exact(mydata, midp=TRUE)$p.value
 fisher.exact(mydata) # mid-P 補正なし確率 fisher.test() と同じ結果を出力する


*** Barnard's exact test [#f19a9e73]
RIGHT:20 July 2019

Unconditional exact test のひとつ~
Fisher exact test は行の合計と列の合計が固定~
Barnard's exact test は一方（列）のみ固定~
２つの比を比較する検定となる。
Fisher exact test はオッズ比が１かどうかの検定。~
データの入力に際して縦横に注意が必要

 library(Exact)
 mydata <- matrix(c(4,6,3,4), nr=2)
 # 4/7 と 6/10 を比較している
 exact.test(mydata, to.plot=FALSE)
 exact.test(mydata, to.plot=FALSE)$p.value



*** 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~


** 二項検定：binom.test() [#if7b43a6]
RIGHT:20 July 2019

確率は $p.value に格納される

 mydata <- c(34, 27)
 binom.test(mydata)
 binom.test(mydata)$p.value


*** 頻度データの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 & 20 July 2019

Pearson の積率相関係数

 corresults <- cor.test(
    mydf3$CourtshipFreq,
    mydf3$Courtshipsec
 )
 cor.test( ~ startmin + Court,
           data=mydfsim,
           method = "pearson",
           conf.level = 0.95)
 )

Kendall のノンパラメトリック相関係数
 cor.test( ~ startmin + Court,
          data=mydfsim,
          method = "kendall",
          continuity = FALSE,
          conf.level = 0.95)


*** 相関を計算：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)



*** 百分率表示する累積折れ線グラフ [#i9777e39]
RIGHT:20 July 2019

参考：http://www.singularpoint.org/blog/r/r-plot-step-function/

パラメータは「...」としているので、
普通の plot()関数で指定できるものは（多分全て）使える

 # http://www.singularpoint.org/blog/r/r-plot-step-function/
 # cumplot <- function(x,...){
 #   plot(sort(x),(1:length(x))/length(x),type="s",...)
 # }
 # 
 mycumplot2 <- function(x, mymaxx, ...){
    mynobs <- length(x)
    nona <- x[!is.na(x)]
    # NA があるときに対応する
    plot(c(0,sort(nona),mymaxx),(c(0,(1:length(nona))/mynobs,length(nona)/mynobs)*100),...)
    # type 指定は関数の外で行うために type="s" を除く＝デフォルト type="p" 
    # 一番右をX軸の最大値（mymaxx）まで延ばす
    # (0,0) を追加し、原点と一番左のデータ点を繋ぐ
 }


** 複数の折れ線グラフ：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, )

*** ギリシャ文字などの対応：expression(paste()) [#c53315a4]
RIGHT:20 July 2019

ギリシャ文字などの指定は、latex の数式モードでの指定と同じ綴りだけど、
「\（バックスラッシュ）」は不要~
「^」 で上付き~
「[]」 で下付き

 text(x=2, y=103, labels=expression(paste(chi[3]^2, "=7.129 ", italic('P'), "=0.067 NS")), cex=1.1, pos=3)


*** グラフタイトルに変数を入れる：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"
                 )
	)


*** 矩形の描画: rect() [#h8d0fa40]
RIGHT:20 July 2019

 rect(-95, 53.2, 60, 48.2) # 色：透過
 rect(0, 0, 60, 48.2, col="gray") # 色：グレー
 # rect(x1,y1,x2,y2)

*** 線分の描画: line() [#cf6bba44]
RIGHT:20 July 2019

 lines(c(0,1), c(4,3))
 # lines(x,y)
 # x, y はベクトル


*** グラデーションの灰色を算出: gray.colors() [#tded589a]
RIGHT:20 July 2019

 gray.colors(4, start =0, end = 1)
 # [1] "#000000" "#9B9B9B" "#D4D4D4" "#FFFFFF"


*** 色の指定 [#k32909d6]
RIGHT:20 July 2019

色々な場面で、色の指定をするときに
  
 col="#000000" # 黒
 col="#9B9B9B" 
 col="#D4D4D4"
 col="#FFFFFF" # 白
 col="yellow" # 黄色


** 描画領域の操作 [#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)


*** 描画領域内のパラメータの取得 [#f52cf791]
RIGHT:20 July 2019

デフォルトのグラフィックパラメタ（沢山！）を保存~
本当はその時のパラメタだけど、R を起動して直ぐならデフォルトパラメタ~
デフォルトの値を使って描画領域内のマージンを指定できる

 defaultpar <- par(no.readonly=TRUE) 
 par(mar=defaultpar$mar + c(4,0,1.28,0)) # マージンの追加  下左上右


*** 分割画面の外側の操作 [#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)


*** グラフのプロットのマーカーの変更 [#ec2a29b3]
RIGHT:20 July 2019

 # pch=0 四角□
 # pch=1 丸○ 
 # pch=2 三角△ 
 # pch=3 十字＋ 
 # pch=4 バツ× 
 # pch=5 ダイヤ◇ 
 # pch=6 逆三角▽ 
 # pch=7 □に× 
 # pch=8 ＋に×  
 # pch=9 ◇に＋  
 # pch=10 ○に＋   
 # まだあるが……



** 『「蜂群＋平均値」に折れ線グラフを重ねる』をふたつ重ねる [#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()

|Today:&counter(today); |Yesterday:&counter(yesterday); |Total:&counter(); since 21 July 2019|

       

Site admin: tomaru

PukiWiki 1.5.1 © 2001-2016 PukiWiki Development Team. Powered by PHP 5.4.16. HTML convert time: 0.024 sec.