---
title: "Kinetic Model Coefficient Analysis"
format: 
  html:
    math: mathjax
  pdf: default
  docx: default
author: Marc Pérez
date: 2025-05-22
---

```{r setup}
#| include: false
#| message: false
#| warning: false


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



```{r}
#| eval: false
#| include: false






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 



```{r}
#|code-fold: true
#|tidy: true


# Wovon hängen Modelparameter ab?

library(lmerTest)

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





r.squaredGLMM(fit.soil.kPS)
r.squaredGLMM(fit.soil.kPS2)

fit.kin.Yrel     <- lmer(Ymain_rel         ~ k * log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D)

fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k + log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
r.squaredGLMM(fit.kin.Ynorm)


fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
r.squaredGLMM(fit.kin.Ynorm)

fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS)  + (1|year), D, subset = Treatment != "P166")
r.squaredGLMM(fit.kin.Ynorm)


fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D)
r.squaredGLMM(fit.kin.Ynorm)

fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS)  + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment == "P0")
r.squaredGLMM(fit.kin.Ynorm)
summary(fit.kin.Ynorm, ddf="Kenward-Roger")

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)




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)
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)
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)
# 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")
r.squaredGLMM(fit.kin.Ynorm)

fit.kin.Ynorm    <- lmer(Ymain_norm       ~ k * log(PS) + (1|year) + (1|Site)  + (1|Site:block), D, subset = Treatment != "P166")
fit.kin.Pexport  <- lmer(annual_P_uptake  ~ k * log(PS) + (1|year) + (1|Site)  + (1|Site:block), D)
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))
car::vif(lm(Ymain_norm       ~ scale(k) * scale(log(PS))   + crop, D))
car::vif(lm(Ymain_norm       ~ k + log(PS)   + crop, D))
car::vif(lm(Ymain_norm       ~ I(k * (PS)) + k + log(PS)    + crop, D))

with(D, hist(I(exp(k) * (PS))))
with(D, hist(log(PS)))

r.squaredGLMM(fit.kin.Ynorm)

anova(fit.soil.PS)
summary(glht(fit.soil.PS))
# Fazit: PS wird von treatment stark beeinfluss, k eher nicht (dafür von site)
anova(fit.soil.k)
summary(glht(fit.soil.k))
anova(fit.soil.kPS)
summary(glht(fit.soil.kPS))
anova(fit.kin.Yrel)
anova(fit.kin.Ynorm)
anova(fit.kin.Pexport)
anova(fit.soil.CO2)
anova(fit.soil.AAE10)
anova(fit.kin.Pbalance)
summary(fit.kin.Pbalance)
fit.kin.Pbalance |> r.squaredGLMM()
fit.kin.Pexport |> r.squaredGLMM()
fit.kin.Yrel |> r.squaredGLMM()
fit.kin.Ynorm |> r.squaredGLMM()
# 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
```{r}
#|code-fold: true
#|tidy: true
# 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)
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)

fit.AAE10.Pexport     <- lmer(annual_P_uptake ~ log(soil_0_20_P_AAE10) + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)

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)


fit.CO2.Pbalance     <- lmer(annual_P_balance      ~ log(soil_0_20_P_CO2)  + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)


fit.AAE10.Pbalance     <- lmer(annual_P_balance      ~ log(soil_0_20_P_AAE10)  + (1|year) + (1|Site)  + (1|Site:block) + (1|Site:Treatment), D)


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)



#summary(glht(fit.PS))
# Fazit: PS wird von treatment stark beeinfluss, k eher nicht (dafür von site)
summary(glht(fit.grud.PS))
summary(glht(fit.grud.k))
summary(glht(fit.grud.kPS))
summary(fit.grud.Yrel)
summary(fit.grud.Pexport)
summary(fit.grud.Pbalance)
fit.grud.Pbalance |> r.squaredGLMM()
fit.grud.Pexport |> r.squaredGLMM()
fit.grud.Yrel |> r.squaredGLMM()


save.image(file = "~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/results_coefficient_analysis")

```



```{r}
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
}
```


```{r}
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")




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



# 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()

# 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:")
print(spearman_capacity)


# --- 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()

# 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:")
print(spearman_kinetic)







save.image(file = "~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/results_coefficient_analysis")


```

