glmer二项式模型：mL与µl计数的过度离散估计差异问题

我们在微生物实验中，针对暴露于不同药物（数据框plate字段）的耐药/易感菌落计数比例数据拟合二项式模型。实验采用20µl培养液平板计数，通常换算为每mL细胞数，导致总计数极高（可达数十亿）。

用glmer()拟合广义线性混合效应模型时，出现严重过度离散；改用每µl细胞数（数值缩小1000倍）后，模型预测结果完全一致，但过度离散值显著降低。我们需要明确：

1. 该现象的核心原因
2. 后续分析的计数单位选择建议

1. 二项式模型的本质与数值缩放的无关性 #

二项式模型的核心是对比例进行建模，对数优势比的参数估计与计数单位完全无关——不管总计数n是放大还是缩小，耐药比例p始终不变，因此模型的预测结果不会改变。但总计数的量级会直接影响模型拟合的数值稳定性和过度离散的计算。

2. 过度离散计算的数值精度问题 #

blmeco::dispersion_glmer()通过皮尔逊卡方/自由度计算过度离散值：

• 当n极大时，计算皮尔逊残差会涉及超大数值，容易触发浮点运算精度误差，导致残差被错误放大，进而让过度离散值异常偏高；
• 二项式模型默认的方差结构n*p*(1-p)在n量级过大时，对真实数据方差的拟合会被数值误差干扰，进一步加剧过度离散的计算偏差。

缩小计数单位后，n的量级回到合理范围，浮点精度误差减少，计算出的过度离散值更接近真实的离散程度（仍存在过度离散，说明二项式模型本身不适合这类数据）。

3. 模型拟合的数值奇异性 #

当n极大时，glmer()优化过程中处理的数值矩阵会出现奇异性（如示例中的“Model is nearly unidentifiable: very large eigenvalue”警告），这会干扰模型收敛稳定性，同时也会影响过度离散的计算准确性。

• 优先匹配原始实验操作的单位：实验用20µl培养液计数，直接使用每µl的计数（或平板原始菌落数换算为每µl），避免不必要的数值放大；
• 选择量级适中的单位：尽量让总计数n处于10²~10⁴的区间，减少浮点精度误差，提升模型拟合的稳定性；
• 以比例为核心，忽略单位的“常规性”：只要耐药比例不变，计数单位不影响预测结果，无需拘泥于“每mL”这类常规单位；
• 结合模型类型调整：如果使用brms拟合β-二项式模型（本身适配过度离散数据），计数单位的影响会更小，但仍建议选择量级适中的单位，提升拟合效率。

# does it matter what plate counts are corrected to for a glmer?

# Load packages ####
library(lme4)
#> Loading required package: Matrix
library(blmeco)
#> Loading required package: MASS
library(tidyverse)

# read in example
d_example <- tibble::tribble(
               ~treatment, ~plate, ~count_cor_ul, ~count_susc_ul, ~total_count_ul,            ~prop_ul, ~count_cor_ml, ~count_susc_ml, ~total_count_ml,            ~prop_ml,
                    "0_1",  "amp",          3600,           8900,           12500,               0.288,       3600000,        8900000,        12500000,               0.288,
                    "0_1",  "azi",           750,          11750,           12500,                0.06,        750000,       11750000,        12500000,                0.06,
                    "0_1",  "cef",             5,          12495,           12500,               4e-04,          5000,       12495000,        12500000,               4e-04,
                    "0_1",  "cip",           150,          12350,           12500,               0.012,        150000,       12350000,        12500000,               0.012,
                    "0_1",  "gen",          5450,           7050,           12500,               0.436,       5450000,        7050000,        12500000,               0.436,
                    "0_1",  "tet",           285,          12215,           12500,              0.0228,        285000,       12215000,        12500000,              0.0228,
                    "0_2",  "amp",           650,          11500,           12150,  0.0534979423868313,        650000,       11500000,        12150000,  0.0534979423868313,
                    "0_2",  "azi",           210,          11940,           12150,  0.0172839506172839,        210000,       11940000,        12150000,  0.0172839506172839,
                    "0_2",  "cef",            40,          12110,           12150, 0.00329218106995885,         40000,       12110000,        12150000, 0.00329218106995885,
                    "0_2",  "cip",           150,          12000,           12150,  0.0123456790123457,        150000,        1.2e+07,        12150000,  0.0123456790123457,
                    "0_2",  "gen",          7900,           4250,           12150,   0.650205761316872,       7900000,        4250000,        12150000,   0.650205761316872,
                    "0_2",  "tet",           110,          12040,           12150, 0.00905349794238683,        110000,       12040000,        12150000, 0.00905349794238683,
                    "0_3",  "amp",          3400,          24100,           27500,   0.123636363636364,       3400000,       24100000,        27500000,   0.123636363636364,
                    "0_3",  "azi",          1750,          25750,           27500,  0.0636363636363636,       1750000,       25750000,        27500000,  0.0636363636363636,
                    "0_3",  "cef",           115,          27385,           27500, 0.00418181818181818,        115000,       27385000,        27500000, 0.00418181818181818,
                    "0_3",  "cip",           455,          27045,           27500,  0.0165454545454545,        455000,       27045000,        27500000,  0.0165454545454545,
                    "0_3",  "gen",         12800,          14700,           27500,   0.465454545454545,      12800000,       14700000,        27500000,   0.465454545454545,
                    "0_3",  "tet",          4000,          23500,           27500,   0.145454545454545,         4e+06,       23500000,        27500000,   0.145454545454545,
                    "0_4",  "amp",          3900,           2700,            6600,   0.590909090909091,       3900000,        2700000,         6600000,   0.590909090909091,
                    "0_4",  "azi",            50,           6550,            6600, 0.00757575757575758,         50000,        6550000,         6600000, 0.00757575757575758,
                    "0_4",  "cef",            10,           6590,            6600, 0.00151515151515152,         10000,        6590000,         6600000, 0.00151515151515152,
                    "0_4",  "cip",            50,           6550,            6600, 0.00757575757575758,         50000,        6550000,         6600000, 0.00757575757575758,
                    "0_4",  "gen",          4500,           2100,            6600,   0.681818181818182,       4500000,        2100000,         6600000,   0.681818181818182,
                    "0_4",  "tet",           300,           6300,            6600,  0.0454545454545455,         3e+05,        6300000,         6600000,  0.0454545454545455,
                    "0_5",  "amp",           850,           7150,            8000,             0.10625,        850000,        7150000,           8e+06,             0.10625,
                    "0_5",  "azi",           350,           7650,            8000,             0.04375,        350000,        7650000,           8e+06,             0.04375,
                    "0_5",  "cef",            80,           7920,            8000,                0.01,         80000,        7920000,           8e+06,                0.01,
                    "0_5",  "cip",           220,           7780,            8000,              0.0275,        220000,        7780000,           8e+06,              0.0275,
                    "0_5",  "gen",          4200,           3800,            8000,               0.525,       4200000,        3800000,           8e+06,               0.525,
                    "0_5",  "tet",           640,           7360,            8000,                0.08,        640000,        7360000,           8e+06,                0.08
               )

# ok want to fit a binomial model to see if resistance changes between drugs and antibiotics

# try and fit model to the ml data
mod1 <- glmer(cbind(count_cor_ml, count_susc_ml) ~ plate + (1|treatment), d_example, family = binomial)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
#>  - Rescale variables?

# check for overdispersion (it should not be over 1.4)
blmeco::dispersion_glmer(mod1)
#> [1] 656.0779

# fit using microlitre data
mod2 <- glmer(cbind(count_cor_ul, count_susc_ul) ~ plate + (1|treatment), d_example, family = binomial)

# check for overdispersion (it should not be over 1.4)
blmeco::dispersion_glmer(mod2)
#> [1] 20.75102

# do the predictions of the model change?

# we can plot the model predictions easy enough
preds1 <- emmeans::emmeans(mod1, ~plate, type = 'response') %>%
  data.frame() %>%
  janitor::clean_names() %>%
  mutate(pred = 'ml')
preds2 <- emmeans::emmeans(mod2, ~plate, type = 'response') %>%
  data.frame() %>%
  janitor::clean_names() %>%
  mutate(pred = 'ul')

preds <- bind_rows(preds1, preds2) %>%
  left_join(., dplyr::select(d_example, plate) %>% distinct())
#> Joining, by = "plate"

# plot data
ggplot(d_example) +
  geom_line(aes(plate, prop_ul, group = treatment), col = 'light grey') +
  geom_point(aes(plate, prob, col = pred), preds, position = position_dodge(0.7), size = 3) +
  geom_linerange(a