From f15ca867f49a65e588abaa70893ca9d0110d6753 Mon Sep 17 00:00:00 2001 From: jaehun shon <2021122006@yonsei.ac.kr> Date: Wed, 14 Aug 2024 10:51:54 +0900 Subject: [PATCH 1/3] list to unlist --- R/forestglm.R | 112 +++++++++++++++++++++++++------------------------- 1 file changed, 57 insertions(+), 55 deletions(-) diff --git a/R/forestglm.R b/R/forestglm.R index 1d26658..70caca2 100644 --- a/R/forestglm.R +++ b/R/forestglm.R @@ -42,7 +42,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, family = "binomial", decimal.estimate = 2, decimal.percent = 1, decimal.pvalue = 3) { . <- NULL - + ### 경고문 ### if (length(formula[[3]]) > 1) stop("Formula must contains only 1 independent variable") if (any(class(data) == "survey.design" & !is.null(var_subgroup))) { @@ -50,7 +50,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, } else if (any(class(data) == "data.frame" & !is.null(var_subgroup))) { if (is.numeric(data[[var_subgroup]])) stop("var_subgroup must categorical.") } - + ## functions with error possible_table <- purrr::possibly(table, NA) possible_prop.table <- purrr::possibly(function(x) { @@ -64,26 +64,26 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, possible_modely <- purrr::possibly(function(x) { purrr::map_dbl(x, .[["y"]], 1) }, 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) + 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)) { ### subgroup 지정 안 한 경우 ### - + # 공변량 있는 경우 formula 변경 if (!is.null(var_cov)) { formula <- as.formula(paste0(deparse(formula), " + ", paste(var_cov, collapse = "+"))) } - + if (any(class(data) == "survey.design")) { model <- survey::svyglm(formula, design = data, x = T, family = family.svyglm) # if (!is.null(model$xlevels) & length(model$xlevels[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") @@ -91,38 +91,38 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, model <- stats::glm(formula, data = data, x = T, family = family) # if (!is.null(model$xlevels) & length(model$xlevels[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") } - + xlev <- ifelse(length(model$xlevels) == 0, NA, model$xlevels) - + # cc, PE, CI, PV 구하기 cc <- summary(model)$coefficients Point.Estimate <- round(stats::coef(model), decimal.estimate)[2:(1 + ncoef)] CI <- round(matrix(c(cc[2:(1 + ncoef), 1] - qnorm(0.975) * cc[2:(1 + ncoef), 2], cc[2:(1 + ncoef), 1] + qnorm(0.975) * cc[2:(1 + ncoef), 2]), - ncol = 2, - dimnames = list(paste0(xlabel, xlev[[1]][-1]), c("2.5 %", "97.5 %")) + ncol = 2, + dimnames = list(paste0(xlabel, xlev[[1]][-1]), c("2.5 %", "97.5 %")) ), decimal.estimate) - + if (family %in% c("binomial", "poisson", "quasipoisson")) { Point.Estimate <- round(exp(stats::coef(model)), decimal.estimate)[2:(1 + ncoef)] CI <- round(exp(matrix(c(cc[2:(1 + ncoef), 1] - qnorm(0.975) * cc[2:(1 + ncoef), 2], cc[2:(1 + ncoef), 1] + qnorm(0.975) * cc[2:(1 + ncoef), 2]), - ncol = 2, - dimnames = list(paste0(xlabel, xlev[[1]][-1]), c("2.5 %", "97.5 %")) + ncol = 2, + dimnames = list(paste0(xlabel, xlev[[1]][-1]), c("2.5 %", "97.5 %")) )), decimal.estimate) } - + # if (length(Point.Estimate) > 1){ # stop("Formula must contain 1 independent variable only.") # } - + # event <- model$y # prop <- round(prop.table(table(event, model$x[, 1]), 2)[2, ] * 100, decimal.percent) pv <- round(summary(model)$coefficients[2:(1 + ncoef), 4], decimal.pvalue) - + # output 만들기 if (ncoef < 2) { data.frame(Variable = "Overall", Count = length(model$y), Percent = 100, `Point Estimate` = Point.Estimate, Lower = CI[1], Upper = CI[2]) %>% dplyr::mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) -> out - + if (family == "binomial") { names(out)[4] <- "OR" } @@ -135,37 +135,38 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, Levels = paste0(xlabel, "=", xlev[[1]]), `Point Estimate` = c("Reference", Point.Estimate), Lower = c("", CI[, 1]), Upper = c("", CI[, 2]) ) %>% dplyr::mutate(`P value` = c("", ifelse(pv >= 0.001, pv, "<0.001")), `P for interaction` = NA) -> out - + + if (family == "binomial") { names(out)[5] <- "OR" } if (family %in% c("poisson", "quasipoisson")) { names(out)[5] <- "RR" } - + rownames(out) <- NULL } - + return(out) } else if (length(var_subgroup) >= 2 | any(grepl(var_subgroup, formula))) { stop("Please input correct subgroup variable.") } else { ### subgroup 지정 한 경우 ### - + # 공변량 있는 경우 formula 변경 if (!is.null(var_cov)) { formula <- as.formula(paste0(deparse(formula), " + ", paste(var_cov, collapse = "+"))) } - + if (any(class(data) == "survey.design")) { ### survey data인 경우 ### - + data$variables[[var_subgroup]] %>% table() %>% names() -> label_val label_val %>% purrr::map(~ possible_svyglm(formula, design = subset(data, get(var_subgroup) == .), x = TRUE, family = family.svyglm)) -> model xlev <- ifelse(length(survey::svyglm(formula, design = data)$xlevels) == 0, NA, survey::svyglm(formula, design = data)$xlevels) - + # pv_int 구하기 pv_int <- tryCatch( { @@ -179,8 +180,9 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, return(NA) } ) + # if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") - + data.design <- data if (family == "binomial") { model.int <- possible_svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep = ""), deparse(formula))), design = data.design, family = quasibinomial()) @@ -191,17 +193,17 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, } else { model.int <- possible_svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep = ""), deparse(formula))), design = data.design, family = quasipoisson()) } - + if (any(is.na(model.int))) { } else if (sum(grepl(":", names(coef(model.int)))) > 1) { pv_anova <- anova(model.int, method = "Wald") pv_int <- round(pv_anova[[length(pv_anova)]][[7]], decimal.pvalue) } - + Count <- as.vector(table(data$variables[[var_subgroup]])) } else { ### survey data가 아닌 경우 ### - + data %>% subset(!is.na(get(var_subgroup))) %>% group_split(get(var_subgroup)) %>% @@ -213,7 +215,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, names() -> label_val xlev <- ifelse(length(stats::glm(formula, data = data)$xlevels) == 0, NA, stats::glm(formula, data = data)$xlevels) model.int <- possible_glm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep = ""), deparse(formula))), data = data, family = family) - + # pv_int 구하기 if (any(is.na(model.int))) { pv_int <- NA @@ -227,10 +229,10 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue) # if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") } - + Count <- as.vector(table(data[[var_subgroup]])) } - + # PE, CI, PV 구하기 if (family %in% c("binomial", "poisson", "quasipoisson")) { Point.Estimate <- model %>% @@ -238,22 +240,22 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, lapply(function(x) { est <- rep(NA, max(length(xlev[[1]]) - 1, 1)) names(est) <- paste0(xlabel, xlev[[1]][-1]) - + for (i in names(est)) { tryCatch(est[i] <- x[i], - error = function(e) est[i] <- NA + error = function(e) est[i] <- NA ) } - + round(exp(est), decimal.estimate) }) - + CI <- model %>% purrr::map(function(model) { cc0 <- tryCatch(summary(model)$coefficients, error = function(e) { return(NA) }) - + ci0 <- matrix(NA, ncol = 2, nrow = max(length(xlev[[1]]) - 1, 1), dimnames = list(paste0(xlabel, xlev[[1]][-1]), c("2.5 %", "97.5 %"))) for (i in rownames(ci0)) { ci0[i, 1] <- tryCatch(cc0[i, 1] - stats::qnorm(0.975) * cc0[i, 2], error = function(e) { @@ -263,7 +265,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, return(NA) }) } - + round(exp(ci0), decimal.estimate) }) } else { @@ -272,22 +274,22 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, lapply(function(x) { est <- rep(NA, max(length(xlev[[1]]) - 1, 1)) names(est) <- paste0(xlabel, xlev[[1]][-1]) - + for (i in names(est)) { tryCatch(est[i] <- x[i], - error = function(e) est[i] <- NA + error = function(e) est[i] <- NA ) } - + round(est, decimal.estimate) }) - + CI <- model %>% purrr::map(function(model) { cc0 <- tryCatch(summary(model)$coefficients, error = function(e) { return(NA) }) - + ci0 <- matrix(NA, ncol = 2, nrow = max(length(xlev[[1]]) - 1, 1), dimnames = list(paste0(xlabel, xlev[[1]][-1]), c("2.5 %", "97.5 %"))) for (i in rownames(ci0)) { ci0[i, 1] <- tryCatch(cc0[i, 1] - stats::qnorm(0.975) * cc0[i, 2], error = function(e) { @@ -297,17 +299,17 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, return(NA) }) } - + round(ci0, decimal.estimate) }) } - + pv <- model %>% purrr::map(function(model) { cc0 <- tryCatch(summary(model)$coefficients, error = function(e) { return(NA) }) - + pvl <- rep(NA, max(length(xlev[[1]]) - 1, 1)) names(pvl) <- paste0(xlabel, xlev[[1]][-1]) for (i in names(pvl)) { @@ -315,15 +317,15 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, return(NA) }) } - + round(pvl, decimal.pvalue) }) - + # output 만들기 if (ncoef < 2) { - data.frame(Variable = paste(" ", label_val), Count = Count, Percent = round(Count / sum(Count) * 100, decimal.percent), `Point Estimate` = unlist(Point.Estimate), Lower = unlist(purrr::map(CI, 1)), Upper = unlist(purrr::map(CI, 2))) %>% - dplyr::mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) -> out - + data.frame(Variable = unlist(paste(" ", label_val)), Count = unlist(Count), Percent = unlist(round(Count / sum(Count) * 100, decimal.percent)), `Point Estimate` = unlist(Point.Estimate), Lower = unlist(purrr::map(CI, 1)), Upper = unlist(purrr::map(CI, 2))) %>% + dplyr::mutate(`P value` = unlist(ifelse(pv >= 0.001, pv, "<0.001")), `P for interaction` = NA) -> out + if (family == "binomial") { names(out)[4] <- "OR" } @@ -336,7 +338,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, Levels = rep(paste0(xlabel, "=", xlev[[1]]), length(label_val)), `Point Estimate` = unlist(lapply(Point.Estimate, function(x) c("Reference", x))), Lower = unlist(lapply(CI, function(x) c("", x[, 1]))), Upper = unlist(lapply(CI, function(x) c("", x[, 2]))) ) %>% dplyr::mutate(`P value` = unlist(lapply(pv, function(x) c("", ifelse(x >= 0.001, x, "<0.001")))), `P for interaction` = NA) -> out - + if (family == "binomial") { names(out)[5] <- "OR" } @@ -344,7 +346,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, names(out)[5] <- "RR" } } - + return(rbind(c(var_subgroup, rep(NA, ncol(out) - 2), ifelse(pv_int >= 0.001, pv_int, "<0.001")), out)) } } @@ -399,7 +401,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, TableSubgroupMultiGLM <- function(formula, var_subgroups = NULL, var_cov = NULL, data, family = "binomial", decimal.estimate = 2, decimal.percent = 1, decimal.pvalue = 3, line = F) { . <- NULL out.all <- TableSubgroupGLM(formula, var_subgroup = NULL, var_cov = var_cov, data = data, family = family, decimal.estimate = decimal.estimate, decimal.percent = decimal.percent, decimal.pvalue = decimal.pvalue) - + if (is.null(var_subgroups)) { return(out.all) } else { From 897124d6f868ee63e3fa8ffaad9e8c197715b780 Mon Sep 17 00:00:00 2001 From: jaehun shon <2021122006@yonsei.ac.kr> Date: Wed, 14 Aug 2024 11:04:56 +0900 Subject: [PATCH 2/3] list datatype to vector --- DESCRIPTION | 2 +- NEWS.md | 5 +++-- R/forestcox.R | 4 ++-- R/forestglm.R | 2 +- 4 files changed, 7 insertions(+), 6 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7337e4a..6bb663b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: jstable Title: Create Tables from Different Types of Regression Version: 1.3.3 -Date: 2024-08-06 +Date: 2024-08-14 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")), diff --git a/NEWS.md b/NEWS.md index 3ba1f6f..b776d81 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,7 @@ #jstable 1.3.3 -*Add cox2.display available in fine-and-gray(competing risk), Multi-State Model (MSM) -*Add TableSubgroupMultiCox available in fine-and-gray(competing risk) +* Update: Add cox2.display available in fine-and-gray(competing risk), Multi-State Model (MSM) +* Update: Add TableSubgroupMultiCox available in fine-and-gray(competing risk) +* Fix: error in `forestcox` and `forestglm` with datatype of P value in table # jstable 1.3.2 diff --git a/R/forestcox.R b/R/forestcox.R index c161e27..0314852 100644 --- a/R/forestcox.R +++ b/R/forestcox.R @@ -467,12 +467,12 @@ TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, # output 만들기 if (ncoef < 2) { out <- data.frame(Variable = paste(" ", label_val), Count = Count, Percent = round(Count / sum(Count) * 100, decimal.percent), `Point Estimate` = Point.Estimate, Lower = CI[, 1], Upper = CI[, 2], check.names = F, row.names = NULL) %>% - mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) + mutate(`P value` = unlist(ifelse(pv >= 0.001, pv, "<0.001")), `P for interaction` = NA) if (!is.null(prop)) { out <- data.frame(Variable = paste(" ", label_val), Count = Count, Percent = round(Count / sum(Count) * 100, decimal.percent), `Point Estimate` = Point.Estimate, Lower = CI[, 1], Upper = CI[, 2], check.names = F, row.names = NULL) %>% cbind(prop) %>% - mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) + mutate(`P value` = unlist(ifelse(pv >= 0.001, pv, "<0.001")), `P for interaction` = NA) } rownames(out) <- NULL diff --git a/R/forestglm.R b/R/forestglm.R index 70caca2..c384bfa 100644 --- a/R/forestglm.R +++ b/R/forestglm.R @@ -323,7 +323,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, # output 만들기 if (ncoef < 2) { - data.frame(Variable = unlist(paste(" ", label_val)), Count = unlist(Count), Percent = unlist(round(Count / sum(Count) * 100, decimal.percent)), `Point Estimate` = unlist(Point.Estimate), Lower = unlist(purrr::map(CI, 1)), Upper = unlist(purrr::map(CI, 2))) %>% + data.frame(Variable = paste(" ", label_val), Count = Count, Percent = round(Count / sum(Count) * 100, decimal.percent), `Point Estimate` = Point.Estimate, Lower = purrr::map(CI, 1), Upper = purrr::map(CI, 2)) %>% dplyr::mutate(`P value` = unlist(ifelse(pv >= 0.001, pv, "<0.001")), `P for interaction` = NA) -> out if (family == "binomial") { From dd21d8a7addb5f074c2dd563d506081a69bfe4f5 Mon Sep 17 00:00:00 2001 From: jaehun shon <2021122006@yonsei.ac.kr> Date: Wed, 14 Aug 2024 13:53:57 +0900 Subject: [PATCH 3/3] . --- R/forestglm.R | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/R/forestglm.R b/R/forestglm.R index c384bfa..f609452 100644 --- a/R/forestglm.R +++ b/R/forestglm.R @@ -136,7 +136,6 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, ) %>% dplyr::mutate(`P value` = c("", ifelse(pv >= 0.001, pv, "<0.001")), `P for interaction` = NA) -> out - if (family == "binomial") { names(out)[5] <- "OR" } @@ -180,7 +179,6 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, return(NA) } ) - # if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") data.design <- data @@ -319,12 +317,12 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, } round(pvl, decimal.pvalue) - }) + }) %>% unlist # output 만들기 if (ncoef < 2) { - data.frame(Variable = paste(" ", label_val), Count = Count, Percent = round(Count / sum(Count) * 100, decimal.percent), `Point Estimate` = Point.Estimate, Lower = purrr::map(CI, 1), Upper = purrr::map(CI, 2)) %>% - dplyr::mutate(`P value` = unlist(ifelse(pv >= 0.001, pv, "<0.001")), `P for interaction` = NA) -> out + data.frame(Variable = paste(" ", label_val), Count = Count, Percent = round(Count / sum(Count) * 100, decimal.percent), `Point Estimate` = unlist(Point.Estimate), Lower = unlist(purrr::map(CI, 1)), Upper = unlist(purrr::map(CI, 2))) %>% + dplyr::mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) -> out if (family == "binomial") { names(out)[4] <- "OR" @@ -413,4 +411,4 @@ TableSubgroupMultiGLM <- function(formula, var_subgroups = NULL, var_cov = NULL, return(rbind(out.all, out.list %>% dplyr::bind_rows())) } } -} +} \ No newline at end of file