4.4
library(languageR)

#ratings の植物と動物の単語とその親密度の数字
ratings[1:5, c("Word", "meanFamiliarity", "Class")]

#ちなみに箱ひげ図で分布を確認してみる
boxplot(meanFamiliarity ~ Class, data=ratings)

#one-way ANOVA を lm 関数で実行     
ratings.lm <- lm(meanFamiliarity ~ Class, data = ratings)
summary(ratings.lm)
summary(lm(meanFamiliarity ~ Class, data = ratings))

#Factor の部分をどうダミー変数を使って処理するかを見るために、ratings から3変数を dummy に格納する
dummy = ratings[, c("Word", "meanFamiliarity", "Class")]
head(dummy)
#新しく Classplant という変数を作りそこに個体数分、１を格納
dummy$Classplant = 1
head(dummy$Classplant)
#Class が animal になっている行だけ、Classplant をゼロに変更
dummy[dummy$Class == "animal", ]$Classplant = 0

#animal =0, plant =1 
dummy[1:5, ]

summary(lm(meanFamiliarity ~ Classplant, data = dummy))

######
#Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   3.5122     0.1386  25.348  < 2e-16 ***
#Classplant    0.8547     0.2108   4.055 0.000117 ***
#######

#ここでの intercept 3.5122 は animal の meanFamiliarity の平均値と同じ
mean(ratings[ratings$Class =="animal", ]$meanFamiliarity)
coef(ratings.lm)[1] #同一の値になることを確認する

# plant の平均は実はこれに CLassplant の 0.8547 を足した値に等しい
# つまり 3.5122 + 0.8547 = 4.3669 となり、この調整結果も t-test で有意となる

#そのまま t.test()という関数を使ってもよい
#テキストの p.94 にあるグルーピングをあらかじめ作っていないと t-test は動かない
plants = ratings[ratings$Class == "plant", ]
animals = ratings[ratings$Class == "animal", ]
#以下が t-test のコマンド
t.test(animals$meanFamiliarity, plants$meanFamiliarity, var.equal = TRUE)

#######結果#######
#	Two Sample t-test
#
# data:  animals$meanFamiliarity and plants$meanFamiliarity 
# t = -4.0548, df = 79, p-value = 0.0001168
# alternative hypothesis: true difference in means is not equal to 0 
# 95 percent confidence interval:
# -1.2742408 -0.4351257 
# sample estimates:
# mean of x mean of y 
# 3.512174  4.366857    <-- x の平均は intercept, y は intercept + Classplant の値と同じ
#################

# 3グループ以上の検定

head(auxiliaries)　# 285のオランダ語の動詞
###
     Verb     Aux VerbalSynsets Regularity
1 blijken    zijn             1  irregular
2 gloeien  hebben             3    regular
3 glimmen zijnheb             2  irregular
4  rijzen    zijn             4  irregular
5  werpen  hebben             3  irregular
6  delven  hebben             2  irregular
###
# Aux は Verb に列挙してある動詞が完了相でどのような助動詞をとるかを示している 
# zijn = be; hebben = have
# VerbalSynsets は Dutch WordNet での当該動詞がいくつの synset で現れるかという数字　→　つまりそれだけ類義語をたくさん持つということ

# Aux に現れる3種類のパターン(zijn, hebben, 両方)の助動詞のとりかたにより synset の数に差があるかを anova で検定

boxplot(VerbalSynsets ~ Aux, data = auxiliaries) # 平均の差を見る前に分布を視覚化しておくと見当がつきやすい
auxiliaries.lm = lm(VerbalSynsets ~ Aux, data = auxiliaries)
anova(auxiliaries.lm)
###
# Analysis of Variance Table
#
# Response: VerbalSynsets
#            Df Sum Sq Mean Sq F value    Pr(>F)    
# Aux         2  117.8  58.901  7.6423 0.0005859 ***
# Residuals 282 2173.4   7.707                      
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
###
# 自由度(2, 282) で F=7.6423 は p=0.000 以下で有意差がある
#次のようにF分布での残り面積を見ても Pr と同じ(値)=0.0005859 がでる
1-pf(7.6423, 2, 282)

# 3グループのどこに差があるかは summary で読み取る
summary(auxiliaries.lm)

### 結果 ###
# Call:
# lm(formula = VerbalSynsets ~ Aux, data = auxiliaries)
# 
# Residuals:
#    Min     1Q Median     3Q    Max 
# -4.069 -1.467 -0.467  1.533 12.533 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   3.4670     0.1907  18.183  < 2e-16 ***  <-- reference level は hebben 3.4670 は hebben の group mean
# Auxzijn       0.5997     0.7417   0.808 0.419488      <-- zijn の平均値は 3.4670 + 0.5997  t-値に差がない
# Auxzijnheb    1.6020     0.4114   3.894 0.000123 ***  <-- 両方とる場合の平均値は 3.4670 + 1.6020　t-値に差がある
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
# 
# Residual standard error: 2.776 on 282 degrees of freedom
# Multiple R-squared: 0.05141,	Adjusted R-squared: 0.04469 
# F-statistic: 7.642 on 2 and 282 DF,  p-value: 0.0005859 
######

# 3つ以上のグループで t検定を繰り返してはいけない理由

# 5% 成功、95%失敗という2項分布、0 = vector of probability, 3 =試行回数、0.05 各試行の成功確率
# pbinom(0,1,0.05)だと1回だけなので、0.05 だが、試行を繰り返すとだんだん危険率が甘くなってしまう
1 - pbinom (0, 3, 0.05)

#多重比較
## Bonferroni Correction: alpha level を比較するグループ数で割る
## TukeyHSD(): 強力だが、それぞれのグループが同数でなければならない
# 現在のデータが不揃いなことを示す：
xtabs(~auxiliaries$Aux)

# そこで Tukey's HSD をやってみるために新しいデータセット warpbreaks(毛糸の織機ごとの縦糸のほつれの統計)を見る
# 線形モデルでデータを格納
warpbreaks.lm = lm(breaks ~ tension, data = warpbreaks) 
# モデルを anova() でテーブルを出す
anova(warpbreaks.lm) 
summary(warpbreaks.lm)

# Tukey HSD()をかけるために aov 関数で aov オブジェクトに結果を格納
warpbreaks.aov = aov (breaks ~ tension, data = warpbreaks) 
summary (warpbreaks.aov)

# F値の有意水準は以下と同じ
1 - pf (7.206, 2, 51) 


# Tukey HSD() 関数
TukeyHSD(warpbreaks.aov)

# Tukey HSD をプロット。confidence interval がゼロにかぶらないペアが有意差がある
plot(TukeyHSD(warpbreaks.aov))

# 先ほどの Aux と Synset の3グループには分散の差が極めて大きいという問題点が・・・
# tapply 関数で Aux の3グループごとの Synset の数の分散を比較する
tapply(auxiliaries$VerbalSynsets, auxiliaries$Aux, var)

# そのため anova を使わず、Kruskal-Wallis Rank Sum Test を用いる 
kruskal.test(auxiliaries$VerbalSynsets, auxiliaries$Aux)

4.4.1 ANCOVA

# ratings のデータで mean sizes ratings を説明するモデルを立てる
# mean familiarity (numeric)、Class（factor) およびその交互作用を取り込んだモデルを組んでみる
ratings.lm = lm(meanSizeRating ~ meanFamiliarity * Class + I(meanFamiliarity^2), data = ratings)
summary (ratings.lm)

# ratings.lm で作った線形モデルで81個の単語に関して fitted()でモデルによる予測値を算出
fitted(ratings.lm)
# それを $fitted というカラムに追加
ratings$fitted = fitted (ratings.lm)

# 値のplotをするが、最初は type="n"でプロットが NULL 状態で出力
plot(ratings$meanFamiliarity, ratings$meanSizeRating, xlab="mean familiarity", ylab="mean size rating", type="n")
# そのうえに P; A というラベルを貼るコマンド
text (ratings$meanFamiliarity, ratings$meanSizeRating, substr(as.character(ratings$Class), 1, 1),col="darkgrey")

plants = ratings[ratings$Class == "plant",]
animals = ratings[ratings$Class == "animal",]
plants = plants[order(plants$meanFamiliarity),]
animals = animals[order(animals$meanFamiliarity),]

#モデルに予測された値をラインを引く
lines(plants$meanFamiliarity, plants$fitted)
lines(animals$meanFamiliarity, animals$fitted)

4.5 P120
# ２つのカテゴリーデータの統計
xt = xtabs(~ Aux + Regularity, data = auxiliaries) 
xt
# 1 = row 行が 1.00 になるように割合を計算
prop.table(xt, 1)
# 2 = column 列が1.00 になるように割合を計算
prop.table(xt, 2) 
# 全体の proportion は全体表を合計で割る
xt/sum(xt)
mosaicplot(xt, col = TRUE)
# 仮想データで比較
x = data.frame(irregular = c(100, 8, 30), regular = c(77, 6, 22)) 
rownames(x) = c("hebben", "zijn", "zijnheb") 
x
chisq.test(xt)
summary(xt) 
chisq.test(x) 
fisher.test(xt) 

4.6
# 統計的有意差が必ずしも役に立つとは限らない
# シミュレーション
n = 100
(x = seq(1, 100, length = n))
y = 0.3 * x + rnorm(n, 0, 80)　# 回帰モデルで得られる x に対応する yを自動生成
model100 = lm(y ~ x)　# その x, y で回帰分析をする
summary(model100)$coef # 内容を見てみる
# n を大きくすると有意にはなる
n = 1000
(x = seq(1, 100, length = n))
y = 0.3 * x + rnorm(n, 0, 80)
model1000 = lm(y ~ x)
summary(model1000) $coef

# しかし plot してみるとほとんど回帰直線の説明力がないことに気づく
plot(x,y)
abline(lm(y ~ x))
confint(model100)
confint(model1000)

# 逆に effects が小さくても、有意味な現象もありえる