Рассчитать прогнозируемые результаты модели, просматривая переменные

У меня есть несколько моделей, подходящих для прогнозирования результата y = x1 + x2 + .....+x22. Это достаточное количество предикторов и значительное количество моделей. Мои клиенты хотят знать, каково предельное влияние каждого X на предполагаемый y. Модели могут включать в себя сплайны и условия взаимодействия. Я могу сделать это, но это громоздко и требует циклов или большого количества копий-вставок, что является медленным или подвержено ошибкам. Могу ли я сделать это лучше, написав свою функцию по-другому и / или используя purrr или *apply функционировать? Воспроизводимый пример приведен ниже. В идеале я мог бы написать одну функцию и применить ее к longdata,

##  create my fake data.

library(tidyverse)
library (rms)
ltrans<- function(l1){ 
  newvar <- exp(l1)/(exp(l1)+1)
  return(newvar)
}

set.seed(123)
mystates <- c("AL","AR","TN")
mydf <- data.frame(idno = seq(1:1500),state = rep(mystates,500))
mydf$x1[mydf$state=='AL'] <- rnorm(500,50,7)
mydf$x1[mydf$state=='AR'] <- rnorm(500,55,8)
mydf$x1[mydf$state=='TN'] <- rnorm(500,48,10)
mydf$x2 <- sample(1:5,500, replace = T)
mydf$x3 <- (abs(rnorm(1500,10,20)))^2
mydf$outcome <- as.numeric(cut2(sample(1:100,1500,replace = T),95))-1
dd<- datadist(mydf)
options(datadist = 'dd')
m1 <- lrm(outcome ~ x1 + x2+ rcs(x3,3), data = mydf)

dothemath <- function(x1 = x1ref,x2 = x2ref,x3 = x3ref) {
  ltrans(-2.1802256-0.01114239*x1+0.050319692*x2-0.00079289232* x3+
             7.6508189e-10*pmax(x3-7.4686271,0)^3-9.0897627e-10*pmax(x3-    217.97865,0)^3+
           1.4389439e-10*pmax(x3-1337.2538,0)^3)}
x1ref <- 51.4
x2ref <- 3
x3ref <- 217.9
dothemath() ## 0.0591
mydf$referent <- dothemath()
mydf$thisobs <- dothemath(x1 = mydf$x1, x2 = mydf$x2, x3 = mydf$x3)
mydf$predicted <- predict(m1,mydf,type = "fitted.ind") ## yes, matches.
mydf$x1_marginaleffect <- dothemath(x1= mydf$x1)/mydf$referent
mydf$x2_marginaleffect <- dothemath(x2 = mydf$x2)/mydf$referent    
mydf$x3_marginaleffect <- dothemath(x3 = mydf$x3)/mydf$referent

## can I do this with long data?
longdata <- mydf %>%
  select(idno,state,referent,thisobs,x1,x2,x3) %>%
  gather(varname,value,x1:x3)

##longdata$marginaleffect <- dothemath(longdata$varname = longdata$value) ##     no, this does not work.
## I need to communicate to the function which variable it is evaluating. 
longdata$marginaleffect[longdata$varname=="x1"] <- dothemath(x1 =         longdata$value[longdata$varname=="x1"])/
                                                longdata$referent[longdata$varname=="x1"]
longdata$marginaleffect[longdata$varname=="x2"] <- dothemath(x2 = longdata$value[longdata$varname=="x2"])/
                                                    longdata$referent[longdata$varname=="x2"]
longdata$marginaleffect[longdata$varname=="x3"] <- dothemath(x3 = longdata$value[longdata$varname=="x3"])/
                                                    longdata$referent[longdata$varname=="x3"]

testing<- inner_join(longdata[longdata$varname=="x1",c(1,7)],mydf[,c(1,10)])
head(testing) ## yes, both methods work.
Источник

1 ответ

Решение

В основном вы говорите о сгруппированных mutateс оговоркой, что dothemath построен так, что вам нужно указать имя переменной, что можно сделать с помощью do.call или же purrr::invoke чтобы вызвать его по именованному списку параметров:

longdata <- longdata %>% 
    group_by(varname) %>% 
    mutate(marginaleffect = invoke(dothemath, setNames(list(value), varname[1])) / referent)

longdata
#> # A tibble: 4,500 x 7
#> # Groups:   varname [3]
#>     idno state referent thisobs varname value marginaleffect
#>    <int> <fct>    <dbl>   <dbl> <chr>   <dbl>          <dbl>
#>  1     1 AL      0.0591  0.0688 x1       46.1          1.06 
#>  2     2 AR      0.0591  0.0516 x1       50.2          1.01 
#>  3     3 TN      0.0591  0.0727 x1       38.0          1.15 
#>  4     4 AL      0.0591  0.0667 x1       48.4          1.03 
#>  5     5 AR      0.0591  0.0515 x1       47.1          1.05 
#>  6     6 TN      0.0591  0.0484 x1       37.6          1.15 
#>  7     7 AL      0.0591  0.0519 x1       60.9          0.905
#>  8     8 AR      0.0591  0.0531 x1       63.2          0.883
#>  9     9 TN      0.0591  0.0780 x1       47.8          1.04 
#> 10    10 AL      0.0591  0.0575 x1       50.5          1.01 
#> # ... with 4,490 more rows

# the first values look similar
inner_join(longdata[longdata$varname == "x1", c(1,7)], mydf[,c(1,10)])
#> Joining, by = "idno"
#> # A tibble: 1,500 x 3
#>     idno marginaleffect x1_marginaleffect
#>    <int>          <dbl>             <dbl>
#>  1     1          1.06              1.06 
#>  2     2          1.01              1.01 
#>  3     3          1.15              1.15 
#>  4     4          1.03              1.03 
#>  5     5          1.05              1.05 
#>  6     6          1.15              1.15 
#>  7     7          0.905             0.905
#>  8     8          0.883             0.883
#>  9     9          1.04              1.04 
#> 10    10          1.01              1.01 
#> # ... with 1,490 more rows

# check everything is the same
mydf %>% 
    gather(varname, marginaleffect, x1_marginaleffect:x3_marginaleffect) %>% 
    select(idno, varname, marginaleffect) %>% 
    mutate(varname = substr(varname, 1, 2)) %>% 
    all_equal(select(longdata, idno, varname, marginaleffect))
#> [1] TRUE

Может быть проще перенастроить dothemath взять дополнительный параметр имени переменной, чтобы избежать гимнастики.

Источник
Другие вопросы по тегам