使用用户自定义的family进行glm回归的模型公式

Question

背景

我正在尝试预测产品线的销售情况（在最后的样本中表示为y_test）。其时间段内的销售量基于另一个产品的所有先前销售情况（x_test），以及这些先前销售情况中有多少仍在使用。然而，直接测量先前销售的产品中有多少仍在使用是不可能的，因此需要推断出生存曲线。

例如，如果您制造特定智能手机型号的配件，那么配件销售量至少部分基于仍在使用的智能手机数量。（顺便说一下，这不是作业。）

细节

我有一些时间序列数据，并希望使用glm或类似方法来拟合回归模型。依赖变量和独立变量之间的关系如下：

其中p是时间期，y_p是依赖变量，x_p是独立变量，c₀和c₁是回归系数，F_t是累积分布函数（例如pgamma），e_p是残差项。

在前三个时间段内，该函数会扩展为以下形式：

#y[1] = c0 + c1*(x[1]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 0, 1)$value))
#y[2] = c0 + c1*(x[1]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 1, 2)$value) + x[2]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 0, 1)$value))
#y[3] = c0 + c1*(x[1]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 2, 3)$value) + x[2]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 1, 2)$value) + x[3]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 0, 1)$value))

因此，我有x_p和y_p的历史数据，并且想要获取最小化残差的系数/参数c₀、c₁、c₂和c₃的值。

我认为解决方案是使用glm并创建自定义family，但我不确定如何做到这一点。 我查看了Gamma family的代码，但没有取得太大进展。我已经能够使用nlminb手动进行优化，但我更喜欢glm或类似函数提供的简单性和有用性（即predict等）。

这里是一些示例数据:

# Survival function (the integral part):
fsurv<-function(q, par) {
  l<-length(q)
  out<-vapply(1:l, function(i) {1-integrate(function(x) {pgamma(x, par[1], par[1]/par[2])}, q[i]-1, q[i])$value}, FUN.VALUE=0)
  return(out)}

# Sum up the products:
frevsumprod <- function(x,y) {
  l <- length(y)
  out <- vapply(1:l, function(i) sum(x[1:i]*rev(y[1:i])), FUN.VALUE=0)
  return(out)}

# Sample data:
p<-1:24 # Number of periods
x_test<-c(1188, 2742, 4132) # Sample data
y_test<-c(82520, 308910, 749395, 801905, 852310, 713935, 624170, 603960, 640660, 553600, 497775, 444140) # Sample data
c<-c(-50.161147,128.787437,0.817085,13.845487) # Coefficients and parameters, from another method that fit the data

# Pad the data to the correct length:
pad<-function(p,v,padval=0) {
  l<-length(p)
  padv<-l-length(v)
  if(padv>0) (v<-c(v,rep(padval,padv)))
  return(v)
}
x_test<-pad(p,x_test)
y_test<-pad(p,y_test,NA)

y_fitted<-c[0+1]+c[1+1]*frevsumprod(x_test,fsurv(p,c[(2:3)+1])) # Fitted values from regression

library(ggplot2)
ggplot(data.frame(p,y_test,y_fitted))+geom_point(aes(p,y_test))+geom_line(aes(p,y_fitted)) # Plot actual and fit

- dnlbrky

---

谢谢@DWin。我编辑了问题，但仍然没有任何尝试来回答它。您有其他更好的提问建议吗？另外，我认为几乎所有的问题陈述都与任何潜在的答案相关，而不仅仅是对我的“失败努力”的叙述。 - dnlbrky

我建议你将此标记为请求迁移到CrossValidated。（或者你可以在那里交叉发布并附上一条注释说明它在SO上没有得到回应。） - IRTFM

Answer 1

使用glm无法完成此操作。在glm中，family指定线性预测器与y的平均值的连接方式。请参见?family和wiki。特别地，您需要能够编写一个family列表，其中包含（某些）函数，例如：

> fam <- poisson()
> str(fam)
List of 12
 $ family    : chr "poisson"
 $ link      : chr "log"
 $ linkfun   :function (mu)  
 $ linkinv   :function (eta)  
 $ variance  :function (mu)  
 $ dev.resids:function (y, mu, wt)  
 $ aic       :function (y, n, mu, wt, dev)  
 $ mu.eta    :function (eta)  
 $ initialize:  expression({  if (any(y < 0))  stop("negative values not allowed for the 'Poisson' family")  n <- rep.int(1, nobs| __truncated__
 $ validmu   :function (mu)  
 $ valideta  :function (eta)  
 $ simulate  :function (object, nsim)  
 - attr(*, "class")= chr "family"
> 
> fam <- Gamma()
> str(fam)
List of 12
 $ family    : chr "Gamma"
 $ link      : chr "inverse"
 $ linkfun   :function (mu)  
 $ linkinv   :function (eta)  
 $ variance  :function (mu)  
 $ dev.resids:function (y, mu, wt)  
 $ aic       :function (y, n, mu, wt, dev)  
 $ mu.eta    :function (eta)  
 $ initialize:  expression({  if (any(y <= 0))  stop("non-positive values not allowed for the 'gamma' family")  n <- rep.int(1, n| __truncated__
 $ validmu   :function (mu)  
 $ valideta  :function (eta)  
 $ simulate  :function (object, nsim)  
 - attr(*, "class")= chr "family"

其中eta是指线性预测器。也就是说，至少需要指定一个反函数链接linkinv，它仅仅通过参数和协变量之间的点积来依赖于协变量。而你的不是这样的，因为它以非线性方式依赖于c_2和c_3。

- Benjamin Christoffersen