回归方程产生的模型超出了所有数据范围。

Question

我很困惑为什么我的回归方程在数据集中超出了所有数据的范围。我有一种感觉，这个方程对于具有大波动的数据非常敏感，但我仍然感到困惑。任何帮助都将不胜感激，统计学肯定不是我的母语！

参考文献：这是一个地球化学热力学问题：我试图将Maier-Kelley方程拟合到一些实验数据上。Maier-Kelley方程描述了平衡常数(K)，在这种情况下是白云石在水中溶解时随温度(T，以开尔文为单位)变化的方式。 log K = A + B.T + C/T + D.logT + E/T^2

简而言之（如果感兴趣，请参阅Hyeong和Capuano,2001），平衡常数(K)与Log_Ca_Mg（钙镁离子活性比）相同。

实验数据使用来自不同位置和不同深度的地下水数据（因此由FIELD和DepthID标识为我的随机变量）。

我包括了3个数据集。

（问题）数据集1：https://pastebin.com/fe2r2ebA

（工作中）数据集2：https://pastebin.com/gFgaJ2c8

（工作中）数据集3：https://pastebin.com/X5USaaNA

使用以下代码，对数据集1进行操作：

> dat1 <- read.csv("PATH_TO_DATASET_1.txt", header = TRUE,sep="\t")
> fm1 <- lmer(Log_Ca_Mg ~ 1 + kelvin + I(kelvin^-1) + I(log10(kelvin)) + I(kelvin^-2) + (1|FIELD) +(1|DepthID),data=dat1)

Warning messages:
1: Some predictor variables are on very different scales: consider rescaling 
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.0196619 (tol = 0.002, component 1)
3: Some predictor variables are on very different

> summary(fm1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Log_Ca_Mg ~ 1 + kelvin + I(kelvin^-1) + I(log10(kelvin)) + I(kelvin^-2) +      (1 | FIELD) + (1 | DepthID)
   Data: dat1

REML criterion at convergence: -774.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5464 -0.4538 -0.0671  0.3736  6.4217 

Random effects:
 Groups   Name        Variance Std.Dev.
 DepthID  (Intercept) 0.01035  0.1017  
 FIELD    (Intercept) 0.01081  0.1040  
 Residual             0.01905  0.1380  
Number of obs: 1175, groups:  DepthID, 675; FIELD, 410

Fixed effects:
                   Estimate Std. Error         df t value Pr(>|t|)
(Intercept)       3.368e+03  1.706e+03  4.582e-02   1.974    0.876
kelvin            4.615e-01  2.375e-01  4.600e-02   1.943    0.876
I(kelvin^-1)     -1.975e+05  9.788e+04  4.591e-02  -2.018    0.875
I(log10(kelvin)) -1.205e+03  6.122e+02  4.582e-02  -1.968    0.876
I(kelvin^-2)      1.230e+07  5.933e+06  4.624e-02   2.073    0.873

Correlation of Fixed Effects:
            (Intr) kelvin I(^-1) I(10()
kelvin       0.999                     
I(kelvn^-1) -1.000 -0.997              
I(lg10(kl)) -1.000 -0.999  0.999       
I(kelvn^-2)  0.998  0.994 -0.999 -0.997
fit warnings:
Some predictor variables are on very different scales: consider rescaling
convergence code: 0
Model failed to converge with max|grad| = 0.0196619 (tol = 0.002, component 1)

对于数据集2

> summary(fm2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Log_Ca_Mg ~ 1 + kelvin + I(kelvin^-1) + I(log10(kelvin)) + I(kelvin^-2) +      (1 | FIELD) + (1 | DepthID)
   Data: dat2

REML criterion at convergence: -1073.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0816 -0.4772 -0.0581  0.3650  5.6209 

Random effects:
 Groups   Name        Variance Std.Dev.
 DepthID  (Intercept) 0.007368 0.08584 
 FIELD    (Intercept) 0.014266 0.11944 
 Residual             0.023048 0.15182 
Number of obs: 1906, groups:  DepthID, 966; FIELD, 537

Fixed effects:
                   Estimate Std. Error         df t value Pr(>|t|)
(Intercept)      -9.366e+01  2.948e+03  1.283e-03  -0.032    0.999
kelvin           -2.798e-02  4.371e-01  1.289e-03  -0.064    0.998
I(kelvin^-1)      2.623e+02  1.627e+05  1.285e-03   0.002    1.000
I(log10(kelvin))  3.965e+01  1.067e+03  1.283e-03   0.037    0.999
I(kelvin^-2)      2.917e+05  9.476e+06  1.294e-03   0.031    0.999

Correlation of Fixed Effects:
            (Intr) kelvin I(^-1) I(10()
kelvin       0.999                     
I(kelvn^-1) -0.999 -0.997              
I(lg10(kl)) -1.000 -0.999  0.999       
I(kelvn^-2)  0.998  0.994 -0.999 -0.997
fit warnings:
Some predictor variables are on very different scales: consider rescaling
convergence code: 0
Model failed to converge with max|grad| = 0.0196967 (tol = 0.002, component 1)

对于数据集3

> summary(fm2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Log_Ca_Mg ~ 1 + kelvin + I(kelvin^-1) + I(log10(kelvin)) + I(kelvin^-2) +      (1 | FIELD) + (1 | DepthID)
   Data: dat3

REML criterion at convergence: -1590.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.2546 -0.4987 -0.0379  0.4313  4.5490 

Random effects:
 Groups   Name        Variance Std.Dev.
 DepthID  (Intercept) 0.01311  0.1145  
 FIELD    (Intercept) 0.01424  0.1193  
 Residual             0.03138  0.1771  
Number of obs: 6674, groups:  DepthID, 3422; FIELD, 1622

Fixed effects:
                   Estimate Std. Error         df t value Pr(>|t|)
(Intercept)       1.260e+03  1.835e+03  9.027e-02   0.687    0.871
kelvin            1.824e-01  2.783e-01  9.059e-02   0.655    0.874
I(kelvin^-1)     -7.289e+04  9.961e+04  9.044e-02  -0.732    0.866
I(log10(kelvin)) -4.529e+02  6.658e+02  9.028e-02  -0.680    0.872
I(kelvin^-2)      4.499e+06  5.690e+06  9.104e-02   0.791    0.860

Correlation of Fixed Effects:
            (Intr) kelvin I(^-1) I(10()
kelvin       0.999                     
I(kelvn^-1) -1.000 -0.997              
I(lg10(kl)) -1.000 -0.999  0.999       
I(kelvn^-2)  0.998  0.994 -0.999 -0.998
fit warnings:
Some predictor variables are on very different scales: consider rescaling
convergence code: 0
unable to evaluate scaled gradient
Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

我已经绘制了“所有数据”，但对于回归分析，红线上方或绿线下方没有数据。 在任何温度下，仅包括在红线和绿线之间具有log_ca_mg值的点进行回归分析。

看着数据集1的回归结果就很离谱，但是由于红线以上没有数据，这让我感到困惑不解。回归分析处在无数据区域，对于其他两个数据集来说并不是问题。即使是样本量更小的数据集（n=200），其回归分析也大致处在相同的区域。三个数据集在单独绘制时看起来相对类似。

我有些迷茫，请帮忙理解此情况，谢谢。

- Hamish Robertson

---

请原谅我之前的评论，它有点草率。但是看了你的代码后，我认为你看到的“模型无法收敛”的警告很值得注意。也许这个链接可以帮到你：https://stats.stackexchange.com/questions/242109/model-failed-to-converge-warning-in-lmer。此外，关于预测变量在非常不同的尺度上的警告也可能是需要考虑的问题。请参考gung在这篇文章中的回答：https://stats.stackexchange.com/questions/29781/when-conducting-multiple-regression-when-should-you-center-your-predictor-varia。 - xilliam

你在编码方面遇到了问题还是实际统计学背后的问题？如果是模型本身，我建议你在交叉验证上提出问题。但是针对你的问题，你是否担心所有预测变量之间的极端相关性？缺乏收敛表明你模型中的预测值将与实际值相差甚远。 - Mike

嗨，哈米什，我认为相关性问题表明可能存在“共线性”问题，这是一种情况，其中一个预测变量本身预测另一个变量。共线性往往会使模型不稳定，这意味着输入的微小变化会导致模型输出的大幅变化。你能描述一下为什么线性模型包含变量开尔文的所有转换吗？我的意思是kelvin^-1、log10(Kelvin)和kelvin^-2？你有一些包括它们的理由吗？如果没有充分的理由来包括它们，那么放弃其中一些可能是可行的。 - xilliam

你尝试过按照警告建议进行重新缩放吗？例如，Log_Ca_Mg ~ 1 + kelvin + I(100 * kelvin^-1) + I(log10(kelvin)) + I(1e4 * kelvin^-2) + (1 | FIELD) + (1 | DepthID)。 - Roland

嗨，罗兰。该模型的五个术语基本上可以立即用于其他热力学参数（如果感兴趣，请参见aqion.de/site/121，它在底部部分）。我想知道重新缩放是否会影响回归方程的有效性。我认为不会，对吧？您只是有效地改变了不同系数的“重点”。 - Hamish Robertson

Answer 1

我认为你的做法有误。听起来你想要估计Maier-Kelley方程中的参数A、B、C、D和E。你可以使用非线性最小二乘法来实现这一点，而不是使用线性混合效应模型。

首先定义一个函数来复制公式：

MK_eq <- function(A, B, C, D, E, Temp)
{
  A + B * Temp + C / Temp + D * log10(Temp) + E / (Temp^2)
}

现在我们使用nls函数来估计A到E的值：

mod1 <- nls(Log_Ca_Mg ~ MK_eq(A, B, C, D, E, kelvin), 
            start = list(A = 1, B = 1, C = 1, D = 1, E = 2), data = dat1)

coef(mod1)
#>             A             B             C             D             E 
#>  4.802008e+03  6.538166e-01 -2.818917e+05 -1.717040e+03  1.755566e+07 

我们可以通过在275到400之间以0.1为增量获取每个Kelvin值的预测，从而创建一个“回归线”：

new_data <- data.frame(kelvin = seq(275, 400, 0.1))
new_data$Log_Ca_Mg <- predict(mod1, newdata = new_data)

我们可以通过在样本上绘制我们的预测来证明这是一个很好的近似：

ggplot(dat1, aes(x = kelvin, y = Log_Ca_Mg)) + 
  geom_point() + 
  geom_line(data = new_data, linetype = 2, colour = "red", size = 2)

注意，为了简单起见，我避免讨论随机效应 - 可以使用 nlme 包进行混合效应非线性最小二乘，但这更复杂，这里的讨论比我在此处能做的更详细。

- Allan Cameron

---

我认为“需要这样做”并不准确。所提供的方程在参数上是线性的，因此可以使用标准线性方法进行估计。 - user2957945

@user2957945 好的，这确实会给出“正确”的结果。我会将其改为“可以”。 - Allan Cameron

OP的数据中似乎有重复的测量。为什么你建议只是忽略这个问题？ - Roland

如果线性混合效应模型无法收敛，非线性混合效应模型能够收敛的话，我会感到惊讶。期待完整的实例演示。 - Roland

@Roland 你可能是对的；在探索数据时，我们发现 DepthID 有 676 个唯一的因子水平，FIELD 有 410 个唯一的因子水平，而总共只有 1175 个数据点。 - Allan Cameron

Answer 2

以下是诊断模型可能出现的问题的尝试。本文将使用数据集1进行讨论：

根据您的问题描述，当使用数据集1运行原始模型时，会收到警告：

# original model
fm1 <- lme4::lmer(Log_Ca_Mg ~ 1 + kelvin + I(kelvin^-1) + I(log10(kelvin)) + I(kelvin^-2) + (1|FIELD) +(1|DepthID),data=dat1)

一些预测变量的规模差异很大：考虑重新调整收敛代码：0 模型无法收敛，最大|grad|=0.0196619 (tol=0.002, component 1)

这些和其他信息表明你的模型存在问题，可能与预测变量的规模有关。

由于 fm1 中有几个预测变量是 'kelvin' 变量的转换，我们还可以使用 car 包中的 vif 函数检查模型的共线性：

# examine collinearity with the vif (variance inflation factors)
> car::vif(fm1)
kelvin     I(kelvin^-1) I(log10(kelvin))     I(kelvin^-2) 
716333          9200929          7688348          1224275 

这些vif值表明fm1模型存在较高的共线性问题。 我们可以尝试删除其中一些预测变量，以检查一个更简单的模型：

fm1_b <- lme4::lmer(Log_Ca_Mg ~ 1 + kelvin + I(kelvin^-1) + (1|FIELD) +(1|DepthID),data=dat1)

当我们运行代码时，仍然会收到有关预测变量处于不同比例的警告：

警告信息：某些预测变量在比例上差异很大：请考虑重新缩放

同时，vif值要小得多：

# examine collinearity with the vif (variance inflation factors)
  > car::vif(fm1_b)
kelvin I(kelvin^-1) 
46.48406     46.48406 

根据gung在评论中提出的建议，我们可以看到当我们居中我们的开尔文变量时会发生什么：

dat1$kelvin_centered <- as.vector(scale(dat1$kelvin, center= TRUE, scale = FALSE ))
# Make a power transformation on the kelvin_centered variable
dat1$kelvin_centered_pwr <- dat1$kelvin_centered^-1

并检查它们是否相关

# check the correlation of the centered vars
cor(dat1$kelvin_centered, dat1$kelvin_centered_pwr)
> cor(dat1$kelvin_centered, dat1$kelvin_centered_pwr)
[1] 0.08056641

使用居中变量构建不同的模型：

# construct a modifed model
fm1_c <- lme4::lmer(Log_Ca_Mg ~ 1 + kelvin_centered + kelvin_centered_pwr + (1|FIELD) +(1|DepthID),data=dat1)

值得注意的是，我们在使用该模型运行代码时并没有看到任何警告。而且vif值非常低：

car::vif(fm1_c)

> car::vif(fm1_c)
    kelvin_centered kelvin_centered_pwr 
           1.005899            1.005899 

结论

原始模型存在高度共线性。共线性可能会使模型不稳定，这可能是为什么fm1无法收敛并且您在图形中看到奇怪预测的原因。模型fm1_c可能是或不是适合您目的的正确模型。它至少提供了一种理解原始模型问题的视角。

- xilliam

---

嗨，哈米什，是的。删除预测变量会改变事情。我主要是想帮助你找到前进的方向。你可以尝试添加kelvin^-2项，例如dat1$kelvin_centered_pwr2 <- dat1$kelvin_centered^-2，然后在模型中包括kelvin_centered_pwr2。至于log10(kelvin)项，我不确定该建议什么。当你应用log10(dat1$kelvin_centered)时，输出中会出现NaN。也许有更有经验的人会提供帮助。你还可以在统计堆栈交换网站上提出新问题。 - xilliam

嗨Xilliam，感谢您在这方面的努力，我非常感激。针对您之前的问题，我可以放弃使用预测器并使用更简单的模型，但是它们与我正在使用的5项模型有所不同（如果感兴趣，请参见https://www.aqion.de/site/121，它在底部部分中）。该模型的五个术语本质上可立即用于其他热力学参数。 - Hamish Robertson

再次感谢Xilliam :) 我会检查这些东西的。 - Hamish Robertson

最后一个想法是，在您的fm1方程式中，参数A-E在哪里？现在，fm1已编码以包括“kelvin”，我相信对应于您的方程式中的“T”变量。 - xilliam