r - 具有用户指定族的 glm 回归模型公式

Question

背景

我正在尝试预测产品线的销售额(最后样本中的 y_test)。它在一段时间内的销售额基于其他产品 (x_test) 的所有先前销售额以及这些先前销售额中有多少仍在使用。但是，无法直接衡量那些以前销售的仍在使用的产品的数量，因此需要推断出生存曲线。

例如，如果您为特定的智能手机型号生产配件，配件销售至少部分取决于仍在使用的智能手机的数量。 (这不是作业，顺便说一句。)

详情

我有一些时间序列数据，想使用 glm 拟合回归模型或类似的东西。因变量和自变量之间的关系是这样的:


其中 p 是时间段，yp 是因变量，xp 是自变量，c0 和 c1 是回归系数，Ft 是累积分布函数(例如 pgamma)，ep 是残差。

通过前三个时间段，该函数将扩展为如下所示:

#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))

因此，我有 xp 和 yp 的历史数据，并且我想获得最小化残差的系数/参数 c0、c1、c2 和 c3 的值。

我认为解决方案是使用 glm并创建一个自定义系列，但我不知道该怎么做。 我查看了 Gamma 的代码家人，但没有走多远。我已经能够使用 nlminb“手动”进行优化，但我更喜欢 predict 提供的简单性和实用性(即 glm 和其他)或类似的功能。

以下是一些示例数据:

# 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

Answer 1

这不能用 glm 来完成. family在 glm指定线性预测变量如何与 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指的是线性预测器。 IE。至少你需要指定一个反向链接函数，linkinv , 其中 只有通过参数和协变量之间的点积取决于协变量。你的不是，因为它以非线性方式依赖于 c_2 和 c_3。