使用spgwr包LMZ.F3GWR.test()计算p值时遇奇异错误求助

你遇到的错误本质是局部子样本的加权设计矩阵出现严重共线性，导致矩阵无法求逆。虽然全局数据无共线性，但自适应带宽下的局部样本可能出现变量线性相关、样本量过小等问题，以下是针对性解决方法：

1. 对齐核函数或调整带宽参数 #

当前模型拟合用的是gwr.Gauss核，但带宽选择用的是gwr.bisquare，核函数不一致可能导致局部权重分布差异，引发奇异矩阵。另外，自适应带宽可能在部分区域选得过小，导致局部样本不足。

# 方案1：更换核函数为带宽选择时的bisquare
gwr.model = gwr(alias(modelo)$Model,
                data = df_mun_spatial,
                adapt=bw,
                gweight = gwr.bisquare,  # 与gwr.sel保持一致
                hatmatrix=TRUE,
                se.fit=TRUE)

# 方案2：手动扩大带宽（比如乘以1.1），增加局部样本量
gwr.model = gwr(alias(modelo)$Model,
                data = df_mun_spatial,
                adapt=bw*1.1,
                gweight = gwr.Gauss,
                hatmatrix=TRUE,
                se.fit=TRUE)

2. 重新选择更大比例的自适应带宽 #

自适应带宽的adapt参数是按比例选取邻居，比如bw=0.3代表选30%的样本作为局部邻居。尝试提高这个比例，重新计算带宽后拟合模型：

# 指定起始比例为0.4，引导选择更大的自适应带宽
bw_new <- gwr.sel(alias(modelo)$Model,
                  data=df_mun_spatial, adapt = TRUE, gweight = gwr.bisquare,
                  start=0.4)

gwr.model_new = gwr(alias(modelo)$Model,
                    data = df_mun_spatial,
                    adapt=bw_new,
                    gweight = gwr.Gauss,
                    hatmatrix=TRUE,
                    se.fit=TRUE)

3. 给局部矩阵添加岭正则化 #

直接修改LMZ检验的内部逻辑，给奇异的局部设计矩阵添加小的岭惩罚项，避免无法求逆。以下是自定义的带正则化的LMZ检验函数：

LMZ.F3GWR.test_reg <- function(model, lambda=1e-8) {
  y <- model$y
  X <- model$X
  n <- length(y)
  p <- ncol(X)
  w <- model$SDF$gw.weights
  trH <- sum(model$hatmatrix)
  res <- y - model$fitted.values
  sigma2 <- sum(res^2)/(n - trH)
  
  t_vals <- matrix(NA, n, p)
  for (i in 1:n) {
    wj <- w[[i]]
    Xj <- X[wj>0,,drop=FALSE]
    Wj <- diag(wj[wj>0])
    XtWX <- t(Xj) %*% Wj %*% Xj
    # 添加岭惩罚项，避免矩阵奇异
    XtWX_reg <- XtWX + lambda * diag(p)
    
    tryCatch({
      invXtWX <- solve(XtWX_reg)
      beta_j <- invXtWX %*% t(Xj) %*% Wj %*% yj
      var_beta_j <- sigma2 * invXtWX %*% t(Xj) %*% Wj %*% Xj %*% invXtWX
      se_beta_j <- sqrt(diag(var_beta_j))
      t_vals[i,] <- beta_j / se_beta_j
    }, error=function(e) {
      warning(paste("局部模型", i, "添加lambda=", lambda, "后仍奇异"))
      t_vals[i,] <- NA
    })
  }
  
  # 计算双侧t检验的p值
  p_vals <- 2 * pt(-abs(t_vals), df = n - trH)
  colnames(p_vals) <- colnames(X)
  rownames(p_vals) <- rownames(model$SDF)
  return(list(t_values=t_vals, p_values=p_vals))
}

# 使用自定义函数计算p值
lmz_result <- LMZ.F3GWR.test_reg(gwr.model)
# 查看结果
head(lmz_result$p_values)

4. 检查变量的局部恒定值 #

全局标准化后无共线性，但某些局部区域可能存在变量取值完全恒定的情况，导致局部设计矩阵列线性相关。用以下代码检查：

library(spdep)
# 根据自适应带宽生成邻居列表
nb <- knn2nb(knearneigh(coordinates(df_mun_spatial), k=round(bw*nrow(df_mun_spatial))))
listw <- nb2listw(nb)

# 计算每个变量在局部邻居内的变异系数（标准差/均值）
vars <- colnames(gwr.model$X)
local_cv <- lapply(vars, function(var) {
  sapply(1:nrow(df_mun_spatial), function(i) {
    # 获取当前点及邻居的变量值
    neighbor_ids <- c(i, nb[[i]])
    vals <- df_mun_spatial[[var]][neighbor_ids]
    # 避免除以0，添加极小值
    sd(vals)/(mean(vals) + 1e-10)
  })
})
names(local_cv) <- vars

# 查看每个变量局部变异系数接近0的数量（即局部恒定的次数）
lapply(local_cv, function(x) sum(x < 1e-10))

如果某个变量存在大量局部恒定的情况，建议移除该变量，或对其进行分组合并（若为分类变量）。

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