Compare 2b0206c ... +0 ... e8779c2

No flags found

Use flags to group coverage reports by test type, project and/or folders.
Then setup custom commit statuses and notifications for each flag.

e.g., #unittest #integration

#production #enterprise

#frontend #backend

Learn more about Codecov Flags here.

Showing 1 of 6 files from the diff.
Newly tracked file
R/forestglm.R created.
Other files ignored by Codecov

@@ -0,0 +1,228 @@
Loading
1 +
#' @title TableSubgroupGLM: Sub-group analysis table for GLM.
2 +
#' @description Sub-group analysis table for GLM.
3 +
#' @param formula formula with survival analysis.
4 +
#' @param var_subgroup 1 sub-group variable for analysis, Default: NULL
5 +
#' @param var_cov Variables for additional adjust, Default: NULL
6 +
#' @param data Data or svydesign in survey package.
7 +
#' @param family family, "gaussian" or "binomial"
8 +
#' @param decimal.estimate Decimal for estimate, Default: 2
9 +
#' @param decimal.percent Decimal for percent, Default: 1
10 +
#' @param decimal.pvalue Decimal for pvalue, Default: 3
11 +
#' @return Sub-group analysis table. 
12 +
#' @details This result is used to make forestplot.
13 +
#' @examples
14 +
#' library(survival);library(dplyr)
15 +
#' lung %>% 
16 +
#'   mutate(status = as.integer(status == 1),
17 +
#'          sex = factor(sex),
18 +
#'          kk = factor(as.integer(pat.karno >= 70))) -> lung
19 +
#' TableSubgroupGLM(status ~ sex, data = lung, family = "binomial")
20 +
#' TableSubgroupGLM(status ~ sex, var_subgroup = "kk", data = lung, family = "binomial")
21 +
#' 
22 +
#' ## survey design
23 +
#' library(survey)
24 +
#' data.design <- svydesign(id = ~1, data = lung)
25 +
#' TableSubgroupGLM(status ~ sex, data = data.design, family = "binomial")
26 +
#' TableSubgroupGLM(status ~ sex, var_subgroup = "kk", data = data.design, family = "binomial")
27 +
#' @seealso 
28 +
#'  \code{\link[purrr]{safely}},\code{\link[purrr]{map}},\code{\link[purrr]{map2}}
29 +
#'  \code{\link[stats]{glm}}
30 +
#'  \code{\link[survey]{svyglm}}
31 +
#'  \code{\link[stats]{confint}}
32 +
#' @rdname TableSubgroupGLM
33 +
#' @export 
34 +
#' @importFrom purrr possibly map_dbl map map2
35 +
#' @importFrom dplyr group_split select filter mutate bind_cols
36 +
#' @importFrom magrittr %>%
37 +
#' @importFrom survey svyglm
38 +
#' @importFrom stats glm confint coefficients anova gaussian quasibinomial
39 +
#' @importFrom utils tail
40 +
41 +
TableSubgroupGLM <- function(formula, var_subgroup = NULL, var_cov = NULL, data, family = "binomial", decimal.estimate = 2, decimal.percent = 1, decimal.pvalue = 3){
42 +
  
43 +
  . <- NULL
44 +
  
45 +
  if (length(formula[[3]]) > 1) stop("Formula must contains only 1 independent variable")
46 +
  if (any(class(data) == "survey.design" & !is.null(var_subgroup))){
47 +
    if (is.numeric(data$variables[[var_subgroup]])) stop("var_subgroup must categorical.")
48 +
    if (length(levels(data$variables[[as.character(formula[[3]])]])) != 2) stop("Independent variable must have 2 levels.")
49 +
  } else if(any(class(data) == "data.frame" & !is.null(var_subgroup))){
50 +
    if (is.numeric(data[[var_subgroup]])) stop("var_subgroup must categorical.")
51 +
    if (length(levels(data[[as.character(formula[[3]])]])) != 2) stop("Independent variable must have 2 levels.")
52 +
  }
53 +
  
54 +
  
55 +
  ## functions with error
56 +
  possible_table <- purrr::possibly(table, NA)
57 +
  possible_prop.table <- purrr::possibly(function(x){prop.table(x, 1)[2, ] * 100}, NA) 
58 +
  possible_pv <- purrr::possibly(function(x){summary(x)[["coefficients"]][1, ] %>% tail(1)}, NA)
59 +
  possible_glm <- purrr::possibly(stats::glm, NA)
60 +
  possible_svyglm <- purrr::possibly(survey::svyglm, NA)
61 +
  possible_confint <- purrr::possibly(stats::confint, NA)
62 +
  possible_modely <- purrr::possibly(function(x){purrr::map_dbl(x, .[["y"]], 1)}, NA)
63 +
  possible_rowone <- purrr::possibly(function(x){x[-1, ]}, NA)
64 +
  
65 +
66 +
  var_cov <- setdiff(var_cov, c(as.character(formula[[3]]), var_subgroup))
67 +
  family.svyglm <- gaussian()
68 +
  if (family  == "binomial") family.svyglm <- quasibinomial()
69 +
  
70 +
  if (is.null(var_subgroup)){
71 +
    if (!is.null(var_cov)){
72 +
      formula <- as.formula(paste0(deparse(formula), " + ", paste(var_cov, collapse = "+")))
73 +
    }
74 +
    
75 +
    if (any(class(data) == "survey.design")){
76 +
      model <- survey::svyglm(formula, design = data, x= T, family = family.svyglm)
77 +
      if (!is.null(model$xlevels) & length(model$xlevels[[1]]) != 2) stop("Categorical independent variable must have 2 levels.")
78 +
    } else{
79 +
      model <- stats::glm(formula, data = data, x= TRUE, family = family)
80 +
      if (!is.null(model$xlevels) & length(model$xlevels[[1]]) != 2) stop("Categorical independent variable must have 2 levels.")
81 +
    }
82 +
    
83 +
    
84 +
    Point.Estimate <- round(coef(model), decimal.estimate)[-1]
85 +
    CI <- round(confint(model)[-1, ], decimal.estimate)
86 +
    if(family == "binomial"){
87 +
      Point.Estimate <- round(exp(coef(model)), decimal.estimate)[-1]
88 +
      CI <- round(exp(confint(model)[-1, ]), decimal.estimate)
89 +
    }
90 +
    
91 +
    
92 +
    
93 +
    #if (length(Point.Estimate) > 1){
94 +
    #  stop("Formula must contain 1 independent variable only.")
95 +
    #}
96 +
    
97 +
    
98 +
    
99 +
    #event <- model$y
100 +
    #prop <- round(prop.table(table(event, model$x[, 1]), 2)[2, ] * 100, decimal.percent)
101 +
    pv <- round(tail(summary(model)$coefficients[-1, ], 1), decimal.pvalue)
102 +
    
103 +
    data.frame(Variable = "Overall", Count = length(model$y), Percent = 100, `Point Estimate` = Point.Estimate, Lower = CI[1], Upper = CI[2]) %>% 
104 +
    dplyr::mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) -> out
105 +
    
106 +
    return(out)
107 +
  } else if (length(var_subgroup) > 1 | any(grepl(var_subgroup, formula))){
108 +
    stop("Please input correct subgroup variable.")
109 +
  } else{
110 +
    if (!is.null(var_cov)){
111 +
      formula <- as.formula(paste0(deparse(formula), " + ", paste(var_cov, collapse = "+")))
112 +
    }
113 +
    if (any(class(data) == "survey.design")){
114 +
      data$variables[[var_subgroup]] %>% table %>% names -> label_val
115 +
      label_val %>% purrr::map(~possible_svyglm(formula, design = subset(data, get(var_subgroup) == .), x = TRUE, family = family.svyglm)) -> model
116 +
      xlev <- survey::svyglm(formula, design = data)$xlevels
117 +
      xlabel <- names(attr(model[[which(!is.na(model))[1]]]$x, "contrast"))[1]
118 +
      pvs_int <- possible_svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), design = data, family = family.svyglm) %>% summary %>% coefficients
119 +
      pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue)
120 +
      if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.")
121 +
      
122 +
      
123 +
      if (length(label_val) > 2){
124 +
        model.int <- survey::svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), design = data, family = family.svyglm)
125 +
        pv_anova <- anova(model.int, method = "Wald")
126 +
        pv_int <- pv_anova[[length(pv_anova)]][[7]]
127 +
      }
128 +
      Count <- as.vector(table(data$variables[[var_subgroup]]))
129 +
      
130 +
    } else{
131 +
      data %>% filter(!is.na(get(var_subgroup))) %>% group_split(get(var_subgroup)) %>% purrr::map(~possible_glm(formula, data = ., x= T, family = family)) -> model
132 +
      data %>% filter(!is.na(get(var_subgroup))) %>% select(var_subgroup) %>% table %>% names -> label_val
133 +
      xlev <- stats::glm(formula, data = data)$xlevels
134 +
      xlabel <- names(attr(model[[which(!is.na(model))[1]]]$x, "contrast"))[1]
135 +
      model.int <- possible_glm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), data = data, family = family)  
136 +
      pvs_int <- model.int %>% summary %>% coefficients
137 +
      pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue)
138 +
      if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.")
139 +
      
140 +
      if (length(label_val) > 2){
141 +
        pv_anova <- anova(model.int, test = "Chisq")
142 +
        pv_int <- round(pv_anova[nrow(pv_anova), 5], decimal.pvalue)
143 +
      }
144 +
      
145 +
      Count <- as.vector(table(data[[var_subgroup]]))
146 +
    }
147 +
    
148 +
    
149 +
    
150 +
    Estimate <- model %>% purrr::map("coefficients", .default = NA) %>% purrr::map_dbl(2, .default = NA)
151 +
    CI0 <- model %>% purrr::map(possible_confint) %>% purrr::map(possible_rowone) %>% Reduce(rbind, .)
152 +
    Point.Estimate <- round(Estimate, decimal.estimate)
153 +
    CI <- round(CI0, decimal.estimate)
154 +
    if (family == "binomial"){
155 +
      Point.Estimate <- round(exp(Estimate), decimal.estimate)
156 +
      CI <- round(exp(CI0), decimal.estimate)
157 +
    }
158 +
    
159 +
    model %>% purrr::map(possible_pv) %>% purrr::map_dbl(~round(., decimal.pvalue)) -> pv
160 +
    
161 +
    
162 +
    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]) %>%
163 +
      dplyr::mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) -> out
164 +
    
165 +
    return(rbind(c(var_subgroup, rep(NA, ncol(out) - 2), ifelse(pv_int >= 0.001, pv_int, "<0.001")), out))
166 +
  }
167 +
}
168 +
169 +
170 +
171 +
#' @title TableSubgroupMultiGLM: Multiple sub-group analysis table for GLM.
172 +
#' @description Multiple sub-group analysis table for GLM.
173 +
#' @param formula formula with survival analysis.
174 +
#' @param var_subgroups Multiple sub-group variables for analysis, Default: NULL
175 +
#' @param var_cov Variables for additional adjust, Default: NULL
176 +
#' @param data Data or svydesign in survey package.
177 +
#' @param family family, "gaussian" or "binomial"
178 +
#' @param decimal.estimate Decimal for estimate, Default: 2
179 +
#' @param decimal.percent Decimal for percent, Default: 1
180 +
#' @param decimal.pvalue Decimal for pvalue, Default: 3
181 +
#' @param line Include new-line between sub-group variables, Default: F
182 +
#' @return Multiple sub-group analysis table. 
183 +
#' @details This result is used to make forestplot.
184 +
#' @examples 
185 +
#' library(survival);library(dplyr)
186 +
#' lung %>% 
187 +
#'   mutate(status = as.integer(status == 1),
188 +
#'          sex = factor(sex),
189 +
#'          kk = factor(as.integer(pat.karno >= 70)),
190 +
#'          kk1 = factor(as.integer(pat.karno >= 60))) -> lung
191 +
#' TableSubgroupMultiGLM(status ~ sex, var_subgroups = c("kk", "kk1"), 
192 +
#'                       data=lung, line = TRUE, family = "binomial")
193 +
#' 
194 +
#' ## survey design
195 +
#' library(survey)
196 +
#' data.design <- svydesign(id = ~1, data = lung)
197 +
#' TableSubgroupMultiGLM(status ~ sex, var_subgroups = c("kk", "kk1"), 
198 +
#'                       data = data.design, family = "binomial")
199 +
#' @seealso 
200 +
#'  \code{\link[purrr]{map}}
201 +
#'  \code{\link[dplyr]{bind}}
202 +
#' @rdname TableSubgroupMultiGLM
203 +
#' @export 
204 +
#' @importFrom purrr map
205 +
#' @importFrom magrittr %>%
206 +
#' @importFrom dplyr bind_rows
207 +
208 +
209 +
TableSubgroupMultiGLM <- function(formula, var_subgroups = NULL, var_cov = NULL, data, family = "binomial", decimal.estimate = 2, decimal.percent = 1, decimal.pvalue = 3, line = F){
210 +
  
211 +
  . <- NULL
212 +
  out.all <- TableSubgroupGLM(formula, var_subgroup = NULL, var_cov = var_cov, data = data, decimal.estimate = decimal.estimate, decimal.percent = decimal.percent, decimal.pvalue = decimal.pvalue)
213 +
  
214 +
  if (is.null(var_subgroups)){
215 +
    return(out.all)
216 +
  } else {
217 +
    out.list <- purrr::map(var_subgroups, ~TableSubgroupGLM(formula, var_subgroup = ., var_cov = var_cov,  data = data, decimal.estimate = decimal.estimate, decimal.percent = decimal.percent, decimal.pvalue = decimal.pvalue))
218 +
    if (line){
219 +
      out.newline <- out.list %>% purrr::map(~rbind(NA, .))
220 +
      return(rbind(out.all, out.newline %>% dplyr::bind_rows()))
221 +
    } else{
222 +
      return(rbind(out.all, out.list %>% dplyr::bind_rows()))
223 +
    }
224 +
  } 
225 +
}
226 +
227 +
228 +

Learn more Showing 1 files with coverage changes found.

New file R/forestglm.R
New
Loading file...
Files Coverage
R -6.18% 79.22%
Project Totals (13 files) 79.22%
Loading