R/ML_estimation.R

In adzemski/netprobitFE: Network probit with fixed effects

## This R package implements the methods proposed in
## Dzemski, Andreas: An empirical model of dyadic link formation in
## a network with unobserved heterogeneity, Review of Economics and Statistics, forthcoming

## Copyright (C) 2018  Andreas Dzemski

## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <https://www.gnu.org/licenses/>.

#' ML estimation of stage 1
#'
#' @param model list describing the model (generated by define_model)
#' @param start optional start values for theta (currently ignored)
#' @param robust flag, if TRUE don't use speedglm
#'
#' @return fitted model
#' @export
#' @import Matrix
#' @import data.table
#' @importFrom speedglm speedglm.wfit
MLE_stage1 <- function(model, start = NULL, robust = FALSE) {

  if (!robust)
    fitted <- speedglm.wfit(model$links[[model$col_outcome]], model$design_matrix_MLE,
                            family = binomial(link = probit))
  else
    fitted <- glm.fit(model$design_matrix_MLE, model$links[[model$col_outcome]],
                      family = binomial(link = probit))

  coefv <- coef(fitted)
  # back out the sender effect of node 1
  gammaS_1 <- setNames(sum(coefv[sprintf("gammaR_%d", seq_len(model$N))]) -
                         sum(coefv[sprintf("gammaS_%d", seq_len(model$N))][-1]), "gammaS_1")
  # add to estimated coefficients
  coefv <- c(coefv, gammaS_1)

  # write to parameter vectors
  theta_hat <- coefv[model$X_names]
  gammaS_hat <- coefv[sprintf("gammaS_%d", seq_len(model$N))]
  names(gammaS_hat) <- gsub('gammaS_(.*)', '\\1', names(gammaS_hat))
  gammaR_hat <- coefv[sprintf("gammaR_%d", seq_len(model$N))]
  names(gammaR_hat) <- gsub('gammaR_(.*)', '\\1', names(gammaR_hat))

  # return results
  list(theta_hat_MLE = theta_hat, gammaS_hat = gammaS_hat, gammaR_hat = gammaR_hat,
       model = model)
}

#' ML estimation stage 2
#'
#' @param fit_MLE_stage1 returned list from stage 1
#' @param gammaS sender effects
#' @param gammaR receiver effects
#' @param method maximize log likelihood ("mle") or search for root of gradient ("gradient")
#' @param rho_start starting value of rho (reciprocity parameter)
#' @param kappa bound rho away from -1 and 1 (see paper)
#'
#' @return fitted model
#' @export
#'
#' @import data.table
MLE_stage2 <- function(fit_MLE_stage1,
                       gammaS = fit_MLE_stage1$gammaS_hat,
                       gammaR = fit_MLE_stage1$gammaR_hat,
                       method = 'mle', rho_start = 0, kappa = 0.02) {

  long <- copy(fit_MLE_stage1$model$links)

  add_ystar(fit_MLE_stage1$model, fit_MLE_stage1$theta_hat_MLE,
            fit_MLE_stage1$gammaS_hat, fit_MLE_stage1$gammaR_hat, long)

  wide <- long_to_wide(long, fit_MLE_stage1$model$col_sender, fit_MLE_stage1$model$col_receiver,
                       nodes = seq_len(fit_MLE_stage1$model$N))
  wide[, Z_ij := Y_ij*Y_ji]

  if (tolower(method) == "mle")
    c(fit_MLE_stage1, MLE_stage2_mle(wide, -(1 - kappa), 1 - kappa, rho_start))
  else if (tolower(method) == "gradient")
    c(fit_MLE_stage1, MLE_stage2_gradient(wide, -(1 - kappa), 1 - kappa, rho_start))
  else
    stop(sprintf("Invalid method to solve log likelihood. You provided: %s", method))
}

#' Stage 2: solve likelihood directly
#'
#' @param wide network data in upper-triangular wide format
#' @param bound_lower lower bound on rho
#' @param bound_upper upper bound on rho
#' @param rho_start starting value (default = 0)
#'
#' @return fitted model
#'
#' @importFrom stats4 mle
MLE_stage2_mle <- function(wide, bound_lower, bound_upper, rho_start = 0) {

  comp_log_likelihood <- function(wide, rho) {

    wide[, r_ij := compute_r(Ystar_ij, Ystar_ji, rho)]
    wide[, lk_ij := Z_ij * log(r_ij) + (1 - Z_ij) * log(1 - r_ij)]

    wide[, sum(lk_ij)]
  }

  fit <- mle(function(rho) - comp_log_likelihood(wide, rho), start = list(rho = rho_start),
             method = 'L-BFGS-B', lower = bound_lower, upper = bound_upper)

  list(rho_hat_MLE = fit@coef, se_rho_MLE = sqrt(as.vector(fit@vcov)))
}

#' Stage 2: solve likelihood directly
#'
#' @param wide network data in upper-triangular wide format
#' @param bound_lower lower bound on rho
#' @param bound_upper upper bound on rho
#' @param rho_start starting value (default = 0)
#' @param inital_search_radius intial search radius around starting value (default = 0.4)
#' @param convergence_tol convergence tolerance
#'
#' @return fitted model
#'
#' @importFrom stats4 mle
MLE_stage2_gradient <- function(wide, bound_lower, bound_upper, rho_start = 0, inital_search_radius = 0.4,
                                convergence_tol = sqrt(.Machine$double.eps)) {

  score_rho <- function(ds, rho) {

    if (rho > bound_upper | rho < bound_lower) {
      rho <- min(bound_upper, max(bound_lower, rho))
    }

    ds[, r_ij := compute_r(Ystar_ij, Ystar_ji, rho)]
    ds[, J_ij := compute_J(Ystar_ij, Ystar_ji, rho, r_ij)]

    ds[, sum(J_ij*(Z_ij - r_ij))]
  }

  search_interval <- c(max(rho_start - inital_search_radius/2, bound_lower),
                       min(rho_start + inital_search_radius/2, bound_upper))

  root_gradient <- uniroot(function(rho) score_rho(wide, rho), interval = search_interval, extendInt = 'downX',
                           trace = 0, tol = convergence_tol, maxiter = 1e2)

  rho_hat <- root_gradient$root

  list(rho_hat_MLE = rho_hat, se_rho_MLE = as.numeric(NaN))
}