R语言空间分析：二阶邻居缺失致lagsarlm报错求助

问题 #

尝试以新泽西州各市镇的二阶邻居特征为工具变量，采用空间Durbin模型分析时，因部分市镇无二阶邻居出现报错：

Error in eigen(similar.listw_Matrix(get("listw", envir = env)), only.values = TRUE) :
infinite or missing values in 'x'

操作背景：

• 通过poly2nb()获取邻居关系，岛屿型市镇用k=2的k近邻法补充
• 已生成二阶邻居的.nb列表，但存在如市镇26与86这类互为唯一邻居、无二阶邻居的情况
• 已设置nb2listw(zero.policy=TRUE)和lagsarlm(zero.policy=TRUE)，但仍触发报错

解决方案 #

1. 修正二阶邻居生成逻辑 #

原代码存在变量名不一致、空邻居处理错误的问题，重新编写二阶邻居生成代码：

n2 = addnbs(geo_2019)
instru.nei = lapply(n2, function(x) c()) # 初始化空的二阶邻居列表

# 遍历每个市镇生成二阶邻居
for(i in 1:length(n2)){  
  k = 1
  # 遍历当前市镇的所有一阶邻居
  for(j in n2[[i]]){
    # 遍历一阶邻居的邻居，筛选符合条件的二阶邻居
    for(xx in n2[[j]]) {
      if (!(xx %in% n2[[i]]) && (xx != i) && (!(xx %in% instru.nei[[i]]))){
        instru.nei[[i]][k] = xx
        k = k + 1
      }
    }  
  }
  # 清理空元素
  instru.nei[[i]] <- instru.nei[[i]][!is.na(instru.nei[[i]])]
}

2. 删除错误的0赋值代码 #

不要给无二阶邻居的市镇赋值0，空邻居列表在zero.policy=TRUE参数下会被正确识别处理，删除以下代码块：

index_zero <- c() 
for (i in 1:324) {
if(length(instru_nofn[[i]])==0)
  index_zero <- append(index_zero, i)
}

for (i in index_zero) {
  instru_nofn[[i]] = as.integer(0)
}

3. 正确生成权重矩阵并运行模型 #

确保变量名一致，指定权重矩阵的标准化方式为行标准化（style="W"），避免数值异常：

# 生成空间权重矩阵，允许空邻居
queen.non <- nb2listw(instru.nei, zero.policy=TRUE, style = "W")

# 运行空间Durbin模型
sar.durbin = lagsarlm(Certified ~ Median_Age + Percent_Children + Percent_Bachelor + Pop_new + 
                         Percent_sixty + White + Black + Asian + Hispanic + Percent_Unemployed + 
                         Mean_trave_time + Percent_Manufacture + Median_income + Percent_poverty, 
                       data = geo_2019, listw = queen.non, type="mixed", zero.policy=TRUE)

summary(sar.durbin)

4. 可选：强制邻居关系对称 #

二阶邻居关系可能存在不对称性，可通过sym.nb()强制对称，进一步避免权重矩阵异常：

# 强制邻居关系对称
instru.nei_sym <- sym.nb(instru.nei)
queen.non_sym <- nb2listw(instru.nei_sym, zero.policy=TRUE, style = "W")

# 使用对称后的权重矩阵运行模型
sar.durbin_sym = lagsarlm(Certified ~ Median_Age + Percent_Children + Percent_Bachelor + Pop_new + 
                           Percent_sixty + White + Black + Asian + Hispanic + Percent_Unemployed + 
                           Mean_trave_time + Percent_Manufacture + Median_income + Percent_poverty, 
                         data = geo_2019, listw = queen.non_sym, type="mixed", zero.policy=TRUE)

summary(sar.durbin_sym)

内容的提问来源于stack exchange，提问作者Charlie Babe