指数回归方法之间的差异

机器算法验证 回归 广义线性模型 非线性回归 指数分布
2022-03-18 06:15:31

可以用许多不同的方式拟合指数。这篇文章建议在lm响应变量的日志上做一些糟糕的事情。这个 SO post建议使用nlswhich 需要一个起始估计。这篇 SO 帖子建议glm使用 gamma/log 链接功能。在这里,杰出的@Glen-b 解释了方法之间的一些潜在差异。

这些不同方法的优点/缺点和适用范围是什么?这些方法在计算置信区间的程度或方式上是否有所不同?

像现在在家中的所有其他数据科学家一样,我正在处理 Covid 19 数据。

我特别注意到的一件事是我可以lmlog,等来做log10log2但必须用 . 从自然对数转换glm

last_14 = data.frame(rbind(
c(3460,  14,    0),
c(3558,  17,    1),
c(3802,  21,    2),
c(3988,  22,    3),
c(4262,  28,    4),
c(4615,  36,    5),
c(4720,  40,    6),
c(5404,  47,    7),
c(5819,  54,    8),
c(6440,  63,    9),
c(7126,  85,   10),
c(7905, 108,   11),
c(8733, 118,   12),
c(9867, 200,   13)))
names(last_14) = c('World', 'US', 'days')

lm(log(World) ~ days, last_14)
#> 
#> Call:
#> lm(formula = log(World) ~ days, data = last_14)
#> 
#> Coefficients:
#> (Intercept)         days  
#>     8.06128      0.08142

glm(formula = World ~ days,  data=last_14, family=gaussian(link='log'))
#> 
#> Call:  glm(formula = World ~ days, family = gaussian(link = "log"), 
#>     data = last_14)
#> 
#> Coefficients:
#> (Intercept)         days  
#>     8.00911      0.08819  
#> 
#> Degrees of Freedom: 13 Total (i.e. Null);  12 Residual
#> Null Deviance:       54450000 
#> Residual Deviance: 816200    AIC: 199.4

nls(World ~ exp(a + b*days), last_14, start=list(a=5, b=0.03))
#> Nonlinear regression model
#>   model: World ~ exp(a + b * days)
#>    data: last_14
#>       a       b 
#> 8.00911 0.08819 
#>  residual sum-of-squares: 816246
#> 
#> Number of iterations to convergence: 8 
#> Achieved convergence tolerance: 1.25e-06

reprex 包(v0.3.0)于 2020-03-20 创建

2个回答

差异之一是每个模型的可能性。万一读者不记得了,可能性包含了关于数据条件分布的假设。对于 COVID-19,这将是特定日期的感染分布(或报告的新病例或死亡病例等)。无论我们希望结果如何,我们都称它为因此,条件分布(例如今天的新病例数)将是(认为这是为条件)。yy|tyt

  • 在获取日志然后执行的情况下lm,这意味着等效地,的对数正态我们对进行线性回归的原因是因为在对数尺度上,条件均值与方差无关,因为对数正态的均值也是方差的函数。所以Pro:我们知道如何进行线性回归,但是Conlog(y)|tN(μ(x),σ2)ytlog(y)这种方法在对数尺度上做出线性回归假设,这些假设总是可以评估的,但在理论上可能很难证明是正确的?的话,人们没有意识到在对数尺度上进行预测然后取指数实际上会使预测产生偏差因此,当您从对数正态模型进行预测时,您需要考虑到这一点。 exp(σ2/2)

  • 据我了解,nls也假设高斯似然,所以在这个模型中除了现在,我们让结果的条件均值是非线性的。这可能很痛苦,因为没有置信区间不低于 0,因此您的模型可能会估计负数的感染数。显然,这不可能发生。当感染数(或其他)较大时,高斯可以证明是合理的。但是当事情刚刚开始时,这可能不是最好的可能性。此外,如果您使用 拟合数据,您会发现它非常适合后期数据,但不适合早期数据。这是因为后期数据的错误拟合会导致更大的损失,而目标是最小化这种损失。y|tN(exp(β0+βt),σ2)nlsnls

  • frees的方法glm有点小,它允许我们通过链接函数来控制条件分布以及条件均值的形式。在这个模型中,我们称为链接,对于日志链接Pro这些模型更具表现力,但我认为力量来自以不正常的可能性进行推理的能力。这解除了许多限制,例如对称置信区间。缺点是你需要更多的理论来理解正在发生的事情y|tGamma(μ(x),ϕ)μ(x)=g1(β0+β1)gμ(x)=exp(β0+β1t)

使用非线性拟合或线性化拟合拟合指数曲线之间的已知差异是不同点的误差/残差相关性的差异。

您可以在下面的图中注意到这一点。

例子

在那个情节中你可以看到

  • 线性化拟合(虚线)更精确地拟合具有小值的点(参见右侧的图,其中虚线更接近开始时的值)。
  • 非线性拟合更接近具有高值的点。

    modnls <- nls(US ~ a*exp(b*days), start=list(a=100, b=0.3))
    modlm <- lm(log(US) ~ days )
    plot(days,US, ylim = c(1,15000))
    lines(days,predict(modnls))
    lines(days,exp(predict(modlm)), lty=2)
    title("linear scale", cex.main=1)
    legend(0,15000,c("lm","nls"),lty=c(2,1))
    
    plot(days,US, log = "y", ylim = c(100,15000))
    lines(days,predict(modnls))
    lines(days,exp(predict(modlm)), lty=2)
    title("log scale", cex.main=1)
    

在实践中正确建模随机噪声并不总是正确的

在实践中,问题并不经常是使用哪种模型来处理随机噪声(是否应该是某种 glm)。

问题更多的是指数模型(确定性部分)不正确,并且选择是否拟合线性化模型是在第一个点之间的强度与拟合最后一个点之间的选择。线性化模型非常适合小尺寸的值,非线性模型更适合高值的值。

当我们绘制增加的比率时,您可以看到指数模型的不正确性。

当我们绘制作为时间函数的世界变量的增加比率时,您可以看到它是一个非常量变量(在此期间它似乎在增加)。您可以为美国绘制相同的图,但它非常嘈杂,这是因为数字仍然很小,区分嘈杂的曲线会使噪声:信号比更大。

(另请注意,误差项将是递增的,如果您真的希望做对,那么您应该对错误使用某种 arima 类型的模型,或使用其他方式使误差项相关)

非恒定增长率


我仍然不明白为什么lmlog 会给我完全不同的系数。如何在两者之间转换?

glm 和 nls 将误差建模为 线性化模型将误差建模为 但是当你取值的对数时,你会改变相对大小。1000.1 和 1000 以及 1.1 和 1 之间的差值都是 0.1。但是在对数尺度上,它不再是相同的差异了。

yymodelN(0,σ2)
log(y)log(ymodel)N(0,σ2)

这实际上是 glm 进行拟合的方式。它使用线性模型,但对错误进行了转换权重(并且迭代了几次)。请参阅以下两个返回相同结果:

last_14 <- list(days <- 0:13,
                World <- c(101784,105821,109795, 113561,118592,125865,128343,145193,156094,167446,181527,197142,214910,242708),
                US <- c(262,402,518,583,959,1281,1663,2179,2727,3499,4632,6421,7783,13677))
days <- last_14[[1]]
US<- last_14[[3]]
World <- last_14[[2]]


Y <- log(US)
X <- cbind(rep(1,14),days)
coef <- lm.fit(x=X, y=Y)$coefficients
yp <- exp(X %*% coef)
for (i in 1:100) {
  # itterating with different
      # weights
      w <- as.numeric(yp^2)          
      # y-values
      Y <- log(US) + (US-yp)/yp
  # solve weighted linear equation
  coef <- solve(crossprod(X,w*X), crossprod(X,w*Y))
  # If am using lm.fit then for some reason you get something different then direct matrix solution
  # lm.wfit(x=X, y=Y, w=w)$coefficients
  yp <- exp(X %*% coef)
}
coef
# > coef
#           [,1]
#      5.2028935
# days 0.3267964

glm(US ~days,  
    family = gaussian(link = "log"), 
    control = list(epsilon = 10^-20, maxit = 100))

# > glm(US ~days,  
# +     family = gaussian(link = "log"), 
# +     control = list(epsilon = 10^-20, maxit = 100))
#
# Call:  glm(formula = US ~ days, family = gaussian(link = "log"), control = list(epsilon = 10^-20, 
#    maxit = 100))
#
# Coefficients:
# (Intercept)         days  
#      5.2029       0.3268  
#
# Degrees of Freedom: 13 Total (i.e. Null);  12 Residual
# Null Deviance:        185900000 
# Residual Deviance: 3533000    AIC: 219.9