R语言aov模型输出缺失Pr(>F)值，如何获取该值？

代码与当前输出 #

原始代码 #

data = data.frame(
  route = rep(c(1, 2, 3, 4), each = 5),
  day = rep(1:5, times = 4),
  time = c(22,26,25,25,31,25,27,28,26,29,26,29,33,30,33,26,28,27,30,30)
)

model = aov(time ~ route + day + Error(route * day), data = data)

summary(model)输出 #

> summary(model)
Error: route
      Df Sum Sq Mean Sq
route  1  27.04   27.04

Error: day
    Df Sum Sq Mean Sq
day  1  60.03   60.03

Error: route:day
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  1  0.845   0.845               

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 16  65.29   4.081

原因分析 #

当前输出缺失Pr(>F)值，是因为你设定的Error(route * day)误差结构将数据变异拆分为多个独立层，但route和day的主效应没有对应的误差项作为F检验的分母，R无法自动计算p值。同时从输出的自由度可以看出，模型对变异的拆分不符合实验设计（比如route有4个水平，理论自由度应为3，但输出中是1），说明误差结构设定有误。

解决方法 #

方法1：根据实验设计调整模型结构 #

如果你的数据是两因素无重复交叉设计（每个route×day组合仅1个观测值），无需使用Error项，直接拟合标准两因素方差分析模型：

model <- aov(time ~ route + day, data = data)
summary(model)

运行后会自动用残差作为误差项，输出完整的F value和Pr(>F)值。

方法2：手动计算p值 #

如果坚持使用当前的误差结构，可以提取各层的均方手动计算：

1. 计算route主效应的F值与p值：

   f_route <- 27.04 / 0.845
   p_route <- pf(f_route, df1 = 1, df2 = 1, lower.tail = FALSE)
   cat("route的F值：", round(f_route, 2), "，p值：", round(p_route, 4), "\n")
   

2. 计算day主效应的F值与p值：

   f_day <- 60.03 / 0.845
   p_day <- pf(f_day, df1 = 1, df2 = 1, lower.tail = FALSE)
   cat("day的F值：", round(f_day, 2), "，p值：", round(p_day, 4), "\n")
   

注意：这种情况下自由度极低（df2=1），检验效能很差，结果参考价值有限。

方法3：修正重复测量设计的误差结构 #

如果你的实验是重复测量设计（比如route是被试，day是重复测量的时间点），应调整误差项为Error(route)，让R用被试内变异计算day效应的p值：

model <- aov(time ~ day + Error(route), data = data)
summary(model)

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