Back to Article
Coefficient Analysis
Download Source

Kinetic Model Coefficient Analysis

Author

Marc Pérez

Published

May 22, 2025

Other Links
Presentation
In [1]:
Code 

suppressPackageStartupMessages({
  library(multcomp)
  library(car)
  library(tidyr)
  library(lme4)
  library(ggplot2)
  library(ggtext)
  library(ggpmisc)
  library(nlme)
  library(latex2exp)
  library(kableExtra)
  library(broom)
  library(dplyr)
  library(MuMIn)
})

options(warn = -1)
RES <- readRDS("~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/RES.rds")
D <- RES$D
d <- RES$data

Model Agroscope \[Y_{rel}\sim A*(1-e^{rate*P_{CO_2}+Env})\]

Wir ersetzen nur rate mit unserer Schätzung k: \[Y_{rel}\sim A*(1-e^{k*P_{CO_2}+Env})\]

Sind unsere Modelparameter gute Prediktoren?? \[Y_{rel}\sim A*(1-e^{k*PS+Env} )\]

Es gibt noch die Kovariaten Niederschlag pro Jahr, Jahresdurchschnittstemperatur und Temperatur in Jugendphase

In [2]:
Code 

library(GGally)

ggpairs(D, 
  aes(col=Site, shape = Treatment,alpha = 0.6), 
  columns = c("soil_0_20_P_AAE10", "soil_0_20_P_CO2", "PS", "k", "kPS"),
  lower = list(continuous = wrap("points", size = 1.3)),
  upper = list(continuous = "blank", combo = "blank", discrete = "blank"))  # Adjust size here
  




p6 <- ggplot(D,aes(y=soil_0_20_P_AAE10, x=soil_0_20_P_CO2, col=Site, size = Treatment)) +
  geom_point(shape = 7) + 
  scale_x_log10() + scale_y_log10() +
  labs(x=TeX("$P_{H_2O10}(mg/kg Soil)$"),
         y=TeX("$P_{AAEDTA}(mg/kg Soil)$")); p6
  
p7 <- ggplot(D,aes(y=PS, x=soil_0_20_P_CO2, col=Site, size = Treatment)) +
  geom_point(shape = 7) +
  scale_x_log10() + scale_y_log10() +
  labs(x=TeX("$P_{CO_2}(mg/kg Soil)$"),
          y=TeX("$PS(mg/kg Soil)$")); p7

p8 <- ggplot(D,aes(y=k, x=soil_0_20_P_CO2, col=Site, size = Treatment)) +
  geom_point(shape = 7) +
  scale_x_log10() +
  labs(x=TeX("$P_{CO_2}(mg/kg Soil)$"),
         y=TeX("$k(1/s)$")); p8

p9 <- ggplot(D,aes(y=k*PS, x=soil_0_20_P_CO2, col=Site, size = Treatment)) +
  geom_point(shape = 7) +
  scale_x_log10() + scale_y_log10() +
  labs(x=TeX("$P_{CO_2}(mg/kg Soil)$"),
         y=TeX("$v=k*PS(mg/s*kg Soil)$"));p9


p11 <- ggplot(D,aes(y=PS, x=soil_0_20_P_AAE10, col=Site, size = Treatment)) +
  scale_x_log10() + scale_y_log10() +
  geom_point(shape = 7) +
  labs(x=TeX("$P_{AAEDTA}(mg/kg Soil)$"),
         y=TeX("$PS(mg/kg Soil)$")); p11

p12 <- ggplot(D,aes(y=k, x=soil_0_20_P_AAE10, col=Site, size = Treatment)) +
  geom_point(shape = 7) +
   scale_x_log10() + scale_y_log10() +
  labs(x=TeX("$P_{AAEDTA}(mg/kg Soil)$"),
         y=TeX("$k(1/s)$"))

p12

p13 <- ggplot(D,aes(y=k*PS, x=soil_0_20_P_AAE10, col=Site, size = Treatment)) +
  scale_x_log10() + scale_y_log10() +
  geom_point(shape = 7) +
  labs(x=TeX("$P_{AAEDTA}(mg/kg Soil)$"),
         y=TeX("$log(v)=log(k*PS)(mg/s*kg Soil)$"))

p13

Nun noch die Linearen Regressionen, die ausstehend sind:

(1|year) + (1|Site) + (1|Site:block) + (Treatment|Site)

Random intercept per year and site, block nested in site. and Treatment nested in site (could also be modelled as a random slope to allow for correlations)

wir sind abe nicht an einem Treatment effekt interesseiert. drum verwerfen wir Treatment als Random UND Fixed effekt.

  1. Vergleiche PS, k und kPS mit
In [3]:
Code 


# Wovon hängen Modelparameter ab?

library(lmerTest)

Attaching package: 'lmerTest'
The following object is masked from 'package:lme4':

    lmer
The following object is masked from 'package:stats':

    step
Code 
fit.soil.PS  <- lmer(log(PS)     ~ soil_0_20_clay+ soil_0_20_pH_H2O + soil_0_20_Corg + soil_0_20_silt + Fed + Ald  + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
# fit.soil.PS2 <- lmer(log(PS)     ~ soil_0_20_clay+ soil_0_20_pH_H2O + soil_0_20_Corg + soil_0_20_silt  + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
fit.soil.k   <- lmer(k           ~ soil_0_20_clay+ soil_0_20_pH_H2O + soil_0_20_Corg + soil_0_20_silt + Fed + Ald + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.soil.kPS <- lmer(I(log(k*PS))~ soil_0_20_clay+ soil_0_20_pH_H2O + soil_0_20_Corg + soil_0_20_silt + Fed + Ald + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.soil.kPS2<- lmer(I(k*log(PS))~ soil_0_20_clay+ soil_0_20_pH_H2O + soil_0_20_Corg + soil_0_20_silt + Fed + Ald + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.soil.CO2 <- lmer(log(soil_0_20_P_CO2)~ soil_0_20_clay+ soil_0_20_pH_H2O + soil_0_20_Corg + soil_0_20_silt + Fed + Ald + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.soil.AAE10<-lmer(log(soil_0_20_P_AAE10)~ soil_0_20_clay+ soil_0_20_pH_H2O + soil_0_20_Corg + soil_0_20_silt + Fed + Ald + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
r.squaredGLMM(fit.soil.kPS)
           R2m       R2c
[1,] 0.3681407 0.9929139
Code 
r.squaredGLMM(fit.soil.kPS2)
           R2m       R2c
[1,] 0.6521416 0.8191785
Code 
fit.kin.Yrel     <- lmer(Ymain_rel         ~ k * log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k + log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
boundary (singular) fit: see help('isSingular')
Code 
r.squaredGLMM(fit.kin.Ynorm)
            R2m       R2c
[1,] 0.01104564 0.3620402
Code 
fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
boundary (singular) fit: see help('isSingular')
Code 
r.squaredGLMM(fit.kin.Ynorm)
            R2m       R2c
[1,] 0.01391747 0.3596733
Code 
fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS)  + (1|year), D, subset = Treatment != "P166")
r.squaredGLMM(fit.kin.Ynorm)
            R2m       R2c
[1,] 0.05114188 0.2399041
Code 
fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D)
boundary (singular) fit: see help('isSingular')
Code 
r.squaredGLMM(fit.kin.Ynorm)
            R2m       R2c
[1,] 0.02721952 0.2448253
Code 
fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment == "P0")
boundary (singular) fit: see help('isSingular')
Code 
r.squaredGLMM(fit.kin.Ynorm)
            R2m      R2c
[1,] 0.01708162 0.484134
Code 
summary(fit.kin.Ynorm, ddf="Kenward-Roger")
Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: Ymain_norm ~ k * log(PS) + (1 | year) + (1 | Site) + (1 | Site:block)
   Data: D
 Subset: Treatment == "P0"

REML criterion at convergence: 209.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5723 -0.2209  0.0187  0.3483  3.8086 

Random effects:
 Groups     Name        Variance Std.Dev.
 Site:block (Intercept) 0.0000   0.0000  
 year       (Intercept) 0.1024   0.3200  
 Site       (Intercept) 0.1634   0.4042  
 Residual               0.2935   0.5418  
Number of obs: 120, groups:  Site:block, 20; year, 6; Site, 5

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)
(Intercept)  1.738343   2.071122 12.765971   0.839    0.417
k           -1.488991  11.262130 12.439272  -0.132    0.897
log(PS)      0.104222   0.606923 12.181857   0.172    0.866
k:log(PS)    0.006624   3.296268 12.185270   0.002    0.998

Correlation of Fixed Effects:
          (Intr) k      lg(PS)
k         -0.944              
log(PS)    0.985 -0.943       
k:log(PS) -0.932  0.992 -0.946
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Code 
fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment == "P100")
r.squaredGLMM(fit.kin.Ynorm)
            R2m       R2c
[1,] 0.03937992 0.3521636
Code 
fit.grud.CO2.Ynorm     <- lmer(Ymain_norm    ~ log(soil_0_20_P_CO2) + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
r.squaredGLMM(fit.grud.CO2.Ynorm)
          R2m      R2c
[1,] 0.218175 0.358219
Code 
fit.grud.AAE10.Ynorm   <- lmer(Ymain_norm    ~ log(soil_0_20_P_AAE10) + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
r.squaredGLMM(fit.grud.AAE10.Ynorm)
           R2m       R2c
[1,] 0.1976951 0.4740521
Code 
fit.grud.CO2.AAE10.Ynorm <- lmer(Ymain_norm  ~ log(soil_0_20_P_CO2) * log(soil_0_20_P_AAE10) + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
r.squaredGLMM(fit.grud.CO2.AAE10.Ynorm)
           R2m       R2c
[1,] 0.2204812 0.3650344
Code 
# compare with k*log(PS)
fit.kin.Ynorm   <- lmer(Ymain_norm       ~ k*log(PS) + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
boundary (singular) fit: see help('isSingular')
Code 
r.squaredGLMM(fit.kin.Ynorm)
            R2m       R2c
[1,] 0.01391747 0.3596733
Code 
fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS) + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
boundary (singular) fit: see help('isSingular')
Code 
fit.kin.Pexport  <- lmer(annual_P_uptake  ~ k * log(PS) + (1|year) + (1|Site)  + (1|Site:block), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.kin.Pbalance <- lmer(annual_P_balance ~ k * log(PS) + (1|year) + (1|Site)  + (1|Site:block), D)

car::vif(lm(Ymain_norm       ~ (k) * log(PS)  + crop, D))
there are higher-order terms (interactions) in this model
consider setting type = 'predictor'; see ?vif
              GVIF Df GVIF^(1/(2*Df))
k         12.64859  1        3.556485
log(PS)   12.32336  1        3.510465
crop       1.05128  4        1.006271
k:log(PS) 29.55436  1        5.436392
Code 
car::vif(lm(Ymain_norm       ~ scale(k) * scale(log(PS))   + crop, D))
there are higher-order terms (interactions) in this model
consider setting type = 'predictor'; see ?vif
                            GVIF Df GVIF^(1/(2*Df))
scale(k)                1.130309  1        1.063160
scale(log(PS))          1.074904  1        1.036776
crop                    1.051280  4        1.006271
scale(k):scale(log(PS)) 1.019256  1        1.009582
Code 
car::vif(lm(Ymain_norm       ~ k + log(PS)   + crop, D))
            GVIF Df GVIF^(1/(2*Df))
k       1.121012  1        1.058778
log(PS) 1.069713  1        1.034270
crop    1.050571  4        1.006186
Code 
car::vif(lm(Ymain_norm       ~ I(k * (PS)) + k + log(PS)    + crop, D))
                GVIF Df GVIF^(1/(2*Df))
I(k * (PS)) 4.642883  1        2.154735
k           2.035796  1        1.426813
log(PS)     4.711359  1        2.170567
crop        1.052981  4        1.006474
Code 
with(D, hist(I(exp(k) * (PS))))

Code 
with(D, hist(log(PS)))

Code 
r.squaredGLMM(fit.kin.Ynorm)
            R2m       R2c
[1,] 0.01391747 0.3596733
Code 
anova(fit.soil.PS)
Type III Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
soil_0_20_clay   0.0017  0.0017     1 39.832   0.0467   0.82993    
soil_0_20_pH_H2O 0.0007  0.0007     1 38.367   0.0202   0.88762    
soil_0_20_Corg   0.1525  0.1525     1 32.811   4.2087   0.04827 *  
soil_0_20_silt   0.0000  0.0000     1 40.483   0.0002   0.98792    
Fed              4.9643  4.9643     1 36.018 137.0324 7.807e-14 ***
Ald              1.8349  1.8349     1 37.406  50.6502 1.865e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
summary(glht(fit.soil.PS))

     Simultaneous Tests for General Linear Hypotheses

Fit: lmer(formula = log(PS) ~ soil_0_20_clay + soil_0_20_pH_H2O + 
    soil_0_20_Corg + soil_0_20_silt + Fed + Ald + (1 | year) + 
    (1 | Site) + (1 | Site:block) + (1 | Site:Treatment), data = D)

Linear Hypotheses:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept) == 0      21.4438130  2.6915204   7.967   <0.001 ***
soil_0_20_clay == 0   -0.0064543  0.0298519  -0.216    1.000    
soil_0_20_pH_H2O == 0 -0.0179677  0.1262969  -0.142    1.000    
soil_0_20_Corg == 0    0.6116152  0.2981307   2.052    0.202    
soil_0_20_silt == 0   -0.0003927  0.0257799  -0.015    1.000    
Fed == 0              -1.0676543  0.0912051 -11.706   <0.001 ***
Ald == 0              -8.7060240  1.2232898  -7.117   <0.001 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
Code 
# Fazit: PS wird von treatment stark beeinfluss, k eher nicht (dafür von site)
anova(fit.soil.k)
Type III Analysis of Variance Table with Satterthwaite's method
                    Sum Sq   Mean Sq NumDF  DenDF F value   Pr(>F)   
soil_0_20_clay   0.0126524 0.0126524     1 32.124 11.0892 0.002190 **
soil_0_20_pH_H2O 0.0010892 0.0010892     1 38.528  0.9546 0.334637   
soil_0_20_Corg   0.0091554 0.0091554     1 31.652  8.0242 0.007963 **
soil_0_20_silt   0.0012806 0.0012806     1 13.842  1.1224 0.307539   
Fed              0.0001832 0.0001832     1  2.350  0.1606 0.722155   
Ald              0.0002817 0.0002817     1  2.262  0.2469 0.663422   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
summary(glht(fit.soil.k))

     Simultaneous Tests for General Linear Hypotheses

Fit: lmer(formula = k ~ soil_0_20_clay + soil_0_20_pH_H2O + soil_0_20_Corg + 
    soil_0_20_silt + Fed + Ald + (1 | year) + (1 | Site) + (1 | 
    Site:block) + (1 | Site:Treatment), data = D)

Linear Hypotheses:
                       Estimate Std. Error z value Pr(>|z|)   
(Intercept) == 0       0.454128   0.358773   1.266  0.70480   
soil_0_20_clay == 0   -0.016164   0.004854  -3.330  0.00542 **
soil_0_20_pH_H2O == 0 -0.020791   0.021279  -0.977  0.87818   
soil_0_20_Corg == 0    0.136661   0.048244   2.833  0.02876 * 
soil_0_20_silt == 0    0.004076   0.003848   1.059  0.83553   
Fed == 0               0.005157   0.012867   0.401  0.99838   
Ald == 0              -0.072344   0.145598  -0.497  0.99477   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
Code 
anova(fit.soil.kPS)
Type III Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
soil_0_20_clay   0.2899  0.2899     1 40.076  3.4177 0.07189 .
soil_0_20_pH_H2O 0.0186  0.0186     1 38.805  0.2198 0.64181  
soil_0_20_Corg   0.6111  0.6111     1 33.275  7.2042 0.01125 *
soil_0_20_silt   0.0039  0.0039     1 30.169  0.0457 0.83217  
Fed              4.9569  4.9569     1  1.276 58.4378 0.05081 .
Ald              1.7919  1.7919     1  2.143 21.1247 0.03873 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
summary(glht(fit.soil.kPS))

     Simultaneous Tests for General Linear Hypotheses

Fit: lmer(formula = I(log(k * PS)) ~ soil_0_20_clay + soil_0_20_pH_H2O + 
    soil_0_20_Corg + soil_0_20_silt + Fed + Ald + (1 | year) + 
    (1 | Site) + (1 | Site:block) + (1 | Site:Treatment), data = D)

Linear Hypotheses:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept) == 0      21.18866    3.96710   5.341   <0.001 ***
soil_0_20_clay == 0   -0.08522    0.04610  -1.849   0.3025    
soil_0_20_pH_H2O == 0 -0.09210    0.19644  -0.469   0.9965    
soil_0_20_Corg == 0    1.25026    0.46581   2.684   0.0423 *  
soil_0_20_silt == 0    0.00845    0.03953   0.214   1.0000    
Fed == 0              -1.08433    0.14184  -7.644   <0.001 ***
Ald == 0              -8.63139    1.87796  -4.596   <0.001 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
Code 
anova(fit.kin.Yrel)
Type III Analysis of Variance Table with Satterthwaite's method
           Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
k         10361.3 10361.3     1 243.49  7.0248 0.008566 **
log(PS)    9132.5  9132.5     1 243.72  6.1918 0.013504 * 
k:log(PS) 12026.0 12026.0     1 243.87  8.1535 0.004668 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
anova(fit.kin.Ynorm)
Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
k         0.182097 0.182097     1 229.90  0.6950 0.4053
log(PS)   0.029224 0.029224     1 229.48  0.1115 0.7387
k:log(PS) 0.232300 0.232300     1 229.46  0.8866 0.3474
Code 
anova(fit.kin.Pexport)
Type III Analysis of Variance Table with Satterthwaite's method
          Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
k         35.305  35.305     1 247.74  0.5252 0.4693
log(PS)   39.586  39.586     1 248.15  0.5889 0.4436
k:log(PS) 53.624  53.624     1 248.29  0.7978 0.3726
Code 
anova(fit.soil.CO2)
Type III Analysis of Variance Table with Satterthwaite's method
                  Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
soil_0_20_clay   0.01031 0.01031     1 39.033  0.3482 0.55854  
soil_0_20_pH_H2O 0.03990 0.03990     1 36.270  1.3477 0.25326  
soil_0_20_Corg   0.03609 0.03609     1 30.050  1.2192 0.27829  
soil_0_20_silt   0.01225 0.01225     1 36.371  0.4139 0.52405  
Fed              2.24907 2.24907     1  1.237 75.9765 0.04625 *
Ald              0.37933 0.37933     1  1.859 12.8143 0.07780 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
anova(fit.soil.AAE10)
Type III Analysis of Variance Table with Satterthwaite's method
                  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
soil_0_20_clay   0.02689 0.02689     1 33.372  1.1447   0.29233    
soil_0_20_pH_H2O 0.00022 0.00022     1 36.721  0.0095   0.92271    
soil_0_20_Corg   0.55693 0.55693     1 38.895 23.7089 1.902e-05 ***
soil_0_20_silt   0.09825 0.09825     1 34.019  4.1824   0.04864 *  
Fed              0.42176 0.42176     1  1.848 17.9545   0.05892 .  
Ald              0.70761 0.70761     1  2.551 30.1233   0.01795 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
anova(fit.kin.Pbalance)
Type III Analysis of Variance Table with Satterthwaite's method
           Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
k          411.53  411.53     1 228.31  2.6679    0.1038    
log(PS)   2451.67 2451.67     1 236.98 15.8941 8.925e-05 ***
k:log(PS)  335.79  335.79     1 232.66  2.1769    0.1414    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
summary(fit.kin.Pbalance)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: annual_P_balance ~ k * log(PS) + (1 | year) + (1 | Site) + (1 |  
    Site:block)
   Data: D

REML criterion at convergence: 2172.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.2416 -0.5978  0.0314  0.5493  2.8712 

Random effects:
 Groups     Name        Variance Std.Dev.
 Site:block (Intercept)  20.33    4.508  
 year       (Intercept)  59.32    7.702  
 Site       (Intercept)  23.97    4.896  
 Residual               154.25   12.420  
Number of obs: 274, groups:  Site:block, 16; year, 6; Site, 4

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   43.833     10.402 146.531   4.214 4.37e-05 ***
k             84.993     52.035 228.313   1.633    0.104    
log(PS)       16.947      4.251 236.979   3.987 8.92e-05 ***
k:log(PS)     33.029     22.386 232.660   1.475    0.141    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr) k      lg(PS)
k         -0.887              
log(PS)    0.858 -0.901       
k:log(PS) -0.843  0.944 -0.971
Code 
fit.kin.Pbalance |> r.squaredGLMM()
           R2m       R2c
[1,] 0.5718185 0.7438699
Code 
fit.kin.Pexport |> r.squaredGLMM()
            R2m       R2c
[1,] 0.06433359 0.6476313
Code 
fit.kin.Yrel |> r.squaredGLMM()
            R2m       R2c
[1,] 0.02182182 0.4385486
Code 
fit.kin.Ynorm |> r.squaredGLMM()
            R2m       R2c
[1,] 0.01391747 0.3596733
Code 
# Verhalten der Modelparameter und Ertragsdaten auf P-CO2 und P-AAE10

Since we now model two measurement methods, we do not expect correlations by Site/year/etc

In [4]:
Code 
# fit.PS       <- lm(PS            ~ soil_0_20_P_CO2 + soil_0_20_P_AAE10, D)
fit.grud.PS       <- lm(log(PS)         ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10), D)
fit.grud.k        <- lm(k               ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10), D)
fit.grud.kPS      <- lm(I(log(k*PS))    ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10), D)

fit.grud.CO2.Yrel     <- lmer(Ymain_rel       ~ log(soil_0_20_P_CO2) + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
fit.grud.AAE10.Yrel     <- lmer(Ymain_rel       ~ log(soil_0_20_P_AAE10) + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.grud.Yrel     <- lmer(Ymain_rel       ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10) + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)


# this is hopeless, since cannot log becaus of 0's
fit.CO2.Pexport     <- lmer(annual_P_uptake     ~ log(soil_0_20_P_CO2) + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.AAE10.Pexport     <- lmer(annual_P_uptake ~ log(soil_0_20_P_AAE10) + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.grud.Pexport     <- lmer(annual_P_uptake     ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10) + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.CO2.Pbalance     <- lmer(annual_P_balance      ~ log(soil_0_20_P_CO2)  + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.AAE10.Pbalance     <- lmer(annual_P_balance      ~ log(soil_0_20_P_AAE10)  + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
fit.grud.Pbalance     <- lmer(annual_P_balance      ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10)  + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)
boundary (singular) fit: see help('isSingular')
Code 
#summary(glht(fit.PS))
# Fazit: PS wird von treatment stark beeinfluss, k eher nicht (dafür von site)
summary(glht(fit.grud.PS))

     Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = log(PS) ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10), 
    data = D)

Linear Hypotheses:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept) == 0            -0.61404    0.20562  -2.986  0.00521 ** 
log(soil_0_20_P_CO2) == 0    1.30106    0.05963  21.819  < 0.001 ***
log(soil_0_20_P_AAE10) == 0 -0.24892    0.05163  -4.821  < 0.001 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
Code 
summary(glht(fit.grud.k))

     Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = k ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10), 
    data = D)

Linear Hypotheses:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept) == 0             0.074843   0.026807   2.792     0.01 *  
log(soil_0_20_P_CO2) == 0   -0.031362   0.007774  -4.034   <0.001 ***
log(soil_0_20_P_AAE10) == 0  0.029663   0.006731   4.407   <0.001 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
Code 
summary(glht(fit.grud.kPS))

     Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = I(log(k * PS)) ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10), 
    data = D)

Linear Hypotheses:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept) == 0            -2.89619    0.20538 -14.102   <0.001 ***
log(soil_0_20_P_CO2) == 0    1.12637    0.05956  18.912   <0.001 ***
log(soil_0_20_P_AAE10) == 0 -0.10426    0.05157  -2.022    0.072 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
Code 
summary(fit.grud.Yrel)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Ymain_rel ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10) + (1 |  
    year) + (1 | Site) + (1 | Site:block) + (1 | Site:Treatment)
   Data: D

REML criterion at convergence: 1740.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.00124 -0.58425 -0.01625  0.57153  2.82374 

Random effects:
 Groups         Name        Variance Std.Dev.
 Site:block     (Intercept)   2.612   1.616  
 Site:Treatment (Intercept)  15.818   3.977  
 year           (Intercept) 161.127  12.694  
 Site           (Intercept)  54.936   7.412  
 Residual                   196.854  14.030  
Number of obs: 212, groups:  
Site:block, 16; Site:Treatment, 12; year, 5; Site, 4

Fixed effects:
                       Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)            108.4942    17.6735  40.7465   6.139 2.81e-07 ***
log(soil_0_20_P_CO2)     9.6918     4.2287  30.9801   2.292   0.0289 *  
log(soil_0_20_P_AAE10)  -0.9326     4.1643  46.0411  -0.224   0.8238    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) l(_0_20_P_C
l(_0_20_P_C  0.830            
l(_0_20_P_A -0.915 -0.858
Code 
summary(fit.grud.Pexport)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: annual_P_uptake ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10) +  
    (1 | year) + (1 | Site) + (1 | Site:block) + (1 | Site:Treatment)
   Data: D

REML criterion at convergence: 1844.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1941 -0.4246 -0.0380  0.3941  4.7946 

Random effects:
 Groups         Name        Variance  Std.Dev.
 Site:block     (Intercept) 0.000e+00 0.000000
 Site:Treatment (Intercept) 1.036e-06 0.001018
 year           (Intercept) 7.372e+01 8.586275
 Site           (Intercept) 2.427e+01 4.926572
 Residual                   6.686e+01 8.177020
Number of obs: 259, groups:  
Site:block, 16; Site:Treatment, 12; year, 6; Site, 4

Fixed effects:
                       Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)             24.7259     8.5160  64.7088   2.903  0.00504 **
log(soil_0_20_P_CO2)     4.5632     1.7896 205.5645   2.550  0.01151 * 
log(soil_0_20_P_AAE10)   0.7243     1.8948 175.6048   0.382  0.70275   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) l(_0_20_P_C
l(_0_20_P_C  0.802            
l(_0_20_P_A -0.859 -0.899     
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Code 
summary(fit.grud.Pbalance)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: annual_P_balance ~ log(soil_0_20_P_CO2) + log(soil_0_20_P_AAE10) +  
    (1 | year) + (1 | Site) + (1 | Site:block) + (1 | Site:Treatment)
   Data: D

REML criterion at convergence: 2129.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.3544 -0.5456  0.0619  0.5653  2.8045 

Random effects:
 Groups         Name        Variance  Std.Dev. 
 Site:block     (Intercept) 0.000e+00  0.000000
 Site:Treatment (Intercept) 4.284e+02 20.698518
 year           (Intercept) 5.796e+01  7.613248
 Site           (Intercept) 1.382e-06  0.001176
 Residual                   1.145e+02 10.701982
Number of obs: 274, groups:  
Site:block, 16; Site:Treatment, 12; year, 6; Site, 4

Fixed effects:
                       Estimate Std. Error       df t value Pr(>|t|)
(Intercept)              6.5755    15.3463 150.9193   0.428    0.669
log(soil_0_20_P_CO2)    -0.4298     4.2050 260.9577  -0.102    0.919
log(soil_0_20_P_AAE10)  -0.5459     3.4874 264.0262  -0.157    0.876

Correlation of Fixed Effects:
            (Intr) l(_0_20_P_C
l(_0_20_P_C  0.754            
l(_0_20_P_A -0.888 -0.754     
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Code 
fit.grud.Pbalance |> r.squaredGLMM()
              R2m       R2c
[1,] 0.0006758751 0.8095345
Code 
fit.grud.Pexport |> r.squaredGLMM()
            R2m       R2c
[1,] 0.06600072 0.6211875
Code 
fit.grud.Yrel |> r.squaredGLMM()
            R2m       R2c
[1,] 0.07256889 0.5767479
Code 
save.image(file = "~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/results_coefficient_analysis")
In [5]:
Code 
create_coef_table <- function(lmer_models, 
                              covariate_order = NULL, 
                              covariate_labels = NULL, # NEU: Benannter Vektor für Zeilennamen
                              model_labels = NULL) {   # NEU: Benannter Vektor für Spaltennamen

  # Extract coefficients and p-values (Ihre Originalfunktion, keine Änderung hier)
  extract_coef_info <- function(model) {
    # ... (keine Änderung, Ihr Code bleibt hier)
    coef_matrix <- summary(model)|> coef()
    estimates <- coef_matrix[, 1]
    p_values <- coef_matrix[, ncol(coef_matrix)]
    formatted_coef <- sapply(seq_along(estimates), function(i) {
      est_str <- sprintf("%.3f", estimates[i])
      stars <- if (p_values[i] < 0.001) "***" else
               if (p_values[i] < 0.01) "** " else
               if (p_values[i] < 0.05) "* " else  ""
      paste0(stars, est_str)
    })
    names(formatted_coef) <- rownames(coef_matrix)
    return(formatted_coef)
  }

  # Extract R-squared values (Ihre Originalfunktion, keine Änderung hier)
  extract_r_squared <- function(model) {
    # ... (keine Änderung, Ihr Code bleibt hier)
    r2_values <- MuMIn::r.squaredGLMM(model) # MuMIn:: hinzugefügt für Klarheit
    return(c(
      R2m = sprintf("%.3f", r2_values[1, "R2m"]),
      R2c = sprintf("%.3f", r2_values[1, "R2c"])
    ))
  }

  # Daten extrahieren (Ihr Originalcode)
  all_coefs <- lapply(lmer_models, extract_coef_info)
  all_r_squared <- lapply(lmer_models, extract_r_squared)
  all_covariate_names <- unique(unlist(lapply(all_coefs, names)))

  if (is.null(covariate_order)) {
    covariate_order <- c("(Intercept)", sort(all_covariate_names[all_covariate_names != "(Intercept)"]))
  }
  covariate_order <- covariate_order[covariate_order %in% all_covariate_names]
  final_order <- c(covariate_order, "R2m", "R2c")

  # Matrix erstellen (Ihr Originalcode)
  results_matrix <- matrix("",
                           nrow = length(final_order),
                           ncol = length(lmer_models),
                           dimnames = list(final_order, names(lmer_models)))

  # Matrix füllen (Ihr Originalcode)
  for (model_name in names(lmer_models)) {
    model_coefs <- all_coefs[[model_name]]
    for (covar in names(model_coefs)) {
      if (covar %in% covariate_order) {
        results_matrix[covar, model_name] <- model_coefs[covar]
      }
    }
    r2_values <- all_r_squared[[model_name]]
    results_matrix["R2m", model_name] <- r2_values["R2m"]
    results_matrix["R2c", model_name] <- r2_values["R2c"]
  }

  # --- NEU: Zeilen- und Spaltennamen ersetzen ---
  
  # Ersetze die Zeilennamen (Kovariaten), falls covariate_labels übergeben wurde
  if (!is.null(covariate_labels)) {
    # Finde die Übereinstimmungen in den aktuellen Zeilennamen
    row_matches <- match(rownames(results_matrix), names(covariate_labels))
    # Ersetze nur die, die gefunden wurden
    new_rownames <- rownames(results_matrix)
    new_rownames[!is.na(row_matches)] <- covariate_labels[row_matches[!is.na(row_matches)]]
    rownames(results_matrix) <- new_rownames
  }
  
  # Ersetze die Spaltennamen (Modelle), falls model_labels übergeben wurde
  if (!is.null(model_labels)) {
    col_matches <- match(colnames(results_matrix), names(model_labels))
    new_colnames <- colnames(results_matrix)
    new_colnames[!is.na(col_matches)] <- model_labels[col_matches[!is.na(col_matches)]]
    colnames(results_matrix) <- new_colnames
  }
  
  # --- Ende der neuen Sektion ---

  # Convert to data frame for kable
  results_df <- data.frame(Covariate = rownames(results_matrix),
                           results_matrix,
                           check.names = FALSE, # Verhindert, dass R Spaltennamen ändert
                           stringsAsFactors = FALSE)
  
  results_df
}
In [6]:
Code 
lmer_models <- list(
  PS = fit.soil.PS,
  k = fit.soil.k,
  'log(k*PS)' = fit.soil.kPS,
  'k*log(PS)' = fit.soil.kPS2,
  CO2 = fit.soil.CO2,
  AAE10 = fit.soil.AAE10
)



coef_table_soil <- create_coef_table(lmer_models)
kable(coef_table_soil,
row.names = FALSE,
align = c("l", rep("r", ncol(coef_table_soil) - 1)),
escape = FALSE,
caption = "Coefficient Table for Soil covariates. 
Significant codes:  0 '\\*\\*\\*' 0.001 '\\*\\*' 0.01 '\\*' 0.05")
Coefficient Table for Soil covariates. Significant codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05
Covariate PS k log(k*PS) k*log(PS) CO2 AAE10
(Intercept) ***21.444 0.454 ***21.189 1.673 ** 14.014 ** 23.126
Ald ***-8.706 -0.072 * -8.631 -1.136 -4.417 * -9.473
Fed ***-1.068 0.005 -1.084 * -0.200 * -0.845 -0.606
soil_0_20_clay -0.006 ** -0.016 -0.085 * 0.054 0.015 -0.029
soil_0_20_Corg * 0.612 ** 0.137 * 1.250 * -0.342 0.269 ***1.454
soil_0_20_pH_H2O -0.018 -0.021 -0.092 0.063 0.124 0.012
soil_0_20_silt -0.000 0.004 0.008 0.011 -0.015 * -0.049
R2m 0.393 0.212 0.368 0.652 0.362 0.494
R2c 0.996 0.915 0.993 0.819 0.995 0.997
Code 
lmer_models_yield <- list(
  "Yn-STP-CO2" = fit.grud.CO2.Ynorm,
  "Yn-STP-AAE10" = fit.grud.AAE10.Ynorm,
  "Yn-STP-GRUD" = fit.grud.CO2.AAE10.Ynorm,
  "Yn-Kinetic" = fit.kin.Ynorm,
  "Yr-STP-CO2" = fit.grud.CO2.Yrel,
  "Yr-STP-AAE10" = fit.grud.AAE10.Yrel,
  "Yr-STP-GRUD" = fit.grud.Yrel,
  "Yr-Kinetic" = fit.kin.Yrel
)

lmer_models_balance <- list(
  CO2_Pbalance = fit.CO2.Pbalance,
  AAE10_Pbalance = fit.AAE10.Pbalance,
  Grud_Pbalance = fit.grud.Pbalance,
  Kin_Pbalance = fit.kin.Pbalance
)

lmer_models_export <- list(
  CO2_Pexport = fit.CO2.Pexport,
  AAE10_Pexport = fit.AAE10.Pexport,
  Grud_Pexport = fit.grud.Pexport,
  Kin_Pexport = fit.kin.Pexport
)


coef_table_yield <- create_coef_table(lmer_models_yield)
kable(coef_table_yield,
row.names = FALSE,
escape = FALSE,
align = c("l", rep("r", ncol(coef_table_yield) - 1)),
caption = "Coefficient Table for Ynorm and Yrel. 
Significant codes:  0 '\\*\\*\\*' 0.001 '\\*\\*' 0.01 '\\*' 0.05")
Coefficient Table for Ynorm and Yrel. Significant codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05
Covariate Yn-STP-CO2 Yn-STP-AAE10 Yn-STP-GRUD Yn-Kinetic Yr-STP-CO2 Yr-STP-AAE10 Yr-STP-GRUD Yr-Kinetic
(Intercept) ***1.059 ***0.532 ***1.096 0.980 ***104.862 ***75.343 ***108.494 56.375
k 2.262 ** 377.498
k:log(PS) 0.931 ** 171.507
log(PS) -0.063 * -27.486
log(soil_0_20_P_AAE10) ***0.120 -0.006 ** 7.111 -0.933
log(soil_0_20_P_CO2) ***0.162 0.137 ** 8.853 * 9.692
log(soil_0_20_P_CO2):log(soil_0_20_P_AAE10) 0.016
R2m 0.218 0.198 0.220 0.014 0.074 0.063 0.073 0.022
R2c 0.358 0.474 0.365 0.360 0.569 0.537 0.577 0.439
Code 
# Step 1: Prepare the data for correlation analysis
# Create a new dataframe containing only the necessary columns and remove rows with NAs
correlation_data <- D %>%
  select(uid, PS, k,E_mod_10080, E_exp_10080,E_exp_1440, n_1440) %>%
  na.omit()

# --- Capacity Comparison (PS vs. E-value) ---

# Step 2: Visualize the relationship
ggplot(correlation_data, aes(x = E_mod_10080, y = PS)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +
  labs(
    title = TeX("Relationship between Desorbable P and Exchangeable P"),
    x = TeX("Isotopically Exchangeable P, $E_{10080}$ (mg kg$^{-1}$)"),
    y = TeX("Desorbable P, $P_{desorb}$ (mg L$^{-1}$)")
  ) +
  theme_bw()
`geom_smooth()` using formula = 'y ~ x'

Code 
# Step 3: Perform the formal correlation test
# Using Spearman's rank correlation is a robust choice
spearman_capacity <- cor.test(correlation_data$PS, correlation_data$E_exp_10080, method = "spearman")
print("Spearman correlation between PS and E_exp_10080:")
[1] "Spearman correlation between PS and E_exp_10080:"
Code 
print(spearman_capacity)

    Spearman's rank correlation rho

data:  correlation_data$PS and correlation_data$E_exp_10080
S = 63135, p-value = 0.0001128
alternative hypothesis: true rho is not equal to 0
sample estimates:
     rho 
0.404355
Code 
# --- Kinetic Comparison (k vs. n-value) ---

# Step 2: Visualize the relationship
ggplot(correlation_data, aes(x = n_1440, y = k)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +
  labs(
    title = TeX("Relationship between Desorption Rate and Exchange Rate"),
    x = TeX("IEK Kinetic Parameter, $n_{60}$"),
    y = TeX("Desorption Rate Constant, $k$ (min$^{-1}$)")
  ) +
  theme_bw()
`geom_smooth()` using formula = 'y ~ x'

Code 
# Step 3: Perform the formal correlation test
spearman_kinetic <- cor.test(correlation_data$k, correlation_data$n_1440, method = "spearman")
print("Spearman correlation between k and n_60:")
[1] "Spearman correlation between k and n_60:"
Code 
print(spearman_kinetic)

    Spearman's rank correlation rho

data:  correlation_data$k and correlation_data$n_1440
S = 67681, p-value = 0.0006275
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.3614673
Code 
save.image(file = "~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/results_coefficient_analysis")