---
title: "Model Validation and coefficient calculation"
format: 
  html:
    math: mathjax
pdf: default
docx: default
author: Marc Pérez
date: 2025-05-22
---
  
  
  

```{r Loading previous results}
#|include: false
#|echo: false

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)
library(sjPlot) # table functions
library(sjmisc) # sample data
library(lme4) # fitting models
library(report)
library(performance)

load("~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/results_coefficient_analysis")
```

<!-- ## Struktur -->

<!-- ## Warum ist die Arbeit wichtig -->

<!-- P ist endlich, Umweltprobleme \## P ist sehr komplex -->

<!-- Siehe oben \## Wie wird bisher P-Ernährung angegangen -->

<!-- GRUD -->

<!-- ## Warum Bodentest -->

<!-- ## Warum kinetischer Bodentest -->

# Results

The results are presented in three main parts. First, the development and validation of the phosphorus (P) desorption kinetic model are detailed, justifying the final modeling approach. Second, the descriptive trends of both agronomic outcomes and soil P parameters in response to long-term fertilization and site differences are explored visually. Finally, the predictive power of the kinetic and standard P parameters is formally evaluated using linear mixed-effects models.

## Establishment of the P-Desorption Kinetic Model

The primary goal was to derive two key parameters for each soil sample: the desorbable P pool ($P_{desorb}$) and the rate constant ($k$). The analysis proceeded in two stages: an initial test of a linearized model, followed by the implementation of a more robust non-linear model.

### Initial Approach: Linearized Model

Following the conceptual framework of Flossmann and Richter (1982), the first-order kinetic equation was linearized to the form $ln(1 - P/P_{desorb}) = -kt$. For this approach, the asymptote ($P_{desorb}$) was not fitted but was calculated beforehand as the difference between Olsen-P and water-soluble P. A linear model was then fitted to the transformed data for each sample.

```{r}
#| label: fig-linearized-model
#| fig-cap: "Test of the linearized first-order kinetic model. Data is faceted by Site and P fertilization Treatment. The solid line represents a linear model fit."
#| echo: false
#| warning: false
#| message: false

# This is a placeholder for your code from pretest.qmd
 ggplot(d[d$Treatment!="P100"&(d$Repetition==1|d$Repetition==2),], aes(y=Y1, x=t.min., col = Repetition)) +
   geom_point() +
   facet_grid(Site ~ Treatment) +
   labs(x=TeX("$Time (min)$"),
        y=TeX("$ln(1-\\frac{P}{P_S})$")) +
   geom_smooth(method="lm", alpha = 0.3)
```

*(Placeholder for the linearized model plot)*

The results of this linearized approach were inconsistent (Figure X.1). While some samples showed a reasonable linear trend, many exhibited significant deviations, particularly a curvilinear pattern. Furthermore, the model required the intercept to be fixed at zero, but the statistical summary of the fitted models (`lmList`) showed that the estimated intercepts were often significantly different from zero (p < 0.05). For example, the sample `Cadenazzo_P166_1` had a highly significant intercept of -0.27 (p < 0.001). This systematic deviation from the model's core assumption indicated that this linearized approach, which relied on an externally estimated asymptote, was not robust for this dataset.

### Final Approach: Non-Linear Model

Given the limitations of the linearized method, a direct non-linear modeling approach was adopted to estimate both $P_{desorb}$ and $k$ simultaneously from the untransformed data. The model $P(t) = P_{desorb} \times (1 - e^{-k \times t'})$ was fitted to each sample's desorption curve.

```{r}
#| label: fig-nonlinear-model
#| fig-cap: "Non-linear first-order kinetic model fits for P desorption over time. Points represent measured data and solid lines represent the fitted model for each replicate."
#| echo: false
#| warning: false
#| message: false

Res <- nlsList(Pv.mg.L. ~ PS * (1 - exp(-k * (t.dt))) | uid, d[, c("Pv.mg.L.", "uid", "t.dt")],  start=list(PS=0.1,k=0.2))
# summary(Res)
# d$nls_pred <- predict(Res)

# Extract coefficients from the nlsList results
nls_coefs <- coef(Res)
nls_coefs$uid <- rownames(nls_coefs)

# Merge coefficients back to the main dataset
d_plot <- merge(d, nls_coefs, by = "uid")

# Most straightforward approach - create curves manually
time_seq <- seq(min(d$t.dt, na.rm = TRUE), max(d$t.dt, na.rm = TRUE), length.out = 100)

# Create prediction data
pred_data <- d_plot %>%
  select(uid, Site, Treatment, Repetition, PS, k) %>%
  distinct() %>%
  crossing(t.dt = time_seq) %>%
  mutate(pred_Pv = PS * (1 - exp(-k * (t.dt))))

# This is a placeholder for your code from pretest.qmd
 p1 <- ggplot() +
   geom_point(data = d_plot, aes(y = Pv.mg.L., x = t.dt, col = Repetition)) +
   geom_line(data = pred_data, aes(x = t.dt, y = pred_Pv, col = Repetition), size = 0.5) +
   facet_grid(Treatment ~ Site,scales = "free_y") +
   labs(x = TeX("$Time (min)$"),
        y = TeX("$P_{V}(\\frac{mg}{L})$"));
 suppressWarnings(print(p1))
```

*(Placeholder for the non-linear model plot)*

This approach proved to be far more successful (Figure X.2). The non-linear model was able to accurately capture the curvilinear shape of the desorption data for nearly all samples. The successful convergence of the models for each unique sample (`uid`) provided robust, sample-specific estimates for both the desorbable P pool ($P_{desorb}$) and the rate constant ($k$). For example, for the sample `Cadenazzo_P166_1`, the model estimated a $P_{desorb}$ of 0.27 mg/L and a $k$ of 0.086 min⁻¹, both with high statistical significance (p < 0.01).

To account for the hierarchical data structure and obtain the most reliable estimates, the final parameters were extracted from a non-linear mixed-effects model (`nlme`), which models the overall fixed effects for $k$ and $P_{desorb}$ while allowing for random variation among individual samples. These final `nlme`-derived coefficients were used for all subsequent analyses.













## Research Questions:

### How well can current GRUD measurements of $C_P$ predict the relative Yield, P-Uptake and P-Balance?

-   Hypothesis I: The measurements of the equlibrium concentrations of Phosphorus in a solvent do not display significant effects on relative Yield and consequently P-Uptake, since it is strongly dependent on yield. $C_P$ relates strongly to the amount of Phosphorus applied, the P-balance might well be siginificantly correlated to $C_P$ but not explain a lot of variance.


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

ggplot(D[D$soil_0_20_P_CO2!=0,],aes(y=Ymain_norm, 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="relative Yield (%)") +
  facet_wrap( ~ Site, nrow = 2)

ggplot(D[D$soil_0_20_P_AAE10!=0,],aes(y=Ymain_norm, x=soil_0_20_P_AAE10, col=Site, size = Treatment)) +
  geom_point(shape = 7) + 
  scale_x_log10() + scale_y_log10() +
  labs(x=TeX("$P_{AAE10}(mg/kg Soil)$"),
         y="relative Yield (%)") +
  facet_wrap( ~ Site, nrow = 2)

ggplot(D[D$soil_0_20_P_CO2!=0&D$annual_P_uptake!=0,],aes(y=annual_P_uptake, 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("Annual P-Uptake $kg~P/ha$")) +
  facet_wrap( ~ Site, nrow = 2)

ggplot(D[D$soil_0_20_P_AAE10!=0&D$annual_P_uptake!=0,],aes(y=annual_P_uptake, x=soil_0_20_P_AAE10, col=Site, size = Treatment)) +
  geom_point(shape = 7) + 
  scale_x_log10() + scale_y_log10() +
  labs(x=TeX("$P_{AAE10}(mg/kg Soil)$"),
         y=TeX("Annual P-Uptake $kg~P/ha$")) +
  facet_wrap( ~ Site, nrow = 2)

ggplot(D[D$soil_0_20_P_CO2!=0&D$annual_P_balance!=0,],aes(y=annual_P_balance, 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("Annual P-Balance $kg~P/ha$")) +
  facet_wrap( ~ Site, nrow = 2)

ggplot(D[D$soil_0_20_P_AAE10!=0&D$annual_P_balance!=0,],aes(y=annual_P_balance, x=soil_0_20_P_AAE10, col=Site, size = Treatment)) +
  geom_point(shape = 7) + 
  scale_x_log10() + scale_y_log10() +
  labs(x=TeX("$P_{AAE10}(mg/kg Soil)$"),
         y=TeX("Annual P-Balance $kg~P/ha$")) +
  facet_wrap( ~ Site, nrow = 2)



```
Now we want to check the strength of the models in terms of $R^2$ and the significance of the effects in terms of p-values:

```{r}
#|code-fold: true
#|tidy: true
#|echo: false
#|output: asis

#tab_model(fit.grud.Yrel,fit.grud.Puptake,fit.grud.Pbalance)
report(fit.grud.Yrel)
report(fit.grud.Pexport)
report(fit.grud.Pbalance)
```
here I also show the non linear mixed models, following the Mitscherlich saturation curve:

```{r}
#|code-fold: true
#|echo: false


library(nlme)

# Make sure grouping variables are factors
D$year  <- as.factor(D$year)
D$Site  <- as.factor(D$Site)
D$block <- as.factor(D$block)
D$crop <- as.factor(D$crop)
# Fit the model
fit.mitscherlich.CO2.Yrel <- nlme(
  Ymain_rel ~ A * (1 - exp(-k * soil_0_20_P_CO2 + E)),                
  fixed = A + k + E ~ soil_0_20_clay + soil_0_20_pH_H2O + ansum_sun + ansum_prec,
  random = A ~ 1 | year/Site/block,
  data = D,
  subset = Treatment != "P166",
  start = c(
    A = 1100, A1 = 0, A2 = 0, A3 = 0, A4 = 0,
    k = 0.05, k1 = 0, k2 = 0, k3 = 0, k4 = 0,
    E = -3, E1 = 0, E2 = 0, E3 = 0, E4 = 0
  ),
  control = nlmeControl(maxIter = 500),
  na.action = na.omit
)

summary(fit.mitscherlich.CO2.Yrel)
anova(fit.mitscherlich.CO2.Yrel)

model_performance(fit.mitscherlich.CO2.Yrel)

r.square.CO2 <- 1-sum(residuals(fit.mitscherlich.CO2.Yrel)^2)/sum((D$Ymain_rel-mean(D$Ymain_rel,na.rm=TRUE))^2,na.rm = TRUE)

# Fit the model
fit.mitscherlich.CO2.Yrel <- nlme(
  Ymain_rel ~ A * (1 - exp(-k * kPS + E)),                
  fixed = A + k + E ~ soil_0_20_clay + soil_0_20_pH_H2O + ansum_sun + ansum_prec,
  random = A ~ 1 | year/Site/block,
  data = D,
  subset = Treatment != "P166",
  start = c(
    A = 600, A1 = 0, A2 = 0, A3 = 0, A4 = 0,
    k = 3, k1 = 0, k2 = 0, k3 = 0, k4 = 0,
    E = -3, E1 = 0, E2 = 0, E3 = 0, E4 = 0
  ),
  control = nlmeControl(maxIter = 1500),
  na.action = na.omit
)

summary(fit.mitscherlich.CO2.Yrel)
anova(fit.mitscherlich.CO2.Yrel)

model_performance(fit.mitscherlich.CO2.Yrel)

r.square.CO2 <- 1-sum(residuals(fit.mitscherlich.CO2.Yrel)^2)/sum((D$Ymain_rel-mean(D$Ymain_rel,na.rm=TRUE))^2,na.rm = TRUE)



fit.mitscherlich.CO2.Ynorm <- nlme(
  Ymain_norm ~ A * (1 - exp(-k * kPS + E)),                
  fixed = A + k + E ~ soil_0_20_clay + soil_0_20_pH_H2O + ansum_sun + ansum_prec,
  random = A ~ 1 | year/Site/block,
  data = D,
  subset = Treatment != "P166",
  start = c(
    A = 6, A1 = 0, A2 = 0, A3 = 0, A4 = 0,
    k = 3, k1 = 0, k2 = 0, k3 = 0, k4 = 0,
    E = -3, E1 = 0, E2 = 0, E3 = 0, E4 = 0
  ),
  control = nlmeControl(maxIter = 1500),
  na.action = na.omit
)

summary(fit.mitscherlich.CO2.Yrel)
anova(fit.mitscherlich.CO2.Yrel)

model_performance(fit.mitscherlich.CO2.Yrel)

r.square.CO2 <- 1-sum(residuals(fit.mitscherlich.CO2.Yrel)^2)/sum((D$Ymain_rel-mean(D$Ymain_rel,na.rm=TRUE))^2,na.rm = TRUE)

# Fit the model
fit.mitscherlich.CO2.Yrel <- nlme(
  Ymain_norm ~ A * (1 - exp(-k * soil_0_20_P_CO2 + E)),                
  fixed = A + k + E ~ soil_0_20_clay + soil_0_20_pH_H2O + ansum_sun + ansum_prec,
  random = A ~ 1 | year/Site/block,
  data = D,
  subset = Treatment != "P166",
  start = c(
    A = 1100, A1 = 0, A2 = 0, A3 = 0, A4 = 0,
    k = 0.05, k1 = 0, k2 = 0, k3 = 0, k4 = 0,
    E = -3, E1 = 0, E2 = 0, E3 = 0, E4 = 0
  ),
  control = nlmeControl(maxIter = 500),
  na.action = na.omit
)

summary(fit.mitscherlich.CO2.Yrel)
anova(fit.mitscherlich.CO2.Yrel)

model_performance(fit.mitscherlich.CO2.Yrel)

r.square.CO2 <- 1-sum(residuals(fit.mitscherlich.CO2.Yrel)^2)/sum((D$Ymain_rel-mean(D$Ymain_rel,na.rm=TRUE))^2,na.rm = TRUE)




```

With the covariate and random effect used as by Juliane Hirte we obtain  $R^2=$ `{r} r.square.CO2`, I don't know how to interpret that, I fear that the model is overfitting data.


### How do GRUD-measurements of $C_P$ relate to the soil properties $C_\text{org}$-content, clay-content, silt-content and pH?

- Hypothesis II:
Given the known capacity of clay and silt compounds to adsorb orthophosphate a positive correlation between $C_P$ (for both $CO_2$ and AAE10) and silt- and clay-content. $C_\text{org}$ has been reported to positively influence the capacity of Phosphorus as well, it is plausible it also shows a positive correlation with $C_P$. AAE10 also deploys $Na_4EDTA$ which is easily captured by $Mg^{2+}$ and $Ca^{2+}$, therefore it is officially by GRUD advised against being used in soils with $\text{pH}>6.8$, therefore $C_P$-AAE10 will presumably be negatively correlated to pH.

```{r}
#|code-fold: true
#|echo: false


anova(fit.soil.CO2)
fit.soil.CO2 |> r.squaredGLMM()
anova(fit.soil.AAE10)
fit.soil.AAE10 |> r.squaredGLMM()
```

### Can the Inclusion of the net-release-kinetic of Orthophosphate improve the model power of predicting relative Yield, P-Uptake and P-Balance?

- Hypothesis III:
Given the comparably low solubility of $PO_4^{3-}$ in the water-soil interface,  most P is transported to the rhizosphere via diffusion. As a consequence the intensity of $PO_4^{3-}$ might not adequately account for the P-uptake in the harvested plant. Since the diffusion process is in its velocity a kinetic and in its finally reached intensity a thermodynamic process, the inclusion of kinetic parameters might well improve the performance. 

```{r}
#|code-fold: true
#|echo: false

fit.mitscherlich.PS.Yrel <- nlme(
  Ymain_rel ~ A * (1 - exp(-r * PS + E)),                
  fixed = A + r + E ~ k + soil_0_20_clay + soil_0_20_pH_H2O + ansum_sun + ansum_prec,
  random = A ~ 1 | year/Site/block,
  data = D,
  start = c(
    A = 220, A1 = 0, A2 = 0, A3 = 0, A4 = 0, A5 = 0,
    r = 1, r1 = 0, r2 = 0, r3 = 0, r4 = 0, r5 = 0,
    E = -1, E1 = 0, E2 = 0, E3 = 0, E4 = 0, E5 = 0
  ),
  control = nlmeControl(maxIter = 500),
  na.action = na.omit
)

summary(fit.mitscherlich.PS.Yrel)
anova(fit.mitscherlich.PS.Yrel)

model_performance(fit.mitscherlich.PS.Yrel)

r.square.PS <- 1-sum(residuals(fit.mitscherlich.PS.Yrel)^2)/sum((D$Ymain_rel-mean(D$Ymain_rel,na.rm=TRUE))^2,na.rm = TRUE)

fit.mitscherlich.kPS.Yrel <- nlme(
  Ymain_rel ~ A * (1 - exp(-k * PS + E)),                
  fixed = A + E ~ soil_0_20_clay + soil_0_20_pH_H2O + ansum_sun + ansum_prec,
  random = A ~ 1 | year/Site/block,
  data = D,
  start = c(
    A = 220, A1 = 0, A2 = 0, A3 = 0, A4 = 0,
    E = -1, E1 = 0, E2 = 0, E3 = 0, E4 = 0
  ),
  control = nlmeControl(maxIter = 500),
  na.action = na.omit
)

summary(fit.mitscherlich.kPS.Yrel)
anova(fit.mitscherlich.kPS.Yrel)

model_performance(fit.mitscherlich.kPS.Yrel)

r.square.kPS <- 1-sum(residuals(fit.mitscherlich.kPS.Yrel)^2)/sum((D$Ymain_rel-mean(D$Ymain_rel,na.rm=TRUE))^2,na.rm = TRUE)


```

With the covariate and random effect used as by Juliane Hirte we obtain  $R^2=$ `{r} r.square.PS`, I don't know how to interpret that, I fear that the model is overfitting data, the same might be true for the model that used $k\times PS$ as a predictor with $R^2=$ `{r} r.square.kPS`.


I also tried more conservative models, where I log-transformed the concentrations and PS, also I was more cautious with random effects. This resulted in coefficients that were not as straight-forward as the mitscherlich coefficients to interpret.



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

# relative Yield
anova(fit.kin.Yrel)
summary(fit.kin.Yrel)
fit.kin.Yrel |> r.squaredGLMM()

# P-Uptake
anova(fit.kin.Pexport)
summary(fit.kin.Pexport)
fit.kin.Pexport |> r.squaredGLMM()


anova(fit.kin.Pbalance)
summary(fit.kin.Pbalance)
fit.kin.Pbalance |> r.squaredGLMM()





```


### Are the kinetic coefficients $k$ and $PS$($k$ can be interpreted as the relative speed of desorption, $PS$ is the equilibrium concentration of $PO_4^{3-}$ of the observed desorption in the dried fine earth-water suspension 1:20 by weight) related to soil properties?

- Hypothesis IV:
Clay particles as well as organic compounds with negative surface charges provide surfaces for P-sorption, especially their structure, but in general their respective concentration in a soil can be expected to significantly influence the kinetic and thermodynamic of the P-desorption reaction. The $pH$ dictates the form of orthophosphate, with $pH<6.5$, the predominant form will be $H_2PO_4^-$, this should reduce electrical interactions and increase the movement- and therefore diffusion-speed.


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

anova(fit.soil.PS)
summary(glht(fit.soil.PS))
#summary(glht(fit.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))
```


### Is the method presented by Flossmann and Richter (1982) with the double extraction replicable with the soils from the STYCS-trial?

- Hypothesis V:
The authors expect the desorption kinetics to follow a 1. order kinetic, with the relation:
$$ \frac{dP}{dt}=PS(1-e^{-kt})$$
where $PS$ is estimated as $PS=[P_\text{Olsen/CAL}]-[P_{H_2O}]$, denoted as the semi-labile P-pool. The Olsen- and CAL-method deploy extractants that increase the solubility by more than order of magnitude. This presents the problem, that the estimation of $PS$ is likely to high. It was chosen by the authors in order to make the equation linearizable, so if the linearization is not well behaved, a non-linear regression might deliver a better estimation of both parameters.

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

res <- lmList(Y1 ~ t.min. | uid, d[d$Repetition==1|d$Repetition==2,],na.action = na.pass)
summary(res)
ggplot(d, aes(y=Y1, x=t.min., col = Repetition)) +
  geom_point() +
  facet_grid(Site ~ Treatment) + 
  labs(x=TeX("$Time (min)$"),
       y=TeX("$ln(1-\\frac{P}{P_S})$")) +
  geom_smooth(method="lm", alpha = 0.3) 
```

The relation can be improved:

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

Res <- nlsList(Pv.mg.L. ~ PS * (1 - exp(-k * (t.dt))) | uid, d[, c("Pv.mg.L.", "uid", "t.dt")],  start=list(PS=0.1,k=0.2))
# summary(Res)
# d$nls_pred <- predict(Res)

# Extract coefficients from the nlsList results
nls_coefs <- coef(Res)
nls_coefs$uid <- rownames(nls_coefs)

# Merge coefficients back to the main dataset
d_plot <- merge(d, nls_coefs, by = "uid")

# Most straightforward approach - create curves manually
time_seq <- seq(min(d$t.dt, na.rm = TRUE), max(d$t.dt, na.rm = TRUE), length.out = 100)

# Create prediction data
pred_data <- d_plot %>%
  select(uid, Site, Treatment, Repetition, PS, k) %>%
  distinct() %>%
  crossing(t.dt = time_seq) %>%
  mutate(pred_Pv = PS * (1 - exp(-k * (t.dt))))

# Final plot
p1 <- ggplot() +
  geom_point(data = d_plot, aes(y = Pv.mg.L., x = t.dt, col = Repetition)) +
  geom_line(data = pred_data, aes(x = t.dt, y = pred_Pv, col = Repetition), size = 0.5) +
  facet_grid(Treatment ~ Site) +
  labs(x = TeX("$Time (min)$"),
       y = TeX("$P_{V}(\\frac{mg}{L})$")); suppressWarnings(print(p1))

```

Now we see how those parameters depend on the tratment:

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

d$ui <- interaction(d$Site, d$Treatment)

nlme.coef.avg <- list()
nlme.coef <- list()
for (lvl in levels(d$ui)){
  d.tmp <- subset(d, ui == lvl)
  # first get nlsList coefs for comparison only (unused)
  temp_nls <- coef(nlsList(Pv.mg.L. ~ PS * (1 - exp(-k * t.dt)) | uid, 
                    d.tmp[, c("Pv.mg.L.", "uid", "t.dt")], 
                    start = list(PS = 0.1, k = 0.2)))
  nlsList_coefs <- c(apply(temp_nls, 2, \(x) c(mean=mean(x), sd=sd(x))))
  names(nlsList_coefs) <- c("PS.mean", "PS.sd", "k.mean", "k.sd")

  # now do the real thing
  model4 <- nlme(Pv.mg.L. ~ PS * (1 - exp(-k * t.dt)),
                fixed = PS + k ~ 1,
                random = PS + k  ~ 1 | uid,
                data = d.tmp[, c("Pv.mg.L.", "uid", "t.dt")],
                start = c(PS = 0.05, k = 0.12),
                control = nlmeControl(maxIter = 200))
  coef(model4)
  fixef <- model4$coefficients$fixed
  ranefs <- ranef(model4)
  colnames(ranefs) <- paste0("ranef_",colnames(ranefs))
  nlme.coef[[lvl]]  <- cbind(coef(model4), ranefs, Rep=1:nrow(ranef(model4)), ui=lvl, Site=d.tmp[1, "Site"], Treatment=d.tmp[1, "Treatment"], uid = rownames(coef(model4)))
  nlme.coef.avg[[lvl]] <- data.frame(PS=fixef["PS"], k=fixef["k"], ui=lvl, Site=d.tmp[1, "Site"], Treatment=d.tmp[1, "Treatment"], uid = d.tmp$uid)
}

nlme.coef.avg <- do.call(rbind, nlme.coef.avg)
# folgendes datenset wollen wir benutzen um ihn mit dem Boden zu kombinieren
nlme.coef <- do.call(rbind, nlme.coef)
points <- geom_point(position=position_dodge(width=0.5), size = 3, alpha = 0.5)

ggplot(nlme.coef, aes(y=PS  , x=Treatment, col=Site, pch=Treatment)) + points + scale_y_log10()
ggplot(nlme.coef, aes(y=k   , x=Treatment, col=Site, pch=Treatment)) + points
ggplot(nlme.coef, aes(y=k*PS, x=Treatment, col=Site, pch=Treatment)) + points + scale_y_log10()


ggplot(nlme.coef, aes(y=PS  , x=Site, col=Treatment)) + points + scale_y_log10()
ggplot(nlme.coef, aes(y=k   , x=Site, col=Treatment)) + points
ggplot(nlme.coef, aes(y=k*PS, x=Site, col=Treatment)) + points + scale_y_log10()

# k PS macht von der interpretation her Sinn
# aber PS ist log-normal verteilt


fit.PS   <- lm(log(PS)      ~ Treatment + Site, nlme.coef)
fit.k    <- lm(k            ~ Treatment + Site, nlme.coef)
fit.kPS  <- lm(I(log(k*PS)) ~ Treatment + Site, nlme.coef)


Anova(fit.PS)
summary(glht(fit.PS, mcp(Treatment = "Tukey")))
# Fazit: PS wird von treatment stark beeinfluss, k eher nicht (dafür von site)
Anova(fit.k)
summary(glht(fit.k, mcp(Treatment = "Tukey")))

Anova(fit.kPS)
summary(glht(fit.kPS, mcp(Treatment = "Tukey")))

anova(fit.soil.PS)
#summary(glht(fit.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))
```





```{r}
library(knitr)
library(kableExtra)

benchmark_results <- readRDS("cache/benchmark-tables.rds")

# Extract individual tables
Tables_yield <- benchmark_results$Tables_yield
Tables_yield_weather <- benchmark_results$Tables_yield_weather
Tables_earth_treatment <- benchmark_results$Tables_earth_treatment
Tables_earth_notreatment <- benchmark_results$Tables_earth_notreatment
Tables_yield_lm <- benchmark_results$Tables_yield_lm
Tables_yield_weather_lm <- benchmark_results$Tables_yield_weather_lm
Tables_earth_lm <- benchmark_results$Tables_earth_lm
```

```{r}
kable(Tables_yield_lm, 
     caption = "Linear model performance for predicting yield and P-balance variables using different P-dynamics variable sets (without weather data). Rows represent different predictor variable sets: 'k' uses only the release rate constant; 'PS' uses only log-transformed semi-labile P; 'kPS' uses both k and log-transformed PS plus their interaction; 'AAE10' uses only log-transformed AAE10-extractable P; 'CO2' uses only log-transformed CO2-extractable P; 'AAE10_CO2' uses both log-transformed AAE10 and CO2 extractable P plus their log-log interaction; 'AAE10_CO2_kPS' combines AAE10, CO2, k, and PS variables with interactions; 'CO2_kPS' combines CO2, k, and PS variables with interactions. Columns show explained variance for different target variables.",
     digits = 3) |>
 kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
               full_width = FALSE)
```

```{r}
kable(Tables_yield_weather_lm, 
      caption = "Linear model performance for predicting yield and P-balance variables including weather data. Rows represent different predictor variable sets: 'onlyweather' uses only weather variables (annual average temperature, annual sum precipitation, juvenile deviation precipitation/sun/temperature, annual sum sun, plus NA weather indicator); 'k' combines weather variables with release rate constant; 'PS' combines weather variables with log-transformed semi-labile P; 'kPS' combines weather variables with k, log-transformed PS, and their interaction; 'AAE10' combines weather variables with log-transformed AAE10-extractable P; 'CO2' combines weather variables with log-transformed CO2-extractable P; 'AAE10_CO2' combines weather variables with both extractable P measures and their interaction; 'AAE10_CO2_kPS' combines weather variables with all P-dynamics parameters; 'CO2_kPS' combines weather variables with CO2, k, and PS parameters.",
      digits = 3) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = FALSE)
```

```{r}
kable(Tables_earth_lm, 
      caption = "Linear model performance comparison for predicting P-dynamics parameters from soil properties. Rows represent different target P-dynamics variables: 'PS_log' is log-transformed semi-labile phosphorus; 'k' is the phosphorus release rate constant; 'kPS_log' is log-transformed product of release rate and semi-labile P; 'P_AAE10_log' is log-transformed AAE10-extractable phosphorus; 'P_CO2_log' is log-transformed CO2-extractable phosphorus. The 'with_treatment' column uses soil variables (clay content, pH, organic carbon, silt content) plus treatment (P0 P100 P166), while 'without_treatment' uses only the soil variables.",
      digits = 3) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = FALSE)
```
