This function accept glm and geeglm/gee objects. It outputs the variable names, variables levels, followed by the odds ratios, confidence intervals, and pvalues. Additionally, the reference levels are included in the output which is not included in the original R glm objects.

model_output.r

model_output <- function(mod_obj){
  # libraries
  require(geepack); require(tibble); require(dplyr)
  xvar <- names(lm_object$model)[-1]
  yvar <- names(lm_object$model)[1] # idky i need this but just in case
  data <- lm_object$data
  
  # xvar levels
  xvar_levels <- sapply(data[, xvar], levels)
  xvar_levels_unlist <- unlist(lapply(seq_along(xvar_levels), 
                                      function(x) paste(names(xvar_levels)[[x]], xvar_levels[[x]], sep = "")))
  xvar_levels_df <- data.frame(vec_levels = xvar_levels_unlist)
  
  # odds ratio, confidence interval, and pvalues
  # first, determine if object is geeglm/gee object or not
  if("geeglm" %in% class(lm_object) | "gee" %in% class(lm_object)){
    confint.geeglm <- function(object, parm, level = 0.95, ...) {
      cc <- coef(summary(object))
      mult <- qnorm((1+level)/2)
      citab <- with(as.data.frame(cc),
                    cbind(lwr=Estimate-mult*Std.err,
                          upr=Estimate+mult*Std.err))
      rownames(citab) <- rownames(cc)
      citab[parm,]
    }
    confint_obj <- confint.geeglm(lm_object)
  } else {
    confint_obj <- confint(lm_object)
    colnames(confint_obj) <- c("lwr", "upr")
  }
  or_ci_obj <- as.data.frame(round(cbind(exp(coef(lm_object)), exp(confint_obj)), 2))
  or_ci_obj1 <- rownames_to_column(or_ci_obj, var = "variable") # surprisingly, it works
  or_ci_obj1$ci <- sprintf("%.2f %s %.2f", or_ci_obj1$lwr, "-", or_ci_obj1$upr)
  
  pval_obj <- as.data.frame(round(summary(lm_object)$coef, 3))
  pval_obj1 <- rownames_to_column(pval_obj, var = "variable")
  colnames(pval_obj1)[5] <- "pvalue"
  
  mod_df <- data.frame(or_ci_obj1[, c("variable", "V1", "ci")], pval_obj1[, c("pvalue")])
  colnames(mod_df) <- c("variable", "or", "ci", "pvalue")
  
  # merge xvar_levels_df with mod_df to create REF categories
  # left_join doesn't need standard evaluation...strange...is it just 5 dplyr verbs? must be
  suppressWarnings(merged_df <- left_join(xvar_levels_df, mod_df, by = c("vec_levels" = "variable")))
  merged_df$pvalue1 <- sprintf("%.3f", merged_df$pvalue)
  
  merged_df1 <- replace(merged_df, is.na(merged_df) == TRUE | merged_df == "NA", "Ref")
  merged_df2 <- select_(merged_df1, ~vec_levels, ~or, ~ci, ~pvalue1) #standard evaluation
  
  # insert function
  insertRow <- function(existingDF, newrows) {
    new_idx <- as.integer(newrows[,1]) # get indices of the new rows
    new_idx <- sort(new_idx) + seq(0, length(new_idx) - 1) # adjust for rows shifting due to other insertions 
    old_idx <- seq(nrow(existingDF) + length(new_idx))[-new_idx] # ge indices for the old rows
    existingDF[old_idx,] <- existingDF # assign old rows
    existingDF[new_idx,] <- newrows[,-1] # assign new rows
    existingDF
  }
  # insert row names into model
  list_length <- unlist(lapply(seq_along(xvar_levels), function(x) length(xvar_levels[[x]])))
  merged_df3 <- insertRow(merged_df2, newrows = cbind(cumsum(list_length) - list_length + 1, 
                                                      xvar, "", "", "")) # ugly but gets the job done
  return(merged_df3)
}