UPDATE: 2023-05-20 20:36:16
ブログからの引っ越し記事。
ここではRstudio community-Quasiquotation inside a formula-で議論されていた内容をもとに自分でも回帰分析を実行する関数を作ってみる。
コミュニティのお題は「Quasiquotation inside a
formula」という内容を議論している。formulaの中で、Tidyevalはどう使うべきか？的な話だと思う。そこで作られていた関数が下記の線形混合モデルをパイプの中で、NSEで使う関数。部分的に変数名を変えていたり、名前空間を追加したりしているが、動きはおなじ。
library(tidyverse)
library(rlang)
library(broom)
library(lme4)
quasi_lmer <- function(data, lhs, ...) {
lhs_enq <- rlang::enquo(lhs)
rhs_enq <- rlang::enquos(...)
rhs_combined_enq <- paste0(purrr::map(.x = rhs_enq,
.f = function(x){rlang::quo_text(x)}),
collapse = " + ")
form <- as.formula(paste0(rlang::quo_text(lhs_enq), " ~ ", rhs_combined_enq))
data %>% lme4::lmer(data = ., formula = form)
}
iris %>%
quasi_lmer(data = .,
lhs = Sepal.Width,
Sepal.Length, Petal.Length, (1 | Species))
Linear mixed model fit by REML ['lmerMod']
Formula: Sepal.Width ~ Sepal.Length + Petal.Length + (1 | Species)
Data: .
REML criterion at convergence: 74.0416
Random effects:
Groups Name Std.Dev.
Species (Intercept) 0.4531
Residual 0.2899
Number of obs: 150, groups: Species, 3
Fixed Effects:
(Intercept) Sepal.Length Petal.Length
1.05585 0.39911 -0.08798
学びがあったのは、...で複数の説明変数や線形混合モデルで使われる(1|x)という独自の記述の渡し方の部分。過去に似たような関数を作ったときにはTidyevalすら使わずに作っていた…実際に、議論の中でもTidyeval使わなくてもできる的な話もある。
paste0(
purrr::map(.x = rhs_enq, .f = function(x){rlang::quo_text(x)}),
collapse = " + ")確かにクオートしたものを文字列にしたければ、quo_text()で変換できるので、これで...で引き受けたものをもとに、回帰分析の説明変数の部分を作ることができる。イメージはこんな感じ。listをpaste()すると、リストを単一の文字列に変換できる。それをもとに単純にmap()で繰り返している。確かにこうやっても作れる、勉強になる。
paste(list("x1", "x2", "x^3", "(1|x)"), collapse = " + ")
[1] "x1 + x2 + x^3 + (1|x)"これをもとに使い方を少し変えて、単回帰分析を連続して行う関数を作ってみた。何らかの目的変数yに対して、一気にx1、x2、x3、、、x_nと回帰分析を実行する関数。...で渡すのではなく、vars()を使って外部クオートさせている。外部クオートすると、クオージャーのリストが得られるので、それを文字列に変換して、回帰式を作っている。
役に立つかどうかは…わからん。あくまでもTidyevalの練習なので。
quasi_lms <- function(data, lhs, rhs) {
lhs_enq <- rlang::enquo(lhs)
# 下記は不要、ここでは変数名を使って明示的に外部クオートしていることを示したもの
rhs_enq <- rhs
form <- purrr::map(.x = rhs_enq,
.f = function(x){
as.formula(paste0(rlang::quo_text(lhs_enq), " ~ ", rlang::quo_text(x)))
})
res <- form %>%
map(.x = .,
.f = function(x){
lm(formula = x, data = data)
})
return(res)
}回帰式を作っている部分はこんなイメージ。
lhs_q <- rlang::quo(Sepal.Width)
rhs_q <- rlang::quos(Sepal.Width, Petal.Length, Petal.Width)
purrr::map(.x = rhs_q,
.f = function(x){
as.formula(paste0(rlang::quo_text(lhs_q), " ~ ", rlang::quo_text(x)))
})
[[1]]
Sepal.Width ~ Sepal.Width
<environment: 0x7fcd3b8650c8>
[[2]]
Sepal.Width ~ Petal.Length
<environment: 0x7fcd3b85a698>
[[3]]
Sepal.Width ~ Petal.Width
<environment: 0x7fcd58748598>あとはこれをパイプの中で使うと、さきほどの回帰式毎に計算ができるし、map_dfr()と{broom}の関数を組みあせれば下記のようにデータフレームで情報を確認できる。
iris %>%
quasi_lms(lhs = Sepal.Length, rhs = vars(Sepal.Width, Petal.Length, Petal.Width)) %>%
map_dfr(.x = ., .f = function(x){tidy(x)})
# A tibble: 6 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 6.53 0.479 13.6 6.47e- 28
2 Sepal.Width -0.223 0.155 -1.44 1.52e- 1
3 (Intercept) 4.31 0.0784 54.9 2.43e-100
4 Petal.Length 0.409 0.0189 21.6 1.04e- 47
5 (Intercept) 4.78 0.0729 65.5 3.34e-111
6 Petal.Width 0.889 0.0514 17.3 2.33e- 37
iris %>%
quasi_lms(lhs = Sepal.Length, rhs = vars(Sepal.Width, Petal.Length, Petal.Width)) %>%
map_dfr(.x = ., .f = function(x){glance(x)})
# A tibble: 3 x 11
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
<dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int>
1 0.0138 0.00716 0.825 2.07 1.52e- 1 2 -183. 372. 381. 101. 148
2 0.760 0.758 0.407 469. 1.04e-47 2 -77.0 160. 169. 24.5 148
3 0.669 0.667 0.478 299. 2.33e-37 2 -101. 208. 217. 33.8 148
以上でおしまい。