R/cm_net.R

In migest: Methods for the Indirect Estimation of Bilateral Migration

#' Conditional maximization routine for the indirect estimation of origin-destination-type migration flow tables with known net migration totals.
#'
#' The \code{cm_net} function finds the maximum likelihood estimates for fitted values in the log-linear model:
#' \deqn{\log y_{ij} = \log \alpha_{i} + \log \alpha_{i}^{-1} + \log m_{ij} }
#' 
#' @param net_tot Vector of net migration totals to constrain the sum of the imputed cell row and columns. Elements must sum to zero.
#' @param m Array of auxiliary data. By default, set to 1 for all origin-destination-migrant typologies combinations. 
#' @param tol Numeric value for the tolerance level used in the parameter estimation.
#' @param maxit Numeric value for the maximum number of iterations used in the parameter estimation.
#' @param verbose Logical value to indicate the print the parameter estimates at each iteration. By default \code{FALSE}.
#' @param alpha0 Vector of initial estimates for alpha
#'
#' @return
#' Conditional maximisation routine set up using the partial likelihood derivatives. The argument \code{net_tot} takes the known net migration totals.
#' The user must ensure that the net migration totals sum globally to zero.
#' 
#' Returns a \code{list} object with
#' \item{mu }{Array of indirect estimates of origin-destination matrices by migrant characteristic}
#' \item{it }{Iteration count}
#' \item{tol }{Tolerance level at final iteration}
#' @author Guy J. Abel, Peter W. F. Smith
#' @export
#'
#' @examples
#' m <- matrix(data = 1:16, nrow = 4)
#' # m[lower.tri(m)] <- t(m)[lower.tri(m)]
#' addmargins(m)
#' sum_net(m)
#' 
#' y <- cm_net(net_tot = c(30, 40, -15, -55), m = m)
#' addmargins(y$n)
#' sum_net(y$n)
#' 
#' m <- matrix(data = c(0, 100, 30, 70, 50, 0, 45, 5, 60, 35, 0, 40, 20, 25, 20, 0),
#'             nrow = 4, ncol = 4, byrow = TRUE,
#'             dimnames = list(orig = LETTERS[1:4], dest = LETTERS[1:4]))
#' addmargins(m)
#' sum_net(m)
#' 
#' y <- cm_net(net_tot = c(-100, 125, -75, 50), m = m)
#' addmargins(y$n)
#' sum_net(y$n)
cm_net <- function(net_tot = NULL, m = NULL, tol = 1e-06, maxit = 500, verbose = TRUE,
                   alpha0 = rep(1, length(net_tot))) {
   # net_tot = c(-100, 125, -75, 50); m = NULL; tol = 1e-06;  maxit = 500; verbose = TRUE
   # m = matrix(data = 1:16, nrow = 4); net_tot = c(30, 40, -15, -55); tol = 1e-06;  maxit = 500; verbose = TRUE; alpha0 = rep(1, length(net_tot))
   R <- unique(c(dim(m), length(net_tot)))
   if (length(R) != 1)
      stop("The m matrix must be square and with the same dimensions as the length of net total vector (net_tot).")
   if (sum(net_tot) %% 1 > tol)
      message("Convergence will not be obtained as net_tot does not sum to zero.")
   
   #set up offset
   if (is.null(m))
      m <- matrix(1, nrow = R, ncol = R)
   if (is.null(dimnames(m)))
      dimnames(m) <- list(orig = LETTERS[1:R], dest = LETTERS[1:R])
   
   alpha <- alpha0
   it <- 0;  
   d_max <- tol * 2
   
   while (d_max > tol & it < maxit) {
      if (verbose == TRUE & (it < 20 | it %% 10 == 0)){
         cat("iteration:", it, "\n")
         cat("alpha parameters:", alpha, "\n")
         cat("\n")
      }
      alpha_old <- alpha
      
      # m <- s1/sum(s1)
      for(i in 1:R){
         # peter's code
         # (-net_tot[h] + sqrt(net_tot[h]^2 + 4 * sum(1/alpha[, (i-1)] * m[h,]) * sum(alpha[,(i-1)] * m[,h]))) / (2 * sum(1/alpha[, (i-1)] * m[h,]))
         # (-net_tot[i] + sqrt(net_tot[i]^2 + 4 * sum(1/alpha * m[i,]) * sum(alpha * m[,i]))) / (2 * sum(1/alpha * m[i,]))
         p <- quadratic_eqn(a = sum(1/alpha_old * m[i,]), 
                            b = net_tot[i], 
                            c = -sum(alpha_old * m[,i]))
         p <- p[p>0]
         if(length(p) == 0)
           p <- 1
         if(is.infinite(p) | is.na(p) | is.nan(p))
           p <- 1
         alpha[i] <- p
      }
      d_max <- max(abs(alpha_old - alpha))
      it <- it + 1
   }
   f <- diag(alpha) %*% m %*% diag(1/alpha)
   dimnames(f) <- dimnames(m)
   return(list(n = f, theta = c(alpha = alpha)))
}