#contents *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]] **データの扱いについてのメモ 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 だけが出てくる。 *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 (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/]] |Today:&counter(today); |Yesterday:&counter(yesterday); |Total:&counter(); since 04 April 2006|