R/fuzzy_regression.R

In flex: Fuzzy Linear Squares Estimation with Explicit Formula (FLEX)

#' Fuzzy Linear Regression
#'
#' Fits a fuzzy linear regression model given fuzzified predictors and response variables.
#'
#' @param X_fuzzy A list of fuzzified predictor values.
#' @param Y_fuzzy A list of fuzzified response values.
#' @param p An integer specifying the number of predictors.
#' @param X_crisp Optional. The original crisp predictor matrix or data frame. Used to retrieve variable names. Default is \code{NULL}.
#'
#' @return A list object of class \code{fuzzy_lm} containing:
#' \item{Coefficients}{A data frame with estimated coefficients, standard errors, t-values, p-values, and significance stars.}
#' \item{Residuals}{The residuals from the fitted model.}
#' \item{Predictions}{The predicted fuzzified response values.}
#' \item{RSS}{The residual sum of squares.}
#' \item{R_squared}{The coefficient of determination (R-squared).}
#' \item{Sigma_squared}{The estimated variance of the residuals.}
#' \item{Degrees_of_Freedom}{The degrees of freedom for the model.}
#'
#' @examples
#' # Simulate complex data for fuzzy linear regression
#' set.seed(123)
#'
#' # Generate a dataset with 100 observations and 4 predictors
#' n <- 100
#' X_crisp <- data.frame(
#'   Age = round(runif(n, 20, 70)),            # Random ages between 20 and 70
#'   Income = round(runif(n, 20000, 120000)),  # Random incomes between 20k and 120k
#'   Education = round(runif(n, 10, 20)),      # Random years of education between 10 and 20
#'   Experience = round(runif(n, 1, 40))       # Random years of work experience between 1 and 40
#' )
#'
#' # Define true coefficients
#' beta <- c(5.0, 1.2, -0.5, 0.8, 0.05)  # Intercept and coefficients for the predictors
#'
#' # Generate the crisp response variable with noise
#' Y_crisp <- round(beta[1] + as.matrix(X_crisp) %*% beta[-1] + rnorm(n, mean = 0, sd = 50))
#'
#' # Fuzzify the predictor and response variables
#' X_fuzzy <- fuzzify_crisp_matrix(as.matrix(X_crisp), spread = 10)  # Larger spread for predictors
#' Y_fuzzy <- fuzzify_crisp_vector(Y_crisp, spread = 20)             # Larger spread for responses
#'
#' # Fit the fuzzy linear model
#' object <- fuzzy_lm(X_fuzzy, Y_fuzzy, p = 4, X_crisp = X_crisp)
#'
#' # Print the coefficients
#' print("Fuzzy Linear Model Coefficients:")
#' print(object$Coefficients)
#'
#' # Example residuals and predictions
#' print("Example Residuals:")
#' print(head(object$Residuals, 6))
#'
#' print("Example Predictions:")
#' print(head(object$Predictions, 6))
#'
#' @export
fuzzy_lm <- function(X_fuzzy, Y_fuzzy, p, X_crisp = NULL) {
  n <- length(Y_fuzzy)
  if (n != length(X_fuzzy)) stop("X_fuzzy and Y_fuzzy must have the same number of observations.")

  # Automatically retrieve variable names if X_crisp is provided
  if (!is.null(X_crisp)) {
    var_names <- c("Intercept", colnames(X_crisp))
  } else {
    var_names <- c("Intercept", paste0("X", 1:p))  # Default names if X_crisp is not provided
  }

  # Embedded Helper: Construct the fuzzy design matrix with sign adjustment
  fuzzy_x_matrix <- function(X_fuzzy, p, beta = NULL) {
    X_design <- matrix(list(), n, p + 1)
    for (i in 1:n) {
      X_design[[i, 1]] <- list(l = 1, x = 1, r = 1)  # Intercept
      for (j in 1:p) {
        l_xij <- X_fuzzy[[i]][[j]]$l
        x_xij <- X_fuzzy[[i]][[j]]$x
        r_xij <- X_fuzzy[[i]][[j]]$r

        # Apply sign-based adjustment if beta is provided
        if (!is.null(beta) && beta[j + 1] < 0) {
          X_design[[i, j + 1]] <- list(
            l = x_xij + (r_xij - x_xij),
            x = x_xij,
            r = x_xij - (x_xij - l_xij)
          )
        } else {
          X_design[[i, j + 1]] <- list(
            l = x_xij - (x_xij - l_xij),
            x = x_xij,
            r = x_xij + (r_xij - x_xij)
          )
        }
      }
    }
    return(X_design)
  }

  # Embedded Helper: Compute normal equations (XtX and XtY)
  fuzzy_normal_eq <- function(X_design, Y_fuzzy) {
    XtX <- matrix(0, p + 1, p + 1)
    XtY <- rep(0, p + 1)
    for (i in 1:n) {
      for (k in 1:(p + 1)) {
        XtY[k] <- XtY[k] + fuzzy_mults(X_design[[i, k]], Y_fuzzy[[i]])
        for (j in 1:(p + 1)) {
          XtX[k, j] <- XtX[k, j] + fuzzy_mults(X_design[[i, k]], X_design[[i, j]])
        }
      }
    }
    return(list(XtX = XtX, XtY = XtY))
  }

  # Initial design matrix without beta sign adjustment
  X_design <- fuzzy_x_matrix(X_fuzzy, p)
  eq <- fuzzy_normal_eq(X_design, Y_fuzzy)
  XtX_inv <- tryCatch(solve(eq$XtX), error = function(e) stop("XtX matrix is singular and cannot be inverted."))
  beta_hat <- XtX_inv %*% eq$XtY

  # Recalculate design matrix with beta sign adjustment
  X_design <- fuzzy_x_matrix(X_fuzzy, p, beta = beta_hat)
  eq <- fuzzy_normal_eq(X_design, Y_fuzzy)
  XtX_inv <- tryCatch(solve(eq$XtX), error = function(e) stop("XtX matrix is singular and cannot be inverted."))
  beta_hat <- XtX_inv %*% eq$XtY

  # Compute residuals
  Y_pred <- compute_pred(list(beta_hat = beta_hat), X_fuzzy)
  residuals <- compute_res(Y_fuzzy, Y_pred)

  # Compute RSS, sigma_squared, and R-squared
  rss <- sum(sapply(1:n, function(i) fuzzy_d_squared(Y_fuzzy[[i]], Y_pred[[i]])))
  df <- n - (p + 1)
  if (df <= 0) stop("Degrees of freedom are non-positive. Check your object or data.")
  sigma_squared <- rss / df
  mean_y <- list(
    l = mean(sapply(Y_fuzzy, `[[`, "l")),
    x = mean(sapply(Y_fuzzy, `[[`, "x")),
    r = mean(sapply(Y_fuzzy, `[[`, "r"))
  )
  total_ss <- sum(sapply(1:n, function(i) fuzzy_d_squared(Y_fuzzy[[i]], mean_y)))
  r_squared <- 1 - (rss / total_ss)

  # Compute standard errors, t-values, p-values, and significance stars
  se_beta <- sqrt(diag(sigma_squared * XtX_inv))
  t_values <- beta_hat / se_beta
  p_values <- 2 * pt(-abs(t_values), df = df)
  significance <- .add_significance_stars(p_values)

  # Create results data frame with variable names in the leftmost column
  Coefficients <- data.frame(
    Variable = var_names,
    Estimate = as.vector(beta_hat),
    `Std. Error` = se_beta,
    `t value` = as.vector(t_values),
    `Pr(>|t|)` = p_values,
    Significance = significance,
    check.names = FALSE
  )

  # Return the `fuzzy_lm` object
  return(structure(list(
    Coefficients = Coefficients,
    Residuals = residuals,
    Predictions = Y_pred,
    RSS = rss,
    R_squared = r_squared,
    Sigma_squared = sigma_squared,
    Degrees_of_Freedom = df
  ), class = "fuzzy_lm"))
}

#' Define generic for coefficients
#'
#' A generic function to retrieve coefficients from model objects.
#'
#' @param object The model object from which to extract coefficients.
#' @param ... Additional arguments (ignored).
#' @return A data frame of coefficients and related statistics.
#' @examples
#' # Simulate data and fit a fuzzy linear model
#' set.seed(123)
#' X_crisp <- matrix(round(runif(300, 2, 10)), nrow = 100, ncol = 3)
#' beta <- c(1.5, -0.8, 2.0)
#' Y_crisp <- round(X_crisp %*% beta + rnorm(100, mean = 0, sd = 1))
#' X_fuzzy <- fuzzify_crisp_matrix(X_crisp, spread = 1)
#' Y_fuzzy <- fuzzify_crisp_vector(Y_crisp, spread = 1)
#' object <- fuzzy_lm(X_fuzzy, Y_fuzzy, p = 3)
#'
#' # Extract coefficients
#' coefficients(object)
#' @export
coefficients <- function(object, ...) {
  UseMethod("coefficients")
}

#' Accessor for Coefficients
#' @param object An object of class `fuzzy_lm`.
#' @param ... Additional arguments (ignored).
#' @return A data frame of coefficients and statistics.
#' @examples
#' # Simulate data and fit a fuzzy linear model
#' set.seed(123)
#' X_crisp <- matrix(round(runif(300, 2, 10)), nrow = 100, ncol = 3)
#' beta <- c(1.5, -0.8, 2.0)
#' Y_crisp <- round(X_crisp %*% beta + rnorm(100, mean = 0, sd = 1))
#' X_fuzzy <- fuzzify_crisp_matrix(X_crisp, spread = 1)
#' Y_fuzzy <- fuzzify_crisp_vector(Y_crisp, spread = 1)
#' object <- fuzzy_lm(X_fuzzy, Y_fuzzy, p = 3)
#'
#' # Extract coefficients
#' coefficients(object)
#' @export
coefficients.fuzzy_lm <- function(object, ...) {
  stopifnot(inherits(object, "fuzzy_lm"))
  object$Coefficients
}

#' Define generic for residuals
#'
#' A generic function to retrieve residuals from model objects.
#'
#' @param object The model object from which to extract residuals.
#' @param ... Additional arguments (ignored).
#' @return A list of fuzzy residuals.
#' @examples
#' # Simulate data and fit a fuzzy linear model
#' set.seed(123)
#' X_crisp <- matrix(round(runif(300, 2, 10)), nrow = 100, ncol = 3)
#' beta <- c(1.5, -0.8, 2.0)
#' Y_crisp <- round(X_crisp %*% beta + rnorm(100, mean = 0, sd = 1))
#' X_fuzzy <- fuzzify_crisp_matrix(X_crisp, spread = 1)
#' Y_fuzzy <- fuzzify_crisp_vector(Y_crisp, spread = 1)
#' object <- fuzzy_lm(X_fuzzy, Y_fuzzy, p = 3)
#'
#' # Extract residuals
#' head(residuals(object))
#' @export
residuals <- function(object, ...) {
  UseMethod("residuals")
}

#' Accessor for Residuals
#'
#' @param object An object of class `fuzzy_lm`. The model object.
#' @param ... Additional arguments (currently ignored).
#' @return A list of fuzzy residuals.
#' @examples
#' # Simulate data and fit a fuzzy linear model
#' set.seed(123)
#' X_crisp <- matrix(round(runif(300, 2, 10)), nrow = 100, ncol = 3)
#' beta <- c(1.5, -0.8, 2.0)
#' Y_crisp <- round(X_crisp %*% beta + rnorm(100, mean = 0, sd = 1))
#' X_fuzzy <- fuzzify_crisp_matrix(X_crisp, spread = 1)
#' Y_fuzzy <- fuzzify_crisp_vector(Y_crisp, spread = 1)
#' object <- fuzzy_lm(X_fuzzy, Y_fuzzy, p = 3)
#'
#' # Extract residuals
#' head(residuals(object))
#' @export
residuals.fuzzy_lm <- function(object, ...) {
  stopifnot(inherits(object, "fuzzy_lm"))
  object$Residuals
}
#' Define generic for predictions
#'
#' @param object An object of class `fuzzy_lm`. The model object.
#' @param ... Additional arguments (currently ignored).
#' @return A list of fuzzy predictions.
#' @examples
#' # Simulate data and fit a fuzzy linear model
#' set.seed(123)
#' X_crisp <- matrix(round(runif(300, 2, 10)), nrow = 100, ncol = 3)
#' beta <- c(1.5, -0.8, 2.0)
#' Y_crisp <- round(X_crisp %*% beta + rnorm(100, mean = 0, sd = 1))
#' X_fuzzy <- fuzzify_crisp_matrix(X_crisp, spread = 1)
#' Y_fuzzy <- fuzzify_crisp_vector(Y_crisp, spread = 1)
#' object <- fuzzy_lm(X_fuzzy, Y_fuzzy, p = 3)
#'
#' # Extract predictions
#' head(predictions(object))
#' @export
predictions <- function(object, ...) {
  UseMethod("predictions")
}

#' Accessor for Predictions
#'
#' @param object An object of class `fuzzy_lm`. The model object.
#' @param ... Additional arguments (currently ignored).
#' @return A list of fuzzy predictions.
#' @examples
#' # Simulate data and fit a fuzzy linear model
#' set.seed(123)
#' X_crisp <- matrix(round(runif(300, 2, 10)), nrow = 100, ncol = 3)
#' beta <- c(1.5, -0.8, 2.0)
#' Y_crisp <- round(X_crisp %*% beta + rnorm(100, mean = 0, sd = 1))
#' X_fuzzy <- fuzzify_crisp_matrix(X_crisp, spread = 1)
#' Y_fuzzy <- fuzzify_crisp_vector(Y_crisp, spread = 1)
#' object <- fuzzy_lm(X_fuzzy, Y_fuzzy, p = 3)
#'
#' # Extract predictions
#' head(predictions(object))
#' @export
predictions.fuzzy_lm <- function(object, ...) {
  stopifnot(inherits(object, "fuzzy_lm"))
  object$Predictions
}

#' Summary for Fuzzy Linear Regression
#'
#' @param object An object of class `fuzzy_lm`. The model object.
#' @param ... Additional arguments (currently ignored).
#' @return Prints a summary of the fuzzy linear regression results.
#' @examples
#' # Simulate data and fit a fuzzy linear model
#' set.seed(123)
#' X_crisp <- matrix(round(runif(300, 2, 10)), nrow = 100, ncol = 3)
#' beta <- c(1.5, -0.8, 2.0)
#' Y_crisp <- round(X_crisp %*% beta + rnorm(100, mean = 0, sd = 1))
#' X_fuzzy <- fuzzify_crisp_matrix(X_crisp, spread = 1)
#' Y_fuzzy <- fuzzify_crisp_vector(Y_crisp, spread = 1)
#' object <- fuzzy_lm(X_fuzzy, Y_fuzzy, p = 3)
#'
#' # Summarize the model
#' summary(object)
#' @export
summary.fuzzy_lm <- function(object, ...) {
  stopifnot(inherits(object, "fuzzy_lm"))

  coefficients <- object$Coefficients
  if (!is.data.frame(coefficients)) {
    stop("Coefficients must be a data frame.")
  }

  # Clean up column names for better readability
  colnames(coefficients) <- c("Variable", "Estimate", "Std. Error", "t-value", "Pr(>|t|)", "Significance")

  cat("\nFuzzy Linear Regression Summary:\n")
  print(coefficients)

  cat("\nModel Statistics:\n")
  cat("Residual sum of squares (RSS):", round(object$RSS, 4), "\n")
  cat("R-squared:", round(object$R_squared, 4), "\n")
  cat("Degrees of freedom:", object$Degrees_of_Freedom, "\n")
}