4

我很困惑为什么我会产生一个超出数据集中所有数据范围的回归方程。我有一种感觉，这个方程对传播范围很大的数据非常敏感，但我仍然感到困惑。任何帮助将不胜感激，统计数据当然不是我的第一语言！

作为参考，这是一个地球化学热力学问题：我试图将 Maier-Kelley 方程拟合到一些实验数据。Maier-Kelley 方程描述了平衡常数 (K)，在这种情况下，白云石溶解在水中，如何随温度变化（在这种情况下，T 为开尔文）。

日志 K = A + BT + 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），它也大致在同一区域。单独绘制时，这三个数据集看起来相对相似。

我有点迷路了。任何有助于理解这一点的帮助将不胜感激。

4

2 回答 2

1

接下来是尝试诊断您的模型可能出现的问题。它将使用数据集 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 模型未能与 max|grad| 收敛 = 0.0196619（tol = 0.002，组件 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可能是也可能不是适合您的目的的正确模型。它至少提供了一个镜头来了解您的原始模型的问题。

于 2020-02-25T10:27:19.077 回答
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

我们可以通过以 0.1 为增量对 275 到 400 之间的每个开尔文值进行预测来创建“回归线”：

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，但它涉及更多，这里的讨论比我在这里更详细地描述了如何做到这一点.

于 2020-03-01T13:18:25.710 回答