我最喜欢的理解两个连续预测变量（例如热量和年份）之间相互作用的方法是在 x 轴上用一个预测变量的全范围和代表另一个预测变量潜在有趣值的几条不同的线来绘制它——我通常选择三个，代表预测变量的低、中和高水平，但您可以根据您的数据来玩弄它。

对于您的具体情况，这里是一些示例代码。首先，由于您没有提供示例数据，我将生成一些（非常粗略地）与您提供的参数估计相匹配的数据：

set.seed(24601) # to get exactly the same results as mine, use this set.seed()

# generate data
year <- 1996:2006
heat <- rnorm(n=20, mean=50, sd=15)
location <- letters[1:10]

library(tidyverse)
my_data <- base::expand.grid(year=year, heat=heat, location=location) %>% 
  mutate(yield = -75411.13 + 38.89*year - 82.56*heat + 0.04*year*heat + rnorm(nrow(.), sd = 50))

基本上，我只是使用了您模型中的固定效应估计值并添加了一些正常噪声。当然，出于您自己的目的，您将使用您的真实数据和真实拟合的模型。

模型 1（相互作用，仅限线性效应）

# run model
library(lme4)
mdl <- lmer(yield ~ year * heat  + (1|location), data = my_data)

从您的问题来看，您似乎主要对整体固定效应感兴趣（而不是对个别位置的估计等），所以我将重点关注这一点。如果您想可视化各个位置的结果，我推荐这篇很棒的博客文章：https ://tjmahr.github.io/plotting-partial-pooling-in-mixed-effects-models/

这是我估计的固定效应（它们与您的估计有些不同，因此您的情节看起来也会有些不同，但我认为它足够接近有用）：

> (fe <- summary(mdl)$coefficients[1:4,1])
  (Intercept)          year          heat     year:heat 
-74219.711477     38.293178   -117.339456      0.057416 

现在让我们绘制它来看看这种交互。我将选择放在heatx 轴上，并为 的低、中和高值设置三个代表线year，但如果您愿意，您可以切换并以其他方式进行。我将使用这些值（ 的全部范围heat和 的三个示例值）根据模型的固定效应估计为每个组合year生成预测分数。yield

# select values for plotting
# full range of heat, and selected low, med and high values for year
plot_df <- expand.grid(heat=heat, year=c(min(year), mean(year), max(year))) %>% 
  mutate(yield = fe[1] + fe[2]*year + fe[3]*heat + fe[4]*year*heat)

现在我将这些预测yield值绘制为线条。

plot_interaction <- ggplot(plot_df, aes(y=yield, x = heat, color = as.factor(year))) + 
  geom_line(size=2) + 
  labs(color = "year") + 
  # tweak plot appearance
  scale_color_brewer(palette = "Dark2") + 
  theme(legend.position = "top") + 
  theme_classic()

plot_interaction

那么我们在这里看到了什么？首先，相对于和的个体效应而言，相互作用非常小year——heat线条看起来几乎是平行的。但是您从参数估计中知道它与零有很大不同，因此即使它很小，它也存在。当= 0heat时有负面影响（理想情况下，您将在估计模型之前将两个预测变量居中，所以这将是“在平均值处”），并且正交互项表示随着增加，负面影响会减少，所以它减弱了。您应该在图中寻找的是，高年份的斜率比低年份的斜率要浅一些。yearyearyearheat

模型 2（线性和二次效应，具有交互作用）

第二个模型看起来要复杂得多，但实际上它几乎一样容易绘制！我们将执行与之前相同的过程，首先从模型中获取固定效应，然后生成一些绘图数据，其中包含 的全范围heat、选定的值和基于固定效应估计的year预测值。yield

首先，由于我没有您的真实数据，因此我正在重新生成另一个数据集，该数据集将（大致）匹配第二个模型的参数估计值，因此我们可以获得或多或少准确的图。

my_data <- base::expand.grid(year=year, heat=heat, location=location) %>% 
  mutate(yield = -97816.61 + 50.07*year + 20499.87*heat - 632.2*heat*heat - 10.22*year*heat +.31*year*heat*heat + rnorm(nrow(.), sd = 50))

运行第二个模型：

mdl <- lmer(yield ~ year*(heat + I(heat^2)) + (1|location), data = my_data)

获得固定效果：

fe <- summary(mdl)$coefficients[1:6,1]

> round(fe, 3)
   (Intercept)           year           heat      I(heat^2)      year:heat year:I(heat^2) 
    -95838.400         49.083      20410.924       -631.207        -10.176          0.310 

足够接近。:) 现在生成绘图数据。

# full range of heat, and selected low, med and high values for year
plot_df <- expand.grid(heat=heat, year=c(min(year), mean(year), max(year))) %>% 
  mutate(yield = fe[1] + fe[2]*year + fe[3]*heat + fe[4]*heat*heat + fe[5]*year*heat + fe[6]*year*heat*heat)

我们可以使用与以前完全相同的绘图代码，只需为其提供更新的数据框：

plot_interaction %+% plot_df

现在我们可以看到，on 的一般影响heat是yield负的，具有加速效应，因此在较高的热量水平下会有一个更陡峭的负斜率。此外，我们可以看到 的影响heat取决于year（即相互作用），因此heat在早年比晚年下降得更快。

您应该能够使用此代码为您自己的数据生成图表。我总是发现在解释交互时具有视觉效果非常有帮助，尤其是在具有多个系数的模型中（例如您的模型 2）。

（请注意，如果您的数据中 和 的范围heat与year我生成这些数据的结果完全不同，那么即使参数估计值相似，您的图也可能看起来完全不同。）

关于绘制和解释交互的最后说明

我从你的参数估计中猜测是你在估计模型之前没有居中year，heat这就是为什么我也没有。但你可能应该。它不仅缓解了一些模型估计问题，例如多重共线性，而且使解释交互更容易，因为 0 值更有意义。例如，您获得的参数估计值是当= 0时每个额外单位heat的预测变化。除非您已经居中（或者您正在研究古代的热量和产量），否则您可能正在谈论值远远超出您的数据范围。如果您首先居中，您对 的影响的估计将更有意义。yieldheatyearyearheat

对热求导，我们得到边际效应，即：

$\frac{\partial yield}{\partial heat}=2500-2*632*heat-10*year+2*0.31*year*heat$

我们可以通过以下方式重新排列 ME 感兴趣的部分：

$year*(-10+0.62*heat)$

这意味着我们有更多的热量，然后它的积极影响逐年下降。

所以假设我们有 10 个热量。那么接下来的每一年，效果都会比上一年下降。但是，假设例如。然后每年热效应将额外增加。$-10+0.62*10\approx-4$$heat>\frac{10}{0.62}$$heat=100$$-10+0.62*100\approx52$

此外，热量可能会逐年变化，因此其效果会有所不同。我们预期的时间趋势是正还是负取决于热值。