偏回归

Partial regression

发布于

2024年11月7日

本文使用R程序包Keng及其中的数据depress，采用四种方式检验在控制了情绪定向的应对方式(em1)、回避型应对方式(am1) 后，问题定向的应对方式(pm1)对抑郁(dm1)的独特效应。这四种方式是:

多元回归及回归系数的t检验。
分层回归及其F检验。
单个参数（pm1的偏回归系数）的PRE。
使用残差检验pm1的独特效应。

首先安装R与RStudio软件。

之后安装Keng的最新版：

install.packages("devtools")
devtools::install_github("qyaozh/Keng", dependencies = TRUE, build_vignettes = TRUE)

在R中加载Keng程序包，并调用depress数据。

library(Keng)
data(depress)

1 多元回归

首先建立两个多元回归模型。Model C （Compact model）做dm1对em1和am1的回归。Model A（Augmented model）做dm1对pm1、 em1和am1的回归。

# multiple regression
fitC <- lm(dm1 ~ em1 + am1, depress)
fitA <- lm(dm1 ~ pm1 + em1 + am1, depress)
summary(fitA)


Call:
lm(formula = dm1 ~ pm1 + em1 + am1, data = depress)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.63018 -0.24748 -0.00681  0.21045  1.01320 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.10497    0.25517   8.249 1.25e-12 ***
pm1         -0.16705    0.04988  -3.349  0.00119 ** 
em1          0.19504    0.05712   3.415  0.00096 ***
am1         -0.06675    0.05992  -1.114  0.26822    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.337 on 90 degrees of freedom
Multiple R-squared:  0.2491,    Adjusted R-squared:  0.224 
F-statistic: 9.949 on 3 and 90 DF,  p-value: 9.928e-06

可见，pm1的偏回归系数是-0.16705, t(90) = -3.349, p = 0.00119。

2 分层回归

第二，通过F检验对上面的分层回归的两个回归模型Model C与Model A进行比较。与Model C相比，Model A增加了 pm1这一个自变量。在SPSS中，该F检验被视为对R² change的检验。

anova(fitC, fitA)

Analysis of Variance Table

Model 1: dm1 ~ em1 + am1
Model 2: dm1 ~ pm1 + em1 + am1
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1     91 11.498                                
2     90 10.224  1    1.2743 11.217 0.001185 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

可见，F (1, 90) = 11.217, p = 0.001185. 该F检验与Model A 多元回归中对pm1的t检验是等价的，它们都检验了控制em1与am1的作用后、pm1的独特效应。此时，F分子上的df为1，F = t²，F分母的df即为t的df。

3 单个参数的PRE

第三，使用PRE检验pm1的独特效应。

format(compare_lm(fitC, fitA), digits = 3, nsmall = 3)

                     Baseline   C          A          A vs. C   
SSE                  "1.36e+01" "1.15e+01" "1.02e+01" "1.27e+00"
n                    "9.40e+01" "9.40e+01" "9.40e+01" "9.40e+01"
Number of parameters "1.00e+00" "3.00e+00" "4.00e+00" "1.00e+00"
df                   "9.30e+01" "9.10e+01" "9.00e+01" "1.00e+00"
R_squared            "      NA" "1.55e-01" "2.49e-01" "9.36e-02"
f_squared            "      NA" "1.84e-01" "3.32e-01" "1.25e-01"
R_squared_adj        "      NA" "1.37e-01" "2.24e-01" "      NA"
PRE                  "      NA" "1.55e-01" "2.49e-01" "1.11e-01"
F(PA-PC,n-PA)        "      NA" "8.38e+00" "9.95e+00" "1.12e+01"
p                    "      NA" "4.58e-04" "9.93e-06" "1.19e-03"
PRE_adj              "      NA" "1.37e-01" "2.24e-01" "1.01e-01"
power_post           "      NA" "9.59e-01" "9.97e-01" "9.12e-01"

可见，F (1, 90) = 11.217, p = 0.00119。此时，PRE的F检验与分层回归的F检验是等价的。

4 使用残差进行回归

第四，使用残差检验pm1的独特效应。

做dm1对em1和am1的回归，得到dm1的残差dm_res，该残差排除了em1和am1对dm1的效应。

再做pm1对em1和am1的回归，得到pm1的残差pm_res，该残差排除了em1和am1对pm1的效应。

之后做dm_res和pm_res的相关，我们就得到了dm1和pm1的偏相关系数/净相关系数。

dm_res <- lm(dm1 ~ em1 + am1, depress)$residuals
pm_res <- lm(pm1 ~ em1 + am1, depress)$residuals
resDat <- data.frame(dm_res, pm_res)
cor(dm_res, pm_res)

[1] -0.3329009

可见，dm1和pm1的偏相关系数为-0.3329009。

做dm_res对pm_res的一元回归，我们就得到了pm1对dm1的独特效应。

summary(lm(dm_res ~ pm_res, data.frame(dm_res, pm_res)))


Call:
lm(formula = dm_res ~ pm_res, data = data.frame(dm_res, pm_res))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.63018 -0.24748 -0.00681  0.21045  1.01320 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  5.153e-17  3.438e-02   0.000  1.00000   
pm_res      -1.670e-01  4.933e-02  -3.386  0.00104 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3334 on 92 degrees of freedom
Multiple R-squared:  0.1108,    Adjusted R-squared:  0.1012 
F-statistic: 11.47 on 1 and 92 DF,  p-value: 0.001044

可见，pm_res的回归系数与前文多元回归Model A中pm1的偏回归系数相同。然而，它们的t与p是不同的。为什么呢？前文用PRE检验了pm1的独特效应，这里我们使用PRE检验pm_res的独特效应，然后将它们进行比较。注意，单个回归系数的PRE的F检验与该回归系数的t检验是等价的；另外，Model A与Model C是相对的关系，统计目的不同，Model C与Model A所指的模型也不同。

fitC <- lm(dm_res ~ 1, resDat)
fitA <- lm(dm_res ~ pm_res, resDat)
format(compare_lm(fitC, fitA), digits = 3, nsmall = 3)

                     Baseline   C          A          A vs. C   
SSE                  "11.49827" "11.49827" "10.22400" " 1.27427"
n                    "94.00000" "94.00000" "94.00000" "94.00000"
Number of parameters " 1.00000" " 1.00000" " 2.00000" " 1.00000"
df                   "93.00000" "93.00000" "92.00000" " 1.00000"
R_squared            "      NA" " 0.00000" " 0.11082" " 0.11082"
f_squared            "      NA" " 0.00000" " 0.12464" " 0.12464"
R_squared_adj        "      NA" " 0.00000" " 0.10116" "      NA"
PRE                  "      NA" " 0.00000" " 0.11082" " 0.11082"
F(PA-PC,n-PA)        "      NA" "      NA" "11.46646" "11.46646"
p                    "      NA" "      NA" " 0.00104" " 0.00104"
PRE_adj              "      NA" " 0.00000" " 0.10116" " 0.10116"
power_post           "      NA" "      NA" " 0.91784" " 0.91784"

比较pm_res和pm1的PRE。可见，PRE是相等的，即pm1的独特效应本质上没有变化。但是，PRE的F检验的df2是不同的，这使得F与p值也发生了变化。换句话说，pm1的独特效应PRE是不变的，PRE基于模型比较检验PRE的显著性，但用来检验PRE的显著性的Model C与Model A是不同的，因而模型比较的结论（i.e., F检验和t检验的结果）也不同。打个比方，一勺盐放在一杯水里感觉齁咸，但一勺盐放在一锅水里感觉咸淡刚好。

回想PRE的F检验的公式，我们将得到以下结论：PRE保持不变，则PRE的显著性（F值与p值）是由Model C的自由度与Model A新增的参数的数量决定的。

因此，我们也将得到以下结论：给定一组自变量的PRE，这一组自变量的PRE达到显著的统计检验力是由样本量n和回归模型中参数的总数决定的。同理，给定一组自变量的PRE、回归模型中参数的总数，以及所需的统计检验力，我们可以计算出所需的样本量——这便是通过回归模型计划样本量的方法。