検定を一気に行うなど

データのまとめ

データのまとめ

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)*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もそうしてみた)

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)

パラメトリック検定

等分散の検定:Bartlett 検定など

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

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 検定がよさそう*3

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

分割表の検定もいろいろあって、困る、かも……*4

Fisher test

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

20 July 2019

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 を推奨
上の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

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

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

二項検定:binom.test()

20 July 2019

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

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

頻度データの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 & 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()

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

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

Today:1Yesterday:1Total:165 since 21 July 2019

*1 Zar, J. H. (2010) Biostatistical Analysis, 5th ed, Pearson Prentice-Hall, Upper Saddle River, NJ
*2 Zar には、P値の補正に Benjamini-Hochberg method を使うべきとは書いていない。実質的には使っているのだと思うけれども。
*3 Zar, J. H. (2010) Biostatistical Analysis, 5th ed, Pearson Prentice-Hall, Upper Saddle River, NJ
*4 ワタシは困っているのかも
*5 Zar, J. H. (2010) Biostatistical Analysis, 5th ed, Pearson Prentice-Hall, Upper Saddle River, NJ

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