r - R中的约束多元线性回归

Question

假设我必须估计回归中的系数 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]范围内。

我也尝试应用 glm 的方法与 offset 一起使用选项: Regression for a Rate variable in R
和
http://en.wikipedia.org/wiki/Poisson_regression#.22Exposure.22_and_offset
但这并不成功。

请注意，在我的数据因变量中:比例 是/x 既是零膨胀又是一膨胀。
任何想法，在 R ('glm'，'lm')中构建模型的合适方法是什么？

Answer 1

你在正确的轨道上:如果 0 ≤ y ≤ x 则 0 ≤ (y/x) ≤ 1。这建议拟合 y/x到 glm(...) 中的逻辑模型.详细信息如下，但考虑到你只有 6 分，这很合适。

主要问题是模型无效，除非(y/x) 中的错误。是具有恒定方差的正态(或者，等效地，y 中的误差随着 x 增加)。如果这是真的，那么我们应该得到一个(或多或少)线性 QQ 图，我们这样做了。

一个细微差别:glm逻辑模型的接口(interface)需要两列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")