教育支出对犯罪率影响的多水平混合模型构建疑问咨询

我正在开展一项为期15年的研究，分析各学区生均教育支出对其服务城镇犯罪率的影响。我的因变量是1,676,191条城镇犯罪数据观测值，数据结构存在几个层级关系：

• 城镇与学区存在交叉结构：一个城镇对应多个学区，意味着单个城镇的生均支出有多个取值
• 学区与县重叠，城镇嵌套于县，层级结构可表示为(1|PLACE_ID)+(1|COUNTY_ID)或(1|PLACE_ID/COUNTY_ID)

我本身熟悉混合效应模型，也参考过相关资料，但还是对自己构建的交叉效应模型合法性存疑，同时有数据处理上的困惑。

模型代码与变量说明 #

模型代码 #

glmer.total <- glmer(CRIME_TOTAL ~ cent.log.pop + cent.log.pop.dens + year + cent.log.unemployment_rate + cent.schooldist.prop5.17.pov + cent.log.per.cap + diff.dem + cent.log.enforcement +cent.EXP_STUDENT + (year|PLACE_ID/COUNTY_ID) + (year|full_district_id) + (1|STATE), family = "poisson", control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), REML = FALSE, total.years, na.action = "na.omit")

变量说明 #

所有变量已做相应预处理（中心化/对数转换）：

• cent.log.pop：中心化后的城镇人口对数
• cent.log.pop.dens：中心化后的城镇人口密度对数
• year：年份
• cent.log.unemployment_rate：中心化后的县失业率对数
• cent.schooldist.prop5.17.pov：中心化后的学区5-17岁儿童贫困率
• cent.log.per.cap：中心化后的县人均收入对数
• diff.dem：县总统选举民主党得票差异
• cent.log.enforcement：中心化后的城镇警力对数
• cent.EXP_STUDENT：中心化后的学区生均支出（单位：/1000）
• 分组变量：PLACE_ID（城镇）、COUNTY_ID（县）、full_district_id（学区）、STATE（州）

核心疑问 #

1. 是否需要按年份对城镇的生均支出取均值？
2. 上述模型代码是否合法？

报错信息 #

运行代码后出现两个错误：

extra argument(s) ‘REML’ disregarded
Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GHrule(0L), compDev = compDev, : (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

数据片段（dput格式） #

structure(list(STATE = c("alabama", "alabama", ...), ...), row.names = c(NA, -68L), class = c("tbl_df", "tbl", "data.frame"))

解答 #

1. 关于生均支出是否需要按年份取城镇均值 #

这个完全取决于你的研究假设：

• 如果你的核心关注是学区层面的生均支出对城镇犯罪率的影响，且同一个城镇内不同学区的支出差异是你想捕捉的关键变异，那么不需要取均值——保留原始的学区-城镇-年份层面的支出数据是合理的，这能反映城镇内部不同学区资源投入的差异对整体犯罪率的潜在影响。
• 但如果你的研究假设是城镇整体的教育投入水平对犯罪率的影响，那按年份对城镇内各学区的生均支出取均值是更合适的，此时你是在聚合学区层面的变量到城镇层面，对应“城镇平均教育投入”的构念。

2. 模型代码的合法性与报错解决 #

代码合法性问题

你的模型结构本身有尝试的空间，但有几个关键问题需要修正：

• (year|PLACE_ID/COUNTY_ID) 这个结构表示的是县嵌套于城镇的随机斜率（随year变化），但根据你描述的层级“城镇嵌套于县”，这个结构是反的！正确的嵌套随机斜率应该是(year|COUNTY_ID/PLACE_ID)，也就是城镇嵌套在县中，每个县内的城镇有各自随年份变化的斜率。
• 你同时加入了(year|full_district_id)（学区层面随年份变化的随机斜率）和城镇-县层级的随机斜率，这确实是在捕捉交叉效应，但模型复杂度极高，很容易导致拟合失败。

报错解决

• 第一个错误：extra argument(s) ‘REML’ disregarded
  这是因为泊松分布的广义线性混合模型（GLMM）不支持REML参数——只有线性混合模型（LMM）才支持REML拟合，GLMM只能用ML（最大似然）。直接去掉REML = FALSE这个参数即可，它在GLMM中是多余的。
• 第二个错误：(maxstephalfit) PIRLS step-halvings failed to reduce deviance
  这个错误通常是模型拟合困难导致的，建议从以下几个方向排查：
  1. 简化模型复杂度：先从只包含随机截距的模型开始拟合，比如(1|STATE) + (1|COUNTY_ID/PLACE_ID) + (1|full_district_id)，确认模型能收敛后，再逐步加入随year变化的随机斜率。
  2. 检查变量尺度：虽然你已经中心化了变量，但再确认所有连续预测变量的尺度是否在同一数量级（比如生均支出除以1000是合理的，可检查其他对数转换变量的尺度）。
  3. 更换优化器：你当前用了nloptwrap，可以尝试bobyqa或Nelder_Mead，修改为control = glmerControl(optimizer = "bobyqa")测试拟合情况。
  4. 排查数据分离：泊松模型容易出现数据分离问题（比如某些分组下犯罪率为0的观测过多），可以检查分组变量的观测数量，比如是否存在某个full_district_id或PLACE_ID的观测数过少。

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