R_statistics
をテンプレートにして作成
開始行:
#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 ...
*** boxplot に使われる summary data の出力:グラフの描画...
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 Biol...
kruskal.test(Courtshipsec ~ EDnum, data = mydf2)
#
library(FSA)
mydt <- dunnTest(Courtshipsec ~ EDnum, data = mydf2, met...
#
cat("\nDunn test\n")
# なぜか、この表示が sink()するとファイルに記録されなか...
cat("Dunn (1964) Kruskal-Wallis multiple comparison\n")
cat(" p-values adjusted with the Benjamini-Hochberg met...
#
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, 5t...
に
「ほぼ」((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 Statistic...
に従って、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$」方式でないと動かないようなので、boxp...
*** 2標本検定 [#f7575e3c]
RIGHT:2017年8月7日
- Wilcoxon-Mann-Whitney test または、Exact Wilcoxon rank ...
- 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("...
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 permut...
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/sta...
に従ってはいるが、出力が~
exactRankTests は、W値~
coin を使った場合は、Z値~
となる。
古くにU検定を教わった身には、W値のほうが、正規化されたZ値...
それに現時点では開発停止と書いてあるけれども、package `ex...
2017-03-01版がある。~
参照:https://cran.r-project.org/web/packages/exactRankTe...
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/sta...
を見ていただくとして。
論文
- Brunner, E., & Munzel, U. (2000). The Nonparametric Beh...
- 名取真人. (2014). マン・ホイットニーのU検定と不等分散時...
*** Kruskal-Wallis test [#acb5ec04]
RIGHT:09 May 2004
3つ以上の$Ageの$Durationについて検定。Kruskal-Wallis検定...
> 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 ...
カテゴリ毎にサマリを出力(追記 22 February 2015)
> by(mydata2, mydata2$Age, summary)
** パラメトリック検定 [#fbd6e1ba]
*** 等分散の検定:Bartlett 検定など [#x1061cc4]
RIGHT:20 July 2019
> なお,事前検定を行うことが不適切であることはだんだん理...
https://oku.edu.mie-u.ac.jp/~okumura/stat/ttest.html
ではあるのだけれど、検定をする必要もあるかもしれない。
また、http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc043...
も参考に。
参考: http://www.cookbook-r.com/Statistical_analysis/Homo...
# 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...
oneway.test(Courtshipsec ~ EDnum, data=mydf2)
更に……~
文句を言われないなら、分散分析はせずに Tukey HSD 検定がよ...
*** n元配置:分散分析、Tukey HSD 検定 [#s26a9cb0]
RIGHT:21 March 2018
分散分析を行ってから、Tukey HSD 検定~
最近は、Tukey HSD 検定だけでもよいみたいですね~
プログラム引用: https://rcompanion.org/rcompanion/d_05.h...
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 ...
** 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 × 2 分割表の検定においては:~
[[Lydersen et al. (2009) Recommended tests for associatio...
は現時点では、unconditional test か mid-p を推奨~
[[Ruxton & Neuhauser (2010) Good practice in testing for ...
は、現時点(2010)では、Yates 補正なしカイ自乗検定か mid-...
上の2つの論文とも、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 は一方(列)のみ固定~
2つの比を比較する検定となる。
Fisher exact test はオッズ比が1かどうかの検定。~
データの入力に際して縦横に注意が必要
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...
> mydata
> fisher.test(mydata, workspace=20000000)
workspaceを大きく取らないとデータによってはエラーになる場...
*** 分割表の独立性の検定 [#nbaf3812]
RIGHT:21 March 2018
dataframe は 2x2 ではなくて、nx2 のデータが格納されている...
まず全体(nx2)の分割表の独立性検定を行い、
続いて、pariwise comparisons を行う~
Fisher Exact Test、カイ二乗検定、G検定で
カイ二乗検定では期待値を表示させて、
検定を行ってよい条件かどうかを見てわかるようにする~
さらに、独立性の検定の期待値が適切かどうかを表示させる((Z...
- 期待値は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) ...
#
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[expected...
}
#
# 独立性の検定の期待値が適切かどうかを表示させる
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 w...
#
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 Willi...
#
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" "wi...
}
#
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/rcompan...
> 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 では、correctio...
オプションを付ければ Yates correction を行うが、
適合性の検定では行わない。~
そこで、手計算と同じようにせざるを得ない。~
差分の絶対値から0.5を引きそれを二乗するように式を変えれば...
> 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$parame...
"\t",
100*mybntest$conf.int[1],
"\t",
100*mybntest$conf.int[2],
"\t"
)
}
** survival analysis [#e1a975ad]
RIGHT:11 April 2004
ショウジョウバエの配偶行動を見るときに、観察時間内に交尾...
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 ...
** 相関 [#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 に3つ以上のデータが含まれている場合
# mycor <- cor(mydf3) # use指定をしないと NA データがあ...
# mycor <- cor(mydf3, use="complete.obs") # これは全体の...
mycor <- cor(mydf3, use="pairwise.complete.obs") # これ...
mycod <- mycor^2 # coeffient of determination
*** 相関を個々に計算 [#i5201e46]
RIGHT:21 March 2018
dataframe に3つ以上のデータが含まれている場合
for(i in 1:7) {
for(j in 1:7) {
mycortest <- cor.test(mydf3[, i], mydf3[, j], me...
cat("data i: ", names(mydf3)[i], "\ndata j: ", n...
print(mycortest)
}
}
#
for(i in 1:7) {
for(j in 1:7) {
mycortest <- cor.test(mydf3[, i], mydf3[, j], me...
cat("data i: ", names(mydf3)[i], "\ndata j: ", n...
print(mycortest)
}
}
#
for(i in 1:7) {
for(j in 1:7) {
mycortest <- cor.test(mydf3[, i], mydf3[, j], me...
cat("data i: ", names(mydf3)[i], "\ndata j: ", n...
print(mycortest)
}
}
** Comparing Regressions [#k65cfd16]
*** 異なる条件下で得られた複数の回帰直線を比較する [#xae7...
RIGHT:01 February 2008
-傾きに差はある?
-切片に差はある?
長くなったのでこちら -> [[複数の回帰直線の比較>R_Comparin...
|Today:&counter(today); |Yesterday:&counter(yesterday); |...
終了行:
#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 ...
*** boxplot に使われる summary data の出力:グラフの描画...
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 Biol...
kruskal.test(Courtshipsec ~ EDnum, data = mydf2)
#
library(FSA)
mydt <- dunnTest(Courtshipsec ~ EDnum, data = mydf2, met...
#
cat("\nDunn test\n")
# なぜか、この表示が sink()するとファイルに記録されなか...
cat("Dunn (1964) Kruskal-Wallis multiple comparison\n")
cat(" p-values adjusted with the Benjamini-Hochberg met...
#
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, 5t...
に
「ほぼ」((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 Statistic...
に従って、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$」方式でないと動かないようなので、boxp...
*** 2標本検定 [#f7575e3c]
RIGHT:2017年8月7日
- Wilcoxon-Mann-Whitney test または、Exact Wilcoxon rank ...
- 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("...
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 permut...
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/sta...
に従ってはいるが、出力が~
exactRankTests は、W値~
coin を使った場合は、Z値~
となる。
古くにU検定を教わった身には、W値のほうが、正規化されたZ値...
それに現時点では開発停止と書いてあるけれども、package `ex...
2017-03-01版がある。~
参照:https://cran.r-project.org/web/packages/exactRankTe...
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/sta...
を見ていただくとして。
論文
- Brunner, E., & Munzel, U. (2000). The Nonparametric Beh...
- 名取真人. (2014). マン・ホイットニーのU検定と不等分散時...
*** Kruskal-Wallis test [#acb5ec04]
RIGHT:09 May 2004
3つ以上の$Ageの$Durationについて検定。Kruskal-Wallis検定...
> 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 ...
カテゴリ毎にサマリを出力(追記 22 February 2015)
> by(mydata2, mydata2$Age, summary)
** パラメトリック検定 [#fbd6e1ba]
*** 等分散の検定:Bartlett 検定など [#x1061cc4]
RIGHT:20 July 2019
> なお,事前検定を行うことが不適切であることはだんだん理...
https://oku.edu.mie-u.ac.jp/~okumura/stat/ttest.html
ではあるのだけれど、検定をする必要もあるかもしれない。
また、http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc043...
も参考に。
参考: http://www.cookbook-r.com/Statistical_analysis/Homo...
# 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...
oneway.test(Courtshipsec ~ EDnum, data=mydf2)
更に……~
文句を言われないなら、分散分析はせずに Tukey HSD 検定がよ...
*** n元配置:分散分析、Tukey HSD 検定 [#s26a9cb0]
RIGHT:21 March 2018
分散分析を行ってから、Tukey HSD 検定~
最近は、Tukey HSD 検定だけでもよいみたいですね~
プログラム引用: https://rcompanion.org/rcompanion/d_05.h...
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 ...
** 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 × 2 分割表の検定においては:~
[[Lydersen et al. (2009) Recommended tests for associatio...
は現時点では、unconditional test か mid-p を推奨~
[[Ruxton & Neuhauser (2010) Good practice in testing for ...
は、現時点(2010)では、Yates 補正なしカイ自乗検定か mid-...
上の2つの論文とも、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 は一方(列)のみ固定~
2つの比を比較する検定となる。
Fisher exact test はオッズ比が1かどうかの検定。~
データの入力に際して縦横に注意が必要
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...
> mydata
> fisher.test(mydata, workspace=20000000)
workspaceを大きく取らないとデータによってはエラーになる場...
*** 分割表の独立性の検定 [#nbaf3812]
RIGHT:21 March 2018
dataframe は 2x2 ではなくて、nx2 のデータが格納されている...
まず全体(nx2)の分割表の独立性検定を行い、
続いて、pariwise comparisons を行う~
Fisher Exact Test、カイ二乗検定、G検定で
カイ二乗検定では期待値を表示させて、
検定を行ってよい条件かどうかを見てわかるようにする~
さらに、独立性の検定の期待値が適切かどうかを表示させる((Z...
- 期待値は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) ...
#
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[expected...
}
#
# 独立性の検定の期待値が適切かどうかを表示させる
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 w...
#
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 Willi...
#
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" "wi...
}
#
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/rcompan...
> 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 では、correctio...
オプションを付ければ Yates correction を行うが、
適合性の検定では行わない。~
そこで、手計算と同じようにせざるを得ない。~
差分の絶対値から0.5を引きそれを二乗するように式を変えれば...
> 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$parame...
"\t",
100*mybntest$conf.int[1],
"\t",
100*mybntest$conf.int[2],
"\t"
)
}
** survival analysis [#e1a975ad]
RIGHT:11 April 2004
ショウジョウバエの配偶行動を見るときに、観察時間内に交尾...
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 ...
** 相関 [#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 に3つ以上のデータが含まれている場合
# mycor <- cor(mydf3) # use指定をしないと NA データがあ...
# mycor <- cor(mydf3, use="complete.obs") # これは全体の...
mycor <- cor(mydf3, use="pairwise.complete.obs") # これ...
mycod <- mycor^2 # coeffient of determination
*** 相関を個々に計算 [#i5201e46]
RIGHT:21 March 2018
dataframe に3つ以上のデータが含まれている場合
for(i in 1:7) {
for(j in 1:7) {
mycortest <- cor.test(mydf3[, i], mydf3[, j], me...
cat("data i: ", names(mydf3)[i], "\ndata j: ", n...
print(mycortest)
}
}
#
for(i in 1:7) {
for(j in 1:7) {
mycortest <- cor.test(mydf3[, i], mydf3[, j], me...
cat("data i: ", names(mydf3)[i], "\ndata j: ", n...
print(mycortest)
}
}
#
for(i in 1:7) {
for(j in 1:7) {
mycortest <- cor.test(mydf3[, i], mydf3[, j], me...
cat("data i: ", names(mydf3)[i], "\ndata j: ", n...
print(mycortest)
}
}
** Comparing Regressions [#k65cfd16]
*** 異なる条件下で得られた複数の回帰直線を比較する [#xae7...
RIGHT:01 February 2008
-傾きに差はある?
-切片に差はある?
長くなったのでこちら -> [[複数の回帰直線の比較>R_Comparin...
|Today:&counter(today); |Yesterday:&counter(yesterday); |...
ページ名: