blin_mle: Estimate the BLIN model using maximum likelihood estimator

Description

This function estimates the bipartite logitudinal influence network (BLIN) model $Y_t = A^T \sum_{k=1}^{lag} Y_{t-k} + \sum_{k=1}^{lag} Y_{t-k} B + X_t \beta + \tau E_t$ using maximum likelihood estimator.

Usage

blin_mle(Y, X = NULL, type = "full", lag = 1, rankA = NULL,
  rankB = rankA, maxit = 1000, tol = 1e-08, init = "I",
  sigma_init = 1, verbose = FALSE, calcses = FALSE, randseed = NA)

Arguments

Response 3-mode array.

Optional 4-mode array of covariates, defaults to no covariates.

type

Optional string specifying BLIN model type: full, reduced_rank, or sparse. Defaults to full.

lag

Optional numeric specifying autoregressive lag in model, defaults to 1.

rankA

Optional numeric rank of influence network matrix $A$ for reduced rank model type, defaults to full rank.

rankB

Optional numeric rank of influence network matrix $B$, defaults to rank of $A$.

maxit

Optional numeric maximum number of iterations for full and reduced rank block coordinate descents, defaults to 1e3.

tol

Optional numeric convergence tolerance for full and reduced rank block coordinate descents, defaults to 1e-8.

init

Optional string specifying initialization type for full and reduced rank block coordinate descents, defaults to "I", identity for $A$ and $B$. Also allows "random" for random initialization of $A$ and $B$.

sigma_init

Optional numeric standard deviation for random initialization of $A$ and $B$ in full and reduced rank block coordinate descents, defaults to 1.

verbose

Optional logical specifying whether progress should be printed out (TRUE) or not (FALSE). Defaults to FALSE.

calcses

Optional logical specifying whether standard errors should be calculated (TRUE) or not (FALSE). Defaults to FALSE. Only standard errors for the full BLIN model are implemented.

randseed

Optional numeric specifying seed for random ininitialization of $A$ and $B$ in full and reduced rank block coordinate descents, defaults to NA (no seed set).

Value

fit

A blin object containing summary information.

Details

This function estimates the continuous BLIN model, $$Y_t = A^T Y_{t-1} + Y_{t-1} B + X_t \beta + \tau E_t$$, where $ \{ Y_t\}_t $ is a set of $S \times L$ matrices representing the bipartite relation data at each observation $t$. The set $\{X_t \}_t$ is a set of $S \times L \times p$ arrays describing the influence of the coefficient vector $beta$. Finally, each matrix $E_t$ is assumed to consist of iid standard normal random variables. The matrices $A$ and $B$ are square matrices respesenting the influence networks among $S$ senders and $L$ receivers, respectively.

This function estimates the BLIN model using maximum likelihood (and related) methods. The "full" model places no restrications on the influence networks $A$ and $B$, and estimates these matrices (along with $\beta$) by block coordinate descent. In addition, if calcses==TRUE, the standard errors for each coefficient will be estimated. Note that the standard error procedure may require large amounts of memory to build the BLIN design matrix; a warning is produced if the estimated size of the desgn is greater than 0.5GB.

The "reduced rank" BLIN model assumes that the matrix $A$ has decomposition $A = UV^T$, where each of $U$ and $V$ is an $S \times \code{rankA}$ matrix, and the matrix $B$ has decomposition $B = WZ^T$, where each of $W$ and $Z$ is an $L \times \code{rankB}$ matrix. This model is also estimated using block coordinate descent.

Finally, the "sparse" BLIN model assumes that $A$ and $B$ matrices have many entries that are small or zero. The cv.glmnet(.) function from the glmnet package is used to estimate the entries in $A$, $B$, and $beta$. The object resuling from cv.glmnet(.) is returned in this case.

Notice that the diagonals of $A$ and $B$ are not identifiable. However, the sum of each diagonal entry in $A$ and $B$, i.e. $a_{ii} + b_{jj}$, is identifiable. Thus, the diagonal sums are broken out as separate estimates under the name diagAB.

If calcses = TRUE and type = full, then standard errors will be returned. These standard errors are based on the assumption that each $E_t$ consists of iid standard normal random variables. In this case, the full design matrix is built, which we call $W$ here. Then, the variance-covariance matrix of the estimated coefficients is formed by $\hat{\tau}^2 (W^T W)^{-1}$, where $\hat{\tau}^2$ is the usual unbiased estimator of the error variance.

Examples

Run this code

# NOT RUN {
S <- 5
L <- 4
tmax <- 10
data <- generate_blin(S,L,tmax, lag=2, sparse=.8, seed=1)

fit <- blin_mle(data$Y, data$X, lag=2, calcses=TRUE)
summary(fit)

# }