このページは分割しました 21 July 2019

以下は分割前のページ


Rの基本的な機能 & tips

R 本体(Windows版)のアップデート

22 October 2017

Updating R from R (on Windows) – using the {installr} package
Windows用Rのバージョンアップ

 install.packages("installr")
 library(installr)
 updateR() 

Object の調査:str()

21 March 2018

str() 関数はいろいろできるらしいのだけれど……

str(obj)

数の扱い

二進数 digitsBase()

23November2013

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) ← ベクトルの長さを指定したいとき

二進数に限らず利用できる

重複しないランダムな数 sample()

23November2013
> sample(0:20,3)
[1]  0 11  6

0から20のうちランダムに3つの数字を重複せずに取り出す

数値の切り捨て

21 March 2018

万能ではないかもしれない
さらに、これは、車輪の再発明になっているかもしれない (が、round down が見つけられなかったのよ)

sprintf を使うと、round されてしまう
統計量は round ではなく切り捨てがよいことも多いので (四捨五入したら有意になってしまうかもしれない*1)、 文字列の操作とすることにして round down を実現する
文字数を数えて、最終的な文字列の長さを決める
元々の小数点以下桁数がわからないものを、小数点以下 dp桁に切り捨てる
マイナス記号が付くと文字列が長くなるので……
関数化しておくとよいのはわかるけれど、 後々のメンテナンスを考え(ブラックボックス化すると面倒なので) 関数化はしないでおく*2
切り捨て桁(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 になる

文字列

文字列の一致を検査:!is.na(match())

21 March 2018

文字列が一致のときは、match が integer の 1 を返す
文字列が不一致のときは、match が integer の NA を返す
NA は FALSE の代りにはならないので、is.na() を利用

if(!is.na(match(myed1, myed2)){
    # 一致したときの処理をココに
}

文字列の表示:cat()

21 March 2018

cat は concatenate (結合、連結)であって猫ではない(← お約束事として書いておく)

cat("\n
All Data: Courtship Latency > 4
\n")

文字を置換:sub() gsub()

23 March 2018

mystring は、list でも、dataframeの要素でもよい
gsub は、同じ文字列に繰り返し適用
正規表現が使えるらしい

mystring <- sub("Re", "", mystring)
mystring <- gsub("-", "/", mystring) # - を / に変える
mystring <- gsub(".", "", mystring) # 全ての文字を消す

出力

ファイル出力:sink()

21 March 2018
sink("out.txt")
# この部分のコマンドの結果がファイル出力される
sink()

pdf出力:pdf()

21 March 2018

サイズを黄金比 \(\left(1:\frac{1+\sqrt{5}}{2}\right)\) にしてみた

pdf("myout.pdf", width=(7*(1+sqrt(5))/2), height=7, pointsize=5)
# この部分の描画の結果が出力される
dev.off()

データの扱い

データの入力

テキストファイルから取り込む:read.table()

21 March 2018

データフレームとして取り込まれる
header=TRUE 1行目がカラム名として取り込まれる
fill=TRUE 空のデータがあるとき、空欄が取り込まれる(事故が起こらないか注意が必要)

mydf <- read.table("data.txt", header=TRUE)
mydf <- read.table("data.txt", header=TRUE, fill=TRUE)

小さなデータなら直接書き込むほうが簡単

21 March 2018

データの部分をスプレッドシートからコピー&ペーストできる利点あり
textConnection() を使って偽装することで実現できる

myinput <- ("
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(myinput) <- "Crosses of D. melanogaster"
#
# dataframe に読み込ませる
# textConnection() を使って偽装する
mydf <- read.table(textConnection(myinput))

複数のデータを書いて、一辺に取り込ませることもできる

21 March 2018

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")
#
}

重み付けデータの扱い

21 March 2018

重み付けデータは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
... 

データの操作

行列の特定の数を置換

23November2013

xという名前の行列の要素のうち、0のものを-1に置換

x[(x==0)]<- -1

names 属性を使ってデータに名前を付ける:names()

21 March 2018

dataframe ではなくても可

names(mydf) <- "Crosses with D. melanogaster males "

カラム名を変更:names()

21 March 2018

dataframe ではなくても可
全部変更(または名付け)

names(mydf2) <- c("cross", "time")

dataframe の操作

dataframe の全てのデータに代入

21 March 2018
mydf$Order <- 0

ここで、全データに同じ値を代入しようとして mydf <- 0 とすると痛い目にあう(笑)
R は自動で型変換をするので、この場合は、各データに代入されず、 同じ名前で dataframe から数値(integer)になってしまうというオチ

> mydf <- 0
> mydf
[1] 0

dataframe に条件に一致するものだけに代入する(置換)

21 March 2018
mydf[mydf<0] <- NA
mydf2[mydf2==0] <- NA

dataframe に条件に一致するものだけに代入する(置換)

23 March 2018

factor(ここでは MyFactor)が、Mylist のいずれかと一致する場合、Cat1 に 0 を代入
for loop を使わずに1行でできてしまう ← すごい
文字列が一致のときは、match が integer の 1 を返す
文字列が不一致のときは、match が integer の NA を返す
NA は FALSE の代りにはならないので、is.na() を利用

mydf$Cat1[is.na(match(mydf$MyFactor, Mylist))==FALSE] <- 0

dataframe から条件に一致するものだけ取り出す

23 March 2018

mydf の Factor1 と同じものを、Factor2 に持つ mydf2 のデータを選び出す
mydf$Factor1 の代わりに上のように Mylist にもできるはず
FALSE の後の コンマ は必須
一致の判断は上述

mydf22 <- mydf2[!is.na(match(mydf2$Factor2, mydf$Factor1)),]       # 行全て
mydf22 <- mydf2[is.na(match(mydf2$Factor2, mydf$Factor1))==FALSE,] # 「!」でなくても
mydf22 <- mydf2$Cat1[!is.na(match(mydf2$Factor2, mydf$Factor1)),]  # あるカラムだけの取り出し

dataframe から条件に一致するものだけ取り出す

25 March 2018

条件は & や | を使ってふたつ以上の使える

mydf$Cat1[mydf$Cat1 > 0] # あるカラムについて正の値のみ取り出し
mydf$Cat2[mydf$Cat1 > 0] # あるカラムについて正の値の場合の、別なカラムの取り出し
mydf$Cat2[!(mydf$Cat1 > 0)] # あるカラムについて正ではない場合の、別なカラムの取り出し
mydf$Cat3[!(mydf$Cat1 > 0) & (mydf$Cat2 > 0)] # Cat1 が0以下 かつ Cat2 が正の場合の、Cat3の取り出し
mydf$Cat3[!(mydf$Cat1 > 0) | (mydf$Cat2 > 0)] # Cat1 が0以下 または Cat2 が正の場合の、Cat3の取り出し

dataframe から条件に一致するものだけ取り出す:subset()

21 March 2018
mydf2 <- subset(mydf, Order>0)

dataframe から条件に一致するものだけ取り出す:subset()

21 March 2018

subset= で行の抽出条件 ← Observed > 5 のデータの行を取り出すとき
select= で列の抽出条件、取り出し順序の変更もできる
as.matrix(): dataframe を matrix に変換

mydata <- as.matrix(subset(
                           mydf,
                           subset=Observed>5,
                           select=c(Courtship, Observed)
                          )
                   )

dataframe の factor を文字列として取り出す:as.character()

21 March 2018

これは for文の中の一部分のコードで i や j を使っての操作

mychar1 <- as.character(mydf1$factor1[i])
mychar2 <- as.character(mydf2$factor2[j])

dataframe から条件に一致する factor を文字列リストとして取り出す:as.character()

23 March 2018

Cat1 が 0 である factor(Myfactor)の取り出し

mycharlist <- as.character(mydf$Myfactor[mydf$Cat1==0])

dataframe にデータを加える(横):cbind()

21 March 2018

ここでは演算していますが、別な dataframe を右に付けることもできる
サイズは一致していないとヘンなことになるので注意が必要

mydf <- cbind(mydf, mydf$N-mydf$Copulated)

dataframe にデータを加える(下):rbind()

21 March 2018

カラム名が一致しないと上手くいかない

mydf <- rbind( mydf, data.frame(columnname1="Dummy", columnname2=NA))

dataframe からカラムを取り出し、結合する:cbind()

21 March 2018
mydf3 <- cbind(mydf2[,3:4], mydf2[,6:10])

カラム名を変更:names()

21 March 2018

全部変更(または名付け)

names(mydf2) <- c("cross","time")

カラム名を変更:names()[]

21 March 2018

何番目のカラムかを意識しないといけません

names(mydf)[3] <- "NotC"

空の因子レベルを取り除く:droplevels()

21 March 2018

操作した dataframe には、 元の dataframe の因子の情報(因子レベル)が含まれている
そこで、空の因子レベルを取り除かないと正しく検定できなくなってしまう
RWiki 因子Tips大全: http://www.okadajp.org/RWiki/?因子Tips大全

mydf2$Factor1 <- droplevels(mydf2$Factor1)

dataframe を表示させる:show()

21 March 2018

for 文の中では表示されないので、show() で対応する
http://d.hatena.ne.jp/yusap/20130514/1368486988 ヘのコメントを参照

show(
    mydf
)

dataframe の行数のカウント:nrow()

21 March 2018
x <- nrow(mydf)

データの扱いについてのメモ

11 July 2011
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 だけが出てくる。

検定を一気に行うなど

データのまとめ

データのまとめ

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
.
.
.

boxplot に使われる summary data の出力:グラフの描画ではない

21 March 2018

plot=FALSE がミソ

boxplot(mydf2$time~mydf2$cross, plot=FALSE)

ノンパラメトリック検定

Kruskal-Wallis 検定

21 March 2018

事後検定として Dunn 検定
プログラム引用:An R Companion for the Handbook of Biological Statistics, Kruskal-Wallis Test

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

2017年10月11日

Kruskal-Wallis test(ノンパラメトリック分散分析)と それに引き続く multiple comaprison(多重比較)を行う。

multiple comaprison は、 Zar (2010)*3 に 「ほぼ」*4 従って、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 に従って、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もそうしてみた)

2標本検定

2017年8月7日
 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 を使うのは 奥村さんの解説 に従ってはいるが、出力が
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検定

奥村さんの解説 を見ていただくとして。

論文

Kruskal-Wallis test

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
.
.
.

カテゴリ毎にサマリを出力(追記 22 February 2015)

> by(mydata2, mydata2$Age, summary)

パラメトリック検定

n元配置:分散分析、Tukey HSD 検定

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

Fisher test

03 Oct 2005
> mydata <- matrix(c(4,6,3,4), nr=2)
> fisher.test(mydata)

Fisher test for 2 by n

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を大きく取らないとデータによってはエラーになる場合がある

分割表の独立性の検定

21 March 2018

dataframe は 2x2 ではなくて、nx2 のデータが格納されているとき
まず全体(nx2)の分割表の独立性検定を行い、 続いて、pariwise comparisons を行う
Fisher Exact Test、カイ二乗検定、G検定で

カイ二乗検定では期待値を表示させて、 検定を行ってよい条件かどうかを見てわかるようにする
さらに、独立性の検定の期待値が適切かどうかを表示させる*5

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

カイ二乗検定

14 January 2016

適合性の検定

メンデルの法則などのようなものの検定

Zar (1981) Example 5.1 のデータ
null hypothesis (H0) は 3:1 (=メンデル分離するかどうかを検定)

参考にしたページ
An R Companion for the Handbook of Biological Statistics by Salvatore S. Mangiafico
G-test of Goodness-of-Fit

> 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検定を行いたいなら

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

頻度データの95%信頼限界の計算:binom.test()

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

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
.
.
.

相関

相関の計算・検定:cor.test()

21 March 2018
corresults <- cor.test(
   mydf3$CourtshipFreq,
   mydf3$Courtshipsec
)

相関を計算:cor()

21 March 2018

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

相関を個々に計算

21 March 2018

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)
    }
}

Comparing Regressions

異なる条件下で得られた複数の回帰直線を比較する

01 February 2008

長くなったのでこちら -> 複数の回帰直線の比較

グラフ

散布図を描く:plot()

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 まで、7分割=目盛を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)"
     )

棒グラフ

棒グラフ:barplot()

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()

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()

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()

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)

複数の折れ線グラフ:matplot()

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()

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 が別にあるとき

ダミーデータ追加

21 March 2018

図を描くとき「空き」を作るためなどに
場合によっていは、名前の文字列は空白がいいのかもしれない (でもひとつしか空きは作れないはず)

mydf2 <- rbind( mydf2, data.frame(cross="Dummy", time=NA))

グラフ領域の操作

グラフの軸ラベルの向き:las=

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=

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()

21 March 2018

substitute と paste のコンビネーション なので、適用できない場面も多々ある
グラフタイトルはおっけい

title(main=
      substitute(paste(
      bold('A. '),  bolditalic('D. melanogaster'), bold(' males'))
      ),
      cex.main=1.5, cex=1.5, )

グラフタイトルに変数を入れる:paste()

21 March 2018

paste を使って文字列を結合すれば、変数がきちんと展開される
italic などが使いたい場合は substitute() paste() list() を利用する

title(main=paste("r", " = ", coefficient, ", ", "P", " = ", pvalue))

グラフタイトルに変数を入れる:太字、斜体の対応も必要な場合 substitute() paste() list()

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
)

軸ラベルの調整:

24 March 2018

参考:http://uncorrelated.hatenablog.com/entry/20130708/1373256551

par(mgp=c(2.5, 1, 0)) 
# デフォルト par(mgp=c(3, 1, 0))

グラフ内に線分を描く:lines()

21 March 2018

lines(X座標のベクトル, Y座標のベクトル)
lines(点1,点2,点3)の形式ではないので注意

lines(c(3.5,5.4),c(95,95))

グラフ描画領域内に文字を書く:text()

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()

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 の操作

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"
                )
	)

描画領域の操作

描画前のおまじない

24 March 2018

ファイルへの書き込みのときはなくてもよい

plot.new()

描画領域の分割 split.screen(c())

24 March 2018

c(横, 縦)

split.screen(c(3,2)) # 3列2行

描画領域を更に分割 split.screen(c(*, *), screen = *)

24 March 2018

分割する画面を指定する

split.screen(c(1,2), screen = 1)
split.screen(c(2,2), screen = 2)

描画面の指定 screen()

24 March 2018

横向きに数が振られていく
再分割した場合は、更に数字が増える

screen(1)
# ココに 1 の描画内容を
screen(6)
# ココに 6 の描画内容を

描画画面の余白 par(oma=c(9,3,0,0))

24 March 2018

split.screen() と screen() の間に指定することで、外側に余白を作れる

split.screen(c(3,5))
par(oma=c(9,3,0,0)) # 余白を作る  下左上右
screen(1)

描画領域内のマージンの指定:下左上右 par(mar=c())

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)

分割画面の外側の操作

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= )

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()

24 March 2018
split.screen(c(3,5))
# 画面分割後の操作をここに
close.screen(all.screens=TRUE)

『「蜂群+平均値」に折れ線グラフを重ねる』をふたつ重ねる

2017年8月7日

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()

ESS

Emacs Speaks Statistics:Emacs で R を使う
Emacs には後からインストールをする必要あり
install-memo5.1を参照

ESS起動

03 Oct 2005
M-x R

コマンド入力後

C-c C-j
C-m (改行)

で、1行ずつコマンドを実行する。

コマンドメモ

30 March 2004
> 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

cygwin shell または、コマンドプロンプトで

rscript command.R
rscript command.R & ← バックグラウンドで実行する場合

※ pathを通しておかなければなりません(MS-DOS!のようですが)。

参考:R の入出力画面を経由せずに実行:他のプログラムからの呼び出し

Rのリンク

Today:5Yesterday:7Total:153392 since 04 April 2006

*1 そもそも、そんな微妙なデータでものを言うのは危ないので避けるべき。サンプルサイズを増やすか、データの取り方を考え直すかなどすべき
*2 コピー&ペーストのミスの可能性は高まるので注意は必要
*3 Zar, J. H. (2010) Biostatistical Analysis, 5th ed, Pearson Prentice-Hall, Upper Saddle River, NJ
*4 Zar には、P値の補正に Benjamini-Hochberg method を使うべきとは書いていない。実質的には使っているのだと思うけれども。
*5 Zar, J. H. (2010) Biostatistical Analysis, 5th ed, Pearson Prentice-Hall, Upper Saddle River, NJ

添付ファイル: filebeeswarmmeanplot.png 198件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 21 Jul 2019 (日) 10:25:02 (360d)