Come ottenere l'inverso di una funzione di collegamento (utilizzando $family$linkinv) su un modello memorizzato in una tabella annidata?
Sto elaborando l'output di un modello generato con glm. L'output del modello viene archiviato in una tabella annidata. Voglio calcolare l'intervallo di confidenza attraverso la trasformazione da type= "link" a inverse-link (usando $family$linkinv). Tuttavia, non riesco a farlo funzionare dplyr::mutatein una tabella nidificata perché il modo di estrarre il $family$linkinvè dall'oggetto del modello usando , che sembra non funzionare come previsto nel formato nidificato.model$family$linkinv(x)
sfondo
Questa domanda attuale si basa su una domanda precedente (e risposta scelta) che ho pubblicato sul test del livello di gradimento dei frutti con diversi predittori utilizzando un modello lineare. Conduco una ricerca per capire quale frutto è più simpatico: mango, banana o mela. A tal fine, vado avanti e assaggio 100 persone a caso. Chiedo loro di valutare, su una scala da 1 a 5, il grado di gradimento di ciascuno dei frutti.
Mentre la domanda precedente aveva a che fare con lm, qui sto cercando di utilizzare il quasibinomiale glm. Il problema è che voglio ottenere intervalli di confidenza ma il mio metodo ( glm %>% predict) restituisce SE in "spazio di collegamento", quindi devo passare attraverso un processo di conversione ( dettagliato in questa risposta SO ) per ottenere ciò che voglio.
Dati
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
Voglio testare i miei dati in una scala percentuale, quindi prima li trasformo sottraendo 1 e poi dividendo per 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
Ora uso una pipe per eseguire il modello glm per ogni frutto, ottenere SE nello spazio dei collegamenti e convertire SE in 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]>
Unnesting the predicted_valuesgets:
> 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
Ecco il problema: ora voglio mutare altre due colonne all'interno predicted_values per la trasformazione del collegamento inverso per lower_ci_linke upper_ci_link, ma questo non riesce
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)
)
#########################################################################
))
E ottengo:
Errore: problema con l'
mutate()inputpredicted_values. x Problema con l'mutate()inputlower_ci_inverse_link. x tentativo di applicare la non funzione i Inputlower_ci_inverse_linkè . i L'errore si è verificato nella rigamodel_fit$family$linkinv(lower_ci_link)
- i Input
predicted_valuesèmap(...). i L'errore si è verificato nella riga 1.
Presumo che il problema sia che sto cercando di mutare nuove colonne all'interno predicted_values, ma usando i riferimenti che si trovano a un livello più alto nella tabella nidificata.model_fit$family$linkinv(lower_ci_link)model_fit
Domanda di conclusione
Come posso mutare le colonne di collegamento inverso all'interno predicted_values dell'uso e alla fine ottenere (scorrere fino alle due colonne più a destra):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
APPENDICE
DIMOSTRAZIONE DI COME SONO IN GRADO DI OTTENERE CI CHE DESIDERO SENZA TUBO O DATAFRAME
Il metodo seguente si basa sull'assegnazione di variabili per diversi passaggi lungo il percorso. A scopo dimostrativo, mostra come eseguire il modello e ottenere il $family$linkinvper un solo frutto.
Dati
Come prima, è fruit_liking_dfdopo aver eseguito la trasformazione aritmetica in decimali, quindi:
> 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
Modello
Mi concentrerò solo sui i_love_appledati della colonna ed eseguirò glmsu di essi.
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
)
Predizione
Ora corro predict()su my_modelutilizzando i dati di previsione da 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
Ora converto dalla misura SE che ho ottenuto all'intervallo prediction_link_typedi confidenza (CI) moltiplicando SE per critval(che è stato assegnato con 1.96). Assegno due vettori separati: uno con CI con limite superiore e un altro con CI con limite inferiore:
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)
Quasi lì! Ho i valori CI ma sono nello spazio "link" (perché predict()usati type = "link"). Per riconvertire i valori CI da "link", utilizzo la funzione di collegamento inverso:
lower_ci_inverse_link <- my_model$family$linkinv(lower_ci_link) upper_ci_inverse_link <- my_model$family$linkinv(upper_ci_link)
In sintesi
Anche se questo metodo "vettori" ottiene il lavoro fatto, è non è quello che sto cercando. Invece, voglio incorporare la conversione di "link -> SE -> CI -> inverselink" tramite il tubo introdotto all'inizio di questa domanda.
Risposte
Per fare riferimento ai dati trasmessi mapè necessario utilizzare .x. Prova la risposta seguente.
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 sembra :
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]>
Per ottenere tutti i valori come colonne separate puoi utilizzare 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>