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 +
Files Coverage
R 79.22%
Project Totals (13 files) 79.22%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading