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 回答