|
| 1 | +# General Consonance Functions Using Profile Likelihood, Wald, or the bootstrap method for linear models. |
| 2 | + |
| 3 | +curve_gen <- function(model, var, method = "default", replicates = 1000, steps = 10000) { |
| 4 | + if (is.list(model) != TRUE) { |
| 5 | + stop("Error: 'model' must be an object with a statistical model") |
| 6 | + } |
| 7 | + if (is.character(method) != TRUE) { |
| 8 | + stop("Error: 'method' must be a character vector") |
| 9 | + } |
| 10 | + if (is.numeric(replicates) != TRUE) { |
| 11 | + stop("Error: 'replicates' must be a numeric vector") |
| 12 | + } |
| 13 | + if (is.numeric(steps) != TRUE) { |
| 14 | + stop("Error: 'steps' must be a numeric vector") |
| 15 | + } |
| 16 | + intrvls <- (0:steps) / steps |
| 17 | + if (method == "default") { |
| 18 | + results <- mclapply(intrvls, FUN = function(i) confint(object = model, level = i)[var, ]) |
| 19 | + } else if (method == "Wald") { |
| 20 | + results <- mclapply(intrvls, FUN = function(i) confint.default(object = model, level = i)[var, ]) |
| 21 | + } else if (method == "lm") { |
| 22 | + results <- mclapply(intrvls, FUN = function(i) confint.lm(object = model, level = i)[var, ]) |
| 23 | + } else if (method == "boot") { |
| 24 | + effect <- coef(model)[[var]] |
| 25 | + boot_dist <- replicate(replicates, |
| 26 | + expr = coef(lm(model$call$formula, |
| 27 | + data = model$model[sample(nrow(model$model), replace = T), ] |
| 28 | + ))[[var]] |
| 29 | + ) - effect |
| 30 | + results <- mclapply(intrvls, FUN = function(i) effect - quantile(boot_dist, probs = (1 + c(i, -i)) / 2)) |
| 31 | + } |
| 32 | + |
| 33 | + df <- data.frame(do.call(rbind, results)) |
| 34 | + intrvl.limit <- c("lower.limit", "upper.limit") |
| 35 | + colnames(df) <- intrvl.limit |
| 36 | + df$intrvl.level <- intrvls |
| 37 | + df$pvalue <- 1 - intrvls |
| 38 | + df$svalue <- -log2(df$pvalue) |
| 39 | + df <- head(df, -1) |
| 40 | + return(df) |
| 41 | +} |
| 42 | + |
| 43 | +# RMD Check |
| 44 | +utils::globalVariables(c("df", "lower.limit", "upper.limit", "intrvl.level", "pvalue", "svalue")) |
0 commit comments