In BENMMI: Benthic Multi-Metric Index

Which MMI correlates best with pressure?

to_log("INFO", "Entering subsection 'Which MMI correlates best with pressure?'...")

number_of_combinations <- c(1L, 3L, 7L)[number_of_indicators]

Multimetric indices are constructed by fitting all r number_of_combinations possible weighted linear combinations of the set of indices r toString(paste0(indicators, "*")) to the pressure.

Our general model is

MMI^* = b₀ + b₁ × PRESSURE

where MMI^* is the r if (tolower(settings$legendtext) == "eqr") {"EQR of the"} else {"normalized"} multimetric index, PRESSURE are the values in the PRESSURE column of the input file, and b₀ and b₁ are the intercept and the slope respectively.

MMI^* is a weighted linear combination of r if (tolower(settings$legendtext) == "eqr") {"the ecological quality ratios of the"} else {"normalized values"} of the selected indices:

r paste(paste0("w<sub>", 1:number_of_indicators, "</sub>"), paste0(indicators, "*"), collapse = " + ", sep = " × ").

where weights r toString(paste0("w<sub>", 1:number_of_indicators, "</sub>")) are non-negative and sum to one.

For each of the r number_of_combinations possible combinations of r title_text_plur, the weights are optimized mathematically in order to maximize the precision of the model mentioned above. Duplicated models (with approximately the same set of weights) have been removed.

# initialize weight matrix and add single metric weights
W <- diag(number_of_indicators)
dimnames(W) <- list(NULL, indicators_EQR)
W <- rbind(W, 
    matrix(
        data = 0, 
        nrow = number_of_combinations - number_of_indicators, 
        ncol = number_of_indicators
    )
)



# bimetric indicators
i <- number_of_indicators
if (number_of_indicators >= 2L) {
    f_obj <- function(w, id) {
            (sprintf(
                "I(%e * %s + %e * %s) ~ PRESSURE",
                w,     indicators_EQR[id[1]],
                1 - w, indicators_EQR[id[2]]
            ) %>%
            as.formula %>%
            lm(data = d_ind) %>%
            summary)$adj.r.squared  # maximum needed
    }
    opt <- optimize(f = f_obj, interval = c(0, 1), maximum = TRUE, id = c(1, 2))
    i <- i + 1L
    W[i, 1] <-     opt$maximum
    W[i, 2] <- 1 - opt$maximum
}


# trimetric indicators
if (number_of_indicators == 3L) {
    opt <- optimize(f = f_obj, interval = c(0, 1), maximum = TRUE, id = c(1, 3))
    i <- i + 1L
    W[i, 1] <-     opt$maximum
    W[i, 3] <- 1 - opt$maximum
    opt <- optimize(f = f_obj, interval = c(0, 1), maximum = TRUE, id = c(2, 3))
    i <- i + 1L
    W[i, 2] <-     opt$maximum
    W[i, 3] <- 1 - opt$maximum
}

if (number_of_indicators == 3L) {
    f_obj <- function(w) {
        if (any(w < 0)) {
            return(Inf)
        }
        w <- w / sum(w)
        (sprintf(
            "I(%e * %s + %e * %s + %e * %s) ~ PRESSURE",
            w[1], indicators_EQR[1],
            w[2], indicators_EQR[2],
            w[3], indicators_EQR[3]
        ) %>% 
        as.formula %>%
        lm(data = d_ind) %>% 
        summary)$adj.r.squared # maximum needed
    }
    opt <- optim(
        par = c(1, 1, 1), 
        fn = f_obj, 
        method = "Nelder-Mead", 
        control = list(maxit = 1000, fnscale = -1) # maximization
    )
    has_converged <- opt$convergence == 0L
    if (!has_converged) {
        stop("Can't find an optimum set of weights", call. = FALSE)
    }
    i <- i + 1L
    W[i, ] <- opt$par / sum(opt$par)
}