修改 lm 系数 - 为什么预测输出像修改截距一样移动？
modifying lm coefficients - why is predict output shifted like it had modified intercept?
我有一个 lm 模型，其死亡率数据取决于每日温度。为了估计对气候变化的可能适应，我想将曲线的斜率降低 10%。
因此，我通过乘以 0.9 修改了 lm 模型的斜率系数。
但是，这个修改后的模型的预测输出是出乎意料的。斜率下降，但曲线不在 x=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")
see plot here
我试图找出输出偏移的原因，我认为它在 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)
see curve plot here
为什么预测输出不同？
为什么示例的截距不为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 值，并且斜率也似乎不正确。
为什么这仍然出乎意料？
我认为问题在于 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' 表示修改后的系数，我们希望以下内容成立：
做一点代数，你会得到：
所以，为了实现这一点，假设您将一阶多项式项的系数乘以 .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))
我认为这就是您要查找的结果：
我有一个 lm 模型，其死亡率数据取决于每日温度。为了估计对气候变化的可能适应，我想将曲线的斜率降低 10%。 因此，我通过乘以 0.9 修改了 lm 模型的斜率系数。
但是，这个修改后的模型的预测输出是出乎意料的。斜率下降，但曲线不在 x=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")
see plot here
我试图找出输出偏移的原因，我认为它在 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)
see curve plot here
为什么预测输出不同？
为什么示例的截距不为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 值，并且斜率也似乎不正确。 为什么这仍然出乎意料？
我认为问题在于 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' 表示修改后的系数，我们希望以下内容成立：
做一点代数，你会得到：
所以，为了实现这一点，假设您将一阶多项式项的系数乘以 .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))
我认为这就是您要查找的结果：