setwd("C:/Documents/StatisticalReform/data")

##### \subsection *{3.7.1 独立な2 群の平均値の比較}
x1 <- c(59 ,48 ,51 ,41 ,39 ,84 ,95 ,56 ,86 ,74) #実験群
x2 <- c(47 ,24 ,38 ,28 ,39 ,74 ,77 ,48 ,40 ,60) #統制群

# 要約統計量
m1 <- mean(x1)
m2 <- mean(x2)
n1 <- length(x1)
n2 <- length(x2)
s1 <- sd(x1)
s2 <- sd(x2)

# Cohenのd
Sp <- sqrt(((n1-1)*s1^2+(n2-1)*s2^2)/(n1+n2))
d <- (m1-m2)/Sp
d

# Hedgesのg
sp <- sqrt(((n1-1)*s1^2+(n2-1)*s2^2)/(n1+n2-2))
g <- (m1-m2)/sp
g

# GlassのDelta
Delta <- (m1-m2)/s2
Delta

# ノンパラメトリックな効果量
library(orddom) # orddomパッケージを利用
xx1<-t(matrix(x1,1))  # 関数orddom()は引数が行列である必要があるので、変換  
colnames(xx1)<-c("Group 1")
xx2<-t(matrix(x2,1))
colnames(xx2)<-c("Group 2")

orddom(xx1,xx2)



##### \subsection *{3.7.2 相関のある2群の平均値の比較}

xd <- x1 - x2
md <- mean(xd)
sd <- sd(xd)

dd <- md/sd
dd


##### \subsection *{3.7.3 1要因被験者間の分散分析}

data <- read.csv("illusionN.csv")
data1 <- data[data$b=="0度",] #これが表3.8のデータ

res <- aov(y~a,data=data1)
sres <- summary(res) 
sres # 分散分析表(表3.9)

SSa <- sres[[1]]$"Sum Sq"[1]
SSe <- sres[[1]]$"Sum Sq"[2]
SSt <- SSa + SSe
eta2 <- SSa/SSt # η^2
eta2

MSa <- sres[[1]]$"Mean Sq"[1]
MSe <- sres[[1]]$"Mean Sq"[2]
dfa <- sres[[1]]$Df[1]
omega2 <- dfa * (MSa - MSe) / (SSt + MSe) # ω^2
omega2


##### \subsection *{3.7.4 1要因被験者内の分散分析}
 
res <- aov(y~a+Error(s/a),data=data1)
sres <- summary(res) 
sres # 分散分析表(表3.9)

SSa <- sres[[2]][[1]]$"Sum Sq"[1]
SSe <- sres[[2]][[1]]$"Sum Sq"[2]
SSs <- sres[[1]][[1]]$"Sum Sq"[1]

SSt <- SSa + SSe + SSs
eta2 <- SSa/SSt # η^2 (3.61)式
eta2p <- SSa/(SSa+SSe) # η^2_p (3.62)式
eta2; eta2p 

dfa <- sres[[2]][[1]]$"Df"[1]
MSa <- sres[[2]][[1]]$"Mean Sq"[1]
MSe <- sres[[2]][[1]]$"Mean Sq"[2]
MSs <- sres[[1]][[1]]$"Mean Sq"[1]
n <- nlevels(data1$s) #群サイズ

omega2 <- dfa*(MSa-MSe)/(SSt+MSs) # ω^2 (3.69)式
omega2p <- dfa*(MSa-MSe)/(dfa*MSa+(n-dfa)*MSe) # ω^2_p  (3.70)式
omega2; omega2p
 

##### \subsection *{3.7.5 2要因被験者間の分散分析}

 
res <- aov(y~a*b,data=data)
sres <- summary(res) 
sres # 分散分析表(表3.12)

SSa <- sres[[1]]$"Sum Sq"[1]
SSb <- sres[[1]]$"Sum Sq"[2]
SSab <- sres[[1]]$"Sum Sq"[3]
SSe <- sres[[1]]$"Sum Sq"[4]
SSt <- SSa + SSb + SSab + SSe

eta2 <- SSa/SSt # 要因Aのη^2 (3.125)式
eta2p <- SSa/(SSa+SSe) # 要因Aのη^2_p (3.126)式
eta2; eta2p

eta2 <- SSab/SSt # 要因Aのη^2 (3.127)式
eta2p <- SSab/(SSa+SSe) # 要因Aのη^2_p (3.128)式
eta2; eta2p

dfa <- sres[[1]]$"Df"[1]
dfb <- sres[[1]]$"Df"[2]
a <- nlevels(data$a) #要因Aの水準数
b <- nlevels(data$b) #要因Bの水準数
MSa <- sres[[1]]$"Mean Sq"[1]
MSb <- sres[[1]]$"Mean Sq"[2]
MSab <- sres[[1]]$"Mean Sq"[3]
MSe <- sres[[1]]$"Mean Sq"[4]
n <- nrow(data)/(a*b) # 群サイズ

sig2a <- dfa/(a*b*n)*(MSa-MSe)
sig2b <- dfb/(a*b*n)*(MSb-MSe)
sig2ab <- (dfa*dfb)/(a*b*n)*(MSab-MSe)
sig2e <- MSe
sig2t <- sig2a + sig2b + sig2ab + sig2e

omega2 <- sig2a / sig2t
omega2p <- sig2a / (sig2a+sig2e)
omega2; omega2p

omega2 <- sig2ab / sig2t
omega2p <- sig2ab / (sig2a+sig2e)
omega2; omega2p
 

\section{第4章 信頼区間 のRプログラム}

##### \subsection *{4.2.1 1標本の場合の母平均の信頼区間}

 
M<-8
s<-0.8
N<-48
se <- s/sqrt(N)

t.crit <- qt(0.975,df=47) # 
CI.upper <- M + se * t.crit # (5.4)式
CI.lower <- M - se * t.crit # (5.4)式
CI.upper; CI.lower
 

##### \subsection *{4.2.4 独立測定における信頼区間}
 
data <- read.csv("shoes.csv") # 履き物データの読み込み
res <- aov(time~cond,data=data)
sres <- summary(res)
sres # 分散分析表 表5.2

n <- 10
MSw <- sres[[1]]$"Mean Sq"[2]  # (5.5)式のMSw
t.crit <- qt(0.975,df=sres[[1]]$Df[2])  # (5.5)式のt.critical

M<- as.vector(by(data$time, data$cond,mean)) #各群の平均ベクトル

CI.upper <- M + t.crit*sqrt(MSw/n) #信頼区間の上限
CI.lower <- M - t.crit*sqrt(MSw/n) #信頼区間の下限
CI.upper; CI.lower
 

##### \subsection {4.2.5 反復測定における信頼区間}
 
data <- read.csv("shoes.csv") # 履き物データの読み込み
res <- aov(time~cond+Error(subj/cond),data=data)
sres <- summary(res)
sres # 分散分析表 表5.3

n <- 10
MSsc <-  sres[[2]][[1]]$"Mean Sq"[2]  # (5.6)式のMSsxc
t.crit <- qt(0.975,df=sres[[1]]$Df[2])  # (5.5)式のt.critical

M<- as.vector(by(data$time, data$cond,mean)) #各群の平均ベクトル

CI.upper <- M + t.crit*sqrt(MSsc/n) #信頼区間の上限
CI.lower <- M - t.crit*sqrt(MSsc/n) #信頼区間の下限
CI.upper; CI.lower
 


##### \subsection *{4.3.1 1標本における頻度の信頼区間}
 
p <- .30
N <- 100
se <- sqrt(p*(1-p)/N)

z.crit <- qnorm(0.975)  # (5.5)式のt.critical

CI.upper <- p + z.crit*se #信頼区間の上限
CI.lower <- p - z.crit*se #信頼区間の下限
CI.upper; CI.lower
 

##### \subsection *{4.3.2 2標本における頻度の信頼区間}
\subsubsection{対応のない場合}
 
p1 <- .30
n1 <- 100
p2 <- .60
n2 <- 100

se <- sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)

z.crit <- qnorm(0.975)  # (5.5)式のt.critical

CI.upper <- (p2 -p1) + z.crit*se #信頼区間の上限
CI.lower <- (p2 -p1) - z.crit*se #信頼区間の下限
CI.upper; CI.lower
 

\subsubsection{対応のある場合}
 
r <- 20
s <- 10
t <- 40
u <- 30
n <- 100
se <- 1/n * sqrt(s + t - (s-t)^2/n)

CI.upper <- (p1 - p2) + z.crit*se #信頼区間の上限
CI.lower <- (p1 - p2) - z.crit*se #信頼区間の下限
 

##### \subsection *{4.4.1 ピアソンの積率相関係数}
 
sandal <- data$time[data$cond=="サンダル"]
barefoot <- data$time[data$cond=="裸足"]
cor.test(sandal,barefoot) 
# 95 percent confidence interval: 以下が信頼区間
 

##### \subsection *{4.5.1 単回帰分析の信頼区間}
 
sandal <- data$time[data$cond=="サンダル"]
barefoot <- data$time[data$cond=="裸足"]

res <- lm(sandal~barefoot) #回帰分析
summary(res) # 回帰分析結果
confint(res) # 回帰係数(と切片)の信頼区間
 

\section{第5章 検定力 のRプログラム}

##### \subsection *{5.6.1 事前の検定力分析 のRプログラム}
 
# dファミリーの効果量についてCohenの基準による「小さな(0.2)」「中程度の(0.5)」
#「大きな(0.8)」差を検出するのに必要な標本サイズを求める場合

library(pwr)
pwr.t.test(d=0.2, power=0.8, sig.level=0.05) #小
pwr.t.test(d=0.5, power=0.8, sig.level=0.05) #中
pwr.t.test(d=0.8, power=0.8, sig.level=0.05) #大
 

##### \subsection *{5.6.2 事後検定力分析 のRプログラム}
 
# 3.7.1節の例において、検定力を分析の事後に求める場合

library(pwr)
pwr.t2n.test(n1=10, n2=10, sig.level=0.05, d=0.88)
 
\section{第6章 さらなる改革に向けて のRプログラム}

##### \subsection *{6.1 メタ分析 のRプログラム}

 
library(meta)
mdata <- read.table("metadata.dat", header=T)
meta.res <- metacont(mdata$n.e, mdata$mean.e, mdata$sd.e, mdata$n.c,
   mdata$mean.c, mdata$sd.c,sm="SMD")
forest(meta.res,comb.random=T)
 

##### \subsection *{6.3 $p_{rep}$ のRプログラム} \label{sec:appprep}

 
x1 <- c(59 ,48 ,51 ,41 ,39 ,84 ,95 ,56 ,86 ,74) #実験群
x2 <- c(47 ,24 ,38 ,28 ,39 ,74 ,77 ,48 ,40 ,60) #統制群

res <- t.test(x1,x2,var.equal=T)

t <- as.numeric(res$statistic) # t値
df <- as.numeric(res$parameter) # df
p <- as.numeric(res$p.value) # p値
p <- p/2 # 両側検定のp値なので半分に

# Lecoutre (2010) のprep
p.rep <- pt(abs(t)/sqrt(2),df=df)
p.rep 

# Kileen (2005)のp_rep
p.rep.old1 <- pnorm(qnorm(1 - p)/sqrt(2))
p.rep .old1

# Kileen (2005)のp_repの近似式
p.rep.old2 <- (1 + (p/(1-p))^(2/3))^(-1)
p.rep.old2