如何将对数回归拟合到 R

如何解决如何将对数回归拟合到 R

我有一个每日降雨量 (x) 和观测值 (y) 的散点图，它看起来像是 x^-2 图的右/正 x 值一半或以 1/2 为底的对数图。基本上，当 x 值非常低时，y 值非常高。 x 值越大，y 变得越低。但 y 值下降的速度变慢，而且 y 永远不会为负。

这是一个代表性示例：

 rain <- c(1,1.2,1.3,2.5,3.2,4.2,5,7,7.5,10.3,11.7,12.9,14.1,15,15.5,17.5,18.3,20,20.2,20.3,25,28,30,34,40)

 obs <- c(42,44,43.9,43.5,35,22,18.4,15.3,10,6.2,5.7,4,3.7,2.3,2,2.7,3.5,3,2.9,1.6,2.2,0.8)

现在我想为这个散点图拟合一个回归模型。我已经尝试了多项式回归直到 x^-4，但我也想尝试对数回归，因为我认为它可能会成为一个更高质量的模型。

到目前为止，我使用多项式模型做了以下工作：

    y <- data$obs
    x <- data$rain
    xsq <- x^-2
    xcub <- x^-3
    xquar <- x^-4
    

    fit4 <- lm(y~x+xsq+xcub+xquar) # I did the same for fit 1-3; until fit 4 it becomes more significant
    xv <- seq(min(x),max(x),0.01)
    yv <- predict(fit5,list(x=xv,xsq=xv^-2,xcub=xv^-3,xquar=xv^-4))
    lines(xv,yv)

这就是我为对数模型尝试的方法，但它只返回与曲线不匹配的直线。我觉得 log() 不是我真正需要的功能。

xlog <- log(x)
fitlogx <- lm(y~xlog)
xv <- seq(min(xlog),max(xlog),0.01)
yv <- predict(fitlogx,list(x=xv))
abline(fitlogx)

ylog <- log(y)
fitlogy <- lm(ylog~x)
xv <- seq(min(x),0.01)
yv <- predict(fitlogy,list(x=xv))
abline(fitlogy)

现在我想知道如何拟合一个有意义的对数函数。如果您知道另一种可能有用的回归模型，我也非常感谢您提供任何建议。

解决方法

您的 obs 变量非常适合 rain 的倒数。例如

dev.new(width=12,height=6)
oldp <- par(mfrow=c(1,2))
plot(obs~rain)
lines(rain,1/rain*40)

曲线需要高一点。我们可以反复猜测，例如试试rain*60，但使用nls函数来获得方程的最佳最小二乘拟合更容易：

obs.nls <- nls(obs~1/rain*k,start=list(k=40))
summary(obs.nls)
# 
# Formula: obs ~ 1/rain * k
# 
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
# k   57.145      4.182   13.66 8.12e-13 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 6.915 on 24 degrees of freedom
# 
# Number of iterations to convergence: 1 
# Achieved convergence tolerance: 2.379e-09
plot(obs~rain)
pred <- predict(obs.nls)
points(rain,pred,col="red",pch=18)
pred.rain <- seq(1,40,length.out=100)
pred.obs <- predict(obs.nls,list(rain=pred.rain))
lines(pred.rain,pred.obs,col="blue",lty=2)

因此，k 的最佳估计值为 57.145。 nls 的主要缺点是您必须为系数提供起始值。它也可能无法收敛，但对于我们在这里使用的简单函数，只要您能估计出合理的起始值，它就可以正常工作。

如果 rain 可以有零值，则可以添加截距：

obs.nls <- nls(obs ~ k / (a + rain),start=list(a=1,k=40))
summary(obs.nls)
# 
# Formula: obs ~ k/(a + rain)
# 
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
# a   1.4169     0.4245   3.337  0.00286 ** 
# k 117.5345    16.6878   7.043 3.55e-07 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 4.638 on 23 degrees of freedom

Number of iterations to convergence: 10 
Achieved convergence tolerance: 6.763e-06

请注意，标准误差较小，但曲线高估了 rain > 10 的实际值。