「検定の解釈」と「ABテスト,検定の統一的な記述」を可能にするinferパッケージ【tidymodels, infer, R】




# 推論に便利なinfer



inferは統計的推論を「簡単に・統一的に」実行できるパッケージ。
今までt.test()やchisq.test()のようにパッケージが分かれていた。
検定だけでなく可視化部分も実行できることから解釈しやすくなった。


inferは以下の4 stepと

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/0432730d-3c6e-41a2-9bbd-f4e44668d4a6.png)

[part of tidymodels 0.5.3　より](https://infer.netlify.app/)

さらに可視化を加えて5 stepで実行する。


- specify()
  - 実行目的を定義
- hypothesize()
  - 帰無仮説の宣言
- generate()
  - 帰無仮説をもとにデータを生成
- calculate()
  - 統計量の分布を作成
- visualize()
  - 統計量分布やp値の可視化


米国の社会的統計調査のデータ「gss」で試しながら紹介


```{r}
library(tidyverse)
library(tidymodels)
```




#step 1 : specify

> specify(x, formula, response = NULL, explanatory = NULL, success = NULL)

データをx
fomulaに式を指定
応答変数をresponse Yにあたる
説明変数をexplanatory Xにあたる
successは割合を検定したい時に指定





gssの調査に回答した者の年齢を調べようと思う場合



```{r}
gss %>%
   specify(response = age)
```



「所属政党」が「年齢」に影響を与えるか調べる場合


```{r}

gss %>%
  specify(age ~ partyid)

```


式で指定せずとも、変数名を名前で入力してもいい
Xとして指定されたものはfactorとして認識され、クラスで認識されるようになる


```{r}

gss %>%
  specify(response = age, explanatory = partyid)
```



# step 2 : hypothesise()

> hypothesise(x, null, p = NULL, mu = NULL, med = NULL, sigma = NULL)

独立性の検定か
点推定か
「null」に指定する

推定したい項目(対象)が何であるか引数に入力

割合をp
平均はmu
中央と分散も検定したければ指定可能


「大学の卒業」と「所属政党」の関係性を独立性の検定で確かめる場合


```{r}
gss %>%
 specify(college ~ partyid, success = "degree") %>%
 hypothesize(null = "independence")

```


一週間の労働時間の検定
政党ごとに平均の点推定をおこなう場合


```{r}
gss %>%
  specify(response = hours) %>%
  hypothesize(null = "point", mu = 40)

```




# step 3 : generate()

> generate(x, reps = 1, type = NULL, ...)

サンプルを疑似生成する場合には、何回生成するかをrepsへ


typeでは生成方法を指定する

bootstrap：置き換え在り
permute：置き換え無し
simulate：理論的な分布を想定している場合、その分布からサンプリング


週40時間労働を検証するために手元のデータを使いbootstrap標本を作成
再標本番号はreplicate列になる



```{r}
gss %>%
 specify(response = hours) %>%
 hypothesize(null = "point", mu = 40) %>%
 generate(reps = 200, type = "bootstrap")
```



```{r}
gss %>%
 specify(partyid ~ age) %>%
 hypothesize(null = "independence") %>%
 generate(reps = 200, type = "permute")

```




# step 4 : calculate()

calculate()で統計量の分布を作り出す

カバーしている統計量は

> stat = c("mean", "median", "sum", "sd", "prop", "count", "diff in means",
  "diff in medians", "diff in props", "Chisq", "F", "slope", "correlation", "t", "z",
  "ratio of props", "odds ratio")


今回はbootstrap標本の各標本から平均を計算している

```{r}
gss %>%
 specify(response = hours) %>%
 hypothesize(null = "point", mu = 40) %>%
 generate(reps = 200, type = "bootstrap") %>%
 calculate(stat = "mean")
```


平均の差の検定では、どちらからどちらを引くかを選ぶ必要があり、これはorderで指定する


```{r}

gss %>%
 specify(age ~ college) %>%
 hypothesize(null = "independence") %>%
 generate(reps = 200, type = "permute") %>%
 calculate("diff in means", order = c("degree", "no degree"))

```





単に手元のデータから統計量を計算する場合、直接calculate()を実行する

F統計量を計算する場合


```{r}
F_hat <- gss %>% 
  specify(age ~ partyid) %>%
  calculate(stat = "F")

F_hat
```


一週間の平均労働時間の点推定(手元の平均値)


```{r}
point_estimate <- gss %>%
  specify(response = hours) %>%
  calculate(stat = "mean")
point_estimate
```








# step 5 : visualize()

visualize()はggplotのhistgramのラッパーであり、
calculate()の計算結果が"stat"になることを認識してヒストグラムをかく。
(このvisualize()は便利なのでヒストグラム書きたい時に変数名を無理やり変更して使ったりしてます)


```{r}
null_distn <- gss %>%
   specify(age ~ partyid) %>%
   hypothesize(null = "independence") %>%
   generate(reps = 1000, type = "permute") %>%
   calculate(stat = "F")

visualize(null_distn) +
  shade_p_value(obs_stat = F_hat, direction = "greater")
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/39fa8b66-9d58-1dd3-c2e2-d04307fdafcf.png)




手元のF統計量を元にした再標本化の分布を元にp値の範囲を可視化する
p値の計算を行う関数は以下


```{r}
null_distn %>%
  get_p_value(obs_stat = F_hat, direction = "greater")
```














































#点推定をしてみる

先ほど計算した一週間の平均労働時間について深堀する

理論的な平均40であった時の統計量の正規分布を求める。
統計量の分布は母集団が分からなければサンプリングを繰り返すこともできないので得られない。
そこで疑似的にbootstrapによって母集団の分散を推定する。
これは別記事を既に書いているのでそちら参照。
[ゼロから始めるbootstrap](https://qiita.com/Ringa_hyj/items/30ed377738c7d25cbeea)






```{r}
null_dist <- gss %>%
  specify(response = hours) %>%
  hypothesize(null = "point", mu = 40) %>%
  generate(reps = 100000, type = "bootstrap") %>%
  calculate(stat = "mean")

null_dist %>%
  visualize()
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/e6f308bf-86b0-6609-7ce8-11101ea2dac6.png)




理論的な平均40のグラフからどれだけ手元の平均値が離れているかを描画してみる。





```{r}
null_dist %>%
  visualize() +
  shade_p_value(obs_stat = point_estimate, direction = "two-sided")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/373ba9a6-90f2-7f2b-656b-7c40ac378aed.png)




p値を確認してみると、約0.03と大変に小さい値であり、
41.38時間とは、平均40時間労働という理想の世界ではほとんど観察されないデータであることが分かった。

もちろん有意水準のαの設定を厳しくして0.01とするならばこの値は平均40時間労働という理想の分布に従っていると考えられる。


```{r}
p_value <- null_dist %>%
  get_p_value(obs_stat = point_estimate, direction = "two-sided")
p_value
```

ブートストラップ法の推定値の信頼区間を求めるには、次のように書きます。



```{r}
null_dist %>%
  get_confidence_interval(point_estimate = point_estimate,
                          level = .95,
                          type = "se")
```








#区間推定してみる

信頼区間を描画することで、推定区間を知る。
理論的な分布がヒストグラムで描かれるのに対して、ブートストラップ法によるサンプリング結果が範囲として描かれる。
平均時間40を理論的な値としているため、この緑色の範囲が40を含んでいなければ、
「今回のgssデータは平均時間40の世界から得られたデータとは考えにくい」
と考察することができる。


```{r}
ci <- null_dist %>%
  get_confidence_interval(point_estimate = point_estimate,
                          level = .95,
                          type = "se")

null_dist %>%
  visualize() +
  shade_confidence_interval(ci)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/6072ee7f-a659-5120-8b1f-4acc6ecdcb1c.png)




こっちのほうが分かりやすい？


```{r}
null_dist %>%
  visualize() +
  shade_confidence_interval(ci)+
shade_p_value(obs_stat = point_estimate, direction = "two-sided")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/0c7bf7aa-745b-7ca0-fa53-da31f308b79a.png)




#検定をしてみる



Chi二乗検定
F検定
t検定
などを統一的なパッケージで実行したい場合に便利。
結果はtibbleで帰ってくるため、デフォルトで入っている検定関数よりもデータとして扱いやすい。


##F検定

先ほどF統計量の計算例を見せた。
この検定の結果やp値について、より直観的に理解したい場合可視化が役に立つ。

generate()を省略してF統計量を計算することを明示しておく。


```{r}
null_f_distn_theoretical <- gss %>%
   specify(age ~ partyid) %>%
   hypothesize(null = "independence") %>%
   calculate(stat = "F")
```


この情報を元に、F分布を描画。
理論的な分布が必要ならば「theoretical」を選択する。


```{r}
visualize(null_f_distn_theoretical, method = "theoretical")
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/68ad6fd7-68c9-cf10-8e6c-5fa55fa55ca0.png)


手元のデータから計算したF統計量を上書きする。


```{r}
F_hat <- gss %>% 
  specify(age ~ partyid) %>%
  calculate(stat = "F")

visualize(null_f_distn_theoretical, method = "theoretical") +
  shade_p_value(obs_stat = F_hat, direction = "greater")
```


理論に基づく帰無分布と無作為化に基づく帰無分布がどのように関連しているかを知るために、無作為化に基づく帰無分布をvisualize()にパイプして、 method = "both" と指定することもできます。



```{r}
null_f_distn <- gss %>%
   specify(age ~ partyid) %>%
   hypothesize(null = "independence") %>%
   generate(reps = 1000, type = "permute") %>%
   calculate(stat = "F")

visualize(null_f_distn, method = "both") +
  shade_p_value(obs_stat = F_hat, direction = "greater")
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/d48a304d-1e1d-f427-56ee-3ff1617aa452.png)



##t検定


```{r}
gss %>%
   drop_na(college) %>%
   t_test(formula = hours ~ college,
      order = c("degree", "no degree"),
      alternative = "two-sided")
```


##chi二乗検定


```{r}
chisq_test(gss, college ~ finrela)
```

##比率の検定



```{r}
prop_test(gss, 
          college ~ sex,  
          order = c("female", "male"))
          
prop_test(gss,
          college ~ NULL,
          p = 0.2)
```












#母集団からサンプリングを繰り返しサンプリング分布を描く場合 : rep_sample_n()


> rep_sample_n(tbl, size, replace = FALSE, reps = 1, prob = NULL)


ひたすらfor文で繰り返すのもセンスが悪い。
結果をtibbleで効率よく返してくれるrep_sample_n()を使う。
リサンプルおよびサンプリング系は「rsample」の領分なので
サンプリングで実現させたいものがあればrsampleパッケージから探してみることをお勧めします。



```{r}
population <- dplyr::storms %>%select(status)

samples <- population %>%
  rep_sample_n(size = 50, reps = 1000)
samples
```



```{r}
samples %>%
  group_by(replicate) %>%
  summarize(prop_hurricane = mean(status == "hurricane"))%>%
  ggplot(aes(x = prop_hurricane)) +
  geom_density() +
  labs(x = "p_hat", y = "Number of samples",
  title = "Sampling distribution of p_hat from 1000 samples of size 50")
```



![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/33537a9d-1c53-a874-49d6-98a7899980c9.png)









































#公式referenceから分析の流れを紹介




アポリポEタンパクが認知機能に影響を与えるか？


```{r}

library(tidymodels)

data(ad_data, package = "modeldata")

cross_table<- ad_data %>%
  select(Genotype, Class)

table(cross_table)
```


アポリポタンパクの組み合わせによって、認知機能に影響が出ているか出ていないかをクロス集計
この割合は発言している場合としていない場合で一緒か？


```{r}
ct<-as.data.frame(table(cross_table))

ct %>%
ggplot()+
  aes(x=Freq, y=Genotype, fill=Class)+ 
  geom_bar(stat = "identity", position = "fill")

```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/248361ea-868d-a46c-2f9c-94bbdb230ac7.png)




カイ二乗統計量をもとめます



```{r}
observed_indep_statistic <- ad_data %>%
  specify(Genotype ~ Class) %>%
  calculate(stat = "Chisq")
observed_indep_statistic
```



帰無分布の二つの作り方をみてみます
一つは無作為に帰無分布をつくるほうほう


```{r}
null_distribution_simulated <- ad_data %>%
  specify(Genotype ~ Class) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 5000, type = "permute") %>%
  calculate(stat = "Chisq")

null_distribution_simulated %>%
  visualize() + 
  shade_p_value(observed_indep_statistic,
                direction = "greater")
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/33a2b040-5c18-f89a-925e-8ac245d8c2ee.png)


もう一つは理論的な帰無分布を作る方法
generarionを消しただけ


```{r}
ad_data %>%
  specify(Genotype ~ Class) %>%
  hypothesize(null = "independence") %>%
  visualize(method = "theoretical") + 
  shade_p_value(observed_indep_statistic,
                direction = "greater")

```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/02e860b6-c94b-97c4-08f7-7e4c525d256a.png)


理論的な帰無分布はシミュレーション(generate())が必要ないので、
generate()したものに上から書く事は簡単に出来る。
bothを指定しましょう


```{r}
null_distribution_simulated %>%
  visualize(method = "both") + 
  shade_p_value(observed_indep_statistic,
                direction = "greater")

```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/60f10746-5721-316b-f70a-83400a063bdb.png)



p値の計算


```{r}
p_value_independence <- null_distribution_simulated %>%
  get_p_value(obs_stat = observed_indep_statistic,
              direction = "greater")
p_value_independence
```


上記のように手順を分かりやすく記述することもできるし、検定用の関数を使えば一発で計算してくれる


```{r}
chisq_test(ad_data, Genotype ~ Class)
```





この遺伝子に関する調査が行われ、論文が2004年に出ました
その結果と、Rに入っているデータを比較してみましょう
論文で調査した対象の遺伝子の割合は


```{r}
meta_rates <- c("E2E2" = 0.71, "E2E3" = 11.4, "E2E4" = 2.32,
                "E3E3" = 61.0, "E3E4" = 22.6, "E4E4" = 2.22)
meta_rates <- meta_rates/sum(meta_rates)

```



今回のデータから見られる遺伝子の割合は


```{r}
obs_rates <- table(ad_data$Genotype)/nrow(ad_data)
round(cbind(obs_rates, meta_rates) * 100, 2)

```




論文の値と、今回のデータのサンプルの値がちがうかどうか
つまり同じような母集団に対して同じようなサンプリングが出来ているのかを調べるため、
比率の検定をカイ二乗検定で行ってみましょう
まずは統計量を計算します


```{r}
observed_gof_statistic <- ad_data %>%
  specify(response = Genotype) %>%
  hypothesize(null = "point", p = meta_rates) %>%
  calculate(stat = "Chisq")
observed_gof_statistic
```




```{r}
null_distribution_gof <- ad_data %>%
  specify(response = Genotype) %>%
  hypothesize(null = "point", p = meta_rates) %>%
  generate(reps = 5000, type = "simulate") %>%
  calculate(stat = "Chisq")


null_distribution_gof %>%
  visualize() + 
  shade_p_value(observed_gof_statistic ,
                direction = "greater")

```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/de63ac4e-96ab-fe26-9604-5916a9d2e3e1.png)




```{r}
p_value_gof <- null_distribution_gof %>%
  get_p_value(observed_gof_statistic,
              direction = "greater")

p_value_gof

```




統計量23から計算したp値は0.002とかなり低く、
可視化した結果を見ても、母比率が同じデータから抽出したサンプルとは思えません。
論文側のデータと同じデータとは言え無さそうです。
今回の調査をinferのラッパーであるchesq_testを使って表現すると、以下の一行です。


```{r}
chisq_test(ad_data, response = Genotype, p = meta_rates)
```







#参考


[Getting to Know infer](https://infer.netlify.app/articles/infer.html)
[infer R Package](https://github.com/tidymodels/infer)
[Tidy Chi-Squared Tests with infer](https://infer.netlify.app/articles/chi_squared.html)
[Tidy t-Tests with infer](https://infer.netlify.app/articles/t_test.html)