如何对差异中的差异分析执行统计检验?

机器算法验证 假设检验 统计学意义 差异中的差异
2022-04-01 10:08:37

我正在尝试分析两个不同时间点的两个区域的结果。在 A 区我们进行了干预,而在 B 区我们没有。然后我们询问客户是否喜欢我们服务的特定部分,我们在干预前后都这样做了。两个时期的客户不同,所以我们没有配对。在接受调查的人中给我们最高分的人数是:

A区,时间=0:130个中有64个

A区,时间=1:118个中有82个

B区,时间=0:100个中有44个

B区,时间=1:100分中的60分

这又给出了以下比例:

A0 = 64/130 = 0.49

A1 = 82/118 = 0.69

B0 = 44/100 = 0.44

B1 = 60/100 = 0.6

根据差异中的差异分析,效果估计为:

E = [(A1-A0) - (B1-B0)] = 0.04

但是我如何进行统计测试以查看这种差异是否显着?

3个回答

差异的差异就是所谓的统计交互作用(正如 Dimitriy Masterov 已经指出的那样)。您想测试干预时与不干预时的时间效应是否不同。

您的数据最自然地建模为二项式,即,假设所有客户独立响应,每次在每个区域接受调查的总人数中得分最高的人数服从二项式分布。测试与二项式数据交互的标准统计方法是运行二项式逻辑回归。

在 R 中,代码如下。首先输入数据:

> NTopScore <- c(64,82,44,60)
> N <- c(130,118,110,100)
> Area <- factor(c("A","A","B","B"))
> Time <- factor(c(0,1,0,1))
> Proportion <- NTopScore / N

然后拟合逻辑回归。在 R 中,这是通过运行广义线性模型并告诉 R 数据应被视为二项式来完成的:

> fit <- glm(Proportion~Area*Time, family=binomial, weights=N)
> summary(fit)

Call:
glm(formula = Proportion ~ Area * Time, family = binomial, weights = N)

Deviance Residuals: 
[1]  0  0  0  0

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.03077    0.17543  -0.175  0.86076   
AreaB       -0.37469    0.26202  -1.430  0.15271   
Time1        0.85397    0.26599   3.211  0.00132 **
AreaB:Time1 -0.04304    0.38768  -0.111  0.91160   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2.3047e+01  on 3  degrees of freedom
Residual deviance: 7.1054e-15  on 0  degrees of freedom
AIC: 28.523

Number of Fisher Scoring iterations: 3

我们看到交互作用的 p 值(差异差异)是,显然不显着。P=0.9116

该模型在对数赔率(logit)尺度上进行拟合。AreaB 参数显示在时间 0 时区域 B 的比例低于区域 A。时间 1 参数显示时间 1 在区域 A 中的比例高于时间 0。AreaB:Time1 参数是差异的差异。

拟合逻辑回归的另一种方法是分别估计区域 A 和 B 的前后时间效应。这表明这两个区域的时间效应几乎相同,无论您是否有干预:

> fit <- glm(Proportion~Area+Area:Time, family=binomial, weights=N)
> summary(fit)

Call:
glm(formula = Proportion ~ Area + Area:Time, family = binomial, 
    weights = N)

Deviance Residuals: 
[1]  0  0  0  0

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.03077    0.17543  -0.175  0.86076   
AreaB       -0.37469    0.26202  -1.430  0.15271   
AreaA:Time1  0.85397    0.26599   3.211  0.00132 **
AreaB:Time1  0.81093    0.28204   2.875  0.00404 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance:  2.3047e+01  on 3  degrees of freedom
Residual deviance: -5.3291e-15  on 0  degrees of freedom
AIC: 28.523

Number of Fisher Scoring iterations: 3    

A区的时间效应为0.86397,B区的时间效应为0.81093。时间效应的差异是,这等于我们在第一个回归中看到的交互作用估计。0.810930.86397=0.04304

最简单的解决方案是使用 DID 的线性回归公式:

  1. 将二元客户评级回归到一个常数、一个后虚拟变量、一个区域 A 虚拟变量以及最后两个变量的交互作用上。添加其他回归量测量在个体水平上随时间变化的特征可能是合适的,但其分布在组水平上随时间而变化。如果它吸收了一些剩余方差,这可能有助于解决下面的显着性问题。
  2. DID 是职位和 A 组相互作用的系数。您可以进行假设检验,即它为零,或者只查看 p 值或 t-stat。
  3. 由于处理在区域内没有变化,通常的标准误差会偏离(通常太小,但有时太大,当集群内误差相关性为负时)。典型的解决方案是按面积或横截面对标准误差进行聚类,但只有 2-4 个聚类,这将无法正常工作,因为这不足以让渐近函数发挥作用。我真的没有很好在这里为您解决。Stata 确实在下面计算了一些东西,但它不太可能是可靠的(即使有少量的集群调整)。我确实怀疑正确的调整不会产生意义,因为 DID 系数的常规标准误差是如此之大。在这种情况下,人们经常会使用模拟来衡量 SE 的距离。

这是对您的数据的分析:

. clear

. input area_a time noobs rating

        area_a       time      noobs     rating
  1. 1 0 64 1 
  2. 1 0 66 0
  3. 1 1 82 1 
  4. 1 1 36 0
  5. 0 0 44 1
  6. 0 0 56 0
  7. 0 1 60 1
  8. 0 1 40 0
  9. end

. egen cs = group(area_a time)

. reg rating i.area_a##i.time [fw=noobs]

      Source |       SS           df       MS      Number of obs   =       448
-------------+----------------------------------   F(3, 444)       =      6.05
       Model |  4.34181458         3  1.44727153   Prob > F        =    0.0005
    Residual |  106.149257       444  .239074903   R-squared       =    0.0393
-------------+----------------------------------   Adj R-squared   =    0.0328
       Total |  110.491071       447  .247183605   Root MSE        =    .48895

------------------------------------------------------------------------------
      rating |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    1.area_a |   .0523077   .0650368     0.80   0.422    -.0755105    .1801259
      1.time |        .16   .0691484     2.31   0.021     .0241012    .2958988
             |
 area_a#time |
        1 1  |   .0426076   .0929871     0.46   0.647    -.1401419     .225357
             |
       _cons |        .44   .0488953     9.00   0.000     .3439051    .5360949
------------------------------------------------------------------------------

. reg rating i.area_a##i.time [fw=noobs], vce(cluster area) // or vce(cluster cs)

Linear regression                               Number of obs     =        448
                                                F(0, 1)           =          .
                                                Prob > F          =          .
                                                R-squared         =     0.0393
                                                Root MSE          =     .48895

                                 (Std. Err. adjusted for 2 clusters in area_a)
------------------------------------------------------------------------------
             |               Robust
      rating |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    1.area_a |   .0523077   1.08e-16  4.8e+14   0.000     .0523077    .0523077
      1.time |        .16   1.01e-16  1.6e+15   0.000          .16         .16
             |
 area_a#time |
        1 1  |   .0426076   1.30e-16  3.3e+14   0.000     .0426076    .0426076
             |
       _cons |        .44   1.01e-16  4.4e+15   0.000          .44         .44
------------------------------------------------------------------------------

两种规范中交互系数的解释都是干预后喜欢您的服务的 4.3 个百分点。

线性模型的优点是更容易解释对概率的加性效应,而不是对数赔率的乘性效应。此外,在非线性 logit 模型中,聚类是有问题的,因为系数只能按比例确定,交互作用的解释和识别假设可能会变得棘手,并且 DID 不再只是四种均值的交叉差异(参见下文引用的关于后一点的 Puhani 论文)。最后,在具有所有交互作用的完全饱和模型中,logit 和线性模型将给出交叉差异边际效应的相同点估计(尽管这不是您关心的参数):

. /* Cross difference to mimic OLS, but wrong */
. logit rating i.area_a##i.time [fw=noobs], nolog

Logistic regression                             Number of obs     =        448
                                                LR chi2(3)        =      17.87
                                                Prob > chi2       =     0.0005
Log likelihood = -298.57102                     Pseudo R2         =     0.0291

------------------------------------------------------------------------------
      rating |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    1.area_a |   .2103904   .2671347     0.79   0.431    -.3131839    .7339647
      1.time |   .6466272   .2867945     2.25   0.024     .0845203    1.208734
             |
 area_a#time |
        1 1  |   .2073448   .3911528     0.53   0.596    -.5593006    .9739902
             |
       _cons |  -.2411621   .2014557    -1.20   0.231    -.6360081    .1536839
------------------------------------------------------------------------------

. margins r.area_a#r.time 

Contrasts of adjusted predictions
Model VCE    : OIM

Expression   : Pr(rating), predict()

------------------------------------------------
             |         df        chi2     P>chi2
-------------+----------------------------------
 area_a#time |          1        0.21     0.6456
------------------------------------------------

--------------------------------------------------------------------
                   |            Delta-method
                   |   Contrast   Std. Err.     [95% Conf. Interval]
-------------------+------------------------------------------------
       area_a#time |
(1 vs 0) (1 vs 0)  |   .0426076   .0926461     -.1389755    .2241906
--------------------------------------------------------------------

我相信 4.6 个百分点的正确边际效应是这样给出的:

/* Puhani's DID Estimator */
gen at = area*time
logit rating i.area_a i.time i.at [fw=noobs], nolog
margins, at(area_a==1 time==1 at ==1) at(area_a==1 time==1 at==0) contrast(atcontrast(a._at) wald)

“非线性“Difference-in-Differences”模型中的治疗效果、交叉差异和交互项,Patrick A. Puhani,经济学快报,2012 年,115 (1),85-87。

感谢您的提问和回复。我有完全相同的问题,只是我有 8 个变量而不是 4 个变量,如下所示:

N_for_KPI <- c(683,538,2225,1458,294,307,922,781)
N <- c(1951,1564,5683,4507,819,862,2479,2511)
Wave <- factor(c("A","A","B","B","C","C"))
Brand <- factor(c(0,1,0,1,0,1))
data = data.frame(N_for_KPI,N)
Proportion <-N_for_KPI / N
Proportion

fit <- glm(Proportion~Wave*Brand, family=binomial, weights=N)
summary(fit)

虽然我得到了显着性结果,但我得到了一个错误:

> N_for_KPI <- c(683,538,2225,1458,294,307,922,781)
> N <- c(1951,1564,5683,4507,819,862,2479,2511)
> Wave <- factor(c("A","A","B","B","C","C"))
> Brand <- factor(c(0,1,0,1,0,1))
> Proportion <-N_for_KPI / N
> Proportion
[1] 0.3500769 0.3439898 0.3915186 0.3234968 0.3589744 0.3561485 
0.3719242 0.3110315
> fit <- glm(Proportion~Wave*Brand, family=binomial, weights=N)
Error in model.frame.default(formula = Proportion ~ Wave * Brand, 
weights = N,  : 
variable lengths differ (found for 'Wave')
> summary(fit)
Call:
glm(formula = Proportion ~ Wave * Brand, family = binomial, weights 
= N)

Deviance Residuals: 
[1]  0  0  0  0

Coefficients:
         Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -2.9422     0.1047 -28.096   <2e-16 ***
WaveB          0.0394     0.1203   0.328    0.743    
Brand1        -0.1574     0.1507  -1.045    0.296    
WaveB:Brand1  -0.4487     0.1786  -2.512    0.012 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance:  4.5383e+01  on 3  degrees of freedom
Residual deviance: -4.9938e-13  on 0  degrees of freedom
AIC: 35.137

Number of Fisher Scoring iterations: 3

除了错误,当我更改为 N_for_KPI 和 N 的新样本大小时,显着结果没有改变

您能否帮助建议如何在该模型中拟合 8 个变量?非常感谢你!