# Down's syndromeのデータ
if (!require(dsrTest)) { install.packages("dsrTest", dep=TRUE); library(dsrTest) }
levels(downs.mi$Age)[7] <- levels(downs.mi$Age)[6] # なぜか40+という変なカテゴリが入っているので
barplot(tapply(downs.mi$Cases, downs.mi$Age, sum)/tapply(downs.mi$Births, downs.mi$Age, sum))
barplot(tapply(downs.mi$Cases, downs.mi$BirthOrder, sum)/tapply(downs.mi$Births, downs.mi$BirthOrder, sum))
downs.mi$Risk <- downs.mi$Cases/downs.mi$Births
barplot(Risk ~ BirthOrder+Age, data=downs.mi, beside=TRUE)

# プール化の例として交絡の例(2)のデータを分析
fisher.test(matrix(c(139, 443, 230, 502), 2)) # 年齢は無視してオッズ比を求める
y <- array(c(2, 53, 1, 61, 3, 121, 5, 152, 14, 95, 7, 114, 27, 103, 12, 66, 
 51, 64, 40, 81, 29, 7, 101, 28, 13, 0, 64, 0), dim=c(2, 2, 7))
for (i in 1:7) {
 print(fisher.test(y[, , i])) 
 # 年齢層別にオッズ比を求めると，喫煙者の方が非喫煙者より死亡リスクが高い年齢層が多い
}
mantelhaen.test(y)

x <- read.delim("https://minato.sip21c.org/advanced-statistics/ecopx.txt") #読み込み
x$Q04c <- ifelse(x$Q04>=2.5, 1, 0) # スコアなので1/0に変換
x$Q10c <- ifelse(x$Q10>=4.8, 1, 0) # スコアなので1/0に変換 
res <- glm(Q04c ~ Q10c+AGE+SEX, family=binomial(logit), data=x) # ロジスティック回帰
summary(res) #回帰分析の結果を出力
exp(coef(res)) #対数オッズ比なので係数のexp()をとってオッズ比を得る
exp(confint(res)) # 信頼区間のexp()をとってオッズ比の95%信頼区間を得る
library(fmsb) # fmsbパッケージを読み込む
NagelkerkeR2(res) # モデルの説明力を示すNagelkerkeのR2の計算
library(stargazer) # きれいに出力するためのstargazerパッケージを読み込む
sink("logistic-res.tex") # このファイルをTeXLiveでコンパイルする
cat("\\documentclass[a4paper,12pt]{article}\n\\begin{document}\n")
stargazer(res, ci=TRUE, apply.coef=exp, type="latex", 
 notes=sprintf("Nagelkerke's R-squared: %4.2f", NagelkerkeR2(res)$R2))
cat("\\end{document}\n")
sink()

y <- predict(res, data.frame(Q10c=x$Q10c, AGE=x$AGE, SEX=x$SEX), "response")
if (require(Epi)==FALSE) { install.packages("Epi", dep=TRUE); library(Epi) }
ROC(y, x$Q04c)

# 効果量
library(effsize)
data(sleep) # 2種類の催眠剤(group)を10人（IDが個人を示す）に与えたときの睡眠時間増加(extra)
t.test(extra ~ group, data=sleep) # 独立2群としてt検定
with(sleep, t.test(extra[group == 1], extra[group == 2], paired = TRUE)) # 個人差考慮
cohen.d(extra ~ group, data=sleep) # Cohenのd
cohen.d(sleep$extra[sleep$group==2], sleep$extra[sleep$group==1], paired=TRUE) # 個人差考慮
cohen.d(extra ~ group, data=sleep, hedges.correction=TRUE) # Hedgesのg
cohen.d(sleep$extra[sleep$group==2], sleep$extra[sleep$group==1], 
 paired=TRUE, hedges.correction=TRUE) # 個人差考慮
library(lsr) # 一元配置分散分析における主効果の大きさを示すetaSquared()を使う
res <- aov(weight ~ feed, data=chickwts) # 鶏体重データのANOVA
summary(res) # 普通のANOVA
etaSquared(res) # η^2を計算

# 鶏体重データを多重比較
pairwise.t.test(chickwts$weight, chickwts$feed, p.adjust.method="holm") # Holm法
pairwise.t.test(chickwts$weight, chickwts$feed, p.adjust.method="fdr") # FDR法
oneway.test(weight ~ feed, data=chickwts) # 等分散を仮定しないWelchの拡張によるANOVA
if (!require(rosetta)) { install.packages("rosetta", dep=TRUE);
 library(rosetta) }
posthocTGH(chickwts$weight, chickwts$feed, method="games-howell") # Games-Howell法

# 回帰分析
library(foreign)
x <- read.spss("./IDIR71SV2/IDIR71FL.SAV") # データはDHSから
# V012: age, V024: Province, V025: Urban/Rural, V106: Highest education
# V113: Drinking water type, V116: Toilet facility type, V201: Parity
# V361: Contraceptive use, M4: Duration of breastfeeding, V481: Health insurance
print(res1 <- lm(V201 ~ V012 + V025 + V106 + V361, data=x))
summary(res1)
library(lmerTest)
print(res2 <- lmer(V201 ~ V012 + V106 + V361 + (1 | V025), data=x, REML=FALSE))
summary(res2)
print(res3 <- glm(V201 ~ V012 + V025 + V106 + V361, data=x, family="poisson"))
summary(res3)
library(stargazer)
sink("IDIR71SV-parity.tex")
cat("\\documentclass[a4paper,12pt]{article}\n\\begin{document}\n")
# lm and glm can be show in the same table
stargazer(res1, res3)
cat("\\end{document}\n")
sink()