R/mnlogit.R

In philswatton/logitPlus: What the Package Does (One Line, Title Case)

#' Multinomial Logistic Regression
#'
#' An R implementation of the multinomial logistic regression model. Allows specification
#' of random effects if desired. Allows use of a multinomial fractional response rather than
#' a nominal variable if desired.
#'
#' @param data Dataframe containing relevant independent variables and optionally dependent variable
#' @param y Optional matrix of responses. Should either be binary or rows should sum to 1
#' @param formula Model formula including lme4 syntax. If y is supplied, LHS should be 'y', i.e. y ~ RHS. Currently, only random intercepts, i.e. y ~ x + (1 | z).
#' @param starts Ignored: not yet implemented
#' @param est Ignored: currently only \code{MSL}.
#' @param control Currently only minimally implemented for \code{est = "MSL"}. Named list containing '\code{draws}' which specifies the number of Halton draws for the random effects. Only used if random effects are specified.
#' @param weights Ignored: not yet implemented
#' @param na.action Ignored: not yet implemented
#' @param contrasts Ignored: not yet implemented
#' @param offset Ignored: not yet implemented
#'
#' @return
#'
#' @references
#' \insertRef{bhat2001}{logitPlus}
#'
#' \insertRef{train2003}{logitPlus}
#'
#' \insertRef{haan2006}{logitPlus}
#'
#' \insertRef{becker2014}{logitPlus}
#'
#' @export
mnlogit <- function(data, y=NULL, formula,
                    starts=NULL, est="MSL", control=list(draws=100),
                    weights=NULL, na.action=NULL, contrasts=NULL, offset=NULL) {

  # Input validation will go here
  # if (!("draws" %in% names(control))) stop("If est=MSL, the number of draws must be specified") #will need to change this once options other than

  # Data preparation
  mdata <- logitPlus:::form_to_data(data=data, y=y, formula=formula)
  y <- mdata$y
  X <- mdata$X
  has_REs <- !is.null(mdata$Z)
  if (has_REs) Z <- mdata$Z

  n <- nrow(X)
  M <- ncol(y)
  K <- ncol(X)
  if (has_REs) A <- length(Z) #number of random effects to estimate

  # Calculate starting values - eventually need to implement the possibility of user-provided starting values
  B <- matrix(rep(numeric(K), M-1), K, M-1) #each column is one outcome's need to work out offsets
  if (has_REs) {
    a <- vector("list", A) #populate this list with unique elements
    nelem <- (M-1)^2 #dimensions for each variance-covariance matrix
    nuniq <- ((M-1) * ((M-1) + 1))/2 #number of *unique* components on each variance-covariance matrix

    for (i in 1:A) {
      var <- rep(exp(0.5), M-1)#numeric(M-1) #the diagonal is M-1 long
      cov <- rep(logitPlus::tanh(0.3), nuniq - (M-1))#numeric(nuniq - (M-1)) #the rest of the unique values are covariances

      a[[i]] <- list(var = var, cov = cov)
    }
  }

  # Halton sequences - for use in MSL estimation of REs
  if (has_REs) {
    R <- control$draws
    e <- vector("list", A)

    for (i in 1:A) {
      Q <- ncol(Z[[i]])
      G <- (R * Q) + 10 #Bhat 2001

      e[[i]] <- logitPlus::halton(draw=G,
                                  dim=(M-1),
                                  scramble=T,
                                  randomise=F,
                                  normal=T)
      e[[i]] <- e[[i]][11:nrow(e[[i]]),] #filter out first 10
      e[[i]] <- lapply(1:Q, function(i) temp2[((i-1)*R + 1:R),])
    }


#
#     sim <- function() {
#       logitPlus::halton(draw=G,
#                         dim=(M-1),
#                         scramble=T,
#                         randomise=F,
#                         normal=T)
#     }
  }

  # Construct the likelihood function to optimise - defaults to ML if no REs
  construct_ll_func <- function(has_REs=has_REs, X=X, y=y, Z=Z, e=e, n=n, K=K, M=M, A=A, draws=control$draws) {

    llfunc <- function(param = list(B, a))

      if (has_REs) {

        VC <- vector("list", A) #this will contain the complete variance-covariance matrices
        W <- vector("list", A)
        U <- vector("list", A) #this will contain the actual values of the random effect



        for (i in 1:A) {
          VC[[i]] <- logitPlus::symmetric(a[[i]]$var, a[[i]]$cov)
          W[[i]] <- t(chol(VC[[i]]))

        }
      }

      for (i in 1:nrow(X)) {
        for (j in 1:M) {
          for (r in 1:draws) {

          }
        }
      }
    return(llfunc)
  }

  llfunc <- construct_ll_func(has_REs)


  chol(VC[[i]])

  # Quasi-random draws to construct simulated REs

  # Estimation


}