R/lba.mle.R

In lba: Latent Budget Analysis for Compositional Data

lba.mle <- function(obj     , 
                    A       , 
                    B       , 
                    K       , 
                    tolG    , 
                    tolA    , 
                    tolB    ,
                    itmax.unide,
                    itmax.ide,
                    trace.lba,
                    toltype , 
                    what,
                    ...) 
{
 # The default toltype is 'all'
 I <- nrow(obj)
 J <- ncol(obj)

 obj[obj==0] <- 1e-4

 P <- obj/rowSums(obj)
 #
 # ------------------------------------------------------------------------------
 # creating random generated values for alpha(k|i)
 if (is.null(A)) {
  A <- matrix(runif(I * K), 
              ncol = K)
  A <- A/rowSums(A)
 }
 #
 # ------------------------------------------------------------------------------
 # creating random generated values for beta(j|k)
 if (is.null(B)) {
  B <- matrix(runif(J * K), 
              ncol = K)
  B <- t(t(B)/colSums(B))
 }

 iter_unide <- 0  # counter of interactions numbers

 mle_func <- function(obj,
                      I,
                      J,
                      K,
                      A,
                      B){

  pip <- rowSums(obj)/sum(obj) #this is pi+
  pjki <- rep(0,I*J*K) #this will become pijk/pi+ see van der Ark page 80
  m <- 0
  # this makes i=1,j=1,k=1; i=2,j=1,k=1; i=3,j=1,k=1...; i=2,j=1,k=1 and so on.
  for(k in 1:K) for(j in 1:J) for(i in 1:I) { 
   m <-  m +1
   pjki[m] <- A[i,k]*B[j,k] 
  }

  mi <- matrix(pjki, I)                                          
  pijk <- mi*pip #this is pijk see van der Ark page 80 
  nijk <- as.vector(pijk)*sum(obj) 
  val_function <- -sum(nijk * log(pjki)) # -loglikelihood function

 } 
 repeat {
  iter_unide <- iter_unide + 1

  if(trace.lba){
   cat('Unidentified iteration: ',iter_unide,'\n')
   cat('fval = ', signif(mle_func(obj,I,J,K,A,B),4), '\n') 
  }  
  #
  # -------------------------------------------------------------------------------
  # previous G2 to be used in |G2 - G2a|<tol
  qij <- A %*% t(B)
  kij <- obj/rowSums(obj)
  G2a <- 2 * sum(obj * log(kij/qij))
  akia <- A
  bjka <- B

  #
  # -------------------------------------------------------------------------------
  ini <- list()
  # ----------------------------------------------------------------------------
  # STARTING THE (expectation)E-STEP
  # -------------------------------------------------------------------------------
  for (k in 1:K) ini[[k]] <- c(A[, k], 
                               B[, k])
  #
  # ------------------------------------------------------------------------------
  # initial values list with vectors a1(I),b1(J);... ak(I),bk(J); ... aK(I),
  # bK(J) where a's are alfa values and b's are beta values. The lenghts of the
  # a's are I and of the b's are J. The elements of the list are c(a1,b1), ...
  # c(ak,bk),... c(aK,bK). 

  pijk <- lapply(ini, function(xx) outer(xx[1:I], 
                                         xx[(I + 1):(I + J)]) * (rowSums(obj)/sum(obj)))
  # The elements of the list pijk are matrices (I x J) for fixed k and the number
  # of elements is K. They are alpha(k|i)beta(j|k)*pi+ = pijk|i*pi+

  mpijk <- sapply(pijk, 
                  function(x) x, simplify = TRUE)
  # Transform pijk into a matrix I*JxK. The columns are the colmuns of each
  # element (matrix) one after the other, of the first then the second and so on.

  sij <- rowSums(mpijk)
  # The sum over k of alpha(k|i)beta(j|k), vector of lenght IxJ

  lnijk <- lapply(pijk, 
                  function(xx) obj * xx)
  # The product of N(ij)*pijk. It is a list of K matrices IxJ each

  nijk <- lapply(lnijk, 
                 function(x) x/sij)
  # The list containing the values of n(hat)(ijk). Each element is a matrix
  # IxJ and the number of elements of the list is K.

  #-------------------------------------------------------------------------------
  #-------------------------------------------------------------------------------
  #Starting computation for maximization step
  #-------------------------------------------------------------------------------        

  sjnijk <- sapply(nijk, 
                   function(x) rowSums(x), simplify = TRUE)
  # sum over j for i and k fixed. sjnijk is a IxK matrix

  sjknijk <- apply(sjnijk, 
                   1, 
                   function(x) sum(x))
  # sum over j and k for i fixed. sjknijk is a vector of lenght I

  sinijk <- sapply(nijk, 
                   function(x) colSums(x), simplify = TRUE)
  # sum over i for j and k fixed. sinijk is a JxK matrix

  sijnijk <- apply(sinijk, 
                   2, 
                   function(x) sum(x))
  # sum over i and j and k fixed. sjknijk is a vector of lenght K
  # -----------------------------------------------------------------------------
  # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  # maximization step
  # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  A <- sjnijk/sjknijk
  # next values alphahat1 to be used in the next step of the iteration

  B <- t(t(sinijk)/sijnijk)
  # next values betahat1 to be used in the next step of the iteration

  # aki; bjk

  pij <- A %*% t(B)
  # Values of pihat(j|i) = the sum over k of alphahat(k|i)*betahat(j|k)

  G2 <- 2 * sum(obj * log(obj/(pij * rowSums(obj))))
  chi2 <- sum(((obj - pij * rowSums(obj))^2)/obj)
  # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  at <- max(abs(A - akia))
  bt <- max(abs(B - bjka))
  Gt <- abs(G2 - G2a)

  #
  # ------------------------------------------------------------------------------
  if(iter_unide > itmax.unide){
   warning("maximum number of the unidentified iteractions exceeded" )
   break 
  }

  if (toltype == "all" & Gt < tolG & at < tolA & bt < tolB) {
   break  #default
  } else 
   if (toltype == "G2" & Gt < tolG) {
    break
   } else 
    if (toltype == "ab" & at < tolA & bt < tolB) {
     break
    }

 }  #closing the repeat

 val_func <- mle_func(obj,I,J,K,A,B) # -loglikelihood function


 if(K==1){

  Aoi <- A

  Boi <- B

  iter_ide  <- 'not applicable'

 } else {

  AB_outerAB_inner <- identif(A,
                              B,
                              P,
                              K,
                              trace.lba,
                              itmax.ide,
                              what)

  Aoi <- AB_outerAB_inner[[1]]

  Boi <- AB_outerAB_inner[[2]]

  iter_ide <- AB_outerAB_inner[[3]]  

 }

 pij <- A %*% t(B)

 residual <- P - pij

 pip <- rowSums(obj)/sum(obj) #this is pi+ 

 aux_pk <- pip %*% Aoi

 pk <- matrix(aux_pk[order(aux_pk,
                       decreasing = TRUE)],
              ncol = dim(aux_pk)[2])

 Aoi <- matrix(Aoi[,order(aux_pk,
                   decreasing = TRUE)],
               ncol = dim(aux_pk)[2])
 Boi <- matrix(Boi[,order(aux_pk,
                   decreasing = TRUE)],
               ncol = dim(aux_pk)[2])

 A <- matrix(A[,order(aux_pk,
               decreasing = TRUE)],
             ncol = dim(aux_pk)[2])
 B <- matrix(B[,order(aux_pk,
               decreasing = TRUE)],
             ncol = dim(aux_pk)[2])

 colnames(pk) <- colnames(Aoi) <- colnames(Boi) <- colnames(A) <- colnames(B) <- paste('LB',1:K,sep='')

 rownames(Aoi) <- rownames(A) <- rownames(P)
 rownames(Boi) <- rownames(B) <- colnames(P)

 rescB <- rescaleB(x=obj,
                   Aoi,
                   Boi)

 colnames(rescB) <- colnames(B)
 rownames(rescB) <- rownames(B) 

 res <- list(P,       
             pij,
             residual,
             A, 
             B,
             Aoi,
             Boi,
             rescB,
             pk,
             val_func,
             iter_unide,
             iter_ide)

 names(res) <- c('P',
                 'pij',
                 'residual',
                 'A',
                 'B',
                 'Aoi',
                 'Boi',
                 'rescB',
                 'pk',
                 'val_func',
                 'iter_unide',
                 'iter_ide')

 class(res) <- "lba.mle" 

 invisible(res)            
}