備忘ログ

チラシの裏的備忘録&メモ

Rの`t.test()`や`cor.test()`などの引数`subset`の使い方

R

t.test()やcor.test()などの引数subsetの使い方

Rのt.test()やcor.test()など{stats}の関数でformula形式で値を指定した場合に使える引数subsetの使い方がヘルプを読んでもいまいちわかりにくかったのでメモしておく。

使い方としては{base}のsbuset()の引数subsetと同じで指定した条件でデータを取り出し、そのデータだけを分析してくれるようにするというもの。

例えば、irisデータの場合にSepal.LengthとSepal.Widthの相関関係をSpeciesで"setosa"となっているものだけでみたい場合は次のようにするとできる。

cor.test(~ Sepal.Length + Sepal.Width, data = iris, subset = Species == "setosa")

## 
##  Pearson's product-moment correlation
## 
## data:  Sepal.Length and Sepal.Width
## t = 7.6807, df = 48, p-value = 6.71e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5851391 0.8460314
## sample estimates:
##       cor 
## 0.7425467

つまり、次のようにデータの一部を取り出したのもを分析するのと同じになる。

# {base}のsubset()を使う場合
cor.test(~ Sepal.Length + Sepal.Width, data = subset(iris, subset = Species == "setosa"))

# {dplyr}のfilter()を使う場合
# `|>`をつかってパイプして書いてみる
# 最後にプレースホルダーで_を使っている
iris |> 
  dplyr::filter(Species == "setosa") |> 
  cor.test(~ Sepal.Length + Sepal.Width, data = _)

また、引数subsetの中の条件の書き方は、{base}のsubset()と同じ書き方をする。

複数条件指定した場合は、&などで条件式を増やしたり%in%で項目の複数条件を指定したら良い。

cor.test(~ Sepal.Length + Sepal.Width, data = iris, subset = Species == "setosa" & Sepal.Length > 5)

## 
##  Pearson's product-moment correlation
## 
## data:  Sepal.Length and Sepal.Width
## t = 2.5447, df = 20, p-value = 0.0193
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.0921590 0.7580924
## sample estimates:
##       cor 
## 0.4945516

これは{base}のsubset()を使った場合は次のようにかける。

cor.test(~ Sepal.Length + Sepal.Width, data = subset(iris, subset = Species == "setosa" & Sepal.Length > 5))

そういう引数になっている。全体としてsubset()使ってもほぼ同じようにかけると思うので、引数のsubsetの優位性はよくわからない。

おまけ

irisのSepal.LengthとSepal.Widthはシンプソンのパラドックス的な状態になっている。

library(ggplot2)
iris |> 
  ggplot() +
  geom_point(aes(x = Sepal.Length, y = Sepal.Width, color = Species)) +
  geom_smooth(aes(x = Sepal.Length, y = Sepal.Width, color = Species), method = "lm" , se = FALSE) +
  geom_smooth(aes(x = Sepal.Length, y = Sepal.Width), method = "lm", se = FALSE)

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

こういう場合に、全体でcor.test()で相関関係をみると

cor.test(~ Sepal.Length + Sepal.Width, data = iris)

## 
##  Pearson's product-moment correlation
## 
## data:  Sepal.Length and Sepal.Width
## t = -1.4403, df = 148, p-value = 0.1519
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.27269325  0.04351158
## sample estimates:
##        cor 
## -0.1175698

となる。

これをsubsetで各種類で投入して分析するとそれぞれの種類ごとの相関関係が見られる。

三種類なのでそんなに手間ではないが、楽をするためにsapply()関数でループを回す。

value <- as.character(unique(iris$Species))
sapply(value, function(value){
  cor.test(~ Sepal.Length + Sepal.Width, 
           data = iris, 
           subset = Species == value)
  }, simplify = FALSE)

## $setosa
## 
##  Pearson's product-moment correlation
## 
## data:  Sepal.Length and Sepal.Width
## t = 7.6807, df = 48, p-value = 6.71e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5851391 0.8460314
## sample estimates:
##       cor 
## 0.7425467 
## 
## 
## $versicolor
## 
##  Pearson's product-moment correlation
## 
## data:  Sepal.Length and Sepal.Width
## t = 4.2839, df = 48, p-value = 8.772e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2900175 0.7015599
## sample estimates:
##       cor 
## 0.5259107 
## 
## 
## $virginica
## 
##  Pearson's product-moment correlation
## 
## data:  Sepal.Length and Sepal.Width
## t = 3.5619, df = 48, p-value = 0.0008435
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2049657 0.6525292
## sample estimates:
##       cor 
## 0.4572278

全体では弱い負の相関関係（しかもp = 0.1519）が、種別にすると正の相関関係がみえてくる。

subsetを使うと簡単に探索できる。が、前述した通り、subset()でもいい。