![[PukiWiki]](image/tomarufront.png "[PukiWiki]")
各$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 . . .
plot=FALSE がミソ
boxplot(mydf2$time~mydf2$cross, plot=FALSE)
事後検定として 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（多重比較）を行う。
multiple comaprison は、
Zar (2010)*1
に
「ほぼ」*2
従って、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もそうしてみた）
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」を指定
奥村さんの解説 を見ていただくとして。
論文
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)
なお，事前検定を行うことが不適切であることはだんだん理解されてきているので，この観点から言えば「等分散検定後に普通の 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)
デフォルトで等分散を仮定しない検定を行う
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 検定がよさそう*3
分散分析を行ってから、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
分割表の検定もいろいろあって、困る、かも……*4
mydata <- matrix(c(4,6,3,4), nr=2) fisher.test(mydata) fisher.test(mydata)$p.value
2 × 2 分割表の検定においては：
Lydersen et al. (2009) Recommended tests for association in 2×2 tables. Statist. Med. 28:1159--1175
は現時点では、unconditional test か mid-p を推奨
Ruxton & Neuhauser (2010) Good practice in testing for an association in contingency tables. Behav. Ecol. Sociobiol 64:1505--1513
は、現時点（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() と同じ結果を出力する
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
> 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を大きく取らないとデータによってはエラーになる場合がある
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")
適合性の検定
メンデルの法則などのようなものの検定
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 = 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
確率は $p.value に格納される
mydata <- c(34, 27) binom.test(mydata) binom.test(mydata)$p.value
二項分布を利用
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"
)
}
ショウジョウバエの配偶行動を見るときに、観察時間内に交尾しないものがある。そこで、交尾したものを「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 . . .
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)
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
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)
}
}
長くなったのでこちら -> 複数の回帰直線の比較
| Today:1 | Yesterday:0 | Total:3231 since 21 July 2019 |
