回归方程在所有数据之外生成模型

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(钙与镁 ionic 活性的比率)相同。

实验数据使用来自不同位置和不同深度的地下水数据(由 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)的数据集，它也大致在相同的区域内。当单独绘制时，这三个数据集看起来相对相似。

我有点失落。任何帮助理解这一点将不胜感激。

Answer 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)

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)

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

自 fm1有几个预测变量是变量“开尔文”的变换，我们还可以使用 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)

当我们运行代码时，我们仍然会收到有关预测变量在不同尺度上的警告:

Warning message: Some predictor variables are on very different scales: consider rescaling

同时，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可能是也可能不是适合您目的的正确模型。它至少提供了一个镜头来了解原始模型的问题。