CODEWITHSUNDEEP
rpurrr
yearntolearn
yearntolearn

Reputation: 1074

obtain means and multiple confidence intervals of multiple models in R

I would love to obtain both 95% and 90% confidence intervals, as well as means, of multiple models in R.

Data 



data <- data.frame(occup = c(2,3,5,4,2,2,6,1,2,0),
                   unoccup = c(1,2,0,3,0,4,1,1,2,2),
                   month = c("feb", "feb", "feb", "feb", "feb", "mar", "mar", "mar", "mar", "mar"))

I have the function

binomNLL_ratio = function(p, k, N) {
  -sum(dbinom(k, prob = p, size = N, log=TRUE))
}

The required libraries 

library(purrr)
library(bbmle)

Run the script

data %>%
   split(.$month) %>% map(~mle2(minuslog = binomNLL_ratio, start = list(p = 0.5), data = list(N = .$occup + .$unoccup, k = .$occup))) %>%
  map(confint, level = 0.95)

This nicely gives me the 95% confidence intervals of each month. I can also replace 0.95 with 0.9 to get the 90% CI, or replace map(confint) with map(coef) to get the mean of each month model.

Nevertheless, ideally, I would love to get 95% CIs, 90% CIs and means of each model in the same data frame. How can I pass the multiple functions and parameters to get the results that I desire?

Thank you for your help.

Upvotes: 0

Views: 372

Answers (1)

dylanjm
dylanjm

Reputation: 2101

This one was a tricky little bugger and I spent too much time trying to make invoke_map() work. In the end, I settled for a slightly hacky way of doing things but it returns a tibble with both confint() and coef() values:

library(tidyverse)
library(bbmle)
#> Loading required package: stats4
#> 
#> Attaching package: 'bbmle'
#> The following object is masked from 'package:dplyr':
#> 
#>     slice

data <- data.frame(occup = c(2,3,5,4,2,2,6,1,2,0),
                   unoccup = c(1,2,0,3,0,4,1,1,2,2),
                   month = c("feb", "feb", "feb", "feb", "feb", "mar", "mar", "mar", "mar", "mar"))

binomNLL_ratio = function(p, k, N) {
  -sum(dbinom(k, prob = p, size = N, log=TRUE))
}

data %>%
  split(.$month) %>% 
  map(~mle2(minuslog = binomNLL_ratio, start = list(p = 0.5), 
            data = list(N = .$occup + .$unoccup, k = .$occup))) %>% 
  imap_dfr(~ {tibble(confints = confint(.x, level = .95), 
                coefs = coef(.x),
                month = .y)})
#> Warning in dbinom(k, prob = p, size = N, log = TRUE): NaNs produced

#> Warning in dbinom(k, prob = p, size = N, log = TRUE): NaNs produced

#> Warning in dbinom(k, prob = p, size = N, log = TRUE): NaNs produced

#> Warning in dbinom(k, prob = p, size = N, log = TRUE): NaNs produced
#> # A tibble: 4 x 3
#>   confints coefs month
#>      <dbl> <dbl> <chr>
#> 1    0.523 0.727 feb  
#> 2    0.881 0.727 feb  
#> 3    0.317 0.524 mar  
#> 4    0.725 0.524 mar

Created on 2019-02-26 by the reprex package (v0.2.1)

Upvotes: 1

Related Questions