1

我有一个 lm 模型，其死亡率数据取决于每日温度。为了估计对气候变化的可能适应，我想将曲线的斜率降低 10%。因此，我通过乘以 0.9 来修改 lm 模型的斜率系数。

然而，这个修改后的模型的预测输出是出乎意料的。斜率下降，但曲线不是在 x=0 处相交，而是在大约 0 处的截距处相交。133.那就是下一个问题，为什么这个截距不是0？

这是一个可重现的示例：

x <- seq(0, 20, 0.1)
y <- seq(0, 20, 0.1)^2

mod <- lm(y ~ poly(x, 2))
mod$coefficients
(Intercept) poly(x, 2)1 poly(x, 2)2 
133.6667   1645.2355    426.9008  

mody <- mod
mody$coefficients[2] <- mody$coefficients[2]*0.9
mody$coefficients[3] <- mody$coefficients[3]*0.9
mody$coefficients
Intercept) poly(x, 2)1 poly(x, 2)2 
133.6667   1480.7120    384.2108 

plot(x, predict(mod), type="l")
lines(x, predict(mody), col="red")

在这里看情节

我试图找出偏移输出的原因，我认为它在 predict() 函数中的某个地方。为了检查修改后的系数，我尝试了这个，它显示了扩展输出。

curve(coef(mod)[1] + coef(mod)[2] * x + coef(mod)[3] * x^2, from=0, to=20, xlab="x", ylab="y")
curve(coef(mody)[1] + coef(mody)[2] * x + coef(mody)[3] * x^2, from=0, to=20,xlab="x", ylab="y", add = T)

在此处查看曲线图

为什么预测输出不同？

为什么示例的 Intercept 不为 0？

或者如何在不使用 predict() 的情况下“手动”“预测”修改后的数据？

谢谢你的帮助！

编辑：DaveArmstrong 的答案通过在 poly() 中使用 raw=TRUE 解决了第一个示例的问题。但是，对于其他数据，它仍然无法正常工作，可能是由于模型中的负系数（？）

这是我的原始数据的一个例子：

x <- seq(15.0, 23.5, 0.1)
y <- c(0.992, 0.998, 1.012, 1.013, 1.015, 1.021, 1.028, 1.027, 1.023, 1.029, 1.032, 1.032, 1.029, 1.036, 1.035, 1.041, 1.043, 1.043, 1.037, 1.037, 1.039, 1.037, 1.041, 1.047, 1.047, 1.048, 1.045, 1.048, 1.044, 1.037, 1.046, 1.042, 1.037, 1.034, 1.032, 1.031, 1.030, 1.034,
1.044, 1.046, 1.036, 1.034, 1.049, 1.050, 1.037, 1.041, 1.046, 1.062, 1.077, 1.084, 1.091, 1.106, 1.114, 1.127, 1.120, 1.122, 1.130,
1.122, 1.135, 1.164, 1.187, 1.186, 1.195, 1.201, 1.197, 1.204, 1.201, 1.205, 1.203, 1.200, 1.205, 1.232, 1.218, 1.218, 1.249, 1.245,
1.253, 1.227, 1.232, 1.252, 1.258, 1.254, 1.248, 1.245, 1.261, 1.289)

org <- lm(y ~ poly(x, 2, raw=TRUE))
coef(org)
(Intercept) poly(x, 2, raw = TRUE)1 poly(x, 2, raw = TRUE)2 
            2.240583377            -0.153426285             0.004822839 
  
orgm <- org
orgm$coefficients[2] <- orgm$coefficients[2]*1.1 #reducing negative coefficient
orgm$coefficients[3] <- orgm$coefficients[3]*0.9

plot(predict(org), ylim=c(0,1.5), type="l")
lines(predict(orgm), col="red")
legend("topleft", legend=c("Original", "Modified"), col=c("black", "red"), lty=c(1,1))

在结果图 ( plot ) 中，修改后的模型以某种方式转移到较低的 y 值，并且斜率也似乎不正确。为什么这仍然出乎意料？

4

1 回答 1

2

我认为问题在于poly()默认情况下该函数正交化多项式回归量。在您的示例中，数据中的平方项之间实际上只有一种关系。如果您改为使用原始多项式执行此操作，它应该可以工作。

x <- seq(0, 20, 0.1)
y <- seq(0, 20, 0.1)^2

mod <- lm(y ~ poly(x, 2, raw=TRUE))
mod$coefficients
# (Intercept) poly(x, 2, raw = TRUE)1 poly(x, 2, raw = TRUE)2 
# -6.961533e-14            1.658415e-14            1.000000e+00 

mody <- mod
mody$coefficients[2] <- mody$coefficients[2]*0.9
mody$coefficients[3] <- mody$coefficients[3]*0.9
mody$coefficients
# (Intercept) poly(x, 2, raw = TRUE)1 poly(x, 2, raw = TRUE)2 
# -6.961533e-14            1.492574e-14            9.000000e-01 

plot(x, predict(mod), type="l")
lines(x, predict(mody), col="red")
legend("topleft", legend=c("Original", "Modified"), col=c("black", "red"), lty=c(1,1))

在此处输入图像描述

对于更多的上下文，这里是在这个例子中正交多项式与原始多项式的关系（第一列给出了将原始多项式与第一个正交多项式回归量相关的系数，第二列给出了将原始多项式与原始多项式相关联的系数。二阶正交多项式回归器）。

p2 <- poly(x, 2)
round(coef(lm(p2 ~ poly(x, 2, raw=TRUE))), 5)
#                                1        2
# (Intercept)             -0.12156  0.15538
# poly(x, 2, raw = TRUE)1  0.01216 -0.04685
# poly(x, 2, raw = TRUE)2  0.00000  0.00234

用正交多项式将这些代入方程，您将得到以下结果（其中在此处输入图像描述表示正交回归量）：

在此处输入图像描述

当您将正交多项式系数乘以 0.9 时，您正在执行以下操作：

在此处输入图像描述

对于原始变量，当您修改正交回归量的系数时，您也在更改截距。

编辑：修改答案以处理真实数据

上面的解决方案有效，因为感兴趣的关系相对简单——一阶项的截距和系数都近似为零。如果不是这种情况，答案就复杂一些。在上面提出的真实数据示例中，x变量的最小值为 15。我的假设是我们希望两条曲线在 15 处相交，但修改后的曲线具有更浅的斜率。为此，我们需要根据原始系数和修改后的系数来考虑这意味着什么。特别是，当 x=15 时，我们需要这两个方程来产生相同的预测。用于b表示原始系数和b'表示修改后的系数，我们希望以下是正确的：

做一点代数，你会得到：

因此，实现这一点，假设您将一阶多项式项的系数乘以 0.9，这将得到：

orgm <- org
orgm$coefficients[2] <- orgm$coefficients[2]*0.9 
orgm$coefficients[2]
# poly(x, 2, raw = TRUE)1 
# -0.1379442

然后我们可以计算原始系数和修改系数之间的差异：

diff <- org$coefficients[2] - orgm$coefficients[2]
diff
# poly(x, 2, raw = TRUE)1 
# -0.01532713

最后，我们可以将这个和二阶多项式回归器的原始系数代入公式，以创建修改后的二阶多项式回归器系数：

orgm$coefficients[3] <- diff/15 + org$coefficients[3] 
orgm$coefficients
# (Intercept) poly(x, 2, raw = TRUE)1 poly(x, 2, raw = TRUE)2 
# 2.239156804            -0.137944190             0.003796868

然后，我们可以制作情节：

plot(x, predict(org), ylim=c(0,1.5), type="l")
lines(x, predict(orgm), col="red")
legend("topleft", legend=c("Original", "Modified"), col=c("black", "red"), lty=c(1,1))

我认为这是您正在寻找的结果：

在此处输入图像描述

于 2021-11-29T16:24:35.173 回答