机器算法验证 - 使用 R 中的回归消除因子对连续比例数据的影响 - 吾爱随笔录

使用 R 中的回归消除因子对连续比例数据的影响

机器算法验证 r 回归物流贝塔分布贝塔回归

2022-03-29 05:29:32

我有一个取决于固定效应因子的连续比例数据集，例如：

df <- data.frame(
   x = c(0.1,0.2,0.12,0.3,0.2,0.35,0.5,0.3,0.45,0.6,0.54,0.77),
   f = as.factor(c(rep("A", 6), rep("B", 6)))
)

我正在寻找一种消除因子影响的方法，例如，如果残差是正态分布的，您可以通过以下方式实现：

fit <- lm(x ~ f, df)
fit$residuals + fit$coef[1]

连续比例数据的合适方法是什么？我希望输出（在去除因子的估计影响后）处于相同的比例（0-1）。

1个回答

我认为，最大的问题是你需要一个数据模型来处理。这些可能是回归模型，也可能只是您的数据适合分布。从理论上讲，您的数据可能与许多模型/分布一致（例如正态分布、由具有适当总试验次数的二项分布产生的成功比例或 beta 分布）。既然你说这些是假数据只是为了说明问题，并且你有连续的比例，我会简单地认为这些比例来自beta 分布。

在其中，您可以使用betareg 包R对 beta 分布式数据进行建模。您应该阅读小插图 ( pdf )，这对于理解有关 beta 回归的想法以及如何使用该包来执行此操作非常有帮助。请特别注意，它们对 beta 分布的参数化与典型的参数化不同。而不是使用 $\alpha$ 和 $\beta$ （例如，标准?pbeta等函数中的shape1and参数），它们使用平均比例，shape2 $\mu$ , 和精度, $\phi$ ，其定义如下：

\begin{aligned} μ & = \frac{α}{(α + β)} & and, & ϕ = α + β \\ thus, \\ α & = μ \times ϕ & and, & β = ϕ - α \end{aligned}

$\begin{align} \mu &= \frac{\alpha}{(\alpha + \beta)} &\text{and,}& &\phi = \alpha + \beta \\ \ &\ &\ & &\ \\ & &\text{thus,}& & \\ \ &\ &\ & &\ \\ \alpha &= \mu\times\phi &\text{and,}& &\beta = \phi - \alpha \end{align}$

我们可以使用这些方程在两个参数化之间移动，从而在两组函数之间移动。然后我们可以将建模链接betareg到分布和分位数函数。操作转换将基于类比qq-plots如何用于比较两个分布。

首先，让我们对这些数据进行建模betareg并确定组之间的精度是否恒定：

library(betareg)
library(lmtest)
  # test of constant precision based on beta distribution:
cp = betareg(x~f, df)
vp = betareg(x~f|f, df)
lrtest(cp, vp)
# Likelihood ratio test
# 
# Model 1: x ~ f
# Model 2: x ~ f | f
#   #Df LogLik Df  Chisq Pr(>Chisq)
# 1   3 9.2136                     
# 2   4 9.4793  1 0.5314      0.466

对于这些（假）数据，几乎没有理由相信精度会有所不同，但对于您的真实数据，它们可能会有所不同。所以我将演示更复杂的版本，如果你愿意，你可以简化它。无论如何，下一步是确定估计的 $\alpha$ ( shape1) 和 $\beta$ ( shape2) 参数，A并B通过将模型的系数转换为分布参数（不要忘记链接函数！）并使用上述公式在两个参数化之间进行转换：

summary(vp)
# ... 
# Coefficients (mean model with logit link):
#             Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  -1.3153     0.2165  -6.074 1.25e-09 ***
# fB            1.4260     0.3178   4.487 7.23e-06 ***
#   
# Phi coefficients (precision model with log link):
#             Estimate Std. Error z value Pr(>|z|)    
# (Intercept)   3.0094     0.5695   5.284 1.26e-07 ***
# fB           -0.5848     0.7943  -0.736    0.462    
# ... 
lo.to.prop = function(x){
  o = exp(x)
  p = o / (o+1)
  return(p)
}
alpha.A = lo.to.prop(coef(vp)[1]            ) * exp(coef(vp)[3]            )
alpha.B = lo.to.prop(coef(vp)[2]+coef(vp)[1]) * exp(coef(vp)[4]+coef(vp)[3])
beta.A  = exp(       coef(vp)[3]            ) - alpha.A
beta.B  = exp(       coef(vp)[4]+coef(vp)[3]) - alpha.B

由于这很复杂，我们可能会对最终值进行常识检查。（如果您不需要对它们进行建模和/或不需要测试精度是否相同，这也建议使用另一种更简单的方法来获取这些参数。）

cbind(c(alpha.A, beta.A), c(alpha.B, beta.B))
#                        [,1]     [,2]
# (Intercept)        4.290073 5.960932
# (phi)_(Intercept) 15.984533 5.336616

  # parameterization for beta distribution:
library(fitdistrplus)
fitdist(with(df, x[f=="A"]), "beta")
# Fitting of the distribution ' beta ' by maximum likelihood 
# Parameters:
#         estimate Std. Error
# shape1  4.290883   2.389871
# shape2 15.987505   9.295580
fitdist(with(df, x[f=="B"]), "beta")
# Fitting of the distribution ' beta ' by maximum likelihood 
# Parameters:
#        estimate Std. Error
# shape1 5.960404   3.379447
# shape2 5.335963   3.010541

现在，我们将原始比例转换B为它们相应的 beta 概率（即，它们在自己的 beta 分布中的位置），并将这些概率转换为A的 beta 分布的分位数：

p.B      = pbeta(df$x[df$f=="B"], shape1=alpha.B, shape2=beta.B)
q.B.as.A = qbeta(p.B,             shape1=alpha.A, shape2=beta.A)

最后一步，让我们看一下新值：

cbind(df$x[df$f=="B"], q.B.as.A)
#             q.B.as.A
# [1,] 0.50 0.18502881
# [2,] 0.30 0.08784683
# [3,] 0.45 0.15784496
# [4,] 0.60 0.24648412
# [5,] 0.54 0.20838604
# [6,] 0.77 0.38246863

在此处输入图像描述

其它你可能感兴趣的问题

上一篇解释 log(y+1) 转换响应的回归系数下一篇为什么我们不应该沉迷于不偏不倚