GLM准泊松模型二次项纳入问题及自由度异常问询

研究背景 #

通过全国所有投喂点同步投喂开展秃鹫种群普查：2004年及2006-2023年期间，每年6月开展2次间隔10-14天的同步普查，监测3种秃鹫（RHV、SBV、WRV）的取食个体数量。

针对每个物种，拟合广义线性模型（GLM），以年份的线性/二次效应、抽样场合（每年2次普查的第1/第2次）的加性固定效应为自变量，建模所有位点的总计数。若首次普查后秃鹫未广泛扩散，更可能出现在第二次普查中，因此普查场合是重要变量。但建模后发现年份二次项未被准确纳入，自由度与预期不符。

数据集子集 #

Date    RHV SBV WRV Total   Census  Month   Year    Season  Occasion
6/10/2004   40  34  88  162 Yes 6   2004    Wet 1
6/20/2004   42  25  90  157 Yes 6   2004    Wet 2
6/10/2005   58  27  149 234 Yes 6   2005    Wet 1
6/20/2005   32  31  83  146 Yes 6   2005    Wet 2
6/10/2007   35  24  160 219 Yes 6   2007    Wet 1
6/20/2007   40  26  150 216 Yes 6   2007    Wet 2
6/10/2008   48  30  113 191 Yes 6   2008    Wet 1
6/20/2008   44  51  191 286 Yes 6   2008    Wet 2
6/10/2009   11  84  50  145 Yes 6   2009    Wet 1
6/20/2009   13  30  209 252 Yes 6   2009    Wet 2

数据摘要（summary(data)结果） #

Date_Num列未用于本次分析，以2004年6月1日为基准计算天数

Date                 RHV             SBV             WRV             Total                      Month  
 Min.   :2004-06-10   Min.   :10.00   Min.   :15.00   Min.   : 42.00   Min.   : 70.0         
 1st Qu.:2009-09-16   1st Qu.:15.25   1st Qu.:28.50   1st Qu.: 71.25   1st Qu.:122.0    
 Median :2014-06-15   Median :21.00   Median :35.00   Median : 89.00   Median :155.0   
 Mean   :2014-05-07   Mean   :25.89   Mean   :37.29   Mean   :104.24   Mean   :167.4                        
 3rd Qu.:2019-03-13   3rd Qu.:38.00   3rd Qu.:44.50   3rd Qu.:136.00   3rd Qu.:214.8                        
 Max.   :2023-06-16   Max.   :58.00   Max.   :84.00   Max.   :209.00   Max.   :289.0                       
      Year        Occasion      Date_Num   
 Min.   :2004   Min.   :1.0   Min.   :   9  
 1st Qu.:2009   1st Qu.:1.0   1st Qu.:1934  
 Median :2014   Median :1.5   Median :3666  
 Mean   :2014   Mean   :1.5   Mean   :3627  
 3rd Qu.:2019   3rd Qu.:2.0   3rd Qu.:5398  
 Max.   :2023   Max.   :2.0   Max.   :6954 

建模代码与结果 #

使用准泊松误差和对数链接函数构建GLM：

modelRHVDate <- glm(Species ~ Year + I(Year^2) + Occasion + Year:Occasion + I(Year^2):Occasion , 
                family = quasipoisson(link = "log"), 
                data = data)
summary(modelRHVDate)

模型输出：

Call:
glm(formula = Species ~ Year + I(Year^2) + Occasion + Year:Occasion + 
    I(Year^2):Occasion, family = quasipoisson(link = "log"), 
    data = data)

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)
(Intercept)         2.560e+04  2.341e+04   1.094    0.282
Year               -2.535e+01  2.327e+01  -1.090    0.284
I(Year^2)           6.274e-03  5.780e-03   1.086    0.286
Occasion           -1.531e+04  1.475e+04  -1.037    0.307
Year:Occasion       1.519e+01  1.466e+01   1.036    0.308
I(Year^2):Occasion -3.769e-03  3.641e-03  -1.035    0.308

(Dispersion parameter for quasipoisson family taken to be 2.874796)

    Null deviance: 228.254  on 37  degrees of freedom
Residual deviance:  99.214  on 32  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

模型选择过程 #

以最低准AICc（QAICc）值选择最优模型，通过QAICc权重展示支持度。计算时先以标准泊松误差重新拟合所有模型，再用R包bbmle（Bolker 2016），以最复杂模型（date2 + occasion）的离散度校正AICc值。

部分模型构建代码：

# 准泊松模型
modelQP_RHV_Year <- glm(Species ~ Year, family = quasipoisson(link = "log"), data = data)
modelQP_RHV_Year2 <- glm(Species ~ I(Year^2), family = quasipoisson(link = "log"), data = data)
# 泊松模型
modelP_RHV_Year <-  update(modelQP_RHV_Year, family=poisson)
modelP_RHV_Year2 <- update(modelQP_RHV_Year2, family=poisson)

模型选择结果：

dqAICc df Weight
modelP_RHV_Year            0.0   2   0.4  
modelP_RHV_Year2           0.0   2   0.4  
modelP_RHV_Year_Occasion   1.9   3   0.1  
modelP_RHV_Year2_Occasion  1.9   3   0.1  
modelP_RHV_Occasion       42.4   2   0.0  

核心问题 #

原本预期：

• 年份二次项模型（~I(Year^2)）自由度应为3（截距+Year+Year²）
• 年份二次项+场合模型自由度应为4

但实际模型选择结果中，modelP_RHV_Year2自由度仅为2，说明二次项未被正确纳入。已用I()包裹二次项（R中^在公式中有特殊含义，非数学幂运算），但问题仍存在，寻求解决方法。

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