diff --git a/DESCRIPTION b/DESCRIPTION index 1b7d341..bba4fad 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: jstable Title: Create Tables from Different Types of Regression -Version: 1.1.6 -Date: 2024-02-15 +Version: 1.1.7 +Date: 2024-02-29 Authors@R: c(person("Jinseob", "Kim", email = "jinseob2kim@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-9403-605X")), person("Zarathu", role = c("cph", "fnd")), person("Yoonkyoung","Jeon", role = c("aut")) @@ -19,3 +19,4 @@ Suggests: knitr, rmarkdown VignetteBuilder: knitr +LazyData: true diff --git a/NEWS.md b/NEWS.md index 935c85f..6c1b49c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# jstable 1.1.7 + +* Add family 'poisson', 'quasipoisson' in `glmshow.display` and `TableSubgroupMultiGLM` +* Add data `mort` + # jstable 1.1.6 * Bugfix `TableSubgroupGLM`: thanks for `weisx2022` diff --git a/R/forestglm.R b/R/forestglm.R index 50c6f47..7776228 100644 --- a/R/forestglm.R +++ b/R/forestglm.R @@ -4,7 +4,7 @@ #' @param var_subgroup 1 sub-group variable for analysis, Default: NULL #' @param var_cov Variables for additional adjust, Default: NULL #' @param data Data or svydesign in survey package. -#' @param family family, "gaussian" or "binomial" +#' @param family family, "gaussian" or "binomial" or 'poisson' or 'quasipoisson' #' @param decimal.estimate Decimal for estimate, Default: 2 #' @param decimal.percent Decimal for percent, Default: 1 #' @param decimal.pvalue Decimal for pvalue, Default: 3 @@ -31,14 +31,13 @@ #' \code{\link[purrr]{safely}},\code{\link[purrr]{map}},\code{\link[purrr]{map2}} #' \code{\link[stats]{glm}} #' \code{\link[survey]{svyglm}} -#' \code{\link[stats]{confint}} #' @rdname TableSubgroupGLM #' @export #' @importFrom purrr possibly map_dbl map map2 #' @importFrom dplyr group_split select filter mutate bind_cols #' @importFrom magrittr %>% #' @importFrom survey svyglm -#' @importFrom stats glm confint coefficients anova gaussian quasibinomial +#' @importFrom stats glm coefficients anova gaussian quasibinomial #' @importFrom utils tail TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, family = "binomial", decimal.estimate = 2, decimal.percent = 1, decimal.pvalue = 3) { @@ -62,14 +61,9 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, }, NA) possible_glm <- purrr::possibly(stats::glm, NA) possible_svyglm <- purrr::possibly(survey::svyglm, NA) - possible_confint <- purrr::possibly(stats::confint, NA) possible_modely <- purrr::possibly(function(x) { purrr::map_dbl(x, .[["y"]], 1) }, NA) - possible_rowone <- purrr::possibly(function(x) { - x[2, ] - }, NA) - xlabel <- setdiff(as.character(formula)[[3]], "+")[1] ncoef <- ifelse(any(class(data) == "survey.design"), ifelse(length(levels(data$variables[[xlabel]])) <= 2, 1, length(levels(data$variables[[xlabel]])) - 1), ifelse(length(levels(data[[xlabel]])) <= 2, 1, length(levels(data[[xlabel]])) - 1) @@ -77,6 +71,8 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, var_cov <- setdiff(var_cov, c(as.character(formula[[3]]), var_subgroup)) family.svyglm <- gaussian() if (family == "binomial") family.svyglm <- quasibinomial() + if (family == "poisson") family.svyglm <- poisson() + if (family == "quasipoisson") family.svyglm <- quasipoisson() if (is.null(var_subgroup)) { if (!is.null(var_cov)) { @@ -91,12 +87,13 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, # if (!is.null(model$xlevels) & length(model$xlevels[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") } - + cc<-summary(model)$coefficients Point.Estimate <- round(stats::coef(model), decimal.estimate)[2] - CI <- round(stats::confint(model)[2, ], decimal.estimate) - if (family == "binomial") { + + CI<-round(c(cc[2, 1] - qnorm(0.975) * cc[2, 2],cc[2, 1] + qnorm(0.975) * cc[2, 2]),decimal.estimate) + if (family %in% c("binomial",'poisson','quasipoisson')) { Point.Estimate <- round(exp(stats::coef(model)), decimal.estimate)[2] - CI <- round(exp(stats::confint(model)[2, ]), decimal.estimate) + CI<-round(exp(c(cc[2, 1] - qnorm(0.975) * cc[2, 2],cc[2, 1] + qnorm(0.975) * cc[2, 2])),decimal.estimate) } @@ -117,6 +114,9 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, if (family == "binomial") { names(out)[4] <- "OR" } + if (family %in% c('poisson','quasipoisson')) { + names(out)[4] <- "RR" + } return(out) } else if (length(var_subgroup) >= 2 | any(grepl(var_subgroup, formula))) { @@ -139,8 +139,13 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, data.design <- data if (family == "binomial") { model.int <- survey::svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep = ""), deparse(formula))), design = data.design, family = quasibinomial()) - } else { + }else if(family=='gaussian'){ model.int <- survey::svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep = ""), deparse(formula))), design = data.design, family = gaussian()) + }else if(family=='poisson'){ + model.int <- survey::svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep = ""), deparse(formula))), design = data.design, family = poisson()) + }else{ + model.int <- survey::svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep = ""), deparse(formula))), design = data.design, family = quasipoisson()) + } if (sum(grepl(":", names(coef(model.int)))) > 1) { pv_anova <- anova(model.int, method = "Wald") @@ -154,7 +159,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, purrr::map(~ possible_glm(formula, data = ., x = T, family = family)) -> model data %>% subset(!is.na(get(var_subgroup))) %>% - select(var_subgroup) %>% + select(all_of(var_subgroup)) %>% table() %>% names() -> label_val xlev <- stats::glm(formula, data = data)$xlevels @@ -178,13 +183,16 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, Estimate <- model %>% purrr::map("coefficients", .default = NA) %>% purrr::map_dbl(2, .default = NA) + CI0 <- model %>% - purrr::map(possible_confint) %>% - purrr::map(possible_rowone) %>% + purrr::map(function(model){ + cc0<-summary(model)$coefficients + c(cc0[2, 1] - qnorm(0.975) * cc0[2, 2],cc0[2, 1] + qnorm(0.975) * cc0[2, 2]) + }) %>% Reduce(rbind, .) Point.Estimate <- round(Estimate, decimal.estimate) CI <- round(CI0, decimal.estimate) - if (family == "binomial") { + if (family %in% c("binomial",'poisson','quasipoisson')) { Point.Estimate <- round(exp(Estimate), decimal.estimate) CI <- round(exp(CI0), decimal.estimate) } @@ -199,6 +207,9 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, if (family == "binomial") { names(out)[4] <- "OR" } + if (family %in% c("poisson",'quasipoisson')) { + names(out)[4] <- "RR" + } return(rbind(c(var_subgroup, rep(NA, ncol(out) - 2), ifelse(pv_int >= 0.001, pv_int, "<0.001")), out)) } @@ -212,7 +223,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, #' @param var_subgroups Multiple sub-group variables for analysis, Default: NULL #' @param var_cov Variables for additional adjust, Default: NULL #' @param data Data or svydesign in survey package. -#' @param family family, "gaussian" or "binomial" +#' @param family family, "gaussian" or "binomial" or 'poisson' or 'quasipoisson' #' @param decimal.estimate Decimal for estimate, Default: 2 #' @param decimal.percent Decimal for percent, Default: 1 #' @param decimal.pvalue Decimal for pvalue, Default: 3 diff --git a/R/glmshow.R b/R/glmshow.R index e26b1ff..781e5e2 100644 --- a/R/glmshow.R +++ b/R/glmshow.R @@ -4,6 +4,7 @@ #' @return coefficient table with NA #' @details DETAILS #' @examples +#' #' coefNA(glm(mpg ~ wt + qsec, data = mtcars)) #' @rdname coefNA #' @export @@ -39,7 +40,9 @@ glmshow.display <- function(glm.object, decimal = 2) { xs <- attr(model$terms, "term.labels") y <- names(model$model)[1] - gaussianT <- ifelse(length(grep("gaussian", model$family)) == 1, T, F) + family<-ifelse(length(grep("gaussian", model$family)) == 1, 1, ifelse(length(grep("binomial", model$family)) >= 1,2,3)) + + data <- model$data ## table @@ -48,14 +51,18 @@ glmshow.display <- function(glm.object, decimal = 2) { } else if (length(xs) == 1) { uni <- data.frame(coefNA(glm.object))[-1, ] rn.uni <- lapply(list(uni), rownames) - if (gaussianT) { + if (family==1) { summ <- paste(round(uni[, 1], decimal), " (", round(uni[, 1] - 1.96 * uni[, 2], decimal), ",", round(uni[, 1] + 1.96 * uni[, 2], decimal), ")", sep = "") uni.res <- matrix(cbind(summ, ifelse(uni[, 4] <= 0.001, "< 0.001", as.character(round(uni[, 4], decimal + 1)))), nrow = nrow(uni)) colnames(uni.res) <- c(paste("Coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value") } else { summ <- paste(round(exp(uni[, 1]), decimal), " (", round(exp(uni[, 1] - 1.96 * uni[, 2]), decimal), ",", round(exp(uni[, 1] + 1.96 * uni[, 2]), decimal), ")", sep = "") uni.res <- matrix(cbind(summ, ifelse(uni[, 4] <= 0.001, "< 0.001", as.character(round(uni[, 4], decimal + 1)))), nrow = nrow(uni)) - colnames(uni.res) <- c(paste("OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value") + if(family==2){ + colnames(uni.res) <- c(paste("OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value") + }else{ + colnames(uni.res) <- c(paste("RR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value") + } } rownames(uni.res) <- rownames(uni) res <- uni.res @@ -68,7 +75,7 @@ glmshow.display <- function(glm.object, decimal = 2) { }) rn.uni <- lapply(uni, rownames) uni <- Reduce(rbind, uni) - if (gaussianT) { + if (family==1) { summ <- paste(round(uni[, 1], decimal), " (", round(uni[, 1] - 1.96 * uni[, 2], decimal), ",", round(uni[, 1] + 1.96 * uni[, 2], decimal), ")", sep = "") uni.res <- t(rbind(summ, ifelse(uni[, 4] <= 0.001, "< 0.001", as.character(round(uni[, 4], decimal + 1))))) colnames(uni.res) <- c(paste("crude coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value") @@ -78,14 +85,16 @@ glmshow.display <- function(glm.object, decimal = 2) { mul.res <- t(rbind(mul.summ, ifelse(mul[, 4] <= 0.001, "< 0.001", as.character(round(mul[, 4], decimal + 1))))) colnames(mul.res) <- c(paste("adj. coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value") } else { + k<-ifelse(family==2,'OR','RR') + summ <- paste(round(exp(uni[, 1]), decimal), " (", round(exp(uni[, 1] - 1.96 * uni[, 2]), decimal), ",", round(exp(uni[, 1] + 1.96 * uni[, 2]), decimal), ")", sep = "") uni.res <- t(rbind(summ, ifelse(uni[, 4] <= 0.001, "< 0.001", as.character(round(uni[, 4], decimal + 1))))) - colnames(uni.res) <- c(paste("crude OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value") + colnames(uni.res) <- c(paste("crude ",k,".(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value") rownames(uni.res) <- rownames(uni) mul <- coefNA(model)[-1, ] mul.summ <- paste(round(exp(mul[, 1]), decimal), " (", round(exp(mul[, 1] - 1.96 * mul[, 2]), decimal), ",", round(exp(mul[, 1] + 1.96 * mul[, 2]), decimal), ")", sep = "") mul.res <- t(rbind(mul.summ, ifelse(mul[, 4] <= 0.001, "< 0.001", as.character(round(mul[, 4], decimal + 1))))) - colnames(mul.res) <- c(paste("adj. OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value") + colnames(mul.res) <- c(paste("adj. ",k,".(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value") } res <- cbind(uni.res[rownames(uni.res) %in% rownames(mul.res), ], mul.res) @@ -129,7 +138,7 @@ glmshow.display <- function(glm.object, decimal = 2) { outcome.name <- y - if (gaussianT) { + if (family==1) { first.line <- paste("Linear regression predicting ", outcome.name, sep = "", "\n") last.lines <- paste("No. of observations = ", length(model$y), "\n", "R-squared = ", round(cor(model$y, predict(model))^2, decimal + 2), "\n", @@ -152,3 +161,35 @@ glmshow.display <- function(glm.object, decimal = 2) { class(results) <- c("display", "list") return(results) } + +#' @title DATASET_TITLE +#' @description DATASET_DESCRIPTION +#' @format A data frame with 17562 rows and 24 variables: +#' \describe{ +#' \item{\code{ccode}}{integer COLUMN_DESCRIPTION} +#' \item{\code{cname}}{character COLUMN_DESCRIPTION} +#' \item{\code{yy}}{integer COLUMN_DESCRIPTION} +#' \item{\code{mm}}{integer COLUMN_DESCRIPTION} +#' \item{\code{dd}}{integer COLUMN_DESCRIPTION} +#' \item{\code{date}}{character COLUMN_DESCRIPTION} +#' \item{\code{nonacc}}{integer COLUMN_DESCRIPTION} +#' \item{\code{cardio}}{integer COLUMN_DESCRIPTION} +#' \item{\code{respir}}{integer COLUMN_DESCRIPTION} +#' \item{\code{influenza}}{integer COLUMN_DESCRIPTION} +#' \item{\code{meanpm10}}{double COLUMN_DESCRIPTION} +#' \item{\code{meanso2}}{double COLUMN_DESCRIPTION} +#' \item{\code{meanno2}}{double COLUMN_DESCRIPTION} +#' \item{\code{meanco}}{double COLUMN_DESCRIPTION} +#' \item{\code{maxco}}{double COLUMN_DESCRIPTION} +#' \item{\code{maxo3}}{double COLUMN_DESCRIPTION} +#' \item{\code{meantemp}}{double COLUMN_DESCRIPTION} +#' \item{\code{maxtemp}}{double COLUMN_DESCRIPTION} +#' \item{\code{mintemp}}{double COLUMN_DESCRIPTION} +#' \item{\code{meanhumi}}{double COLUMN_DESCRIPTION} +#' \item{\code{meanpress}}{double COLUMN_DESCRIPTION} +#' \item{\code{season}}{integer COLUMN_DESCRIPTION} +#' \item{\code{dow}}{integer COLUMN_DESCRIPTION} +#' \item{\code{sn}}{integer COLUMN_DESCRIPTION} +#'} +#' @details DETAILS +"mort" diff --git a/data/mort.rda b/data/mort.rda new file mode 100644 index 0000000..db15e86 Binary files /dev/null and b/data/mort.rda differ diff --git a/man/TableSubgroupGLM.Rd b/man/TableSubgroupGLM.Rd index 191ccfb..f5be0d9 100644 --- a/man/TableSubgroupGLM.Rd +++ b/man/TableSubgroupGLM.Rd @@ -24,7 +24,7 @@ TableSubgroupGLM( \item{data}{Data or svydesign in survey package.} -\item{family}{family, "gaussian" or "binomial"} +\item{family}{family, "gaussian" or "binomial" or 'poisson' or 'quasipoisson'} \item{decimal.estimate}{Decimal for estimate, Default: 2} diff --git a/man/TableSubgroupMultiGLM.Rd b/man/TableSubgroupMultiGLM.Rd index 0e324fd..0a3b1ee 100644 --- a/man/TableSubgroupMultiGLM.Rd +++ b/man/TableSubgroupMultiGLM.Rd @@ -25,7 +25,7 @@ TableSubgroupMultiGLM( \item{data}{Data or svydesign in survey package.} -\item{family}{family, "gaussian" or "binomial"} +\item{family}{family, "gaussian" or "binomial" or 'poisson' or 'quasipoisson'} \item{decimal.estimate}{Decimal for estimate, Default: 2} diff --git a/man/mort.Rd b/man/mort.Rd new file mode 100644 index 0000000..06e46a6 --- /dev/null +++ b/man/mort.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/glmshow.R +\docType{data} +\name{mort} +\alias{mort} +\title{DATASET_TITLE} +\format{ +A data frame with 17562 rows and 24 variables: +\describe{ + \item{\code{ccode}}{integer COLUMN_DESCRIPTION} + \item{\code{cname}}{character COLUMN_DESCRIPTION} + \item{\code{yy}}{integer COLUMN_DESCRIPTION} + \item{\code{mm}}{integer COLUMN_DESCRIPTION} + \item{\code{dd}}{integer COLUMN_DESCRIPTION} + \item{\code{date}}{character COLUMN_DESCRIPTION} + \item{\code{nonacc}}{integer COLUMN_DESCRIPTION} + \item{\code{cardio}}{integer COLUMN_DESCRIPTION} + \item{\code{respir}}{integer COLUMN_DESCRIPTION} + \item{\code{influenza}}{integer COLUMN_DESCRIPTION} + \item{\code{meanpm10}}{double COLUMN_DESCRIPTION} + \item{\code{meanso2}}{double COLUMN_DESCRIPTION} + \item{\code{meanno2}}{double COLUMN_DESCRIPTION} + \item{\code{meanco}}{double COLUMN_DESCRIPTION} + \item{\code{maxco}}{double COLUMN_DESCRIPTION} + \item{\code{maxo3}}{double COLUMN_DESCRIPTION} + \item{\code{meantemp}}{double COLUMN_DESCRIPTION} + \item{\code{maxtemp}}{double COLUMN_DESCRIPTION} + \item{\code{mintemp}}{double COLUMN_DESCRIPTION} + \item{\code{meanhumi}}{double COLUMN_DESCRIPTION} + \item{\code{meanpress}}{double COLUMN_DESCRIPTION} + \item{\code{season}}{integer COLUMN_DESCRIPTION} + \item{\code{dow}}{integer COLUMN_DESCRIPTION} + \item{\code{sn}}{integer COLUMN_DESCRIPTION} +} +} +\usage{ +mort +} +\description{ +DATASET_DESCRIPTION +} +\details{ +DETAILS +} +\keyword{datasets}