在R中运行百分比最小二乘回归答案

【问题标题】：Running percentage least squares regression in R在R中运行百分比最小二乘回归
【发布时间】：2013-09-18 15:35:30
【问题描述】：

我有兴趣在 R 中运行百分比最小二乘回归，而不是普通的最小二乘回归。这也可以称为具有乘法误差的线性模型。之前有人问过这个网站上关于最小二乘百分比的问题，响应者建议研究加权回归，一种可能性是通过 X 值的平方反比对每个观察值进行加权。

stackoverflow.com/questions/15275236/least-square-percentage-regression

但是，这假设我知道每个观察应该先验地加权多少。我不。我不知道百分比误差是 1%、10%、15% 等。我想要的是一个适合的模型

y= b1*x + e

误差项被建模为：

e= b2*x

b2 是回归模型中需要最小化的百分比误差。我还没有找到任何包或任何代码来适合 R 的这种类型的模型。任何关于如何做到这一点的反馈将不胜感激。

【问题讨论】：

这不等同于使用 y 的对数变换并进行普通最小二乘回归吗？要返回未转换的 y，您需要对 RHS 取幂，从而导致乘法项和错误。
@zkurtz 对数变换假定关系采用 y=e^x 的形式，因此如果关系是真正的线性关系，则反向变换的影响大小将减小。对数转换可能适用于确定显着性，但不适用于影响大小。我更愿意按照实际情况对数据进行建模，而不是以估计适当的效应大小为代价进行转换以处理非正态残差。
更明确地说，日志转换将使错误分布正常并处理此模式的异方差性，以便可以在 OLS 框架中处理。此外，在这种情况下，从未转换的模型中获取效应大小可能是合理的。但是，我更愿意运行一个模型来处理这一切。这将通过运行百分比最小二乘回归而不是普通的最小二乘回归来实现。

标签： r regression

【解决方案1】：

我假设您的意思是 Tofallis (2009) 定义的百分比回归。

用他的例子：

Sales <- c(6375,11626,14655,21869,26408,32406,35108,40295,70762,80553,95294,101314,116141,122316,141650,175026,230614,293543)
Expenses <- c(62.5,92.9,178.3,258.4,494.7,1083,1620.6,421.7,509.2,6620.1,3918.6,1595.3,6107.5,4454.1,3163.8,13210.7,1703.8,9528.2)

如果我们应用普通最小二乘法并将销售额作为依赖项变量我们得到模型 Sales = 43942 + 15.00 R&D 截距和斜率的 p 值分别为 0.03 和 0.0015。

fit1 <- lm(Sales ~ Expenses)
summary(fit1)
#                Estimate Std. Error t value Pr(>|t|)   
# (Intercept)   43941.705  18493.079   2.376  0.03033 * 
#   Expenses       14.994      3.915   3.830  0.00148 **

如果我们这样做并执行普通最小二乘，我们会得到模型：Ln(Sales) = 10.341 + 0.000198 R&D，p 值为 0.002 斜率和截距基本上为零。

fit2 <- lm(log(Sales) ~ Expenses)  
summary(fit2)
#                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   1.034e+01  2.535e-01  40.793  < 2e-16 ***
#   Expenses    1.982e-04  5.366e-05   3.694  0.00197 **

最后，我们转向本文提出的方法，最小化平方百分比残差。结果模型被发现是，转换回来后：销售额 = 8817 + 17.88 R&D，p 值为斜率和截距分别为 0.002 和 5×10-5。

fit3 <- lm(Sales ~ Expenses, weights = 1/Sales^2)
summary(fit3)
#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   8816.553   2421.644   3.641   0.0022 ** 
#   Expenses      17.880      3.236   5.525 4.61e-05 ***

所以说到底，这是加权回归。

为了确认这一点，我们还可以使用数值优化：

resfun <- function(par) {
  sum((Sales - par[[1]]*Expenses - par[[2]])^2 / Sales^2)
}

optim(c(10,1000), resfun)
# $par
# [1]   17.87838 8816.44304

optim(c(10,1000), resfun, method="BFGS")
# $par
# [1]   17.97975 8575.71156

（不同的优化器会给出稍微不同的结果。）

【讨论】：

我的意思是 Tofallis 2009 定义的百分比回归。感谢您附上参考资料。你的回复和原论文都很清楚。不过，我确实有一个问题-您链接的 Tofallis 论文建议按 (1/y) 加权，而您建议按 (1/y^2) 加权。 y 在您的示例中是“销售”。为什么您选择按 (1/y^2) 而不是源中推荐的 (1/y) 加权？对我来说，通过平方反比加权而不是逆平方加权会给更大的观察结果甚至更小的权重，这不是我想做的事情！
@colin Weighting with 1/y^2（使用函数lm）意味着平方残差乘以权重（即残差乘以1/y）。一眼看去，这与参考文献一致。如果您不同意，请指出论文的相关部分。如果可以改进模型，您可以随意使用任何权重。
谢谢。我决定用 1/y 加权，而不是 1/y^2。在手稿中，描述位于“系数公式的推导”部分，从第 3 页底部开始，一直到第 4 页。接下来，在检查正态性时，我相信在这种情况下查看是正确的模型残差除以 y 值的平方根的 qqnorm 图，残差与拟合值的关系也应如此？
@colin AFAIU 他们将残差除以y。由于我的代码将残差平方除以y^2，它应该是相同的。但是，此处适用于加权回归的一般建议。研究残差图和杠杆图并比较 AIC。
@colin：我认为 weights =(1/y^2) 在这个例子中起作用。拟合值（残差 = 拟合的 40%）：10 (4)，r^2[1/y] 16/10 = 1.6，r^2[1/y^2] 16/100 = 0.16； 20 (8)，r^2[1/y] 64/20 = 3.2，r^2[1/y^2] 64/400 = 0.16。如果您的模型最小化 Residuals^2，并且您希望 y 的 z% 的所有 Residuals 具有相同的权重（在我的示例中为 40%），那么 weights=(1/y^2)。

【解决方案2】：

查看nlme 包中的gls 函数，以及varClasses 之一，例如varIdent 或varPower。

可能是这样的模型：

gls( y ~ x, data=mydata, weights=varPower(form= ~x) )

【讨论】：