Merge pull request #20 from cyk0315/master

add family poisson, quasipoisson
jinseob2kim · Feb 29, 2024 · 61fe4c6 · 61fe4c6
2 parents ed24d3f + eab97e5
commit 61fe4c6
Show file tree

Hide file tree

Showing 8 changed files with 132 additions and 29 deletions.
diff --git 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 = "[email protected]", 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
@@ -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
@@ -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,21 +61,18 @@ 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)
   )
   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
@@ -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
diff --git a/man/TableSubgroupGLM.Rd b/man/TableSubgroupGLM.Rd
diff --git a/man/TableSubgroupMultiGLM.Rd b/man/TableSubgroupMultiGLM.Rd
diff --git a/man/mort.Rd b/man/mort.Rd