Network autocorrelation model

network_autocorrelation.R

# simulate and estimate a network autocorrelation model
set.seed(45)
N <- 200
A <- matrix(rbinom(N*N, 1, 0.2), N)
diag(A) <- 0
A2 <- matrix(rbinom(N*N, 1, 0.2), N)
diag(A2) <- 0
An <- A / rowSums(A) # row-normalized ????

# params
rho <- 0.2
phi <- 0
beta <- c(0.1, -.5)
sigma <- 0.2

# draw data
nu <- rnorm(N, sd = sigma)
X <- matrix(rnorm(N*2), N)
e <- solve(diag(N) - phi*A, nu)
y <- solve(diag(N) - rho*A, X%*%beta + e) # see https://search.r-project.org/CRAN/refmans/sna/html/lnam.html

x <- A%*%y
solve(crossprod(x), crossprod(x, y)) # rho is estimated correctly!

#' Network linear model
#' 
#' Model structure:
#' y = X beta + rho W_ar y + e
#' e = phi W_ma e + nu
#' 
#' @param y the outcome variable
#' @param X the linear model predictors
#' @param W_ar AR-type weights matrix (spatial autoregressive model, network autoregression)
#' @param W_ma MA-type weights matrix (spatial error model)
#' @param maxit maximum number of iterations
#' @param verbose whether to print the parameter values in the intermediate iterations
#' 
#' @return list of parameter estimates
#' 
#' @export
networklm <- function(y, X, W_ar, W_ma, maxit = 100L, verbose = FALSE) {
  beta_hat <- if (missing(X)) NA else numeric(ncol(X))
  rho_hat  <- if (missing(W_ar)) NA else 0
  phi_hat  <- if (missing(W_ma)) NA else 0
  reg_part <- ar_part <- ma_part <- rep(0, length(y))
  for (i in 1:maxit) {
    
    # compute regression
    if (!missing(X)) {
      beta_hat <- compute_beta(y, X, c(ar_part + ma_part))
      reg_part <- X %*% beta_hat
      if (verbose) cat("iter", i, "beta_hat:", beta_hat, "\n")
    }
    
    # compute ar parameter
    if (!missing(W_ar)) {
      rho_hat  <- compute_rho(y, A_ar, c(reg_part + ma_part))
      ar_part  <- (rho_hat * W_ar) %*% y
      if (verbose) cat("iter", i, "rho_hat:", rho_hat, "\n")
    }
    
    # compute ma parameter
    if (!missing(W_ma)) {
      # NB: some problem here because now we assume res_cur is 
      # observed, which it is not. We may need to add noise? 
      # look at what the Bayesian posterior is?
      # compute partial likelihood and optim it?
      res_cur  <- y - reg_part - ar_part
      phi_hat  <- compute_phi(res_cur, W_ma)
      ma_part  <- (phi_hat * W_ma) %*% res_cur
      if (verbose) cat("iter", i, "phi_hat:", phi_hat, "\n")
    }
  }
  sigma_hat <-  sd(y - reg_part - ar_part - ma_part)
  return(list(beta = beta_hat, rho = rho_hat, phi = phi_hat, sigma = sigma_hat))
}

compute_beta <- function(y, X, offset = 0)    lm.fit(X, y, offset)$coefficients
compute_rho  <- function(y, W_ar, offset = 0) lm.fit(W_ar%*%y, y, offset)$coefficients
compute_phi  <- function(e, W_ma)             lm.fit(W_ma%*%e, e)$coefficients
  
coefs <- networklm(y, X, A, verbose = TRUE)
res   <- summary(sna::lnam(y, X, A, A))


coefs
coef(res)

Maybe it's not in the valid parameter space!