偏回归

Partial regression

发布于

2024年11月7日

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

首先安装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)做dm1em1am1的回归。Model A(Augmented model)做dm1pm1em1am1的回归。

# 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检验被视为对R2 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 多元回归中对pm1t检验是等价的,它们都检验了控制em1与am1的作用后pm1的独特效应。此时,F分子上的df为1,F = t2F分母的df即为tdf

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。此时,PREF检验与分层回归的F检验是等价的。

4 使用残差进行回归

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

dm1em1am1的回归,得到dm1的残差dm_res,该残差排除了em1am1dm1的效应。

再做pm1em1am1的回归,得到pm1的残差pm_res,该残差排除了em1am1pm1的效应。

之后做dm_respm_res的相关,我们就得到了dm1pm1的偏相关系数/净相关系数。

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

可见,dm1pm1的偏相关系数为-0.3329009。

dm_respm_res的一元回归,我们就得到了pm1dm1的独特效应。

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的偏回归系数相同。然而,它们的tp是不同的。为什么呢?前文用PRE检验了pm1的独特效应,这里我们使用PRE检验pm_res的独特效应,然后将它们进行比较。注意,单个回归系数的PREF检验与该回归系数的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_respm1PRE。可见,PRE是相等的,即pm1的独特效应本质上没有变化。但是,PREF检验的df2是不同的,这使得Fp值也发生了变化。换句话说,pm1的独特效应PRE是不变的,PRE基于模型比较检验PRE的显著性,但用来检验PRE的显著性的Model C与Model A是不同的,因而模型比较的结论(i.e., F检验和t检验的结果)也不同。打个比方,一勺盐放在一杯水里感觉齁咸,但一勺盐放在一锅水里感觉咸淡刚好。

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

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