具有嵌套效应且系数为NA的模型

编程语言 2026-07-10

我正在运行一个前后对照影响（before-after-control-impact，BACI）模型，包含1个对照点和1个受影响点，Before期有若干年，After期有2年。我想用Year（作为因子）在Before/After效应内嵌来指定。起初我使用了生存模型（因为我的数据因未检测值而被删失）。该模型运行良好，除了最后一个采样年份（2024）的系数始终是NA。我以为是生存模型的问题，直到我用相同结构运行了一个简单的 lm()，也出现了同样的问题。也就是说，那个最后的系数一定是不可估计的，原因可能是什么？可能是参数过多（overparameterization），结合Ben Bolker的评论这里。原因是什么，我应该怎么做才能得到所有年份的系数？

With 2年 Before与 2年 After的示例：

library(dplyr)
library(survival)
library(ggplot2)

md <- structure(list(Year = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L), levels = c("2021", "2022", "2023", "2024"), class = "factor"), 
   CI = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L), levels = c("Control", "Impact"), class = "factor"), 
    value = c(0.00022, 0.00021, 5e-04, 5e-04, 2e-04, 2e-04, 0.00016, 
    0.00016, 0.00021, 2e-04, 2e-04, 0.00019, 2e-04, 0.00023, 
    0.00024, 0.00018, 0.00022, 0.00023, 5e-04, 5e-04, 5e-04, 
    2e-04, 0.00025, 2e-04, 0.00017, 2e-04, 2e-04, 2e-04, 2e-04, 
    0.00022, 0.00019, 2e-04, 0.00015, 2e-04, 2e-04, 2e-04, 0.00017, 
    0.00016, 0.00021, 2e-04, 0.00023, 2e-04, 5e-04, 5e-04, 0.00023, 
    0.00017, 0.00016, 2e-04, 2e-04, 0.00018, 0.00016, 0.00016, 
    0.00018, 0.00021, 2e-04, 2e-04, 2e-04, 2e-04, 0.00017, 0.00021, 
    2e-04, 5e-04, 2e-04, 2e-04, 0.00023, 2e-04, 2e-04, 2e-04, 
    2e-04, 0.00016, 0.00019, 0.00016, 0.00018, 0.00018, 0.00018, 
    2e-04, 5e-04, 2e-04, 2e-04, 0.00019, 2e-04, 2e-04, 2e-04, 
    0.00017, 0.00016, 2e-04, 0.00016, 2e-04, 2e-04, 2e-04, 2e-04, 
    2e-04, 2e-04, 2e-04, 2e-04, 0.00015, 0.00018, 0.00016, 0.00017, 
    2e-04, 2e-04, 2e-04, 2e-04, 2e-04, 2e-04), censor = c(FALSE, 
    FALSE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, 
    TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, 
    TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, 
    TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, 
    FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, 
    FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, 
    FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, 
    FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, 
    FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, 
    TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, 
    TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, 
    TRUE, TRUE, TRUE, TRUE, TRUE)), row.names = c(NA, -105L), class =  
    "data.frame") %>%
mutate(BA = ifelse(Year %in% c(2020:2022), "Before", "After"))

ggplot(md, aes(x = Year, y = value, fill = CI)) +
    geom_boxplot() + facet_wrap(~CI) + theme_bw()
table(md$Year, md$BA, md$CI)

## survival model
srv1 <- survreg(Surv(time = value, event = (1 - censor), type = "left") ~ (BA/Year) + CI + (BA/Year):CI, data = md, dist = "lognormal")
summary(srv1) 
# shows NA for 2024 in both before and after category (whereas I'd have expected a value under the After category)

# same issue happening with lm()
mod <- lm(value ~ (BA/Year) + CI + (BA/Year):CI, data = md)
summary(mod)

解决方案

第一个模型

我们可以尝试如下代码：

library(emmeans)
fm <- lm(value ~ CI * Year, md)
em <- emmeans(fm, ~ Year | CI)
contrast(em, list(BA = c(0.5, 0.5, -0.5, -0.5)))

这里 fm 是线性模型对象，从中我们得到 em 的emmeans对象。这些emmeans包含预测值（以及其他信息）。例如，第一个如下。它与 em 对象第一行的emmean列中的值相同。

predict(fm, list(CI = "Control", Year = "2021"))

共有8 个这样的预测，因为有2 个CI水平和4 个Year水平，因此共有2 × 4 = 8种组合。

上面 constrast 语句输出的第一行给出的结果与以下相同：

library(multcomp)
fmControl <- lm(value ~ Year + 0, md, subset = CI == "Control")
glht(fmControl, "0.5*Year2022 - 0.5*Year2023 - 0.5*Year2024 = 0")

第二个模型

我们可以把 joint_tests 解释为对一组对比的 Wald检验。在一个满因子设计2^k的情况下，我们可以使用通过Sylvester/Walsh构造生成的Hadamard矩阵，并由 pracma::hadamard 适当缩放以实现正交。该构造总是产生一个对称矩阵，因此尽管通常对比存放在行中，我们也可以在列上使用 apply，这有消除转置的优点。

我们通过一个经过标准化以实现正交的 Hadamard矩阵来实现，方法是除以2。对每一行执行所指示的Wald检验。

library(emmeans)
library(aod) # wald.test
library(pracma) # hadamard

fm2 <- lm(value ~ CI * BA, md)
rg  <- ref_grid(fm2)
b   <- predict(rg)
V   <- vcov(rg)
df  <- df.residual(fm2)

L <- hadamard(4) / sqrt(4)

results <- apply(L, 2, function(L) {
  wald.test(Sigma = V, b = b, L = matrix(L, 1), df = df)$result$Ftest
})

rownames(results) <- c("Int", "BA", "CI", "BA:CI")
print(results)
##             Fstat df1 df2         P
## Int   636.0566952   1 101 0.0000000
## BA      0.1288346   1 101 0.7203926
## CI      2.5700865   1 101 0.1120246
## BA:CI   0.8741506   1 101 0.3520388

joint_tests(fm2)
##  model term df1 df2 F.ratio p.value
##  CI           1 101   0.129  0.7204
##  BA           1 101   2.570  0.1120
##  CI:BA        1 101   0.874  0.3520

注意，我们对其进行了适当简化，但实际上我们展示的Hadamard矩阵构造仅适用于2^k因子设计。这里是2x2的情况，因此可行且简单，这也是我们使用它的原因；然而，实际上这并非 joint_tests 实际使用的对比矩阵。请认识到，如果 A 是任意满秩方阵，且下面的矩阵运算在该矩阵下成立，那么 wald.test 不依赖于我们使用的 A。
如果 L 是我们从上述 pracma::hadamard 推导出的对比矩阵，那么联合检验使用的 A <- attr(rg, "linefct")，事实证明在数值上更稳定。

# A cancels making this invariant to the value of square full rank matrix A
wald.test(A %*% V %*% t(A), A %*% b, L = L %*% solve(A), df = res_df)

旧方法

我们也可以尝试：

library(emmeans)
library(ggplot2)
fm2 <- lm(value ~ CI * BA, md)
em2 <- emmeans(fm2, ~ BA | CI)
cc <- contrast(em2, list(BA = c(1, -1))); cc
plot(cc)
pairs(emmeans(fm2, ~ BA | CI))
joint_tests(fm2)

在这种情况下，BA有两个水平，CI也有两个水平，因此在类似的方式中，在 em2 中将有4 个预测，同样我们也可以对每个CI水平下的指示对比进行检验（并绘制它们）。在这种情况下，对比是成对比较，emmeans提供了一个更简单的函数 pairs，它为你提供对比。

从概念上讲，joint_tests 的工作原理是通过 ref_grid 来平衡数据，然后对其应用方差分析（anova）函数。这实际上并不奏效，因为anova无法处理零残差，但希望这能给出思路。

rg2 <- ref_grid(fm2)
md_bal <- as.data.frame(rg2)
fm3 <- lm(prediction ~ CI * BA, md_bal, weights = 1 / diag(vcov(rg2)),
  contrasts = list(CI = "contr.sum", BA = "contr.sum"))
anova(fm3, df = df.residual(fm2)) # does not work but this is the idea
joint_tests(fm2) # same

不幸的是我找不到一个包来完成最后这一块，但从原理出发：

library(emmeans)
fm2 <- lm(value ~ CI * BA, md)
rg <- ref_grid(fm2)
means <- predict(rg)      # The 4 means: [μ11, μ12, μ21, μ22]
V_diag <- diag(vcov(rg))  # The 4 variances

L_A  <- c(1, 1, -1, -1) / 2
L_B  <- c(1, -1, 1, -1) / 2
L_AB <- c(1, -1, -1, 1) / 2

# F-statistics: (Estimate^2) / Variance_of_Estimate
calc_f <- function(L, means, V_diag, res_df) {
  estimate <- sum(L * means)
  # Variance of a weighted sum (since V is diagonal) is sum(L^2 * V)
  var_est  <- sum(L^2 * V_diag)
  f_val    <- (estimate^2) / var_est
  p_val    <- pf(f_val, df1 = 1, df2 = res_df, lower.tail = FALSE)
  return(c(F = f_val, p = p_val))
}

res_df <- df.residual(fm2)
results <- rbind(
  CI    = calc_f(L_B,  means, V_diag, res_df),
  BA    = calc_f(L_A,  means, V_diag, res_df),
  "CI:BA"   = calc_f(L_AB, means, V_diag, res_df)
)

print(results)

joint_tests(fm2) # same

请注意，这给出了等价的t 检验（p值相同，F = t^2）。也就是说，从定义 calc_f 开始的所有内容都可以用这个来替代，除了它展示的是等价的t 检验而不是F 检验。

library(glht)
calc_t <- function(L, rg) summary(as.glht(contrast(rg, list(L))))
results2 <- rbind(
  CI    = calc_t(L_B,  rg),
  BA    = calc_t(L_A,  rg),
  "CI:BA" = calc_t(L_AB, rg)
)
print(results2)

或使用car包：

library(car)
calc_f2 <- function(L, fm) {
  linearHypothesis(fm, attr(contrast(rg, list(custom = L)), "linfct"))[2, ]
}
results3 <- rbind(
  CI    = calc_f2(L_B,  fm2),
  BA    = calc_f2(L_A,  fm2),
  "CI:BA" = calc_f2(L_AB, fm2)
)
print(results3)

站内所有文章版权归属LeftHeroAI导航站，无授权禁止任何主体转载、抄袭、复制内容，亦不得私自架设镜像站点。一经侵权，本站将通过法律途径追责。