1#' Compute the Statistical Mode of a Vector 2#' @aliases Mode mode 3#' @param x a vector of numeric, factor, or ordered values 4#' @return the statistical mode of the vector. If more than one mode exists, 5#' the last one in the factor order is arbitrarily chosen (by design) 6#' @export 7#' @author Christopher Gandrud and Matt Owen 8 9Mode <- function (x) { 10 # build a table of values of x 11 tab <- table(as.factor(x)) 12 # find the mode, if there is more than one arbitrarily pick the last 13 max_tab <- names(which(tab == max(tab))) 14 v <- max_tab[length(max_tab)] 15 # if it came in as a factor, we need to re-cast it as a factor, with the same exact levels 16 if (is.factor(x)) 17 return(factor(v, levels = levels(x))) 18 # re-cast as any other data-type 19 as(v, class(x)) 20} 21 22## Zelig 3 and 4 backward compatibility 23## This enables backward compatibility, but results in a warning when library attached 24# mode <- Mode 25 26#' Compute the Statistical Median of a Vector 27#' @param x a vector of numeric or ordered values 28#' @param na.rm ignored 29#' @return the median of the vector 30#' @export 31#' @author Matt Owen 32 33Median <- function (x, na.rm=NULL) { 34 v <- ifelse(is.numeric(x), 35 median(x), 36 levels(x)[ceiling(median(as.numeric(x)))] 37 ) 38 if (is.ordered(x)) 39 v <- factor(v, levels(x)) 40 v 41} 42 43#' Create a table, but ensure that the correct 44#' columns exist. In particular, this allows for 45#' entires with zero as a value, which is not 46#' the default for standard tables 47#' @param x a vector 48#' @param levels a vector of levels 49#' @param ... parameters for table 50#' @return a table 51#' @author Matt Owen 52 53table.levels <- function (x, levels, ...) { 54 # if levels are not explicitly set, then 55 # search inside of x 56 if (missing(levels)) { 57 levels <- attr(x, 'levels') 58 table(factor(x, levels=levels), ...) 59 } 60 # otherwise just do the normal thing 61 else { 62 table(factor(x, levels=levels), ...) 63 } 64} 65 66#' Compute central tendancy as approrpriate to data type 67#' @param val a vector of values 68#' @return a mean (if numeric) or a median (if ordered) or mode (otherwise) 69#' @export 70 71avg <- function(val) { 72 if (is.numeric(val)) 73 mean(val) 74 else if (is.ordered(val)) 75 Median(val) 76 else 77 Mode(val) 78} 79 80#' Set new value of a factor variable, checking for existing levels 81#' @param fv factor variable 82#' @param v value 83#' @return a factor variable with a value \code{val} and the same levels 84#' @keywords internal 85setfactor <- function (fv, v) { 86 lev <- levels(fv) 87 if (!v %in% lev) 88 stop("Wrong factor") 89 return(factor(v, levels = lev)) 90} 91 92#' Set new value of a variable as approrpriate to data type 93#' @param val old value 94#' @param newval new value 95#' @return a variable of the same type with a value \code{val} 96#' @keywords internal 97setval <- function(val, newval) { 98 if (is.numeric(val)) 99 newval 100 else if (is.ordered(val)) 101 newval 102 else if (is.logical(val)) 103 newval 104 else { 105 lev <- levels(val) 106 if (!newval %in% lev) 107 stop("Wrong factor", call. = FALSE) 108 return(factor(newval, levels = lev)) 109 } 110} 111 112 113#' Calculate the reduced dataset to be used in \code{\link{setx}} 114#' 115#' #' This method is used internally 116#' 117#' @param dataset Zelig object data, possibly split to deal with \code{by} 118#' argument 119#' @param s list of variables and their tentative \code{setx} values 120#' @param formula a simplified version of the Zelig object formula (typically 121#' with 1 on the lhs) 122#' @param data Zelig object data 123#' @param avg function of data transformations 124#' @return a list of all the model variables either at their central tendancy or 125#' their \code{setx} value 126#' 127#' @keywords internal 128#' @author Christine Choirat and Christopher Gandrud 129#' @export 130 131reduce = function(dataset, s, formula, data, avg = avg) { 132 pred <- try(terms(fit <- lm(formula, data), "predvars"), silent = TRUE) 133 if ("try-error" %in% class(pred)) # exp and weibull 134 pred <- try(terms(fit <- survreg(formula, data), "predvars"), 135 silent = TRUE) 136 137 dataset <- model.frame(fit) 138 139 ldata <- lapply(dataset, avg) 140 if (length(s) > 0) { 141 n <- union(as.character(attr(pred, "predvars"))[-1], names(dataset)) 142 if (is.list(s[[1]])) s <- s[[1]] 143 m <- match(names(s), n) 144 ma <- m[!is.na(m)] 145 if (!all(complete.cases(m))) { 146 w <- paste("Variable '", names(s[is.na(m)]), "' not in data set.\n", 147 sep = "") 148 stop(w, call. = FALSE) 149 } 150 for (i in seq(n[ma])) { 151 ldata[n[ma]][i][[1]] <- setval(dataset[n[ma]][i][[1]], 152 s[n[ma]][i][[1]]) 153 } 154 } 155 return(ldata) 156} 157 158 159 160#' Create QI summary matrix 161#' @param qi quantity of interest in the discrete case 162#' @return a formatted qi 163#' @keywords internal 164#' @author Christine Choirat 165statmat <- function(qi) { 166 if (!is.matrix(qi)) 167 qi <- as.matrix(qi, ncol = 1) 168 m <- t(apply(qi, 2, quantile, c(.5, .025, .975), na.rm = TRUE)) 169 n <- matrix(apply(qi, 2, mean, na.rm = TRUE)) 170 colnames(n) <- "mean" 171 o <- matrix(apply(qi, 2, sd, na.rm = TRUE)) 172 colnames(o) <- "sd" 173 p <- cbind(n, o, m) 174 return(p) 175} 176 177#' Describe Here 178#' @param qi quantity of interest in the discrete case 179#' @param num number of simulations 180#' @return a formatted quantity of interest 181#' @keywords internal 182#' @author Christine Choirat 183statlevel <- function(qi, num) { 184 if (is.matrix(qi)){ 185 #m <- t(apply(qi, 2, table)) / num 186 all.levels <- levels(qi) 187 m <- t(apply(qi, 2, function(x) 188 table(factor(x, levels=all.levels)))) / num 189 } else { 190 m <- table(qi) / num 191 } 192 return(m) 193} 194 195#' Pass Quantities of Interest to Appropriate Summary Function 196#' 197#' @param qi quantity of interest (e.g., estimated value or predicted value) 198#' @param num number of simulations 199#' @return a formatted qi 200#' @keywords internal 201#' @author Christine Choirat 202stat <- function(qi, num) { 203 if (is.null(attr(qi, "levels"))) 204 return(statmat(qi)) 205 else 206 return(statlevel(qi, num)) 207} 208 209#' Generate Formulae that Consider Clustering 210#' 211#' This method is used internally by the "Zelig" Package to interpret 212#' clustering in GEE models. 213#' @param formula a formula object 214#' @param cluster a vector 215#' @return a formula object describing clustering 216cluster.formula <- function (formula, cluster) { 217 # Convert LHS of formula to a string 218 lhs <- deparse(formula[[2]]) 219 cluster.part <- if (is.null(cluster)) 220 # NULL values require 221 sprintf("cluster(1:nrow(%s))", lhs) 222 else 223 # Otherwise we trust user input 224 sprintf("cluster(%s)", cluster) 225 update(formula, paste(". ~ .", cluster.part, sep = " + ")) 226} 227 228 229#' Zelig Copy of plyr::mutate to avoid namespace conflict with dplyr 230#' 231#' @source Hadley Wickham (2011). The Split-Apply-Combine Strategy for Data 232#' Analysis. Journal of Statistical Software, 40(1), 1-29. URL 233#' \url{https://www.jstatsoft.org/v40/i01/}. 234#' @keywords internal 235zelig_mutate <- function (.data, ...) 236{ 237 stopifnot(is.data.frame(.data) || is.list(.data) || is.environment(.data)) 238 cols <- as.list(substitute(list(...))[-1]) 239 cols <- cols[names(cols) != ""] 240 for (col in names(cols)) { 241 .data[[col]] <- eval(cols[[col]], .data, parent.frame()) 242 } 243 .data 244} 245 246#' Convenience function for setrange and setrange1 247#' 248#' @param x data passed to setrange or setrange1 249#' @keywords internal 250 251expand_grid_setrange <- function(x) { 252 # m <- expand.grid(x) 253 set_lengths <- unlist(lapply(x, length)) 254 unique_set_lengths <- unique(as.vector(set_lengths)) 255 256 m <- data.frame() 257 for (i in unique_set_lengths) { 258 temp_df <- data.frame(row.names = 1:i) 259 for (u in 1:length(x)) { 260 if (length(x[[u]]) == i) { 261 temp_df <- cbind(temp_df, x[[u]]) 262 names(temp_df)[ncol(temp_df)] <- names(x)[u] 263 } 264 } 265 if (nrow(m) == 0) m <- temp_df 266 else m <- merge(m, temp_df) 267 } 268 if (nrow(m) == 1) 269 warning('Only one fitted observation provided to setrange.\nConsider using setx instead.', 270 call. = FALSE) 271 return(m) 272} 273 274#' Bundle Multiply Imputed Data Sets into an Object for Zelig 275#' 276#' This object prepares multiply imputed data sets so they can be used by 277#' \code{zelig}. 278#' @note This function creates a list of \code{data.frame} objects, which 279#' resembles the storage of imputed data sets in the \code{amelia} object. 280#' @param ... a set of \code{data.frame}'s or a single list of \code{data.frame}'s 281#' @return an \code{mi} object composed of a list of data frames. 282#' 283#' @author Matt Owen, James Honaker, and Christopher Gandrud 284#' 285#' @examples 286#' # create datasets 287#' n <- 100 288#' x1 <- runif(n) 289#' x2 <- runif(n) 290#' y <- rnorm(n) 291#' data.1 <- data.frame(y = y, x = x1) 292#' data.2 <- data.frame(y = y, x = x2) 293#' 294#' # merge datasets into one object as if imputed datasets 295#' 296#' mi.out <- to_zelig_mi(data.1, data.2) 297#' 298#' # pass object in place of data argument 299#' z.out <- zelig(y ~ x, model = "ls", data = mi.out) 300#' @export 301 302to_zelig_mi <- function (...) { 303 304 # Get arguments as list 305 imputations <- list(...) 306 307 # If user passes a list of data.frames rather than several data.frames as separate arguments 308 if((class(imputations[[1]]) == 'list') & (length(imputations) == 1)){ 309 imputations = imputations[[1]] 310 } 311 312 # Labelling 313 names(imputations) <- paste0("imp", 1:length(imputations)) 314 315 # Ensure that everything is a data.frame 316 for (k in length(imputations):1) { 317 if (!is.data.frame(imputations[[k]])){ 318 imputations[[k]] <- NULL 319 warning("Item ", k, " of the provided objects is not a data.frame and will be ignored.\n") 320 } 321 } 322 323 if(length(imputations) < 1){ 324 stop("The resulting object contains no data.frames, and as such is not a valid multiple imputation object.", 325 call. = FALSE) 326 } 327 if(length(imputations) < 2){ 328 stop("The resulting object contains only one data.frame, and as such is not a valid multiple imputation object.", 329 call. = FALSE) 330 } 331 class(imputations) <-c("mi", "list") 332 333 return(imputations) 334} 335 336#' Enables backwards compatability for preparing non-amelia imputed data sets 337#' for \code{zelig}. 338#' 339#' See \code{\link{to_zelig_mi}} 340#' 341#' @param ... a set of \code{data.frame}'s 342#' @return an \code{mi} object composed of a list of data frames. 343mi <- to_zelig_mi 344 345#' Conduct variable transformations called inside a \code{zelig} call 346#' 347#' @param formula model formulae 348#' @param data data frame used in \code{formula} 349#' @param FUN character string of the transformation function. Currently 350#' supports \code{factor} and \code{log}. 351#' @param check logical whether to just check if a formula contains an 352#' internally called transformation and return \code{TRUE} or \code{FALSE} 353#' @param f_out logical whether to return the converted formula 354#' @param d_out logical whether to return the converted data frame. Note: 355#' \code{f_out} must be missing 356#' 357#' @author Christopher Gandrud 358#' @keywords internal 359 360transformer <- function(formula, data, FUN = 'log', check, f_out, d_out) { 361 362 if (!missing(data)) { 363 if (is.data.frame(data)) 364 is_df <- TRUE 365 else if (!is.data.frame(data) & is.list(data)) 366 is_df <- FALSE 367 else 368 stop('data must be either a data.frame or a list', call. = FALSE) 369 } 370 371 if (FUN == 'as.factor') FUN_temp <- 'as\\.factor' 372 else FUN_temp <- FUN 373 FUN_str <- sprintf('%s.*\\(', FUN_temp) 374 375 f <- as.character(formula)[3] 376 f_split <- unlist(strsplit(f, split = '\\+')) 377 to_transform <- grep(pattern = FUN_str, f_split) 378 379 if (!missing(check)) { 380 if (length(to_transform) > 0) return(TRUE) 381 else return(FALSE) 382 } 383 384 if (length(to_transform) > 0) { 385 to_transform_raw <- trimws(f_split[to_transform]) 386 if (FUN == 'factor') 387 to_transform_raw <- gsub('^as\\.', '', to_transform_raw) 388 to_transform_plain_args <- gsub(FUN_str, '', to_transform_raw) 389 to_transform_plain <- gsub(',\\(.*)', '', to_transform_plain_args) 390 to_transform_plain <- gsub('\\)', '', to_transform_plain) 391 to_transform_plain <- trimws(gsub(',.*', '', to_transform_plain)) 392 393 if (is_df) 394 not_in_data <- !all(to_transform_plain %in% names(data)) 395 else if (!isTRUE(is_df)) 396 not_in_data <- !all(to_transform_plain %in% names(data[[1]])) 397 if (not_in_data) stop('Unable to find variable to transform.') 398 399 if (!missing(f_out)) { 400 f_split[to_transform] <- to_transform_plain 401 rhs <- paste(f_split, collapse = ' + ') 402 lhs <- gsub('\\(\\)', '', formula[2]) 403 f_new <- paste(lhs, '~', rhs) 404 f_out <- as.Formula(f_new) 405 return(f_out) 406 } 407 else if (!missing(d_out)) { 408 409 transformer_fun <- trimws(gsub('\\(.*', '', to_transform_raw)) 410 411 transformer_args_str <- gsub('\\)', '', to_transform_plain_args) 412 transformer_args_list <- list() 413 for (i in seq_along(transformer_args_str)) { 414 args_temp <- unlist(strsplit(gsub(' ', '' , 415 transformer_args_str[i]), ',')) 416 if (is_df) 417 args_temp[1] <- sprintf('data[, "%s"]', args_temp[1]) 418 else if (!isTRUE(is_df)) 419 args_temp[1] <- sprintf('data[[h]][, "%s"]', args_temp[1]) 420 arg_names <- gsub('\\=.*', '', args_temp) 421 arg_names[1] <- 'x' 422 args_temp <- gsub('.*\\=', '', args_temp) 423 424 args_temp_list <- list() 425 if (is_df) { 426 for (u in seq_along(args_temp)) 427 args_temp_list[[u]] <- eval(parse(text = args_temp[u])) 428 } 429 else if (!isTRUE(is_df)) { 430 for (h in seq_along(data)) { 431 temp_list <- list() 432 for (u in seq_along(args_temp)) { 433 temp_list[[u]] <- eval(parse(text = args_temp[u])) 434 names(temp_list)[u] <- arg_names[u] 435 } 436 args_temp_list[[h]] <- temp_list 437 } 438 } 439 if (is_df) { 440 names(args_temp_list) <- arg_names 441 data[, to_transform_plain[i]] <- do.call( 442 what = transformer_fun[i], 443 args = args_temp_list) 444 } 445 else if (!isTRUE(is_df)) { 446 for (j in seq_along(data)) { 447 data[[j]][, to_transform_plain[i]] <- do.call( 448 what = transformer_fun[i], 449 args = args_temp_list[[j]]) 450 } 451 } 452 } 453 return(data) 454 } 455 } 456 else if (length(to_transform) == 0) { 457 if (!missing(f_out)) return(formula) 458 else if (d_out) return(data) 459 } 460} 461 462 463#' Remove package names from fitted model object calls. 464#' 465#' Enables \code{\link{from_zelig_model}} output to work with stargazer. 466#' @param x a fitted model object result 467#' @keywords internal 468 469strip_package_name <- function(x) { 470 if ("vglm" %in% class(x)) # maybe generalise to all s4? 471 call_temp <- gsub('^.*(?=(::))', '', x@call[1], perl = TRUE) 472 else 473 call_temp <- gsub('^.*(?=(::))', '', x$call[1], perl = TRUE) 474 call_temp <- gsub('::', '', call_temp, perl = TRUE) 475 if ("vglm" %in% class(x)) 476 x@call[1] <- as.call(list(as.symbol(call_temp))) 477 else 478 x$call[1] <- as.call(list(as.symbol(call_temp))) 479 return(x) 480} 481 482#' Extract p-values from a fitted model object 483#' @param x a fitted Zelig object 484#' @keywords internal 485 486p_pull <- function(x) { 487 if ("vglm" %in% class(x)) { # maybe generalise to all s4? 488 p_values <- summary(x)@coef3[, 'Pr(>|z|)'] 489 } 490 else { 491 p_values <- summary(x)$coefficients 492 if ('Pr(>|t|)' %in% colnames(p_values)) { 493 p_values <- p_values[, 'Pr(>|t|)'] 494 } else { 495 p_values <- p_values[, 'Pr(>|z|)'] 496 } 497 } 498 499 return(p_values) 500} 501 502#' Extract standard errors from a fitted model object 503#' @param x a fitted Zelig object 504#' @keywords internal 505 506se_pull <- function(x) { 507 if ("vglm" %in% class(x)) # maybe generalise to all s4? 508 se <- summary(x)@coef3[, "Std. Error"] 509 else 510 se <- summary(x)$coefficients[, "Std. Error"] 511 return(se) 512} 513 514#' Drop intercept columns or values from a data frame or named vector, 515#' respectively 516#' 517#' @param x a data frame or named vector 518#' @keywords internal 519 520rm_intercept <- function(x) { 521 intercept_names <- c('(Intercept)', 'X.Intercept.', '(Intercept).*') 522 names_x <- names(x) 523 if (any(intercept_names %in% names(x))) { 524 keep <- !(names(x) %in% intercept_names) 525 if (is.data.frame(x)) 526 x <- data.frame(x[, names_x[keep]]) 527 else if (is.vector(x)) 528 x <- x[keep] 529 names(x) <- names_x[keep] 530 } 531 return(x) 532} 533 534 535#' Combines estimated coefficients and associated statistics 536#' from models estimated with multiply imputed data sets or bootstrapped 537#' 538#' @param obj a zelig object with an estimated model 539#' @param out_type either \code{"matrix"} or \code{"list"} specifying 540#' whether the results should be returned as a matrix or a list. 541#' @param bagging logical whether or not to bag the bootstrapped coefficients 542#' @param messages logical whether or not to return messages for what is being 543#' returned 544#' 545#' @return If the model uses multiply imputed or bootstrapped data then a 546#' matrix (default) or list of combined coefficients (\code{coef}), standard 547#' errors (\code{se}), z values (\code{zvalue}), p-values (\code{p}) is 548#' returned. Rubin's Rules are used to combine output from multiply imputed 549#' data. An error is returned if no imputations were included or there wasn't 550#' bootstrapping. Please use \code{get_coef}, \code{get_se}, and 551#' \code{get_pvalue} methods instead in cases where there are no imputations or 552#' bootstrap. 553#' 554#' @examples 555#' set.seed(123) 556#' 557#' ## Multiple imputation example 558#' # Create fake imputed data 559#' n <- 100 560#' x1 <- runif(n) 561#' x2 <- runif(n) 562#' y <- rnorm(n) 563#' data.1 <- data.frame(y = y, x = x1) 564#' data.2 <- data.frame(y = y, x = x2) 565#' 566#' # Estimate model 567#' mi.out <- to_zelig_mi(data.1, data.2) 568#' z.out.mi <- zelig(y ~ x, model = "ls", data = mi.out) 569#' 570#' # Combine and extract coefficients and standard errors 571#' combine_coef_se(z.out.mi) 572#' 573#' ## Bootstrap example 574#' z.out.boot <- zelig(y ~ x, model = "ls", data = data.1, bootstrap = 20) 575#' combine_coef_se(z.out.boot) 576#' 577#' @author Christopher Gandrud and James Honaker 578#' @source Partially based on \code{\link{mi.meld}} from Amelia. 579#' 580#' @export 581 582combine_coef_se <- function(obj, out_type = 'matrix', bagging = FALSE, 583 messages = TRUE) 584{ 585 is_zelig(obj) 586 is_uninitializedField(obj$zelig.out) 587 if (!(out_type %in% c('matrix', 'list'))) 588 stop('out_type must be either "matrix" or "list"', call. = FALSE) 589 590 if (obj$mi || obj$bootstrap) { 591 coeflist <- obj$get_coef() 592 vcovlist <- obj$get_vcov() 593 coef_names <- names(coeflist[[1]]) 594 595 am.m <- length(coeflist) 596 if (obj$bootstrap & !obj$mi) am.m <- am.m - 1 597 am.k <- length(coeflist[[1]]) 598 if (obj$bootstrap & !obj$mi) 599 q <- matrix(unlist(coeflist[-(am.m + 1)]), nrow = am.m, 600 ncol = am.k, byrow = TRUE) 601 else if (obj$mi) { 602 q <- matrix(unlist(coeflist), nrow = am.m, ncol = am.k, 603 byrow = TRUE) 604 se <- matrix(NA, nrow = am.m, ncol = am.k) 605 for(i in 1:am.m){ 606 se[i, ] <- sqrt(diag(vcovlist[[i]])) 607 } 608 } 609 ones <- matrix(1, nrow = 1, ncol = am.m) 610 comb_q <- (ones %*% q)/am.m 611 if (obj$mi) ave.se2 <- (ones %*% (se^2)) / am.m 612 diff <- q - matrix(1, nrow = am.m, ncol = 1) %*% comb_q 613 sq2 <- (ones %*% (diff^2))/(am.m - 1) 614 if (obj$mi) { 615 if (messages) message('Combining imputations. . .') 616 comb_se <- sqrt(ave.se2 + sq2 * (1 + 1/am.m)) 617 coef <- as.vector(comb_q) 618 se <- as.vector(comb_se) 619 } 620 621 else if (obj$bootstrap & !obj$mi) { 622 if (messages) message('Combining bootstraps . . .') 623 comb_se <- sqrt(sq2 * (1 + 1/am.m)) 624 if (bagging) { 625 coef <- as.vector(comb_q) 626 } else { 627 coef <- coeflist[[am.m + 1]] 628 } 629 se <- as.vector(comb_se) 630 } 631 632 zvalue <- coef / se 633 pr_z <- 2 * (1 - pnorm(abs(zvalue))) 634 635 if (out_type == 'matrix') { 636 out <- cbind(coef, se, zvalue, pr_z) 637 colnames(out) <- c("Estimate", "Std.Error", "z value", "Pr(>|z|)") 638 rownames(out) <- coef_names 639 } 640 else if (out_type == 'list') { 641 out <- list(coef = coef, se = se, zvalue = zvalue, p = pr_z) 642 for (i in seq(out)) names(out[[i]]) <- coef_names 643 } 644 return(out) 645 } 646 else if (!(obj$mi || obj$bootstrap)) { 647 message('No multiply imputed or bootstrapped estimates found.\nReturning untransformed list of coefficients and standard errors.') 648 out <- list(coef = coef(obj), 649 se = get_se(obj), 650 pvalue = get_pvalue(obj) 651 ) 652 653 return(out) 654 } 655} 656 657#' Find vcov for GEE models 658#' 659#' @param obj a \code{geeglm} class object. 660 661vcov_gee <- function(obj) { 662 if (!("geeglm" %in% class(obj))) 663 stop('Not a geeglm class object', call. = FALSE) 664 out <- obj$geese$vbeta 665 return(out) 666} 667 668#' Find vcov for quantile regression models 669#' 670#' @param obj a \code{rq} class object. 671 672vcov_rq <- function(obj) { 673 if (!("rq" %in% class(obj))) 674 stop('Not an rq class object', call. = FALSE) 675 out <- summary(obj, cov = TRUE)$cov 676 return(out) 677} 678 679#' Find odds ratios for coefficients and standard errors 680#' for glm.summary class objects 681#' 682#' @param obj a \code{glm.summary} class object 683#' @param label_mod_coef character string for how to modify the coefficient 684#' label. 685#' @param label_mod_se character string for how to modify the standard error 686#' label. 687 688or_summary <- function(obj, label_mod_coef = "(OR)", 689 label_mod_se = "(OR)"){ 690 if (class(obj) != "summary.glm") 691 stop("obj must be of summary.glm class.", 692 call. = FALSE) 693 694 obj$coefficients[, 1] <- exp(obj$coefficients[, 1]) 695 696 var_diag = diag(vcov(obj)) 697 obj$coefficients[, 2] <- sqrt(obj$coefficients[, 1] ^ 2 * var_diag) 698 699 colnames(obj$coefficients)[c(1, 2)] <- paste( 700 colnames(obj$coefficients)[c(1, 2)], 701 c(label_mod_coef, label_mod_se)) 702 return(obj) 703} 704