リンク関数の逆関数を取得する方法（を使用して $family$linkinv）ネストされたtibbleに保存されたモデルで？
Aug 24 2020
で生成されたモデルの出力を作成していglmます。モデルの出力は、ネストされたティブルに保存されます。type= "link"からinverse-link（を使用$family$linkinv）への変換を通じて信頼区間を計算したいと思います。ただし、を使用してモデルオブジェクトからdplyr::mutateプルする方法であるため、ネストされたティブルで動作させることはできません。これは、ネストされた形式では意図したとおりに機能しないようです。$family$linkinvmodel$family$linkinv(x)
バックグラウンド
この現在の質問は、線形モデルを使用してさまざまな予測子による果物の好みのレベルをテストすることについて投稿した以前の質問（および選択した回答）に基づいています。私は、マンゴー、バナナ、リンゴのどれがより好感が持てるのかを調べるために調査を行っています。この目的のために、私は先に進み、ランダムに100人をサンプリングします。それぞれの果物の好みの程度を1〜5のスケールで評価してもらいます。
前の質問はと関係がありましたがlm、ここでは準二項を利用しようとしていますglm。問題は、信頼区間を取得したいのですが、私のメソッド（glm %>% predict）は「リンク空間」にSEを出力するため、必要なものを取得するために変換プロセス（このSO回答で詳しく説明）を実行する必要があります。
データ
library(tidyverse)
library(magrittr)
set.seed(123)
fruit_liking_df <-
data.frame(
id = 1:100,
i_love_apple = sample(c(1:5), 100, replace = TRUE),
i_love_banana = sample(c(1:5), 100, replace = TRUE),
i_love_mango = sample(c(1:5), 100, replace = TRUE),
age = sample(c(20:70), 100, replace = TRUE),
is_male = sample(c(0, 1), 100, prob = c(0.2, 0.8), replace = TRUE),
education_level = sample(c(1:4), 100, replace = TRUE),
is_colorblinded = sample(c(0, 1), 100, replace = TRUE)
)
> as_tibble(fruit_liking_df)
## # A tibble: 100 x 8
## id i_love_apple i_love_banana i_love_mango age is_male education_level is_colorblinded
## <int> <int> <int> <int> <int> <dbl> <int> <dbl>
## 1 1 3 5 2 50 1 2 0
## 2 2 3 3 1 49 1 1 0
## 3 3 2 1 5 70 1 1 1
## 4 4 2 2 5 41 1 3 1
## 5 5 3 1 1 49 1 4 0
## 6 6 5 2 1 29 0 1 0
## 7 7 4 5 5 35 1 3 0
## 8 8 1 3 5 24 0 3 0
## 9 9 2 4 2 55 1 2 0
## 10 10 3 4 2 69 1 4 0
## # ... with 90 more rows
データをパーセントスケールでテストしたいので、最初に1を引き、次に4で割って変換します。
fruit_liking_df %<>%
mutate_at(vars(starts_with("i_love_")), ~ subtract(., 1) %>% divide_by(., 4))
> as_tibble(fruit_liking_df)
## # A tibble: 100 x 8
## id i_love_apple i_love_banana i_love_mango age is_male education_level is_colorblinded
## <int> <dbl> <dbl> <dbl> <int> <dbl> <int> <dbl>
## 1 1 0.5 1 0.25 50 1 2 0
## 2 2 0.5 0.5 0 49 1 1 0
## 3 3 0.25 0 1 70 1 1 1
## 4 4 0.25 0.25 1 41 1 3 1
## 5 5 0.5 0 0 49 1 4 0
## 6 6 1 0.25 0 29 0 1 0
## 7 7 0.75 1 1 35 1 3 0
## 8 8 0 0.5 1 24 0 3 0
## 9 9 0.25 0.75 0.25 55 1 2 0
## 10 10 0.5 0.75 0.25 69 1 4 0
## # ... with 90 more rows
ここで、パイプを使用して各フルーツのglmモデルを実行し、リンクスペースでSEを取得し、SEをCIに変換します。
## will be needed later
my_new_data_for_pred <- expand_grid(
age = 45,
is_male = .5,
education_level = 2.5,
is_colorblinded = 0.5
)
## will be needed later
critval <- 1.96
model_fits_grouped <-
fruit_liking_df %>%
pivot_longer(starts_with("i_love"), values_to = "fruit") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(
data,
~ glm(
data = .x,
fruit ~ I(age - 45) +
I((age - 45) ^ 2) +
I(is_male - .5) +
I(education_level - 2) +
is_colorblinded,
family = quasibinomial
)
)) %>%
mutate(predicted_values = map(
model_fit,
~ bind_cols(my_new_data_for_pred,
as.data.frame(
predict(
newdata = my_new_data_for_pred,
.x,
type = "link",
interval = "confidence",
level = 0.95,
se.fit = T
)
)) %>%
rowwise() %>%
mutate(
estimate = fit,
lower_ci_link = fit - critval * se.fit,
upper_ci_link = fit + critval * se.fit
)
))
> model_fits_grouped
## # A tibble: 3 x 4
## # Groups: name [3]
## name data model_fit predicted_values
## <chr> <list> <list> <list>
## 1 i_love_apple <tibble [100 x 6]> <glm> <tibble [1 x 10]>
## 2 i_love_banana <tibble [100 x 6]> <glm> <tibble [1 x 10]>
## 3 i_love_mango <tibble [100 x 6]> <glm> <tibble [1 x 10]>
取得のネストを解除するpredicted_values：
> model_fits_grouped %>% unnest(predicted_values)
## # A tibble: 3 x 13
## # Groups: name [3]
## name data model_fit age is_male education_level is_colorblinded fit se.fit residual.scale estimate lower_ci_link upper_ci_link
## <chr> <list> <list> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 i_love_apple <tibble [100 x 6~ <glm> 45 0.5 2.5 0.5 0.0843 0.261 0.709 0.0843 -0.427 0.595
## 2 i_love_banana <tibble [100 x 6~ <glm> 45 0.5 2.5 0.5 -0.0718 0.286 0.781 -0.0718 -0.633 0.489
## 3 i_love_mango <tibble [100 x 6~ <glm> 45 0.5 2.5 0.5 -0.140 0.279 0.762 -0.140 -0.687 0.407
ここで問題です。今、私は2つの以上の列を変異させたい範囲内 predicted_valuesのため、逆リンクのための変革lower_ci_linkとupper_ci_link、これは失敗しました
model_fits_grouped <-
fruit_liking_df %>%
pivot_longer(starts_with("i_love"), values_to = "fruit") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(
data,
~ glm(
data = .x,
fruit ~ I(age - 45) +
I((age - 45) ^ 2) +
I(is_male - .5) +
I(education_level - 2) +
is_colorblinded,
family = quasibinomial
)
)) %>%
mutate(predicted_values = map(
model_fit,
~ bind_cols(my_new_data_for_pred,
as.data.frame(
predict(
newdata = my_new_data_for_pred,
.x,
type = "link",
interval = "confidence",
level = 0.95,
se.fit = T
)
)) %>%
rowwise() %>%
mutate(
estimate = fit,
lower_ci_link = fit - critval * se.fit,
upper_ci_link = fit + critval * se.fit
) %>%
######################### this addition fails ###########################
mutate(
lower_ci_inverse_link = model_fit$family$linkinv(lower_ci_link),
upper_ci_inverse_link = model_fit$family$linkinv(upper_ci_link)
)
#########################################################################
))
そして私は得る：
エラー：
mutate()入力に問題がありますpredicted_values。xmutate()入力に問題がありますlower_ci_inverse_link。x非機能の適用を試みますi入力lower_ci_inverse_linkはです。i行でエラーが発生しましたmodel_fit$family$linkinv(lower_ci_link)
- i入力
predicted_valuesはmap(...)です。i行1でエラーが発生しました。
問題は、内の新しい列を変更しようとしていることだと思いますがpredicted_values、参照を使用すると、ネストされたティブルの上位レベルにあります。model_fit$family$linkinv(lower_ci_link)model_fit
結論の質問
使用中に 逆リンク列を変更predicted_valuesし、最終的に取得するにはどうすればよいですか（右端の2つの列までスクロールします）。model_fit$family$linkinv(lower_ci_link)model_fit$family$linkinv(upper_ci_link)
> model_fits_grouped %>% unnest(predicted_values)
## # A tibble: 3 x 15
## # Groups: name [3]
## name data model_fit age is_male education_level is_colorblinded fit se.fit residual.scale estimate lower_ci_link upper_ci_link lower_ci_inverse_link_*DEMO* upper_ci_inverse_link_*DEMO*
## <chr> <list> <list> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 i_love_apple <tibble [100 x 6]> <glm> 45 0.5 2.5 0.5 0.521 0.0632 0.349 0.521 0.397 0.645 0.111 0.111
## 2 i_love_banana <tibble [100 x 6]> <glm> 45 0.5 2.5 0.5 0.482 0.0701 0.387 0.482 0.345 0.620 0.222 0.222
## 3 i_love_mango <tibble [100 x 6]> <glm> 45 0.5 2.5 0.5 0.465 0.0683 0.377 0.465 0.331 0.599 0.333 0.333
付録
私は私が欲しいものを得ることができるよHOWの実証なしA PIPE OR DATAFRAME
次の方法は、途中のいくつかのステップに変数を割り当てることに依存しています。デモンストレーションのために、モデルを実行して$family$linkinv1つのフルーツを取得する方法を示します。
データ
前と同じように、fruit_liking_df小数への算術変換を行った後なので、次のようになります。
> as_tibble(fruit_liking_df)
## # A tibble: 100 x 8
## id i_love_apple i_love_banana i_love_mango age is_male education_level is_colorblinded
## <int> <dbl> <dbl> <dbl> <int> <dbl> <int> <dbl>
## 1 1 0.5 1 0.25 50 1 2 0
## 2 2 0.5 0.5 0 49 1 1 0
## 3 3 0.25 0 1 70 1 1 1
## 4 4 0.25 0.25 1 41 1 3 1
## 5 5 0.5 0 0 49 1 4 0
## 6 6 1 0.25 0 29 0 1 0
## 7 7 0.75 1 1 35 1 3 0
## 8 8 0 0.5 1 24 0 3 0
## 9 9 0.25 0.75 0.25 55 1 2 0
## 10 10 0.5 0.75 0.25 69 1 4 0
## # ... with 90 more rows
モデル
i_love_apple列データのみに焦点を当てて実行glmします。
my_model <-
glm(
i_love_apple ~
I(age - 45) +
I((age - 45) ^ 2) +
I(is_male - 0.5) +
I(education_level - 2) +
I(is_colorblinded - 0.5),
family = quasibinomial,
data = fruit_liking_df
)
予測
今私はからの予測データpredict()をmy_model使用して実行しますmy_new_data_for_pred：
prediction_link_type <-
predict(object = my_model,
newdata = my_new_data_for_pred,
type = "link", ## <------------ type = "link" is crucial to note
interval = "confidence",
level = 0.95,
se.fit = TRUE)
> prediction_link_type
## $fit ## 1 ## 0.08427577 ## $se.fit
## [1] 0.2606326
## $residual.scale
## [1] 0.7090294
ここでprediction_link_type、SEにcritval（で割り当てられている）を掛けることにより、取得したSEメジャーから信頼区間（CI）に変換し1.96ます。2つの別々のベクトルを割り当てます。1つは上限CIを持ち、もう1つは下限CIを持ちます。
lower_ci_link <- prediction_link_type$fit - (critval * prediction_link_type$se.fit) upper_ci_link <- prediction_link_type$fit + (critval * prediction_link_type$se.fit)
もうすぐです！CI値を取得しましたが、それらは「リンク」スペースにあります（predict()使用されているためtype = "link"）。CI値を「リンク」から元に戻すには、逆リンク関数を使用します。
lower_ci_inverse_link <- my_model$family$linkinv(lower_ci_link) upper_ci_inverse_link <- my_model$family$linkinv(upper_ci_link)
要約すれば
この「ベクトル」メソッドは仕事を成し遂げますが、それは私が探しているものではありません。代わりに、この質問の冒頭で紹介したパイプを介して、「リンク-> SE-> CI->逆リンク」の変換を組み込みたいと思います。
回答
1 RonakShah Aug 25 2020 at 16:46
渡されたデータを参照するには、mapを使用する必要があります.x。以下の答えを試してください。
library(tidyverse)
result <- fruit_liking_df %>%
pivot_longer(starts_with("i_love"), values_to = "fruit") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(
data,
~ glm(
data = .x,
fruit ~ I(age - 45) +
I((age - 45) ^ 2) +
I(is_male - .5) +
I(education_level - 2) +
is_colorblinded,
family = quasibinomial
)
)) %>%
mutate(predicted_values = map(
model_fit,
~ bind_cols(my_new_data_for_pred,
as.data.frame(
predict(
newdata = my_new_data_for_pred,
.x,
type = "link",
interval = "confidence",
level = 0.95,
se.fit = T
)
)) %>%
rowwise() %>%
mutate(
estimate = fit,
lower_ci_link = fit - critval * se.fit,
upper_ci_link = fit + critval * se.fit,
lower_ci_inverse_link = .x$family$linkinv(lower_ci_link),
upper_ci_inverse_link = .x$family$linkinv(upper_ci_link)
)))
result のように見えます：
result
# name data model_fit predicted_values
# <chr> <list> <list> <list>
#1 i_love_apple <tibble [100 × 6]> <glm> <tibble [1 × 12]>
#2 i_love_banana <tibble [100 × 6]> <glm> <tibble [1 × 12]>
#3 i_love_mango <tibble [100 × 6]> <glm> <tibble [1 × 12]>
すべての値を個別の列として取得するには、次を使用できますunnest_wider。
result %>% unnest_wider(predicted_values)
# name data model_fit age is_male education_level is_colorblinded fit se.fit
# <chr> <lis> <list> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#1 i_lo… <tib… <glm> 45 0.5 2.5 0.5 0.0843 0.261
#2 i_lo… <tib… <glm> 45 0.5 2.5 0.5 -0.0718 0.286
#3 i_lo… <tib… <glm> 45 0.5 2.5 0.5 -0.140 0.279
# … with 6 more variables: residual.scale <dbl>, estimate <dbl>, lower_ci_link <dbl>,
# upper_ci_link <dbl>, lower_ci_inverse_link <dbl>, upper_ci_inverse_link <dbl>