在R中进行受限制的多元线性回归

Question

在R中进行受限制的多元线性回归

4

假设我需要在回归中估计系数a和b：

y=a*x+b*z+c

我提前知道y的范围总是y>=0和y<=x，但是回归模型有时会产生超出此范围的y值。

样本数据：

mydata<-data.frame(y=c(0,1,3,4,9,11),x=c(1,3,4,7,10,11),z=c(1,1,1,9,6,7))
round(predict(lm(y~x+z,data=mydata)),2) 
    1     2     3     4     5     6 
-0.87  1.79  3.12  4.30  9.34 10.32

第一个预测值为<0。

我尝试了没有截距的模型：所有的预测值都>0，但是y的第三个预测值>x（4.03>3）

round(predict(lm(y~x+z-1,data=mydata)),2)
   1    2    3    4    5    6 
0.76 2.94 4.03 4.67 8.92 9.68

我也考虑过对比例y/x进行建模，而不是只对y进行建模：

mydata$y2x<-mydata$y/mydata$x
round(predict(lm(y2x~x+z,data=mydata)),2)
   1    2    3    4    5    6 
0.15 0.39 0.50 0.49 0.97 1.04 
round(predict(lm(y2x~x+z-1,data=mydata)),2)
   1    2    3    4    5    6 
0.08 0.33 0.46 0.47 0.99 1.07

但是现在第六个预测值>1，但比例应该在[0,1]范围内。

我还尝试了使用带有offset选项的glm方法： R中用于速率变量的回归和http://en.wikipedia.org/wiki/Poisson_regression#.22Exposure.22_and_offset，但这并不成功。

请注意，在我的数据中，依赖变量proportion y/x既具有零膨胀又具有一膨胀。

有什么适合在R中构建模型的方法吗（'glm'、'lm'）？

- ipj

你尝试了什么，为什么不起作用，我们是否也能得到作业积分？ - Andy Clifton

1

如果您提供数据或代表子集，并像@AndyClifton所说的那样展示您尝试过的内容，那么您更有可能获得帮助。此外，在您的模型中，y出现在LHS和RHS上。这是有意为之吗？ - jlhoward

我会检查产生系数超出此范围的数据集。为什么坚信它们总是在你的界限内？如果数据是模拟的，那么有几个越过边缘也没有问题。 - Roman Luštrik

我的错误已经修正：y仅出现在左侧。数据并非模拟的，而且了解生成数据的过程规则0<=y<=x不能被违反。我还提供了详细的描述和代码示例。 - ipj

现在代码应该可以正常运行。 - ipj

显示剩余3条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jlhoward · Accepted Answer

您正在走上正确的道路：如果0 ≤ y ≤ x，则0 ≤ (y/x) ≤ 1。这表明需要在glm(...)中将y/x拟合到logistic模型中。具体细节如下，但考虑到您只有6个点，这是一个相当好的拟合。

主要问题是，除非(y/x)中的误差为常数方差的正态分布（或等价地，y中的误差随x增加而增加），否则该模型无效。如果这是真的，则我们应该得到一个（或多或少）线性Q-Q图，我们确实得到了这个图。

一个微妙之处：glm logistic模型的接口需要两列y：“成功次数（S）”和“失败次数（F）”。然后它计算概率为S/(S+F)。因此，我们必须提供两列来模仿这个模型：y和x-y。然后glm(...)将计算y/(y+(x-y)) = y/x。

最后，拟合摘要表明x很重要，z可能重要也可能不重要。您可能希望尝试一种排除z并查看是否改善AIC的模型。

fit = glm(cbind(y,x-y)~x+z, data=mydata, family=binomial(logit))
summary(fit)
# Call:
# glm(formula = cbind(y, x - y) ~ x + z, family = binomial(logit), 
#     data = mydata)

# Deviance Residuals: 
#        1         2         3         4         5         6  
# -0.59942  -0.35394   0.62705   0.08405  -0.75590   0.81160  

# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)  
# (Intercept)  -2.0264     1.2177  -1.664   0.0961 .
# x             0.6786     0.2695   2.518   0.0118 *
# z            -0.2778     0.1933  -1.437   0.1507  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# (Dispersion parameter for binomial family taken to be 1)

#     Null deviance: 13.7587  on 5  degrees of freedom
# Residual deviance:  2.1149  on 3  degrees of freedom
# AIC: 15.809

par(mfrow=c(2,2))
plot(fit)         # residuals, Q-Q, Scale-Location, and Leverage Plots

mydata$pred <- predict(fit, type="response")
par(mfrow=c(1,1))
plot(mydata$y/mydata$x,mydata$pred,xlim=c(0,1),ylim=c(0,1), xlab="Actual", ylab="Predicted")
abline(0,1, lty=2, col="blue")