tiidy pglm estimate and get simulations

tidy_pglm.R

tidy_pglm <- function(fit){
  
  summary(fit)$estimate %>%
    tibble::as_tibble() %>%
    dplyr::mutate(param = rownames(summary(fit)$estimate)) %>%
    dplyr::select(param, dplyr::everything()) %>%
    janitor::clean_names(.)
}

predict_glm <- function(fit, new_data, form ){
  coef <- tidy_pglm(fit)
  
  k <- length(fit$estimate)
  
  req_cols <- str_split(as.character(form), "~|\\*|\\+") %>%
    reduce(c) %>%
    discard(~.x == "") %>%
    str_trim
  
  fit_data <- new_data %>%
    dplyr::select(req_cols) %>%
    mutate(.id = 1:nrow(.)) %>%
    drop_na
  
  beta <- coef$estimate %>%
    imap_dfc(~tibble(.x)) %>%
    set_names(coef$param)
  
  sim_est <- new("sim",
                 coef = as.matrix(beta),
                 sigma = rep (sqrt(1), 1)) %>%
    coef %>%
    as_tibble %>%
    dplyr::select(colnames(modelr::model_matrix(fit_data, form))) %>%
    as.matrix()
  
  fit_data %>%
    modelr::model_matrix(form) %>%
    as.matrix %>%
    magrittr::multiply_by_matrix(t(sim_est)) %>%
    as_tibble %>%
    mutate_all(exp) %>%
    rename_all(~str_replace(.x, "V", "sim_")) %>%
    mutate(.id = fit_data$.id) %>%
    right_join(tibble(.id = 1:nrow(new_data))) %>%
    arrange(.id) %>%
    dplyr::select(-.id)
  
}


sim_pglm <- function(fit, n_sim, new_data, form ){
  
  coef <- summary(fit)$estimate %>%
    as_tibble() %>%
    select_at(1:2) %>%
    set_names(c("coef", "sd"))
  
  k <- length(fit$estimate)
  
  V_beta <- vcov(fit) * array(coef$sd,c(k,k)) * t(array(coef$sd,c(k,k)))
  
  beta <- 1:n_sim %>%
    map_dfr(~{
      MASS::mvrnorm(1, coef$coef, V_beta) %>%
        imap_dfc(~tibble(.x) %>% set_names(.y))
    })
  
  req_cols <- str_split(as.character(form), "~|\\*|\\+") %>%
    reduce(c) %>%
    discard(~.x == "") %>%
    str_trim
  
  fit_data <- new_data %>%
    dplyr::select(req_cols) %>%
    mutate(.id = 1:nrow(.)) %>%
    drop_na
  
  sim_est <- new("sim",
                 coef = as.matrix(beta),
                 sigma = rep (sqrt(1), n_sim)) %>%
    coef %>%
    as_tibble %>%
    dplyr::select(colnames(modelr::model_matrix(fit_data, form))) %>%
    as.matrix()
  
  
  
  fit_data %>%
    modelr::model_matrix(form) %>%
    as.matrix %>%
    magrittr::multiply_by_matrix(t(sim_est)) %>%
    as_tibble %>%
    mutate_all(exp) %>%
    rename_all(~str_replace(.x, "V", "sim_")) %>%
    mutate(.id = fit_data$.id) %>%
    right_join(tibble(.id = 1:nrow(new_data))) %>%
    arrange(.id) %>%
    dplyr::select(-.id)
  
}