R/forestglm.R
changed.
R/forestcox.R
changed.
Other files ignored by Codecov
DESCRIPTION
has changed.
60 | 60 | possible_modely <- purrr::possibly(function(x){purrr::map_dbl(x, .[["y"]], 1)}, NA) |
|
61 | 61 | possible_rowone <- purrr::possibly(function(x){x[2, ]}, NA) |
|
62 | 62 | ||
63 | - | ||
63 | + | xlabel <- setdiff(as.character(formula)[[3]], "+")[1] |
|
64 | + | ncoef <- ifelse(any(class(data) == "survey.design"), ifelse(length(levels(data$variables[[xlabel]])) <= 2, 1, length(levels(data$variables[[xlabel]])) - 1), |
|
65 | + | ifelse(length(levels(data[[xlabel]])) <= 2, 1, length(levels(data[[xlabel]])) - 1)) |
|
64 | 66 | var_cov <- setdiff(var_cov, c(as.character(formula[[3]]), var_subgroup)) |
|
65 | 67 | family.svyglm <- gaussian() |
|
66 | 68 | if (family == "binomial") family.svyglm <- quasibinomial() |
116 | 118 | data$variables[[var_subgroup]] %>% table %>% names -> label_val |
|
117 | 119 | label_val %>% purrr::map(~possible_svyglm(formula, design = subset(data, get(var_subgroup) == .), x = TRUE, family = family.svyglm)) -> model |
|
118 | 120 | xlev <- survey::svyglm(formula, design = data)$xlevels |
|
119 | - | xlabel <- setdiff(as.character(formula)[[3]], "+")[1] |
|
120 | 121 | pvs_int <- possible_svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), design = data, family = family.svyglm) %>% summary %>% coefficients |
|
121 | 122 | pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue) |
|
122 | 123 | #if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") |
|
124 | + | model.int <- survey::svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), design = data, family = family.svyglm) |
|
123 | 125 | ||
124 | 126 | ||
125 | - | if (length(label_val) > 2 | length(xlev) > 2){ |
|
126 | - | model.int <- survey::svyglm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), design = data, family = family.svyglm) |
|
127 | + | if (sum(grepl(":", names(coef(model.int)))) > 1){ |
|
127 | 128 | pv_anova <- anova(model.int, method = "Wald") |
|
128 | 129 | pv_int <- round(pv_anova[[length(pv_anova)]][[7]], decimal.pvalue) |
|
129 | 130 | } |
133 | 134 | data %>% subset(!is.na(get(var_subgroup))) %>% group_split(get(var_subgroup)) %>% purrr::map(~possible_glm(formula, data = ., x= T, family = family)) -> model |
|
134 | 135 | data %>% subset(!is.na(get(var_subgroup))) %>% select(var_subgroup) %>% table %>% names -> label_val |
|
135 | 136 | xlev <- stats::glm(formula, data = data)$xlevels |
|
136 | - | xlabel <- setdiff(as.character(formula)[[3]], "+")[1] |
|
137 | 137 | model.int <- possible_glm(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), data = data, family = family) |
|
138 | 138 | pvs_int <- model.int %>% summary %>% coefficients |
|
139 | 139 | pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue) |
|
140 | 140 | #if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") |
|
141 | 141 | ||
142 | - | if (length(label_val) > 2 | length(xlev) > 2){ |
|
142 | + | if (sum(grepl(":", names(coef(model.int)))) > 1){ |
|
143 | 143 | pv_anova <- anova(model.int, test = "Chisq") |
|
144 | 144 | pv_int <- round(pv_anova[nrow(pv_anova), 5], decimal.pvalue) |
|
145 | 145 | } |
67 | 67 | ||
68 | 68 | formula.km <- formula |
|
69 | 69 | var_cov <- setdiff(var_cov, c(as.character(formula[[3]]), var_subgroup)) |
|
70 | + | xlabel <- setdiff(as.character(formula)[[3]], "+")[1] |
|
71 | + | ||
72 | + | ncoef <- ifelse(any(class(data) == "survey.design"), ifelse(length(levels(data$variables[[xlabel]])) <= 2, 1, length(levels(data$variables[[xlabel]])) - 1), |
|
73 | + | ifelse(length(levels(data[[xlabel]])) <= 2, 1, length(levels(data[[xlabel]])) - 1)) |
|
74 | + | ||
70 | 75 | if (is.null(var_subgroup)){ |
|
71 | 76 | if (!is.null(var_cov)){ |
|
72 | 77 | formula <- as.formula(paste0(deparse(formula), " + ", paste(var_cov, collapse = "+"))) |
87 | 92 | #out.kap <- paste(res.kap.times[["n.event"]], " (", round(100 * (1 - res.kap.times[["surv"]]), decimal.percent), ")", sep = "") |
|
88 | 93 | } |
|
89 | 94 | ||
90 | - | Point.Estimate <- round(exp(coef(model)), decimal.hr)[1] |
|
91 | 95 | ||
92 | 96 | ||
97 | + | Point.Estimate <- round(exp(coef(model)), decimal.hr)[1:ncoef] |
|
98 | + | ||
93 | 99 | #if (length(Point.Estimate) > 1){ |
|
94 | 100 | # stop("Formula must contain 1 independent variable only.") |
|
95 | 101 | #} |
|
96 | 102 | ||
97 | 103 | ||
98 | - | CI <- round(exp(confint(model)[1, ]), decimal.hr) |
|
104 | + | CI <- round(exp(confint(model)[1:ncoef, ]), decimal.hr) |
|
99 | 105 | event <- purrr::map_dbl(model$y, 1) %>% tail(model$n) |
|
100 | 106 | #prop <- round(prop.table(table(event, model$x[, 1]), 2)[2, ] * 100, decimal.percent) |
|
101 | - | pv <- round(summary(model)$coefficients[1, "Pr(>|z|)"], decimal.pvalue) |
|
107 | + | pv <- round(summary(model)$coefficients[1:ncoef, "Pr(>|z|)"], decimal.pvalue) |
|
102 | 108 | ||
103 | - | out <- data.frame(Variable = "Overall", Count = model$n, Percent = 100, `Point Estimate` = Point.Estimate, Lower = CI[1], Upper = CI[2], check.names = F) %>% |
|
104 | - | mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) |
|
105 | - | ||
106 | - | if (!is.null(names(prop))){ |
|
109 | + | if (ncoef <= 2){ |
|
107 | 110 | out <- data.frame(Variable = "Overall", Count = model$n, Percent = 100, `Point Estimate` = Point.Estimate, Lower = CI[1], Upper = CI[2], check.names = F) %>% |
|
108 | - | cbind(t(prop)) %>% |
|
109 | - | mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) |
|
110 | - | } |
|
111 | + | mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) |
|
112 | + | ||
113 | + | if (!is.null(names(prop))){ |
|
114 | + | out <- data.frame(Variable = "Overall", Count = model$n, Percent = 100, `Point Estimate` = Point.Estimate, Lower = CI[1], Upper = CI[2], check.names = F) %>% |
|
115 | + | cbind(t(prop)) %>% |
|
116 | + | mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) |
|
117 | + | } |
|
118 | + | } else{ |
|
119 | + | out <- data.frame(Variable = c("Overall", rep("", length(Point.Estimate))), Subgroup = rep("", length(Point.Estimate) + 1), Levels = levels(data[[xlabel]]), `Point Estimate` = c("Reference", Point.Estimate), Lower = c("", CI[, 1]), Upper = c("", CI[, 2]), check.names = F) %>% |
|
120 | + | mutate(`P value` = c("", ifelse(pv >= 0.001, pv, "<0.001")), `P for interaction` = NA) |
|
121 | + | rownames(out) <- NULL |
|
122 | + | } |
|
111 | 123 | ||
112 | 124 | return(out) |
|
113 | 125 | } else if (length(var_subgroup) > 1 | any(grepl(var_subgroup, formula))){ |
120 | 132 | data$variables[[var_subgroup]] %>% table %>% names -> label_val |
|
121 | 133 | label_val %>% purrr::map(~possible_svycoxph(formula, design = subset(data, get(var_subgroup) == .), x = TRUE)) -> model |
|
122 | 134 | xlev <- survey::svycoxph(formula, design = data)$xlevels |
|
123 | - | xlabel <- setdiff(as.character(formula)[[3]], "+")[1] |
|
124 | 135 | pvs_int <- possible_svycoxph(as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))), design = data) %>% summary %>% coefficients |
|
125 | 136 | pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue) |
|
126 | 137 | #if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") |
|
138 | + | model.int <- survey::svycoxph(as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))), design = data) |
|
127 | 139 | ||
128 | - | ||
129 | - | if (length(label_val) > 2 | length(xlev) > 2){ |
|
130 | - | model.int <- survey::svycoxph(as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))), design = data) |
|
140 | + | if (sum(grepl(":", names(coef(model.int)))) > 1){ |
|
141 | + | model.int$call$formula <- as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))) |
|
131 | 142 | pv_anova <- survey::regTermTest(model.int, as.formula(paste0("~", xlabel, ":", var_subgroup))) |
|
132 | 143 | pv_int <- round(pv_anova$p[1], decimal.pvalue) |
|
133 | 144 | } |
149 | 160 | data %>% filter(!is.na(get(var_subgroup))) %>% split(.[[var_subgroup]]) %>% purrr::map(~possible_coxph(formula, data = ., x= T)) -> model |
|
150 | 161 | data %>% filter(!is.na(get(var_subgroup))) %>% select(var_subgroup) %>% table %>% names -> label_val |
|
151 | 162 | xlev <- survival::coxph(formula, data = data)$xlevels |
|
152 | - | xlabel <- setdiff(as.character(formula)[[3]], "+")[1] |
|
153 | 163 | model.int <- possible_coxph(as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))), data = data) |
|
154 | - | pvs_int <- model.int %>% summary %>% coefficients |
|
155 | - | pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue) |
|
156 | - | #if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") |
|
157 | - | ||
158 | - | if (length(label_val) > 2 | length(xlev) > 2){ |
|
164 | + | if (sum(grepl(":", names(coef(model.int)))) == 1){ |
|
165 | + | pvs_int <- model.int %>% summary %>% coefficients |
|
166 | + | pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue) |
|
167 | + | #if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") |
|
168 | + | if (!is.numeric(data[[xlabel]])){ |
|
169 | + | res.kap.times <- data %>% filter(!is.na(get(var_subgroup))) %>% split(.[[var_subgroup]]) %>% purrr::map(~survival::survfit(formula.km, data = .)) %>% purrr::map(~summary(., times = time_eventrate, extend = T)) |
|
170 | + | prop <- res.kap.times %>% purrr::map(~round(100 * (1 - .[["surv"]]), decimal.percent)) %>% dplyr::bind_cols() %>% t |
|
171 | + | colnames(prop) <- xlev[[1]] |
|
172 | + | } else{ |
|
173 | + | prop <- NULL |
|
174 | + | } |
|
175 | + | } else{ |
|
159 | 176 | model.int$call$formula <- as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))) |
|
160 | 177 | model.int$call$data <- as.name("data") |
|
161 | 178 | pv_anova <- anova(model.int) |
|
162 | 179 | pv_int <- round(pv_anova[nrow(pv_anova), 4], decimal.pvalue) |
|
163 | - | } |
|
164 | - | ||
165 | - | if (!is.numeric(data[[xlabel]])){ |
|
166 | - | res.kap.times <- data %>% filter(!is.na(get(var_subgroup))) %>% split(.[[var_subgroup]]) %>% purrr::map(~survival::survfit(formula.km, data = .)) %>% purrr::map(~summary(., times = time_eventrate, extend = T)) |
|
167 | - | prop <- res.kap.times %>% purrr::map(~round(100 * (1 - .[["surv"]]), decimal.percent)) %>% dplyr::bind_cols() %>% t |
|
168 | - | colnames(prop) <- xlev[[1]] |
|
169 | - | } else{ |
|
170 | 180 | prop <- NULL |
|
171 | 181 | } |
|
182 | + | ||
183 | + | ||
172 | 184 | } |
|
173 | 185 | ||
174 | 186 | model %>% purrr::map_dbl("n", .default = NA) -> Count |
|
175 | - | model %>% purrr::map("coefficients", .default = NA) %>% purrr::map_dbl(1) %>% exp %>% round(decimal.hr) -> Point.Estimate |
|
176 | - | model %>% purrr::map(possible_confint) %>% purrr::map(possible_rowone) %>% Reduce(rbind, .) %>% exp %>% round(decimal.hr) -> CI |
|
187 | + | model %>% purrr::map("coefficients", .default = NA) %>% lapply(function(x){round(exp(x[1:ncoef]), decimal.hr)}) %>% unlist -> Point.Estimate |
|
188 | + | #model %>% purrr::map("coefficients", .default = NA) %>% purrr::map_dbl(1) %>% exp %>% round(decimal.hr) -> Point.Estimate |
|
189 | + | model %>% purrr::map(possible_confint) %>% Reduce(rbind, .) %>% exp %>% round(decimal.hr) -> CI |
|
190 | + | #model %>% purrr::map(possible_confint) %>% purrr::map(possible_rowone) %>% Reduce(rbind, .) %>% exp %>% round(decimal.hr) -> CI |
|
177 | 191 | #model %>% purrr::map("y") %>% purrr::map(~purrr::map_dbl(., 1)) %>% purrr::map(~tail(., length(.)/2)) -> event |
|
178 | 192 | #purrr::map2(event, model, ~possible_table(.x, .y[["x"]][, 1])) %>% purrr::map(possible_prop.table) %>% purrr::map(~round(., decimal.percent)) %>% Reduce(rbind, .) -> prop |
|
179 | 193 | model %>% purrr::map(possible_pv) %>% purrr::map_dbl(~round(., decimal.pvalue)) -> pv |
Files | Coverage |
---|---|
R | 85.40% |
Project Totals (13 files) | 85.40% |