如何在4×4网格的GAM拟合图内有序添加edf、调整R方与P值

需求说明 #

• 对16个种群的时序数据分别拟合GAM模型
• 所有种群的原始观测点、拟合趋势线、置信区间按4×4网格布局展示
• 每个子图内分三行依次展示edf、调整R平方、P值三类模型统计量

你的数据集结构 #

> dput(bb)
structure(list(Year = c(1970, 1971, 1972, 1973, 1974, 1975, 1976, 
1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 
1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 
1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 
2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018), F = c(28, 
31, 27, 34, 24, 16, 18, 33, 31, 25, 33, 34, 46, 45, 45, 40, 18, 
18, 25, 14, 21, 16, 21, 20, 21, 20, 18, 24, 26, 17, 24, 16, 18, 
15, 21, 35, 53, 64, 51, 52, 59, 52, 56, 52, 50, 51, 44, 20, 30
), M = c(24, 39, 53, 39, 41, 18, 17, 33, 45, 48, 48, 60, 47, 
34, 36, 35, 28, 27, 24, 32, 29, 30, 35, 35, 38, 40, 37, 34, 43, 
31, 34, 44, 45, 46, 51, 61, 112, 109, 116, 119, 89, 106, 103, 
82, 87, 81, 67, 44, 35), U = c(7, 4, 3, 8, 7, 42, 29, 6, 21, 
14, 18, 7, 13, 13, 12, 40, 12, 51, 24, 0, 0, 0, 14, 79, 60, 64, 
67, 81, 96, 98, 68, 97, 118, 206, 156, 136, 1, 0, 0, 1, 0, 1, 
5, 1, 1, 4, 1, 71, 82), Tot = c(59, 74, 83, 81, 72, 76, 64, 72, 
97, 87, 99, 101, 106, 92, 93, 115, 58, 96, 73, 46, 50, 46, 70, 
134, 119, 124, 122, 139, 165, 146, 126, 157, 181, 267, 228, 232, 
166, 173, 167, 172, 148, 159, 164, 135, 138, 136, 112, 135, 147
), ratio = c(0.857142857142857, 1.25806451612903, 1.96296296296296, 
1.14705882352941, 1.70833333333333, 1.125, 0.944444444444444, 
1, 1.45161290322581, 1.92, 1.45454545454545, 1.76470588235294, 
1.02173913043478, 0.755555555555556, 0.8, 0.875, 1.55555555555556, 
1.5, 0.96, 2.28571428571429, 1.38095238095238, 1.875, 1.66666666666667, 
1.75, 1.80952380952381, 2, 2.05555555555556, 1.41666666666667, 
1.65384615384615, 1.82352941176471, 1.41666666666667, 2.75, 2.5, 
3.06666666666667, 2.42857142857143, 1.74285714285714, 2.11320754716981, 
1.703125, 2.27450980392157, 2.28846153846154, 1.50847457627119, 
2.03846153846154, 1.83928571428571, 1.57692307692308, 1.74, 1.58823529411765, 
1.52272727272727, 2.2, 1.16666666666667), popsize = c(NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, 2074, 2074, 2074, 2074, 2074, 
2074, 2074, 1546, 1546, 1546, 1546, 1546, 1546, 1546, 1546, 1546, 
2826, 2826, 2826, 2826, 2826, 2826), rate = c(NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, 6.70202507232401, 7.9556412729026, 
7.03953712632594, 6.07521697203472, 7.56991321118611, 8.72709739633558, 
12.8736740597878, 14.7477360931436, 15.006468305304, 10.7373868046572, 
11.1901681759379, 10.8020698576973, 11.1254851228978, 9.57309184993532, 
10.2846054333765, 10.608020698577, 4.77707006369427, 4.88322717622081, 
4.81245576786978, 3.96319886765747, 4.77707006369427, 5.20169851380042
)), row.names = c(NA, -49L), class = c("tbl_df", "tbl", "data.frame"
))

你当前的参考脚本 #

library(mgcv)

bb$rate <- (bb$Tot/bb$popsize)*100

m.bb <- gam(Tot~s(Year,k=6),family=poisson,data=bb)
mbb <- summary(m.bb)

YearP=seq(1970,2018,by=1)
mbb.pred=predict(m.bb,newdata=data.frame(Year=YearP),type="response",se.fit=T)

par(cex.lab=1.3,cex.axis=1.3)

plot(YearP, mbb.pred$fit, type="l", ylim=c(0,300),xlim=c(1970,2020),
     xlab="", ylab=" ",lwd=1.6,las=1)
points(bb$Year,bb$Tot,col=ifelse(bb$rate<4.5|is.na(bb$rate),'black','red'),pch=19,cex=1.7)
lines(YearP,mbb.pred$fit+2*mbb.pred$se.fit,lty="dotted",col="black")
lines(YearP,mbb.pred$fit-2*mbb.pred$se.fit,lty="dotted",col="black")
mtext('BB', side=3, line=0.8, at=1969,adj=0,cex=1)

r2 <- mbb$r.sq
rp = vector('expression',2)
rp[1] = substitute(expression(italic(R)^2 == r2), 
                   list(r2 = format(r2,dig=2)))[2]
legend('topright',adj=c(0.2,0.2),legend=rp,bty='n',cex=1.5)

可实现需求的完整代码 #

library(mgcv)

# 计算率指标
bb$rate <- (bb$Tot/bb$popsize)*100

# --------------------------
# 1. 初始化4×4网格布局
# --------------------------
par(mfrow = c(4,4), # 4行4列网格
    cex.lab = 1.1, 
    cex.axis = 1.1,
    mar = c(3,3,2,1)) # 调整子图边距避免内容重叠

# --------------------------
# 单种群绘图逻辑（16个种群套入循环即可自动填充网格）
# --------------------------
# 拟合GAM模型
m.bb <- gam(Tot~s(Year,k=6),family=poisson,data=bb)
mbb <- summary(m.bb)

# 生成预测序列
YearP <- seq(1970,2018,by=1)
mbb.pred <- predict(m.bb, newdata = data.frame(Year=YearP), type = "response", se.fit = T)

# 绘制基础图层
plot(YearP, mbb.pred$fit, type="l", ylim=c(0,300),xlim=c(1970,2020),
     xlab="年份", ylab="种群总数",lwd=1.6,las=1)
# 添加观测点
points(bb$Year,bb$Tot,col=ifelse(bb$rate<4.5|is.na(bb$rate),'black','red'),pch=19,cex=1.3)
# 添加95%置信区间
lines(YearP,mbb.pred$fit+2*mbb.pred$se.fit,lty="dotted",col="black")
lines(YearP,mbb.pred$fit-2*mbb.pred$se.fit,lty="dotted",col="black")
# 添加种群标识
mtext('BB种群', side=3, line=0.5, adj=0, cex=0.9)

# --------------------------
# 新增：三个统计量分三行展示
# --------------------------
# 提取模型统计量
edf_val <- round(mbb$s.table["s(Year)", "edf"], 2)
r2_val <- round(mbb$r.sq, 2)
p_val <- mbb$s.table["s(Year)", "p-value"]
p_text <- ifelse(p_val < 0.001, "< 0.001", paste0("= ", round(p_val, 3)))

# 构造格式化标签
stat_labels <- c(
  substitute(italic(EDF) == x, list(x = edf_val)),
  substitute(italic(R)^2 == x, list(x = r2_val)),
  substitute(italic(P) ~ x, list(x = p_text))
)

# 添加到图右上角（无边框）
legend("topright", legend = stat_labels, bty = "n", cex = 1)

# --------------------------
# 16个种群扩展说明
# --------------------------
# 若你有16个种群数据，可将上述逻辑套入循环：
# 1. 将16个种群的数据集存储为列表，每个元素对应单个种群的data.frame
# 2. 执行for(i in 1:16) { 上述单种群绘图逻辑 } 即可自动填充4×4网格

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