Уравнение регрессии создает модель вне всех данных - PullRequest
Новая
       37

Уравнение регрессии создает модель вне всех данных

3 голосов
/ 21 февраля 2020

Я совершенно не понимаю, почему я создаю уравнение регрессии, которое выходит за пределы диапазона всех данных в наборе данных. У меня такое чувство, что уравнение очень чувствительно к данным с большим разбросом, но я все еще в замешательстве. Любая помощь будет принята с благодарностью, статистика, безусловно, не мой родной язык!

Для справки это проблема геохимической термодинамики: я пытаюсь согласовать уравнение Майера-Келли с некоторыми экспериментальными данными. Уравнение Майера-Келли описывает, как константа равновесия (K), в данном случае доломит, растворяющийся в воде, изменяется с температурой (в данном случае T в Кельвинах).

log K = A + BT + C / T + D.logT + E / T ^ 2

Короче говоря (см. Hyeong and 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 между красной и зеленой линией при любой температуре включены в регрессионный анализ.

enter image description here

Итак, рассмотрим регрессии на Набор данных plot 1 еще далеко, но, поскольку над красной линией нет никаких данных, это просто не смущает меня. Регрессия находится в области, где нет данных. Для двух других наборов данных это не проблема. Даже для наборов данных с меньшими размерами (n = 200) это примерно в той же области. Три набора данных выглядят относительно похожими, когда отображаются отдельно.

Я немного потерян. Любая помощь в понимании этого будет принята.

Ответы [ 2 ]

1 голос
/ 01 марта 2020

Я думаю, вы поступаете неправильно. Похоже, вы пытаетесь оценить параметры 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: 

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)

enter image description here

Обратите внимание, что для простоты я избежал обсуждения случайных эффектов - возможно сделать смешанные эффекты нелинейными наименьшими квадратами, используя пакет nlme, но это более сложный процесс, и обсуждение здесь описывает, как сделать это более подробно, чем я могу здесь.

1 голос
/ 25 февраля 2020

Далее следует попытка диагностировать, что может быть не так с вашей моделью. Для этого обсуждения будет использоваться набор данных 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 (толь = 0,002, компонент 1)

Эта и другая информация указывает на наличие проблем в вашей модели, возможно, связанных с тем, что предикторы находятся в другом масштабе.

Поскольку fm1 имеет несколько предикторов, которые являются преобразованиями переменной 'kelvin', мы также можем проверить модель на коллинеарность с помощью функции car package 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

Следуя предложению Ганга, которое я упомянул в комментариях, мы можем видеть, что происходит, когда мы центрируем наши переменные Кельвина: 

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

A и проверьте, не коррелированы ли они 

# 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 может быть или не быть подходящей моделью для ваших целей. Это, по крайней мере, обеспечивает объектив, чтобы понять проблему с вашей оригинальной моделью.

...