R 中的分数响应回归

Fractional Response Regression in R

我正在尝试对响应变量介于 0 和 1 之间的数据建模,因此我决定在 R 中使用分数响应模型。根据我目前的理解,分数响应模型类似于逻辑回归,但它使用 qausi-likelihood 方法来确定参数。我不确定我是否理解正确。

到目前为止,我尝试的是 frm 包中的 frmglm 中的以下数据,与此相同 OP

library(foreign)
mydata <- read.dta("k401.dta")

此外,我遵循了这个 OP 中使用 glm 的过程。但是,对于具有 frm 的相同数据集,它 returns 不同的 SE

library(frm)
y <- mydata$prate
x <- mydata[,c('mrate', 'age', 'sole', 'totemp1')]
myfrm <- frm(y, x, linkfrac = 'logit')

frm returns,

*** Fractional logit regression model ***

           Estimate Std. Error t value Pr(>|t|)    
INTERCEPT  1.074062   0.048902  21.963    0.000 ***
mrate      0.573443   0.079917   7.175    0.000 ***
age        0.030895   0.002788  11.082    0.000 ***
sole       0.363596   0.047595   7.639    0.000 ***
totemp1   -0.057799   0.011466  -5.041    0.000 ***

Note: robust standard errors

Number of observations: 4734 
R-squared: 0.124 

glm,我用

myglm <- glm(prate ~ mrate + totemp1 + age + sole, data = mydata, family = quasibinomial('logit'))
summary(myglm)

Call:
glm(formula = prate ~ mrate + totemp1 + age + sole, family = quasibinomial("logit"), 
    data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.1214  -0.1979   0.2059   0.4486   0.9146  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.074062   0.047875  22.435  < 2e-16 ***
mrate        0.573443   0.048642  11.789  < 2e-16 ***
totemp1     -0.057799   0.011912  -4.852 1.26e-06 ***
age          0.030895   0.003148   9.814  < 2e-16 ***
sole         0.363596   0.051233   7.097 1.46e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.2913876)

    Null deviance: 1166.6  on 4733  degrees of freedom
Residual deviance: 1023.7  on 4729  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

我应该依靠哪一个?使用 glm 而不是 frm 是否更好,因为我已经看到 SE 估计的 OP 可能不同

您需要将 prate 变量除以 100。您可能还需要升级 frm 的版本。

这两种方法的差异源于稳健标准误差计算中不同的自由度校正。使用类似的默认值,结果将是相同的。请参阅以下示例:

library(foreign)
library(frm)
library(sandwich)
library(lmtest)

df <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/401k.dta")
df$prate <- df$prate/100

y <- df$prate
x <- df[,c('mrate', 'age', 'sole', 'totemp')]

myfrm <- frm(y, x, linkfrac = 'logit')

*** Fractional logit regression model ***

           Estimate Std. Error t value Pr(>|t|)    
INTERCEPT  0.931699   0.084077  11.081    0.000 ***
mrate      0.952872   0.137079   6.951    0.000 ***
age        0.027934   0.004879   5.726    0.000 ***
sole       0.340332   0.080658   4.219    0.000 ***
totemp    -0.000008   0.000003  -2.701    0.007 ***

现在 GLM:

myglm <- glm(prate ~ mrate + totemp + age + sole, 
             data = df, family = quasibinomial('logit'))
coeftest(myglm, vcov.=vcovHC(myglm, type="HC0"))

z test of coefficients:

                 Estimate    Std. Error z value              Pr(>|z|)    
(Intercept)  0.9316994257  0.0840772572 11.0815 < 0.00000000000000022 ***
mrate        0.9528723652  0.1370808798  6.9512     0.000000000003623 ***
totemp      -0.0000082352  0.0000030489 -2.7011              0.006912 ** 
age          0.0279338963  0.0048785491  5.7259     0.000000010291017 ***
sole         0.3403324262  0.0806576852  4.2195     0.000024488075931 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

对于HC0,标准误差是相同的。也就是说,frm默认使用HC0。有关广泛的讨论,请参阅 sandwich 使用的默认值在某些情况下可能更好,但我怀疑一般情况下它并不重要。您已经可以从结果中看出这一点:差异在数值上非常小。