permutations

Updates so number of permutations flows through table generation properly

Author: kdyson

AICc_PERMANOVA.R

@@ -53,7 +53,15 @@ AICc.PERMANOVA2 <- function(adonis2.model) {
 
     # Ok, now extract appropriate terms from the adonis model Calculating AICc
     # using residual sum of squares (RSS or SSE) since I don't think that adonis
-    # returns something I can use as a liklihood function... o wait maximum likelihood may be MSE.
+    # returns something I can use as a likelihood function... maximum likelihood
+    # and MSE estimates are the same when distribution is gaussian See e.g.
+    # https://www.jessicayung.com/mse-as-maximum-likelihood/;
+    # https://towardsdatascience.com/probability-concepts-explained-maximum-likelihood-estimation-c7b4342fdbb1
+    # So using RSS or MSE estimates is fine as long as the residuals are
+    # Gaussian https://robjhyndman.com/hyndsight/aic/ If models have different
+    # conditional likelihoods then AIC is not valid. However, comparing models
+    # with different error distributions is ok (above link).
+    
 
     RSS <- adonis2.model$SumOfSqs[ length(adonis2.model$SumOfSqs) - 1 ]
     MSE <- RSS / adonis2.model$Df[ length(adonis2.model$Df) - 1 ]
@@ -63,27 +71,29 @@ AICc.PERMANOVA2 <- function(adonis2.model) {
     k <- nn - adonis2.model$Df[ length(adonis2.model$Df) - 1 ]
 
 
-    # AIC : 2*k + n*ln(RSS)
+    # AIC : 2*k + n*ln(RSS/n)
     # AICc: AIC + [2k(k+1)]/(n-k-1)
 
     # based on https://en.wikipedia.org/wiki/Akaike_information_criterion;
+    # https://www.statisticshowto.datasciencecentral.com/akaikes-information-criterion/ ;
     # https://www.researchgate.net/post/What_is_the_AIC_formula;
-    # http://avesbiodiv.mncn.csic.es/estadistica/ejemploaic.pdf
+    # http://avesbiodiv.mncn.csic.es/estadistica/ejemploaic.pdf;
+    # https://medium.com/better-programming/data-science-modeling-how-to-use-linear-regression-with-python-fdf6ca5481be 
 
     # AIC.g is generalized version of AIC = 2k + n [Ln( 2(pi) RSS/n ) + 1]
     # AIC.pi = k + n [Ln( 2(pi) RSS/(n-k) ) +1],
 
-    AIC <- 2*k + nn*log(RSS)
+    AIC <- 2*k + nn*log(RSS/nn)
     AIC.g <- 2*k + nn * (1 + log( 2 * pi * RSS / nn))
     AIC.MSE <- 2*k + nn * log(MSE)
     AIC.pi <- k + nn*(1 + log( 2*pi*RSS/(nn-k) )   )
     AICc <- AIC + (2*k*(k + 1))/(nn - k - 1)
     AICc.MSE <- AIC.MSE + (2*k*(k + 1))/(nn - k - 1)
     AICc.pi <- AIC.pi + (2*k*(k + 1))/(nn - k - 1)
 
-    output <- list("AIC" = AIC, "AIC.g" = AIC.g, "AICc" = AICc,
+    output <- list("AIC" = AIC, "AICc" = AICc, "AIC.g" = AIC.g, 
                    "AIC.MSE" = AIC.MSE, "AICc.MSE" = AICc.MSE,
-                   "AIC.pi" = AIC.pi, "AICc.pi" = AICc.pi, "k" = k)
+                   "AIC.pi" = AIC.pi, "AICc.pi" = AICc.pi, "k" = k, "N" = nn)
 
     return(output)   
 

AICc_table_generation.R

@@ -4,7 +4,7 @@ require(tibble)
 
 ## -- For two variables: ------------------------------------
 
-AICc.table.2var <- function(sig.vars, control.var.char = NULL, c.var = 0, matrix.char, perm = 999, type = "AICc", method = "bray") {
+AICc.table.2var <- function(sig.vars, control.var.char = NULL, c.var = 0, matrix.char, perm, type = "AICc", method = "bray") {
 
     varcomb.2.AICc <- tibble(variables = rep("var.name", choose(length(sig.vars),2)),
                              AICc.values = rep(0),
@@ -13,7 +13,7 @@ AICc.table.2var <- function(sig.vars, control.var.char = NULL, c.var = 0, matrix
                              `Var Explnd` = rep(0),
                              Model = rep("model"))
 
-    if (is.character(control.var.char) == TRUE & c.var == 0) {c.var = 1} 
+    if (is.character(control.var.char) == TRUE & c.var == 0) {c.var = length(control.var.char)} 
 
     combo.list <- combn(x = sig.vars, m = 2, simplify = FALSE)
 
@@ -52,11 +52,12 @@ AICc.table.2var <- function(sig.vars, control.var.char = NULL, c.var = 0, matrix
 
     varcomb.2.AICc$`Delta AICc` <- varcomb.2.AICc$AICc.values -
         min(varcomb.2.AICc$AICc.values)
-    varcomb.2.AICc$`Relative Likelihood` <- exp((min(varcomb.2.AICc$AICc.values) -
-                                                     varcomb.2.AICc$AICc.values)/2)
+    varcomb.2.AICc$`Relative Likelihood` <-
+        exp(-.5 * (varcomb.2.AICc$AICc.values - min(varcomb.2.AICc$AICc.values)))
 
     # Relative likelihood compared with best model; see
     # https://en.wikipedia.org/wiki/Likelihood_function
+    # https://www.rdocumentation.org/packages/qpcR/versions/1.4-1/topics/akaike.weights 
 
 
 return(varcomb.2.AICc)
@@ -68,7 +69,7 @@ return(varcomb.2.AICc)
 ## -- For N variables: ---------------------------------------------------
 
 
-AICc.table.Nvar <- function(sig.vars, control.var.char = NULL, c.var = 0, matrix.char, perm = 999, n.var = 1, composite = FALSE, type = "AICc", method = "bray") {
+AICc.table.Nvar <- function(sig.vars, control.var.char = NULL, c.var = 0, matrix.char, perm, n.var = 1, composite = FALSE, type = "AICc", method = "bray") {
 
     if (n.var > length(sig.vars)) { stop("n.var greater than number of significant variables")}
 
@@ -123,8 +124,8 @@ AICc.table.Nvar <- function(sig.vars, control.var.char = NULL, c.var = 0, matrix
 if (composite == FALSE) {
     varcomb.N.AICc$`Delta AICc` <- varcomb.N.AICc$AICc.values -
         min(varcomb.N.AICc$AICc.values)
-    varcomb.N.AICc$`Relative Likelihood` <- exp((min(varcomb.N.AICc$AICc.values) -
-                                                     varcomb.N.AICc$AICc.values)/2)
+    varcomb.2.AICc$`Relative Likelihood` <-
+        exp(-.5 * (varcomb.2.AICc$AICc.values - min(varcomb.2.AICc$AICc.values)))
     }
 
 
@@ -149,7 +150,7 @@ AICc.table.all <- function(sig.vars, control.var.char = NULL, matrix.char, perm
 
         temp <- AICc.table.Nvar(sig.vars = control.var.char, control.var.char = NULL,
                                 matrix.char = matrix.char, n.var = 1, composite = TRUE,
-                                type = type, method = method)
+                                type = type, method = method, perm = perm)
 
         varcomb.all <- rbind(varcomb.all, temp)
 
@@ -161,7 +162,7 @@ AICc.table.all <- function(sig.vars, control.var.char = NULL, matrix.char, perm
 
         temp <- AICc.table.Nvar(sig.vars = sig.vars, control.var.char = control.var.char,
                                 matrix.char = matrix.char, n.var = i, composite = TRUE,
-                                type = type, method = method)
+                                type = type, method = method, perm = perm)
 
         varcomb.all <- rbind(varcomb.all, temp)