R/my_linDA.R

In nclSLR: The nclSLR package.

my_linDA <-
  function(X, y, learn, test, prior, prob)
  {
    # Perform a lienar discriminant analysis
    # X: matrix or data.frame with explanatory variables
    # y: vector or factor with group membership
    # learn: vector of learning observations
    # test: vector of testing observations
    # prior: vector of prior proportions
    # prob: logical indicating results in proability terms

    # how many observations
    n = nrow(X[learn,,drop=FALSE])
    ntest = length(test)
    # how many groups
    ng = nlevels(y[learn])
    glevs = levels(y[learn])
    # how many obs in each group
    nobs_group = as.vector(table(y[learn]))
    # group means
    GM = my_groupMeans(X[learn,,drop=FALSE], y[learn])
    # within-class covariance matrix
    W = my_withinCov(X[learn,,drop=FALSE], y[learn])
    # inverse of Within cov matrix
    W_inv = solve(W)

    # constant terms of classification functions
    cons = rep(0, ng)
    # coefficients of classification functions
    Betas = matrix(0, nrow(W_inv), ng)
    for (k in 1:ng)
    {
      cons[k] = -(1/2) * GM[k,] %*% W_inv %*% GM[k,] + log(prior[k])
      Betas[,k] = t(GM[k,]) %*% W_inv
    }
    # Fisher's Discriminant Functions
    FDF = rbind(cons, Betas)
    rownames(FDF) = c("constant", colnames(X))
    colnames(FDF) = glevs

    # matrix of constant terms
    A = matrix(rep(cons,ntest), ntest, ng, byrow=TRUE)
    # apply discrim functions
    Disc = X[test,,drop=FALSE] %*% Betas + A

    # probability values
    if (prob) {
      # exponential
      Disc <- 1 - exp( -(Disc - apply(Disc, 1, min, na.rm=TRUE)))
      # predicting classes
      pred = Disc / drop(Disc %*% rep(1, ng))
      # predicted class
      pred_class = factor(max.col(pred), levels=seq_along(glevs), labels=glevs)
    } else {
      # predicted class
      pred = apply(Disc, 1, function(u) which(u == max(u)))
      names(pred) = NULL
      # assign class values
      pred_class = factor(pred, levels=seq_along(glevs), labels=glevs)
    }
    dimnames(Disc) = list(rownames(X[test,,drop=FALSE]), glevs)
    # confusion matrix
    conf = table(original=y[test], predicted=pred_class)
    # results
    res = list(FDF=FDF, conf=conf, Disc=Disc, pred_class=pred_class)
    res
  }