Merge branch 'master' of https://github.com/tylermorganwall/skpr

Author: tylermorganwall

R/add_ci_bounds_mc_power.R

@@ -5,10 +5,16 @@
 #'@keywords internal
 add_ci_bounds_mc_power = function(power_results, nsim, conf = 0.95) {
   get_lcbs = function(power) {
-    return(unlist(lapply(power, \(x) (binom.test(x*nsim, n = nsim, conf.level = conf)[["conf.int"]][1]))))
+    return(unlist(lapply(
+      power,
+      \(x) (binom.test(x * nsim, n = nsim, conf.level = conf)[["conf.int"]][1])
+    )))
   }
   get_ucbs = function(power) {
-    return(unlist(lapply(power, \(x) (binom.test(x*nsim, n = nsim, conf.level = conf)[["conf.int"]][2]))))
+    return(unlist(lapply(
+      power,
+      \(x) (binom.test(x * nsim, n = nsim, conf.level = conf)[["conf.int"]][2])
+    )))
   }
   lcb = get_lcbs(power_results[["power"]])
   ucb = get_ucbs(power_results[["power"]])

R/anticoef_from_delta.R

@@ -20,46 +20,45 @@ anticoef_from_delta = function(default_coef, delta, glmfamily) {
   if (glmfamily == "gaussian") {
     if (length(delta) == 1) {
       anticoef = default_coef * delta / 2
-    }
-    else {
+    } else {
       anticoef = default_coef * (delta[2] - delta[1]) / 2
     }
-  }
-  else if (glmfamily == "exponential") {
+  } else if (glmfamily == "exponential") {
     if (length(delta) == 1) {
       anticoef = default_coef * delta / 2
-      warning("default or length 1 delta used with glmfamily == 'exponential'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate.")
-    }
-    else {
+      warning(
+        "default or length 1 delta used with glmfamily == 'exponential'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate."
+      )
+    } else {
       anticoef = gen_exponential_anticoef(default_coef, delta[1], delta[2])
     }
-  }
-  else if (glmfamily == "poisson") {
+  } else if (glmfamily == "poisson") {
     if (length(delta) == 1) {
       anticoef = default_coef * delta / 2
-      warning("default or length 1 delta used with glmfamily == 'poisson'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate.")
-    }
-    else {
+      warning(
+        "default or length 1 delta used with glmfamily == 'poisson'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate."
+      )
+    } else {
       anticoef = gen_exponential_anticoef(default_coef, delta[1], delta[2])
     }
-  }
-  else if (glmfamily == "binomial") {
+  } else if (glmfamily == "binomial") {
     if (length(delta) == 1) {
       anticoef = default_coef * delta / 2
-      warning("default or length 1 delta used with glmfamily == 'binomial'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate.")
-    }
-    else {
+      warning(
+        "default or length 1 delta used with glmfamily == 'binomial'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate."
+      )
+    } else {
       anticoef = gen_binomial_anticoef(default_coef, delta[1], delta[2])
     }
-  }
-  else {
+  } else {
     if (length(delta) == 1) {
       anticoef = default_coef * delta / 2
-    }
-    else {
+    } else {
       anticoef = gen_exponential_anticoef(default_coef, delta[1], delta[2])
     }
-    warning("delta parameter used with unsupported glmfamily. Make sure the generated anticipated coefficients are appropriate.")
+    warning(
+      "delta parameter used with unsupported glmfamily. Make sure the generated anticipated coefficients are appropriate."
+    )
   }
   anticoef
 }

R/anticoef_from_delta_surv.R

@@ -20,37 +20,36 @@ anticoef_from_delta_surv = function(default_coef, delta, distribution) {
   if (distribution == "gaussian") {
     if (length(delta) == 1) {
       anticoef = default_coef * delta / 2
-    }
-    else {
+    } else {
       anticoef = default_coef * (delta[2] - delta[1]) / 2
     }
-  }
-  else if (distribution == "exponential") {
+  } else if (distribution == "exponential") {
     if (length(delta) == 1) {
       anticoef = default_coef * delta / 2
-      warning("default or length 1 delta used with distribution == 'exponential'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate.")
-    }
-    else {
+      warning(
+        "default or length 1 delta used with distribution == 'exponential'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate."
+      )
+    } else {
       anticoef = gen_exponential_anticoef(default_coef, delta[1], delta[2])
     }
-  }
-  else if (distribution == "lognormal") {
+  } else if (distribution == "lognormal") {
     if (length(delta) == 1) {
       anticoef = default_coef * delta / 2
-      warning("default or length 1 delta used with distribution == 'lognormal'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate.")
-    }
-    else {
+      warning(
+        "default or length 1 delta used with distribution == 'lognormal'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate."
+      )
+    } else {
       anticoef = gen_exponential_anticoef(default_coef, delta[1], delta[2])
     }
-  }
-  else {
+  } else {
     if (length(delta) == 1) {
       anticoef = default_coef * delta / 2
-    }
-    else {
+    } else {
       anticoef = gen_exponential_anticoef(default_coef, delta[1], delta[2])
     }
-    warning("delta parameter used with unsupported distribution. Make sure the generated anticipated coefficients are appropriate.")
+    warning(
+      "delta parameter used with unsupported distribution. Make sure the generated anticipated coefficients are appropriate."
+    )
   }
   anticoef
 }

R/calc_conservative_anticoef.R

@@ -27,7 +27,7 @@ calc_conservative_anticoef = function(results, effectsize) {
       if (length(which(abs(powers - min(powers)) < 1E-10)) > 1) {
         numberofequal = length(which(abs(powers - min(powers)) < 1E-10))
         exponents = 1:numberofequal + 1
-        values = rep(-1, numberofequal) ^ exponents
+        values = rep(-1, numberofequal)^exponents
         if (numberofequal > 2) {
           values[3:length(values)] = 0
         }

R/calc_interaction_degrees.R

@@ -9,10 +9,19 @@
 #'@param split_degrees The number of degrees of freedom for each main effect term
 #'@return Degrees of freedom vector
 #'@keywords internal
-calc_interaction_degrees = function(design, model, contrast, nointercept,
-                                    split_layers, split_degrees, split_plot_structure) {
+calc_interaction_degrees = function(
+  design,
+  model,
+  contrast,
+  nointercept,
+  split_layers,
+  split_degrees,
+  split_plot_structure
+) {
   contrastslistspd = list()
-  for (x in names(design[lapply(design, class) %in% c("factor", "character")])) {
+  for (x in names(design[
+    lapply(design, class) %in% c("factor", "character")
+  ])) {
     contrastslistspd[[x]] = contrast
   }
 
@@ -22,82 +31,116 @@ calc_interaction_degrees = function(design, model, contrast, nointercept,
     contrastslistspd = NULL
   }
 
-  model = as.formula(paste0("~", paste(attr(terms.formula(model), "term.labels"), collapse = " + ")))
-  if(!nointercept) {
+  model = as.formula(paste0(
+    "~",
+    paste(attr(terms.formula(model), "term.labels"), collapse = " + ")
+  ))
+  if (!nointercept) {
     character_model = as.character(model)[-1]
   } else {
-    character_model = paste(c(as.character(model)[-1]," + -1"), collapse="")
+    character_model = paste(c(as.character(model)[-1], " + -1"), collapse = "")
   }
   splitterms = unlist(strsplit(character_model, split = " + ", fixed = TRUE))
   ismaineffect = rep(FALSE, length(splitterms))
   ismaineffect[1:length(split_layers)] = TRUE
   names(ismaineffect) = splitterms
   interactions = list()
-  if(max(split_layers) > 0) {
-    for(i in 1:max(split_layers, na.rm = TRUE)) {
+  if (max(split_layers) > 0) {
+    for (i in 1:max(split_layers, na.rm = TRUE)) {
       wholeplotterms = colnames(design)[split_layers == i]
       wholeorwholeinteraction = rep(FALSE, length(splitterms))
       for (term in wholeplotterms) {
-        regex = paste0("(\\b", term, "\\b)|(\\b", term, ":)|(:", term, "\\b)|(\\b", term, "\\s\\*)|(\\*\\s", term, "\\b)")
-        wholeorwholeinteraction = wholeorwholeinteraction | grepl(regex, splitterms, perl = TRUE)
+        regex = paste0(
+          "(\\b",
+          term,
+          "\\b)|(\\b",
+          term,
+          ":)|(:",
+          term,
+          "\\b)|(\\b",
+          term,
+          "\\s\\*)|(\\*\\s",
+          term,
+          "\\b)"
+        )
+        wholeorwholeinteraction = wholeorwholeinteraction |
+          grepl(regex, splitterms, perl = TRUE)
       }
       interactions[[i]] = wholeorwholeinteraction
     }
   } else {
-    for(i in seq_along(1:max(split_layers, na.rm = TRUE))) {
+    for (i in seq_along(1:max(split_layers, na.rm = TRUE))) {
       interactions[[i]] = FALSE
     }
   }
-  if(!nointercept) {
-    degrees_of_freedom = rep(min(split_degrees), length(splitterms)+1)
-    for(i in 1:length(split_layers)) {
-      degrees_of_freedom[i+1] = split_degrees[split_layers[i]+1]
+  if (!nointercept) {
+    degrees_of_freedom = rep(min(split_degrees), length(splitterms) + 1)
+    for (i in 1:length(split_layers)) {
+      degrees_of_freedom[i + 1] = split_degrees[split_layers[i] + 1]
     }
     names(degrees_of_freedom) = c("(Intercept)", splitterms)
   } else {
     degrees_of_freedom = rep(min(split_degrees), length(splitterms))
-    for(i in 1:length(split_layers)) {
+    for (i in 1:length(split_layers)) {
       degrees_of_freedom[i] = split_degrees[split_layers[i]]
     }
     names(degrees_of_freedom) = c(splitterms)
   }
   #get higher order terms
-  for(term in colnames(design)) {
-    arith_terms = splitterms[grepl(paste0("(I\\(", term, "\\^.+\\))"), splitterms, perl = TRUE)]
-    for(arith_term in arith_terms) {
+  for (term in colnames(design)) {
+    arith_terms = splitterms[grepl(
+      paste0("(I\\(", term, "\\^.+\\))"),
+      splitterms,
+      perl = TRUE
+    )]
+    for (arith_term in arith_terms) {
       degrees_of_freedom[arith_term] = degrees_of_freedom[term]
     }
   }
   #Effect terms
-  for(i in seq_along(1:max(split_layers, na.rm = TRUE))) {
+  for (i in seq_along(1:max(split_layers, na.rm = TRUE))) {
     subplotterms = colnames(design)
-    subplotterms = subplotterms[!(subplotterms %in% colnames(design)[split_layers >= i])]
+    subplotterms = subplotterms[
+      !(subplotterms %in% colnames(design)[split_layers >= i])
+    ]
     wholeplotterms = colnames(design)[split_layers < i]
 
     regularmodel = rep(FALSE, length(splitterms))
     for (term in subplotterms) {
-      regex = paste0("(\\b", term, "\\b)|(\\b", term, ":)|(:", term, "\\b)|(\\b", term, "\\s\\*)|(\\*\\s", term, "\\b)")
+      regex = paste0(
+        "(\\b",
+        term,
+        "\\b)|(\\b",
+        term,
+        ":)|(:",
+        term,
+        "\\b)|(\\b",
+        term,
+        "\\s\\*)|(\\*\\s",
+        term,
+        "\\b)"
+      )
       regularmodel = regularmodel | grepl(regex, splitterms, perl = TRUE)
     }
     #get whole:whole interaction terms
     wholewholeinteractionterms = splitterms[!ismaineffect & interactions[[i]]]
-    for(term in wholewholeinteractionterms) {
+    for (term in wholewholeinteractionterms) {
       subterms = unlist(strsplit(term, split = ":", fixed = TRUE))
       maxdegree = c()
-      for(subterm in subterms) {
-        if(!subterm %in% wholeplotterms) {
+      for (subterm in subterms) {
+        if (!subterm %in% wholeplotterms) {
           maxdegree = max(c(maxdegree, degrees_of_freedom[subterm]))
         }
       }
       degrees_of_freedom[term] = maxdegree
     }
     #get whole:non-whole interaction terms
     wholeinteractionterms = splitterms[regularmodel & interactions[[i]]]
-    for(term in wholeinteractionterms) {
+    for (term in wholeinteractionterms) {
       subterms = unlist(strsplit(term, split = ":", fixed = TRUE))
       maxdegree = c()
-      for(subterm in subterms) {
-        if(!subterm %in% wholeplotterms) {
+      for (subterm in subterms) {
+        if (!subterm %in% wholeplotterms) {
           maxdegree = max(c(maxdegree, degrees_of_freedom[subterm]))
         }
       }

R/calcnoncentralparam.R

@@ -8,14 +8,21 @@
 #'@param vinv variance-covariance matrix
 #'@return The non-centrality parameter for the F test
 #'@keywords internal
-calcnoncentralparam = function(X, L, b, vinv=NULL) {
+calcnoncentralparam = function(X, L, b, vinv = NULL) {
   if (is.null(vinv)) {
     vinv = diag(nrow(X))
   }
-  matrix_solve = tryCatch({solve(t(X) %*% vinv %*% X)}, error = function(e) e)
-  if(inherits(matrix_solve,"error")) {
-    stop(sprintf("skpr: Can't calculate non-centrality parameter with given design matrix (X) and variance-covariance (Vinv) matrix due to `t(X) %*% Vinv %*% X` being singular. `solve()` error message: '%s'",
-                 as.character(matrix_solve)))
+  matrix_solve = tryCatch(
+    {
+      solve(t(X) %*% vinv %*% X)
+    },
+    error = function(e) e
+  )
+  if (inherits(matrix_solve, "error")) {
+    stop(sprintf(
+      "skpr: Can't calculate non-centrality parameter with given design matrix (X) and variance-covariance (Vinv) matrix due to `t(X) %*% Vinv %*% X` being singular. `solve()` error message: '%s'",
+      as.character(matrix_solve)
+    ))
   }
   return(t(L %*% b) %*% solve(L %*% matrix_solve %*% t(L)) %*% L %*% b)
 }

R/calculate_block_nesting.R

@@ -9,22 +9,26 @@
 calculate_block_nesting = function(blockgroups, blockstructure) {
   #`nested_level` is a label for which is the highest blocking level (excluding the lowest level)
   # the column is nested within
-  nested_level = rep(1,length(blockgroups))
+  nested_level = rep(1, length(blockgroups))
   for (i in seq_len(length(blockgroups))) {
-    isblockcol = rep(TRUE,ncol(blockstructure))
+    isblockcol = rep(TRUE, ncol(blockstructure))
     temp_block_str = blockgroups[[i]]
     temp_block_str_names = names(temp_block_str)
     for (j in seq_len(length(temp_block_str))) {
       block_name = temp_block_str_names[j]
       if (temp_block_str[j] - 1 > 0) {
-        block_idx = which(blockstructure[,i] == block_name)
-        is_constant = lapply(apply(blockstructure[block_idx, ], 2, unique), length) == 1
+        block_idx = which(blockstructure[, i] == block_name)
+        is_constant = lapply(
+          apply(blockstructure[block_idx, ], 2, unique),
+          length
+        ) ==
+          1
         isblockcol = isblockcol & is_constant
       }
     }
     isblockcol[i] = FALSE
     block_nested_within = which(isblockcol)
-    if(length(block_nested_within) != 0) {
+    if (length(block_nested_within) != 0) {
       nested_level[i] = max(block_nested_within) + 1
     }
   }