使用的R语言版本号为：4.3.3

背景描述

在研究国家财政收入时，我们把财政收入按收入形式分为：各项税收收入、企业收入、债务收入、国家能源交通重点建设基金收入、基本建设贷款归还收入、国家预算调节基金收入、其他收入等。为了建立国家财政收入回归模型，我们以财政收入y（亿元）为因变量，自变量如变量说明中所示。从《中国统计年鉴》获得1978-2023年的统计数据。

产业结构变迁：1978 - 2023 年，我国从计划经济向市场经济转型，农业在经济占比中逐步让位于工业、服务业，工业经历从轻工业为主到重化工业、高端制造升级，建筑业随城镇化快速发展，产业结构动态变化深刻影响财政收入来源结构，如工业税收、建筑业相关税费的演变 。

消费与人口维度：人口规模增长与结构老龄化并行，社会消费总额从物资短缺到消费升级、多元化，消费作为经济增长 “三驾马车” 之一，其规模与结构变化，以及人口因素对公共服务、就业、税源的影响，都与财政收入的生成和支出需求紧密交织 。

灾害影响常态：此期间，自然灾害（如旱灾、水灾）频发，农业、基础设施等易受冲击，受灾面积数据反映生态与经济风险，关联财政救灾投入、农业产业恢复对财政的影响，是研究财政收入稳定性不可忽视的背景因素 。

数据选择

自变量中，农业、工业、建筑业增加值反映产业经济规模，产业发展对财政收入（如税收、产业贡献）影响直接；人口数关联税源、消费与公共服务支出对应的财政收支；社会消费总额体现经济活力与内需，影响相关税收及经济循环对财政的支撑；受灾面积关乎农业等产业损失及财政救灾支出，与财政收入存在间接关联 ，这些数据与财政收入因变量高度契合研究逻辑。

研究价值

政策制定参考：通过明确各因素对财政收入的影响程度，助力政府精准施策。如明晰工业增加值对财政收入的高贡献，可针对性扶持工业创新升级；了解受灾面积的影响，能优化财政救灾预算与农业保险等配套政策 。

经济预测应用：构建的回归模型可基于产业、人口、消费等数据变化，预测财政收入趋势，为财政预算安排、债务管理、公共支出规划提供前瞻性依据，保障财政可持续性 。

数据描述

变量说明

变量 单位 说明
year 时间变量
aav 亿元 农业增加值
iav 亿元 工业增加值
cav 亿元 建筑业增加值
pnum 万人 人口数
tsrc 亿元 社会零售消费总额
da 千公顷 受灾面积
y 亿元 财政收入

导入数据和数据概览

data_frame <- read.csv("../data/project.csv", header = TRUE)
knitr::kable(head(data_frame), caption = "数据前10行概览")
数据前10行概览
year aav iav cav pnum tsrc da y
1978年 1018.5 1621.4 138.9 96259 1558.6 50810 1132.26
1979年 1259.0 1786.5 144.6 97542 1800.0 39370 1146.38
1980年 1359.5 2014.8 196.3 98705 2140.0 50030 1159.93
1981年 1545.7 2067.7 208.0 100072 2350.0 39790 1175.79
1982年 1761.7 2183.0 221.6 101654 2570.0 33130 1212.33
1983年 1960.9 2399.0 271.7 103008 2849.4 34710 1366.95

描述分析

在做回归分析之前，我们首先对数据做一个描述性分析，以了解数据的分布情况。

library(skimr)
skim(data_frame)
Data summary
Name data_frame
Number of rows 46
Number of columns 8
_______________________
Column type frequency:
character 1
numeric 7
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
year 0 1 5 5 0 46 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
aav 0 1 27132.42 27256.26 1018.50 4425.45 15109.95 48008.82 89755.2 ▇▂▁▂▁
iav 0 1 107297.53 124587.81 1621.40 6620.25 42056.40 205460.82 399103.1 ▇▁▂▂▁
cav 0 1 20530.37 26728.69 138.90 824.27 5739.75 35903.70 85691.1 ▇▁▁▁▁
pnum 0 1 123848.22 14119.85 96259.00 113111.25 127185.00 135670.50 141260.0 ▃▂▃▆▇
tsrc 0 1 118827.32 148947.61 1558.60 8151.07 40343.75 199088.92 467098.0 ▇▁▁▁▂
da 0 1 38700.63 13042.64 10539.00 31485.00 40540.00 49732.50 55470.0 ▂▂▃▃▇
y 0 1 58027.37 73076.58 1132.26 2732.95 14890.64 113908.75 216795.4 ▇▁▁▁▂

从数据分析结果来看，所有数值型变量呈现出显著的右偏分布特征（均值大于中位数），其中iav、cav、tsrc和y的极端右偏现象尤为突出，反映出经济指标在后期可能出现爆发式增长或存在特殊事件驱动。变量间对比显示，aav（均值27,132）与iav（均值107,297）的量级差异暗示农业与工业领域的经济活动规模存在显著差距，而cav在均值20,530的基础上最大值达到85,691，表明该指标在2020年后呈现加速增长态势。时间序列分析揭示了明显的长期增长趋势，特别是2000年后多数指标进入指数增长阶段，但pnum（人口数量指标）标准差仅14,119.85，其稳定增长特征与经济类变量形成鲜明对比。值得注意的是，y在2020-2023年突破21万大关的同时，tsrc在2021年出现43.8万的峰值，反映出社会零售消费总额对财政收入贡献显著增加。

数据相关性

knitr::kable(cor(data_frame[, -1]), caption = "数据相关性")
数据相关性
aav iav cav pnum tsrc da y
aav 1.0000000 0.9964996 0.9881188 0.8495634 0.9878964 -0.8449471 0.9881722
iav 0.9964996 1.0000000 0.9925025 0.8198523 0.9915681 -0.8638416 0.9928630
cav 0.9881188 0.9925025 1.0000000 0.7756141 0.9986561 -0.8780677 0.9930390
pnum 0.8495634 0.8198523 0.7756141 1.0000000 0.7863979 -0.5380266 0.7938737
tsrc 0.9878964 0.9915681 0.9986561 0.7863979 1.0000000 -0.8758818 0.9954968
da -0.8449471 -0.8638416 -0.8780677 -0.5380266 -0.8758818 1.0000000 -0.8798825
y 0.9881722 0.9928630 0.9930390 0.7938737 0.9954968 -0.8798825 1.0000000

从相关系数矩阵中，我们可以看到各变量之间的相关关系。例如，财政收入（y）与农业增加值（aav）、工业增加值（iav）、建筑业增加值（cav）和税收收入（tsrc）之间存在高度正相关关系，而与受灾面积（da）之间存在负相关关系。这些信息对于理解各变量之间的相互关系和建立回归模型非常重要。

模型建立

一元线性回归模型

首先，我们使用一元线性回归模型来分析财政收入与建筑业增加值之间的关系。

lm1 <- lm(y ~ cav, data = data_frame)
summary(lm1)
## 
## Call:
## lm(formula = y ~ cav, data = data_frame)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -18142  -3864  -1660   3112  20310 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.288e+03  1.625e+03   1.408    0.166    
## cav         2.715e+00  4.855e-02  55.924   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8705 on 44 degrees of freedom
## Multiple R-squared:  0.9861, Adjusted R-squared:  0.9858 
## F-statistic:  3127 on 1 and 44 DF,  p-value: < 2.2e-16

根据回归分析结果，财政收入（y）与建筑业增加值（cav）之间存在高度显著的正向线性关系。模型显示，cav的回归系数为2.715（p<0.001），表明在样本期间，建筑业增加值每增加1亿元，财政收入平均增加2.715亿元。该模型的拟合优度\(R^2\)达到0.9861，说明建筑业增加值能够解释财政收入变动的98.6%。同时，残差分析显示标准误为8705，最大正向偏离约20310亿元，表明在个别年份可能存在其他未纳入模型的因素（如政策调整或非税收入变化）导致实际财政收入偏离趋势线。从统计显著性（F统计量3127，p<2.2e-16）和经济意义双重维度看，建筑业发展对财政收入具有决定性作用，这一发现对财政规划和产业结构调整具有重要参考价值。

残差分析

根据上面的分析，我们进行残差分析，以检查模型的假设是否成立。

e <- rstandard(lm1)
plot(data_frame$cav, e, xlab = "cav", ylab = "Standardized Residuals", main = "Residuals plot")
abline(h = 0, col = "red")

该残差图显示标准化残差与自变量 cav 之间存在潜在的异方差性问题，表现为低值区域（cav 接近 0）的残差聚集且部分点偏离零线较远，而高值区域（如 cav > 60000）的残差波动范围扩大，误差方差随自变量变化而变化；同时，低值区密集的残差点可能反映模型对局部数据的拟合不足，需警惕异常值或未被捕捉的非线性关系。为改进模型，使用多元回归模型，并引入其他可能影响财政收入的变量。

异方差检验

在多元回归之前，对模型做异方差检验。

abse <- abs(e)
cor.test(data_frame$cav, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$cav and abse
## S = 4564, p-value = 1.047e-07
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.7185322

尝试运用最优权函数处理异方差

s <- seq(-2, 2, 0.5)
result1 <- vector(length = 9, mode = "list")
result2 <- vector(length = 9, mode = "list")
result3 <- vector(length = 9, mode = "list")
for (j in 1:9)
{
  w <- data_frame$cav^(-s[j])
  lm4 <- lm(y ~ cav, weights = w, data_frame)
  result1[[j]] <- logLik(lm4)
  result2[[j]] <- summary(lm4)
  result3[[j]] <- abs(rstandard(lm4))
}
result1
## [[1]]
## 'log Lik.' -565.4661 (df=3)
## 
## [[2]]
## 'log Lik.' -545.1545 (df=3)
## 
## [[3]]
## 'log Lik.' -525.2086 (df=3)
## 
## [[4]]
## 'log Lik.' -504.3782 (df=3)
## 
## [[5]]
## 'log Lik.' -481.5431 (df=3)
## 
## [[6]]
## 'log Lik.' -457.869 (df=3)
## 
## [[7]]
## 'log Lik.' -438.674 (df=3)
## 
## [[8]]
## 'log Lik.' -429.1203 (df=3)
## 
## [[9]]
## 'log Lik.' -427.2758 (df=3)
result2[which.max(result1)]
## [[1]]
## 
## Call:
## lm(formula = y ~ cav, data = data_frame, weights = w)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0982 -0.3147  0.0759  0.4387  0.7180 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 711.80212   41.62660   17.10   <2e-16 ***
## cav           2.57946    0.08908   28.96   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5214 on 44 degrees of freedom
## Multiple R-squared:  0.9501, Adjusted R-squared:  0.949 
## F-statistic: 838.5 on 1 and 44 DF,  p-value: < 2.2e-16
cor.test(data_frame$cav, result3[[which.max(result1)]], method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$cav and result3[[which.max(result1)]]
## S = 17974, p-value = 0.4717
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.1084798

使用最优权函数成功消除了异方差性。

自相关性检验

library(lmtest)
## 载入需要的程辑包：zoo
## 
## 载入程辑包：'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
dwtest(lm1, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  lm1
## DW = 0.1425, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is not 0

由DW值为0.1425知，说明残差存在自相关性，尝试使用差分消除自相关性。

dlm1 <- lm(diff(y)~diff(cav)-1,data=data_frame)
summary(dlm1)
## 
## Call:
## lm(formula = diff(y) ~ diff(cav) - 1, data = data_frame)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11892.7   -455.0    110.9   1198.5   6841.4 
## 
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)    
## diff(cav)   2.4583     0.1666   14.76   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3201 on 44 degrees of freedom
## Multiple R-squared:  0.832,  Adjusted R-squared:  0.8281 
## F-statistic: 217.8 on 1 and 44 DF,  p-value: < 2.2e-16
dwtest(dlm1, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  dlm1
## DW = 1.3952, p-value = 0.03788
## alternative hypothesis: true autocorrelation is not 0

由DW值显著接近2，说明差分后残差不再存在自相关性。

多元线性回归模型

lm_m <- lm(y ~ . - year, data = data_frame)
summary(lm_m)
## 
## Call:
## lm(formula = y ~ . - year, data = data_frame)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11308.7  -3094.1   -429.1   2717.1   9092.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.995e+04  2.010e+04   2.486  0.01733 *  
## aav          4.237e-01  5.328e-01   0.795  0.43132    
## iav          3.406e-01  9.616e-02   3.542  0.00105 ** 
## cav         -3.674e+00  8.094e-01  -4.539 5.30e-05 ***
## pnum        -4.541e-01  2.166e-01  -2.096  0.04263 *  
## tsrc         8.121e-01  1.165e-01   6.972 2.34e-08 ***
## da          -1.240e-01  1.477e-01  -0.839  0.40642    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5075 on 39 degrees of freedom
## Multiple R-squared:  0.9958, Adjusted R-squared:  0.9952 
## F-statistic:  1548 on 6 and 39 DF,  p-value: < 2.2e-16

从上面的结果可以看到，模型拟合优度\(R^2\)为0.9958，说明模型能够解释财政收入变动的99.58%。同时，模型中由很多变量不显著。

# 回归方程显著性检验
library(car)
## 载入需要的程辑包：carData
Anova(lm_m, type = 3)
## Anova Table (Type III tests)
## 
## Response: y
##                 Sum Sq Df F value    Pr(>F)    
## (Intercept)  159147743  1  6.1780  0.017329 *  
## aav           16288915  1  0.6323  0.431317    
## iav          323149250  1 12.5444  0.001048 ** 
## cav          530716300  1 20.6020 5.300e-05 ***
## pnum         113163056  1  4.3929  0.042628 *  
## tsrc        1252073996  1 48.6046 2.344e-08 ***
## da            18146177  1  0.7044  0.406417    
## Residuals   1004655278 39                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

逐步回归

首先选择最优子集

library(leaps)
exps <- regsubsets(y ~ . - year, data = data_frame, nbest = 1, really.big = TRUE) # nebst = 1表示只输出单变量的最优子集
expres <- summary(exps)

res <- data.frame(expres$outmat, 调整R平方 = expres$adjr2)
knitr::kable(res, caption = "调整R平方")
调整R平方
aav iav cav pnum tsrc da 调整R平方
1 ( 1 ) * 0.9908097
2 ( 1 ) * * 0.9926637
3 ( 1 ) * * * 0.9941103
4 ( 1 ) * * * * 0.9952353
5 ( 1 ) * * * * * 0.9952205
6 ( 1 ) * * * * * * 0.9951761 
res <- data.frame(expres$outmat, Cp = expres$cp)
knitr::kable(res, caption = "CP准则")
CP准则
aav iav cav pnum tsrc da Cp
1 ( 1 ) * 41.827053
2 ( 1 ) * * 25.395630
3 ( 1 ) * * * 13.279877
4 ( 1 ) * * * * 4.496959
5 ( 1 ) * * * * * 5.632324
6 ( 1 ) * * * * * * 7.000000 
lm_step <- step(lm_m, direction = "both")
## Start:  AIC=791.37
## y ~ (year + aav + iav + cav + pnum + tsrc + da) - year
## 
##        Df  Sum of Sq        RSS    AIC
## - aav   1   16288915 1020944193 790.11
## - da    1   18146177 1022801455 790.19
## <none>               1004655278 791.37
## - pnum  1  113163056 1117818334 794.28
## - iav   1  323149250 1327804528 802.19
## - cav   1  530716300 1535371578 808.88
## - tsrc  1 1252073996 2256729274 826.59
## 
## Step:  AIC=790.11
## y ~ iav + cav + pnum + tsrc + da
## 
##        Df  Sum of Sq        RSS    AIC
## - da    1   22273335 1043217528 789.10
## <none>               1020944193 790.11
## + aav   1   16288915 1004655278 791.37
## - pnum  1  127007228 1147951421 793.50
## - cav   1  558299481 1579243674 808.17
## - iav   1  788344583 1809288776 814.43
## - tsrc  1 1277151144 2298095337 825.43
## 
## Step:  AIC=789.1
## y ~ iav + cav + pnum + tsrc
## 
##        Df  Sum of Sq        RSS    AIC
## <none>               1043217528 789.10
## + da    1   22273335 1020944193 790.11
## + aav   1   20416073 1022801455 790.19
## - pnum  1  277772181 1320989709 797.96
## - cav   1  635122400 1678339929 808.97
## - iav   1  970339994 2013557522 817.35
## - tsrc  1 1440169027 2483386555 827.00
summary(lm_step)
## 
## Call:
## lm(formula = y ~ iav + cav + pnum + tsrc, data = data_frame)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11686.0  -2962.8   -300.7   1879.0  10584.9 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.014e+04  1.341e+04   2.994  0.00466 ** 
## iav          4.130e-01  6.688e-02   6.175 2.44e-07 ***
## cav         -3.538e+00  7.081e-01  -4.996 1.14e-05 ***
## pnum        -4.041e-01  1.223e-01  -3.304  0.00198 ** 
## tsrc         8.100e-01  1.077e-01   7.523 3.04e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5044 on 41 degrees of freedom
## Multiple R-squared:  0.9957, Adjusted R-squared:  0.9952 
## F-statistic:  2351 on 4 and 41 DF,  p-value: < 2.2e-16

结论与应用

根据多元线性回归模型的结果，我们可以得出以下结论：

根据提供的回归结果，我们得到了以下回归方程：

y = 40140 + 0.413 * iav - 3.538 * cav - 0.404 * pnum + 0.810 * tsrc

模型诊断和检验

多重共线性检验

library(car)
vif(lm_step) # 方差扩大因子>10，说明存在多重共线性
##        iav        cav       pnum       tsrc 
## 122.773342 633.605303   5.274357 454.865127

尝试使用box-cox变换来处理，从结果看，多重共线性的并没有消除

library(MASS)
bc3.2 <- boxcox(y ~ iav + cav + pnum + tsrc, data = data_frame, lambda = seq(-2, 2, 0.01))

summary(bc3.2)
##   Length Class  Mode   
## x 401    -none- numeric
## y 401    -none- numeric
lambda <- bc3.2$x[which.max(bc3.2$y)]
lambda
## [1] -0.4
y_bc <- (data_frame$y^lambda - 1) / lambda
lm3.2_bc <- lm(y_bc ~ iav + cav + pnum + tsrc, data = data_frame)
summary(lm3.2_bc)
## 
## Call:
## lm(formula = y_bc ~ iav + cav + pnum + tsrc, data = data_frame)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.006254 -0.001567  0.000422  0.001028  0.006996 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.035e+00  8.068e-03 252.243  < 2e-16 ***
## iav          2.477e-07  4.024e-08   6.156  2.6e-07 ***
## cav         -1.426e-06  4.261e-07  -3.347  0.00176 ** 
## pnum         3.195e-06  7.359e-08  43.412  < 2e-16 ***
## tsrc         3.955e-08  6.478e-08   0.610  0.54495    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.003035 on 41 degrees of freedom
## Multiple R-squared:  0.996,  Adjusted R-squared:  0.9956 
## F-statistic:  2547 on 4 and 41 DF,  p-value: < 2.2e-16
vif(lm3.2_bc)
##        iav        cav       pnum       tsrc 
## 122.773342 633.605303   5.274357 454.865127

多元线性回归的异方差

abse <- abs(rstandard(lm_step))
cor.test(data_frame$iav, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$iav and abse
## S = 6558, p-value = 1.876e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5955597
cor.test(data_frame$cav, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$cav and abse
## S = 6562, p-value = 1.894e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.595313
cor.test(data_frame$pnum, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$pnum and abse
## S = 6548, p-value = 1.831e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5961764
cor.test(data_frame$tsrc, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$tsrc and abse
## S = 6610, p-value = 2.122e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5923528

这些变量都具有一致的异方差性，使用对数变换来尝试消除异方差性。

lm_step2 <- lm(log(y) ~ iav + cav + pnum + tsrc, data = data_frame)
summary(lm_step2)
## 
## Call:
## lm(formula = log(y) ~ iav + cav + pnum + tsrc, data = data_frame)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47870 -0.10694  0.03154  0.08516  0.63707 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.452e-01  4.998e-01  -0.491 0.626353    
## iav          2.241e-05  2.492e-06   8.990 3.03e-11 ***
## cav         -1.588e-04  2.639e-05  -6.018 4.08e-07 ***
## pnum         7.230e-05  4.558e-06  15.862  < 2e-16 ***
## tsrc         1.539e-05  4.013e-06   3.835 0.000425 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.188 on 41 degrees of freedom
## Multiple R-squared:  0.9908, Adjusted R-squared:  0.9899 
## F-statistic:  1105 on 4 and 41 DF,  p-value: < 2.2e-16

显著性由很明显的提高，然后检验异方差性

abse <- abs(rstandard(lm_step2))
cor.test(data_frame$iav, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$iav and abse
## S = 13846, p-value = 0.3315
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1460993
cor.test(data_frame$cav, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$cav and abse
## S = 13854, p-value = 0.3331
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1456059
cor.test(data_frame$pnum, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$pnum and abse
## S = 13814, p-value = 0.325
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1480728
cor.test(data_frame$tsrc, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data_frame$tsrc and abse
## S = 13866, p-value = 0.3356
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1448659

异方差性已经消除，因此可以认为该模型是有效的。

多元线性回归的自相关性

library(lmtest)
dwtest(lm_step2, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  lm_step2
## DW = 1.1149, p-value = 0.0001341
## alternative hypothesis: true autocorrelation is not 0

接近0还是存在自相关性。

主成分分析

使用主成分分析来尝试解决多重共线性问题。

com1 <- prcomp(data_frame[, -c(1, ncol(data_frame))], center = TRUE, scale. = TRUE)
summary(com1)
## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5     PC6
## Standard deviation     2.3284 0.68861 0.30428 0.09503 0.04528 0.02698
## Proportion of Variance 0.9036 0.07903 0.01543 0.00151 0.00034 0.00012
## Cumulative Proportion  0.9036 0.98260 0.99803 0.99954 0.99988 1.00000
# 查看主成分载荷（即旋转后的系数）
knitr::kable(com1$rotation, caption = "主成分载荷")
主成分载荷
PC1 PC2 PC3 PC4 PC5 PC6
aav 0.4278204 -0.0736293 0.1543286 0.4503279 0.7189035 0.2609859
iav 0.4278084 0.0014473 0.2027416 0.5833026 -0.6599580 -0.0093445
cav 0.4257717 0.0959679 0.3540469 -0.2905008 0.1390044 -0.7618704
pnum 0.3580190 -0.7718343 -0.4915463 -0.1585524 -0.0596769 -0.0760013
tsrc 0.4262403 0.0755326 0.3099598 -0.5887409 -0.1569869 0.5876038
da -0.3779017 -0.6196236 0.6870491 0.0285902 -0.0102393 0.0172666

重新计算VIF值

pc_scores <- com1$x  # 主成分得分
pc_data <- as.data.frame(pc_scores[, 1:2])
y <- data_frame$y
model_pcr <- lm(y ~ PC1 + PC2, data = pc_data)
summary(model_pcr)
## 
## Call:
## lm(formula = y ~ PC1 + PC2, data = pc_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17600.3  -6617.2    171.1   6366.4  16244.2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  58027.4     1230.1  47.173  < 2e-16 ***
## PC1          31155.6      534.1  58.328  < 2e-16 ***
## PC2           4874.1     1806.1   2.699  0.00991 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8343 on 43 degrees of freedom
## Multiple R-squared:  0.9875, Adjusted R-squared:  0.987 
## F-statistic:  1705 on 2 and 43 DF,  p-value: < 2.2e-16
vif(model_pcr)
## PC1 PC2 
##   1   1

可以看出多重共线性已经消除，并且可以看到主成分因子1贡献最大，占了绝大多数的方差，因此可以认为主成分因子1已经包含了大部分信息。查看主成分的载荷，可以认为除了受灾面积其他对于国家财政收入均具有正向作用。

abse <- abs(rstandard(model_pcr))
cor.test(pc_data$PC1, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  pc_data$PC1 and abse
## S = 12830, p-value = 0.1634
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2087573
cor.test(pc_data$PC2, abse, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  pc_data$PC2 and abse
## S = 12294, p-value = 0.1054
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2418131
library(lmtest)
dwtest(model_pcr, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  model_pcr
## DW = 1.0283, p-value = 0.0001048
## alternative hypothesis: true autocorrelation is not 0

回归分析结论

综合以上分析结果，我们可以从以下几个方面进行解读：

一、回归模型公式

这是一个多元线性回归模型，因变量是财政收入 y 的对数形式，自变量包括：工业增加值（iav）、建筑业增加值（cav）、人口数（pnum）和社会零售消费总额（tsrc）。

回归模型表达式为：

\[ \log(y) = -0.2452 + 2.241 \times 10^{-5} \cdot iav - 1.588 \times 10^{-4} \cdot cav + 7.230 \times 10^{-5} \cdot pnum + 1.539 \times 10^{-5} \cdot tsrc + \varepsilon \]

其中： - \(y\)：财政收入（亿元） - \(iav\)：工业增加值（亿元） - \(cav\)：建筑业增加值（亿元） - \(pnum\)：人口数（万人） - \(tsrc\)：社会零售消费总额（亿元） - \(\varepsilon\)：误差项

二、系数解释（注意是对数变换）

由于因变量是 \(\log(y)\)，因此系数的解释应基于百分比变化：

变量 系数估计值 解释（每单位增加的影响）
iav 2.241e-05 工业增加值每增加1亿元，财政收入平均增长约 0.002241%
cav -1.588e-04 建筑业增加值每增加1亿元，财政收入平均下降约 0.01588%
pnum 7.230e-05 人口每增加1万人，财政收入平均增长约 0.00723%
tsrc 1.539e-05 社会零售消费总额每增加1亿元，财政收入平均增长约 0.001539%

三、统计显著性分析

所有自变量在 0.05 水平上均显著，尤其是：

  • iav（工业增加值）：t = 8.99, p < 0.001
  • cav（建筑业增加值）：t = -6.02, p < 0.001
  • pnum（人口数）：t = 15.86, p < 0.001
  • tsrc（社会零售消费总额）：t = 3.84, p = 0.0004

这表明它们对财政收入具有统计意义上的显著影响。

四、拟合优度与整体显著性

  • R² = 0.9908：模型可以解释财政收入对数的 99.08% 的变异，说明模型拟合效果非常好。
  • 调整 R² = 0.9899：考虑了变量数量后的拟合优度，依然非常高。
  • F 统计量 = 1105，p < 2.2e-16：整个模型在统计上高度显著。

五、残差分析简要判断

  • 最小值：-0.47870
  • 第一四分位数：-0.10694
  • 中位数：0.03154
  • 第三四分位数：0.08516
  • 最大值：0.63707

残差分布基本对称，没有极端偏移。

六、经济意义分析

自变量 影响 说明
iav 正 工业发展通常带动税收增长
cav 负 建筑业对财政贡献可能受政策或成本影响
pnum 正 人口多意味着税源广
tsrc 正 消费活跃带来增值税等税收增长
  1. 模型整体显著且拟合良好（R² 接近 0.99），说明所选变量能很好地解释财政收入的变化。
  2. 主要驱动因素为：
    • 工业增加值（正向）
    • 人口数量（正向）
    • 社会零售总额（正向）
  3. 建筑业增加值系数为负，这是因为建筑业通常是资本密集型行业，企业在购置设备、建筑材料等方面的投入较大，折旧和成本抵扣较多，从而导致其税收贡献相对较低。此外，建筑项目周期较长，资金回笼慢，企业盈利能力有限，进一步限制了其对财政收入的直接贡献。同时，建筑业的成本结构复杂，涉及大量非标准化合同和隐蔽工程，容易引发偷税漏税行为，使得实际税收征管难度加大。再者，在某些时期可能存在政策扶持或减免税优惠，也降低了该行业的税收贡献。最后，建筑业的发展往往伴随着政府的基础设施投资，这些投资本身可能带来财政支出的增加，从而在短期内对财政赤字产生压力，间接表现出与财政收入的负相关性。

本研究采用多元线性回归方法，建立财政收入 \(\log(y)\) 与工业增加值（iav）、建筑业增加值（cav）、人口数（pnum）和社会零售消费总额（tsrc）之间的关系模型。结果显示，模型具有高度显著性（F = 1105, p < 0.001），并能解释财政收入 99.08% 的变异。其中，工业增加值、人口数量和社会零售总额对财政收入具有显著正向作用，而建筑业增加值呈显著负相关。该结果提示，在制定财政政策时应重点关注工业和消费领域的发展，同时关注建筑业对财政的潜在影响。