#======================================================
#	3 統計的推測
#======================================================
#======================================================
#	3.1 検定に関連する主要な分布
#======================================================

#------------<正規分布>------------

#図3.1
x <- seq(-4, 6, length=200)
plot(x, dnorm(x), type="l")
#既存のプロットに線を書き加えるには、linesを使用する。
#線種は、ltyで指定する。
lines(x, dnorm(x, mean=2), lty=2) 			#平均を２とする。
lines(x, dnorm(x, mean=2, sd=1.5), lty=3, lwd=1.5)#標準偏差を1.5とする。
#凡例を書き加える。
legend(x=-4, y=0.4, lty=1:3, lwd=c(1, 1, 1.5),
	legend=c("N(0, 1)", "N(2, 1)", "N(2, 2.25)"))

#------------<カイ二乗分布>------------

#ガンマ関数
factorial(5) #5!
gamma(5+1) 　#同じ

#カイ二乗分布
#図3.2
x <- 0:30
plot(x, dchisq(x, df=5), type="l", ylab="dchisq")
lines(x, dchisq(x, df=10), lty=2)
lines(x, dchisq(x, df=20), lty=3, lwd=1.5)
#凡例に数式を表示する。
legend(x=25, y=0.15, lty=1:3, lwd=c(1, 1, 1.5),
	legend=c(expression(chi[5]^2), expression(chi[10]^2),
	expression(chi[20]^2)))

#数式
demo(plotmath)

#------------<t分布>------------

#図3.3
x <- seq(-5, 5, length=200)
plot(x, dt(x, 1), type="l", ylim=c(0, 0.45), ylab="dt")
lines(x, dt(x, 5), lty=2)
lines(x, dt(x, 20), lty=3, lwd=1.5)
legend(x=-4, y=0.44, lty=1:3, lwd=c(1, 1, 1.5), 
	legend=c(expression(t[1]), expression(t[5]),
	expression(t[20])))
lines(x, dt(x, 1063), lty=3, lwd=3, col="red")
abline(h=0)

#------------<F分布>------------

#図3.4
x <- seq(0, 10, length=200)
plot(x, df(x, df1=1, df2=10), type="l", ylab="df", ylim=c(0, 1))
lines(x, df(x, df1=10, df2=10), lty=2)
legend(x=8, y=1.0, lty=1:3,
	legend=c(expression(F["1,10"]), expression(F["10, 10"])))

#======================================================
#	3.2 １標本問題
#======================================================
#------------<データセットの整備>------------
肝臓病 <- read.csv("C:/r/users/datasets/ILPD/ILPD.csv", header=F)
head(肝臓病)　#確認
#変数名を付ける。
names(肝臓病) <- c("年齢", "性別", "総ビリルビン", "直接ビリルビン",
                  "ALP", "GPT_ALT", "GOT_AST", "TP", "ALB", 
                  "A/G比","肝臓病")　
table(肝臓病$肝臓病) #確認
肝臓病$肝臓病[肝臓病$肝臓病==2] <- 0
肝臓病$肝臓病 <- as.factor(肝臓病$肝臓病)
summary(肝臓病)      #要約

#A/G比の計算
肝臓病$"A/G比" <- 肝臓病$ALB/(肝臓病$TP-肝臓病$ALB)

#データの部分集合
肝臓病男性 <- 肝臓病[肝臓病$肝臓病==1 & 肝臓病$性別=="Male", ]
正常者男性 <- 肝臓病[肝臓病$肝臓病==0 & 肝臓病$性別=="Male", ]
肝臓病女性 <- 肝臓病[肝臓病$肝臓病==1 & 肝臓病$性別=="Female", ]
正常者女性 <- 肝臓病[肝臓病$肝臓病==0 & 肝臓病$性別=="Female", ]

#------------<QQプロット>------------

#図3.5
qqnorm(肝臓病男性$"A/G比", main="A/G比：肝臓病（男性）")
qqline(肝臓病男性$"A/G比")

#図3.6
qqnorm(肝臓病女性$"A/G比", main="A/G比：肝臓病（女性）")
qqline(肝臓病女性$"A/G比")

#------------<Kolmogorov-Smirnov検定>------------

#図3.7
#経験累積分布関数 A/G比（男性）
x <- ecdf(肝臓病男性$"A/G比")        
plot(x, xlab="A/G比", col="steelblue", 
	ylab="累積分布関数", main="経験分布と理論分布（A/G比：男性患者）")
#理論累積分布関数 A/G比（男性）
xx <- seq(0, 2, length=200)
mu <- mean(肝臓病男性$"A/G比")
sigma <- sd(肝臓病男性$"A/G比")
lines(xx, pnorm(xx, mean=mu, sd=sigma), lty=2, col="red", lwd=2)
arrows(x0=1.3, x1=1.05, y0=0.65, y1=0.65, length=0.15, code=2)   #矢印	
text(1.35, 0.65, "理論分布", adj=0)                              #テキスト

#図3.8
#経験累積分布関数A/G比（女性）
x <- ecdf(肝臓病女性$"A/G比")        
plot(x, xlab="A/G比", col="steelblue", 
	ylab="累積分布関数", main="経験分布と理論分布（A/G比：女性患者）")
#理論累積分布関数 A/G比（女性）
xx <- seq(0, 3, length=200)
mu <- mean(肝臓病女性$"A/G比")
sigma <- sd(肝臓病女性$"A/G比")
lines(xx, pnorm(xx, mean=mu, sd=sigma), lty=2, col="red", lwd=2)
arrows(x0=1.35, x1=1.05, y0=0.65, y1=0.65, length=0.15, code=2)
text(1.4, 0.65, "理論分布", adj=0)

#------------<Kolmogorov-Smirnov検定>------------

mu <- mean(肝臓病男性$"A/G比")
sigma <- sd(肝臓病男性$"A/G比")
ks.test(肝臓病男性$"A/G比", "pnorm", mean=mu, sd=sigma)
set.seed(123)
ks.test(jitter(肝臓病男性$"A/G比"), "pnorm", mean=mu, sd=sigma)

mu <- mean(肝臓病女性$"A/G比")
sigma <- sd(肝臓病女性$"A/G比")
ks.test(jitter(肝臓病女性$"A/G比"), "pnorm", mean=mu, sd=sigma)

#------------<t検定>------------

mu.0 <- mean(正常者男性$"A/G比")    #男性正常者のA/G比平均
#男性患者のA/G比の平均が、正常者のA/G比平均と異なるかどうかのt検定
h <- t.test(肝臓病男性$"A/G比", mu=mu.0)
h

#検算
mu.a <- mean(肝臓病男性$"A/G比")        #標本平均
sigma.2 <- var(肝臓病男性$"A/G比")      #標本分散
n <- length(肝臓病男性$"A/G比")         #標本数
t.value <- (mu.a-mu.0)/sqrt(sigma.2/n)  #t値
t.value
n - 1                            #自由度（df）
pt(-abs(t.value), df=n-1) * 2    #p値

#女性
mu.0 <- mean(正常者女性$"A/G比")
(h2 <- t.test(肝臓病女性$"A/G比", mu=mu.0))

#戻り値	
names(h) 				#リストhの各成分名
(male.ci <- h$conf.int)           #男性患者の95％信頼区間
(female.ci <- h2$conf.int)        #女性患者の95％信頼区間

#図3.9
op <- par(mfrow=c(2, 1), mar=par()$mar - 1)
plot(c(0.85, 1.1), c(0, 2), type="n",      #type="n"：何も描かない
	main="男性患者のA/G比平均の信頼区間", 
	xlab="A/G比", ylab="", axes=F)      #y軸のラベルも、xy軸も描かない
axis(side=1)                               #x軸のみを描く
arrows(x0=male.ci[1], x1=male.ci[2], y0=1, y1=1,
	length=0.15, code=3, col="steelblue", lwd=3) #男性患者の信頼区間の矢印
points(h$est, 1, pch=16, col="steelblue", cex=2)    #男性患者の標本平均
mu.0 <- mean(正常者男性$"A/G比")
abline(v=mu.0, lty=2, col="darkgreen")              #男性正常者平均の垂直線
text(mu.0+0.008, 1, "男性正常者平均", srt=-90)      #正常者の注記。srtは文字列回転

plot(c(0.85, 1.1), c(0, 2), type="n",
	main="女性患者のA/G比平均の信頼区間", 
	xlab="A/G比", ylab="", axes=F)
axis(side=1)
arrows(x0=female.ci[1], x1=female.ci[2], y0=1, y1=1,
	length=0.15, code=3, col="steelblue", lwd=3)  #女性患者の信頼区間の矢印
points(h2$est, 1, pch=16, col="steelblue", cex=2)    #女性患者の標本平均
mu.0 <- mean(正常者女性$"A/G比")
abline(v=mu.0, lty=2, col="darkgreen")
text(mu.0+0.008, 1, "女性正常者平均", srt=-90)
par(op)

#信頼区間の検算
mu.a + sqrt(sigma.2/n) * qt(c(0.025, 0.975), df=n-1)

#------------<Wilcoxon符号順位検定>------------

set.seed(123)     #乱数を再現可能とするため
(x <- rnorm(20, mean=50, sd=10))    #正規乱数を20個作成　　 　
(rnk <- rank(abs(x-50)))   #差の絶対値が最も小さい値からの昇順の順位
(V <- sum(rnk[(x-50) > 0]))  #差の符号が正である観測値の順位の合計

wilcox.test(x, mu=50)

#A/G比
mu.0 <- mean(正常者男性$"A/G比")
wilcox.test(肝臓病男性$"A/G比", mu=mu.0, conf.int=T)

mu.0 <- mean(正常者女性$"A/G比")
wilcox.test(肝臓病女性$"A/G比", mu=mu.0, conf.int=T) 

#------------<比率に関する検定>------------

#男性
n.0 <- length(肝臓病$年齢[肝臓病$性別=="Male"])   #男性観測数
r.0 <- length(肝臓病$肝臓病[肝臓病$性別=="Male"
		& 肝臓病$肝臓病 == 1])             #男性肝臓病患者数　　　　　　　　　　　　　 
(prop.0 <- r.0/n.0)    #男性肝臓病患者比率

#女性
n.a <- length(肝臓病$年齢[肝臓病$性別=="Female"])   #女性観測数
r.a <- length(肝臓病$肝臓病[肝臓病$性別=="Female"
		& 肝臓病$肝臓病 == 1])             #女性肝臓病患者数　　　　　　　　　　　　　 
(prop.a <- r.a/n.a)    #女性肝臓病患者比率

#検定
h <- binom.test(r.a, n.a, p=prop.0)
h

#検算（p値：観測比率が理論値以下の場合）
relErr <- 1 + 1e-07
d <- dbinom(r.a, n.a, p=prop.0)
m <- n.a * prop.0
i <- seq.int(from=ceiling(m), to=n.a)
y <- sum(dbinom(i, n.a, prop.0) <= d * relErr)
pbinom(r.a, n.a, prop.0) + pbinom(n.a - y, n.a, prop.0, lower=F)

#検算（信頼区間）
女性.ci <- numeric(2) #初期化
(女性.ci[1] <- qbeta(0.025, r.a, n.a - r.a + 1)) #下限
(女性.ci[2] <- qbeta(0.975, r.a + 1, n.a - r.a)) #上限

#逆転バージョン
h2 <- binom.test(r.0, n.0, p=prop.a)
h2

#検算（p値：観測比率が理論値よりも大きい場合）
relErr <- 1 + 1e-07
d <- dbinom(r.0, n.0, p=prop.a)
m <- n.0 * prop.a
i <- seq.int(from=0, to=floor(m))
y <- sum(dbinom(i, n.0, prop.a) <= d * relErr)
pbinom(y - 1, n.0, prop.a) + pbinom(r.0 - 1, n.0, prop.a, lower=F)

#検算（信頼区間）
男性.ci <- numeric(2) #初期化
(男性.ci[1] <- qbeta(0.025, r.0, n.0 - r.0 + 1)) #下限
(男性.ci[2] <- qbeta(0.975, r.0 + 1, n.0 - r.0)) #上限

#図3.10
op <- par(mfrow=c(2, 1), mar=par()$mar - 1)
plot(c(0.55, 0.8), c(0, 2), type="n",      #type="n"：何も描かない
	main="女性の患者比率の信頼区間", 
	xlab="患者比率", ylab="", axes=F)   #y軸のラベルも、xy軸も描かない
axis(side=1)                               #x軸のみを描く
arrows(x0=女性.ci[1], x1=女性.ci[2], y0=1, y1=1,
	length=0.15, code=3, col="steelblue", lwd=3) #女性の信頼区間の矢印
points(h$est, 1, pch=16, col="steelblue", cex=2)    #女性患者比率の標本平均
abline(v=prop.0, lty=2, col="darkgreen")            #男性患者比率の垂直線
text(prop.0+0.008, 1, "男性患者比率", srt=-90)      #注記。srtは文字列回転

plot(c(0.55, 0.8), c(0, 2), type="n",
	main="男性の患者比率の信頼区間", 
	xlab="患者比率", ylab="", axes=F)
axis(side=1)
arrows(x0=男性.ci[1], x1=男性.ci[2], y0=1, y1=1,
	length=0.15, code=3, col="steelblue", lwd=3)  #男性患者の信頼区間の矢印
points(h2$est, 1, pch=16, col="steelblue", cex=2)    #男性患者の標本平均
abline(v=prop.a, lty=2, col="darkgreen")
text(prop.a+0.008, 1, "女性患者比率", srt=-90)
par(op)

#======================================================
#	3.3 ２標本問題
#======================================================

#------------<等分散に関する検定：肝臓病>------------

#男性
(h <- var.test(肝臓病男性$ALP, 正常者男性$ALP))

#検算
var.a <- var(肝臓病男性$ALP)    #肝臓病男性の分散
var.b <- var(正常者男性$ALP)    #正常者男性の分散
(f.value <- var.a/var.b)            #F値
(m <- length(肝臓病男性$ALP))    #肝臓病男性の件数
(n <- length(正常者男性$ALP))    #正常者男性の件数
p.value <- pf(f.value, m-1, n-1)
2 * min(p.value, 1 - p.value)            #両側検定のため２倍する。
(male.ci <- f.value/qf(c(0.975, 0.025), m-1, n-1)) #95％信頼区間

#女性
(h2 <- var.test(肝臓病女性$ALP, 正常者女性$ALP))

summary(肝臓病$ALP)
head(rev(sort(肝臓病$ALP)), 20)

#------------<等分散に関する検定：自動車>------------

data(car.test.frame, package="rpart")
summary(car.test.frame)

#燃費
日本車.燃費 <- car.test.frame$Mileage[car.test.frame$Country=="Japan"]
それ以外.燃費 <- car.test.frame$Mileage[car.test.frame$Country!="Japan"]
(h <- var.test(日本車.燃費, それ以外.燃費))

#価格
日本車.価格 <- car.test.frame$Price[car.test.frame$Country=="Japan"]
それ以外.価格 <- car.test.frame$Price[car.test.frame$Country!="Japan"]
(h2 <- var.test(日本車.価格, それ以外.価格))

#図3.11
plot(c(0, 7), c(0, 3), type="n", 
	main="燃費と価格の分散比の95％信頼区間", 
	xlab="分散比", ylab="", axes=F)
axis(side=1)
arrows(x0=h$conf.int[1], x1=h$conf.int[2], y0=2, y1=2,
	length=0.15, code=3, col="steelblue", lwd=3)	#燃費の95％信頼区間
arrows(x0=h2$conf.int[1], x1=h2$conf.int[2], y0=1, y1=1,
	length=0.15, code=3, col="tomato4", lwd=3)	#価格の95％信頼区間
points(h$est, 2, pch=16, col="steelblue", cex=2)	#燃費の標本分散
points(h2$est, 1, pch=16, col="tomato4", cex=2)		#価格の標本分散
abline(v=1, lty=2, col="darkgreen")			#等分散の垂直線
text(0.3, 2, "燃費", cex=1.5)
text(0.45, 1, "価格", cex=1.5)
text(1+0.3, 2.5, "等分散", srt=-90, cex=1.5)		#等分散の注記

#------------<２標本検定>------------

#等分散
(h <- t.test(日本車.燃費, それ以外.燃費, var.equal=T))

#検算
var.a <- var(日本車.燃費)     #日本車.燃費の分散
var.b <- var(それ以外.燃費)   #それ以外.燃費の分散
n.a <- length(日本車.燃費)    #日本車.燃費の件数
n.b <- length(それ以外.燃費)  #それ以外.燃費の件数
V <- ((n.a-1)*var.a + (n.b-1)*var.b)/((n.a-1) + (n.b-1)) #併合分散
(mu.a <- mean(日本車.燃費))   #日本車.燃費の平均
(mu.b <- mean(それ以外.燃費)) #それ以外.燃費の平均
mileage.diff <- mu.a - mu.b
(t.value <- (mu.a - mu.b)/sqrt(V*(1/n.a + 1/n.b))) #t値
(df.all <- (n.a-1) + (n.b-1))      #自由度
pt(-abs(t.value), df=df.all) * 2   #p値。両側検定のため２倍する。

#信頼区間（参考までに）
a <- 日本車.燃費
b <- それ以外.燃費
ci.s <- sqrt(((1/n.a + 1/n.b) * 
	(sum((a - mu.a)^2) + sum((b - mu.b)^2)))/df.all)
(mileage.ci <- (mu.a - mu.b) + qt(c(0.025, 0.975), df.all) * ci.s)

#Welchの方法
(h2 <- t.test(日本車.価格, それ以外.価格))

#検算（参考までに）
var.a <- var(日本車.価格)    #日本車.価格の分散
var.b <- var(それ以外.価格)  #それ以外.価格の分散
n.a <- length(日本車.価格)   #日本車.価格の件数
n.b <- length(それ以外.価格) #それ以外.価格の件数
V <- var.a/n.a + var.b/n.b   #併合分散
(mu.a <- mean(日本車.価格))  #日本車.価格の平均
(mu.b <- mean(それ以外.価格))#それ以外.価格の平均
price.diff <- mu.a - mu.b
(t.value <- (mu.a - mu.b)/sqrt(V)) #t値
const <- 1/(1+(n.a*var.b)/(n.b*var.a))
(df.all <- 1/((const^2)/(n.a-1)+((1-const)^2)/(n.b-1))) #自由度
pt(-abs(t.value), df=df.all) * 2   #p値。両側検定のため２倍する。

#信頼区間（参考までに）
ci.s <- sqrt(V)
(price.ci <- (mu.a - mu.b) + qt(c(0.025, 0.975), df.all) * ci.s)

#------------<対応のある２標本検定>------------

#データ
wear.A <- c(14.0, 8.8, 11.2, 14.2, 11.8,
        6.4, 9.8, 11.3, 9.3, 13.6)
wear.B <- c(13.2, 8.2, 10.9, 14.3, 10.7, 
        6.6, 9.5, 10.8, 8.8, 13.3)

#通常のt検定
t.test(wear.A, wear.B)

#対応のあるt検定
t.test(wear.A, wear.B, paired=T)

#------------<Wilcoxon順位和検定>------------

wilcox.test(肝臓病男性$ALP, 正常者男性$ALP, conf.int=T)

#------------<Fisher正確検定>------------

#肝臓病
n.a <- length(肝臓病$年齢[肝臓病$性別=="Male"])   #男性観測数
n.b <- length(肝臓病$年齢[肝臓病$性別=="Female"]) #女性観測数
s.a <- nrow(肝臓病男性)  #肝臓病男性の数
s.b <- nrow(肝臓病女性)  #肝臓病女性の数
f.a <- n.a - s.a  #正常者男性の数
f.b <- n.b - s.b  #正常者女性の数
mat.肝臓病 <- matrix(c(s.a, f.a, s.b, f.b), nrow=2, byrow=T)
fisher.test(mat.肝臓病)

#------------<prop.test>------------

prop.test(mat.肝臓病)

#======================================================
#	3.4 適合度検定
#======================================================

#------------<２標本Kolmogorov-Smirnov検定>------------

set.seed(123)
ks.test(jitter(肝臓病男性$ALP), jitter(肝臓病女性$ALP))
ks.test(jitter(肝臓病男性$"A/G比"), jitter(肝臓病女性$"A/G比"))

#------------<独立性に関するカイ二乗検定>------------

quantile(肝臓病$ALB)
bks1 <- c(0, 2.5, 3.0, 3.5, 6)   #ALB区分
quantile(肝臓病$"A/G比")
bks2 <- c(0, 0.75, 1, 1.15, 3.1) #A/G比区分

mat <- table(cut(肝臓病$ALB, breaks=bks1),
        cut(肝臓病$"A/G比", breaks=bks2))  #分割表
mat
chisq.test(mat) #カイ二乗検定