One and two-stage least squares fits for spatial Durbin models

one-two-ls-durbin-fits.md

# pak::pak("gpiras/sphet")

library(fastRG)
#> Loading required package: Matrix
library(sphet)

# simulate from y = alpha + trt * gamma + G trt * delta + lambda G y

set.seed(26)

n <- 500
k <- 5

B <- matrix(0.1, nrow = k, ncol = k)
diag(B) <- 0.5

model <- sbm(
  n = n,
  B = B,
  pi = rep(1 / k, k),
  allow_self_loops = FALSE,
  poisson_edges = FALSE,
  expected_degree = n^0.6
)

A <- sample_sparse(model)
G <- Matrix::rowScale(A, 1 / rowSums(A))

trt <- rbinom(n, size = 1, prob = 0.5)
Gtrt <- drop(G %*% trt)
varepsilon <- rnorm(n)

alpha <- 3
gamma <- 5
lambda <- 0.4
delta <- 7

undiffused <- alpha + trt * gamma + Gtrt * delta + varepsilon

I <- Matrix::Diagonal(n, 1)
y <- qr.solve(I - lambda * G, undiffused)
Gy <- drop(G %*% y)

node_data <- data.frame(
  y = y,
  lambda = Gy,
  trt = trt,
  lag_trt = Gtrt
)

# fit with OLS. consistent and ASN under degree restrictions on A, but
# inefficiency

ols_fit <- lm(y ~ trt + lambda + lag_trt, data = node_data)
summary(ols_fit)
#> 
#> Call:
#> lm(formula = y ~ trt + lambda + lag_trt, data = node_data)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -2.70292 -0.66418 -0.03877  0.68882  3.04257 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -0.58494    3.59725  -0.163  0.87089    
#> trt          4.99302    0.09532  52.381  < 2e-16 ***
#> lambda       0.70037    0.29874   2.344  0.01945 *  
#> lag_trt      5.05672    1.94054   2.606  0.00944 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.013 on 496 degrees of freedom
#> Multiple R-squared:  0.866,  Adjusted R-squared:  0.8652 
#> F-statistic:  1068 on 3 and 496 DF,  p-value: < 2.2e-16

# fit with two-stage least squares. consistent and ASN without requiring
# minimum degree on A to be omega(sqrt(n)), and also efficiency in the class
# of linear estimators

tsls_fit <- spreg(
  y ~ trt,
  data = node_data,
  listw = G,
  model = "lag",
  Durbin = TRUE
)

summary(tsls_fit)
#> ======================================================
#> ======================================================
#>                Spatial Durbin Model
#> ======================================================
#> ======================================================
#> 
#> Call:
#> spreg(formula = y ~ trt, data = node_data, listw = G, model = "lag", 
#>     Durbin = TRUE)
#> 
#> Residuals:
#>     Min.  1st Qu.   Median  3rd Qu.     Max. 
#> -2.67964 -0.66634 -0.04018  0.68730  3.03068 
#> 
#> Coefficients:
#>             Estimate Std. Error t-value Pr(>|t|)    
#> (Intercept) 0.225817   4.692484  0.0481   0.9616    
#> trt         4.999734   0.098538 50.7392   <2e-16 ***
#> lag_trt     5.451183   2.432059  2.2414   0.0250 *  
#> lambda      0.632528   0.390924  1.6180   0.1057    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2024-02-28 with reprex v2.0.2}