Skip to content

Commit

Permalink
Apply automatic stylistic changes
Browse files Browse the repository at this point in the history
  • Loading branch information
github-actions[bot] committed Aug 14, 2024
1 parent 843397f commit 0212a49
Showing 1 changed file with 56 additions and 55 deletions.
111 changes: 56 additions & 55 deletions R/forestglm.R
Original file line number Diff line number Diff line change
Expand Up @@ -42,15 +42,15 @@

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))) {
if (is.numeric(data$variables[[var_subgroup]])) stop("var_subgroup must categorical.")
} 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) {
Expand All @@ -64,65 +64,65 @@ 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.")
} else {
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"
}
Expand All @@ -135,37 +135,37 @@ 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(
{
Expand All @@ -180,7 +180,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data,
}
)
# 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())
Expand All @@ -191,17 +191,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)) %>%
Expand All @@ -213,7 +213,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
Expand All @@ -227,33 +227,33 @@ 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 %>%
purrr::map("coefficients", default = NA) %>%
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) {
Expand All @@ -263,7 +263,7 @@ TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data,
return(NA)
})
}

round(exp(ci0), decimal.estimate)
})
} else {
Expand All @@ -272,22 +272,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) {
Expand All @@ -297,33 +297,34 @@ 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)) {
pvl[i] <- tryCatch(cc0[i, 4], error = function(e) {
return(NA)
})
}

round(pvl, decimal.pvalue)
}) %>% unlist

}) %>%
unlist()

# 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

if (family == "binomial") {
names(out)[4] <- "OR"
}
Expand All @@ -336,15 +337,15 @@ 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"
}
if (family %in% c("poisson", "quasipoisson")) {
names(out)[5] <- "RR"
}
}

return(rbind(c(var_subgroup, rep(NA, ncol(out) - 2), ifelse(pv_int >= 0.001, pv_int, "<0.001")), out))
}
}
Expand Down Expand Up @@ -399,7 +400,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 {
Expand All @@ -411,4 +412,4 @@ TableSubgroupMultiGLM <- function(formula, var_subgroups = NULL, var_cov = NULL,
return(rbind(out.all, out.list %>% dplyr::bind_rows()))
}
}
}
}

0 comments on commit 0212a49

Please sign in to comment.