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

## Context

To explore, whether the proposed mechanisms and experiments to assess their dynamic, in a first step the Treatment levels $0P$ and $166P$ for all sites were analyzed. The experiments were conducted as displayed in the original paper of Flossmann & Richter with adjustments according to developments in technique and available equipment of the soil laboratory. Instead of the CAL-method, the Olsen-method was used to measure and estimate the quantity of P.

## Model of P-release after Flossman & Richter
$$\frac{dP}{dt}=k\times(P^S-P)$$
The constant $P^S$ denotes the amount of semi-labile P and was originally estimated as $P_\text{Olsen}-P_{H_2O}$. Subsequently the DE is solved exactly, since the soil is as $t=0$ mixed with deionized water, it was assumed that $P(0)=0$

## Exact Solution
$$P(t)=P^{\text{S}}-C\times e^{-kt}$$

for $P(0)=P_0$ we receive:

$$P(t)=P^S-(P^S-P_0)\times e^{-kt}$$

If we set $P(0)=0$ we receive:

$$P(t)=P^S\times(1-e^{-kt})q$$

## Linearization
Now we linearize the DE, so that a linear model can be employed to test the relation and estimate the parameters of interest:

$$P(t)=P^S-(P^S-P_0)\times e^{-kt}$$

$$P(t)-P^S=-(P^S-P_0)\times e^{-kt}$$

$$P^S-P(t)=(P^S-P_0)\times e^{-kt}$$

$$1-\frac{P(t)}{P^S}=(1-\frac{P_0}{P^S})\times e^{-kt}$$

Given $P_0=0$,

$$log(1-\frac{P(t)}{P^S})=-kt$$

## Setup and preparation of dataset


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

RES <- list()

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


try(setwd("~/Documents/Master Thesis/Master-Thesis-P-kinetics/"), silent=TRUE)

options(warn = -1)
# Dataset preparation
d <- tryCatch(read.csv(file = "./data/release-Data.csv", na.strings = c("", "n.a.")))
if (inherits(d, "try-error")){
  d <- read.csv(file = "~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/release-Data.csv", na.strings = c("", "n.a."))
}
d$Treatment <- d$Treatment |> factor()
d$Pv.mg.L. <- as.numeric(d$Pv.mg.L.)
d$Pv_labile.mg.L. <- as.numeric(d$Pv_labile.mg.L.)
d$Pv_Olsen.mg.L. <- as.numeric(d$Pv_Olsen.mg.L.)
d$uid <- paste(d$Site, d$Treatment, as.character(d$Repetition), sep="_")
d$Repetition <- as.factor(d$Repetition)
d$t.dt <- d$t.min. + 3
d[d$Treatment=="P0"&d$Site=="Ellighausen"&d$Repetition==4,"Pv_Olsen.mg.L."] <- NA
d$flos <- as.numeric(d$Pv_Olsen.mg.L.)-as.numeric(d$Pv_labile.mg.L.)
d$flos[d$flos <= 0] <- quantile(d$flos[d$flos > 0], 0.025, na.rm=TRUE)/2
d$Y1 <- log(1 - d$Pv.mg.L. / d$flos) # one timeseries removed, since flos too low
d$Y1[is.nan(d$Y1)] <- NA


str(d)
#> data.frame:	353 obs. of  21 variables:
#>  $ Site              : chr  "Cadenazzo" "Cadenazzo" "Cadenazzo" "Cadenazzo" ...
#>  $ Treatment         : Factor w/ 3 levels "0P","100P","166P": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ Repetition        : Factor w/ 4 levels "1","2","3","4": 2 2 2 2 2 2 1 1 1 1 ...
#>  $ Nummer            : num  1 1 1 1 1 1 2 2 2 2 ...
#>  $ mSoil_Olsen.g.    : chr  "10.05" "10.05" "10.05" "10.05" ...
#>  $ mSoil_H2O.g.      : num  10.3 10.3 10.3 10.3 10.3 ...
#>  $ Abs_labile.1.     : chr  "0.26" "0.26" "0.26" "0.26" ...
#>  $ Abs_Olsen.1.      : chr  "1.02" "1.02" "1.02" "1.02" ...
#>  $ Abs_Olsen1_10.1.  : chr  "0.28" "0.28" "0.28" "0.28" ...
#>  $ t.min.            : num  5 10 20 30 45 60 5 10 20 30 ...
#>  $ Abs.1.            : num  0.21 0.22 0.23 0.24 0.26 0.25 0.22 0.24 0.25 0.27 ...
#>  $ Pv_labile.mg.L.   : num  0.05 0.05 0.05 0.05 0.05 0.05 0.06 0.06 0.06 0.06 ...
#>  $ Pv_Olsen.mg.L.    : num  0.35 0.35 0.35 0.35 0.35 0.35 0.42 0.42 0.42 0.42 ...
#>  $ Pv_Olsen1_10.mg.L.: chr  "0.46" "0.46" "0.46" "0.46" ...
#>  $ Pv.mg.L.          : num  0.03 0.04 0.04 0.05 0.05 0.05 0.04 0.05 0.05 0.06 ...
#>  $ Dilution.1.       : num  20 20 20 20 20 20 20 20 20 20 ...
#>  $ V_H2O.L.          : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ...
#>  $ uid               : chr  "Cadenazzo_0P_2" "Cadenazzo_0P_2" "Cadenazzo_0P_2" "Cadenazzo_0P_2" ...
#>  $ flos              : num  0.3 0.3 0.3 0.3 0.3 0.3 0.36 0.36 0.36 0.36 ...
#>  $ Y1                : num  -0.105 -0.143 -0.143 -0.182 -0.182 ...
#>  $ nls_pred          : num  0.0283 0.04 0.0468 0.048 0.0482 ...
#>   ..- attr(*, "label")= chr "Fitted values"
```

Now we can see, whether our linearized model displays a linear relation.

```{r}
#| echo: false
#| warning: false
#| message: 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) 

```

If the parameter for the plateau could be estimated directly by using a non-linear non-least-squares model, we could omit the Olsen-measurement in the future.

LG: our nls is very sensitive to moderatly high Pv.mg.L at small time points. Since the ... disolves already before we start measureing, we will add 3 min to our time-measurement.

```{r}
#| echo: true
#| 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))))

# 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,scales = "free_y") +
  labs(x = TeX("$Time (min)$"),
       y = TeX("$P_{V}(\\frac{mg}{L})$")); suppressWarnings(print(p1))

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)

```


LG: hier machen wir folgendes:

1. Visualisiere Daten
2. for k*PS use sqrt-scale
3. Erkenne, dass keine offenslichtichen verletzuungen für ein lineares modell vorhanden sind
4. fitte ordinary linear squares model, with Treatment as the factor of interest and Site as covariate (analougous to paired t-test and equivalent to anova with Site as block factor)
5. Perform a classical Type II anova (using the car::Anova function)
6. Perform (post-hoc) TukeyHSD test (using multcomp package)

```{r}

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

Results:

1. for PS Treatment explains a lot, and site not so much. c.f. plot for a monotonous relationship
2. for k, the Treatment seems to be little relevant



```{r}
# new Data set, that gives info about Soil
allP <- tryCatch(readRDS("./data/all_P.rds"))
if (inherits(d, "try-error")){
  allP <- tryCatch(readRDS("~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/all_P.rds"))
}
allP$rep <- allP$rep %>% as.roman() %>% as.integer()
allP$uid <- paste(allP$location,allP$treatment_ID,as.character(allP$rep),sep = "_")
```


```{r}
# 1. merge this with nlme.coef
soil_data_2022 <- allP %>%
  filter(year == 2022) %>%
  select(soil_0_20_P_AAE10, soil_0_20_P_CO2,uid) %>%
  # Wir benennen die Spalten um, um später Konflikte zu vermeiden.
  rename(
    soil_AAE10_from_2022 = soil_0_20_P_AAE10,
    soil_CO2_from_2022 = soil_0_20_P_CO2
  )

 
Stycs_data <- allP %>% filter(year>=2017 & year<=2022) %>% 
  left_join(soil_data_2022, by = "uid")


Stycs_data_2022 <- Stycs_data %>%
  mutate(
    soil_0_20_P_AAE10 = ifelse(year != 2022, soil_AAE10_from_2022, soil_0_20_P_AAE10),
    soil_0_20_P_CO2   = ifelse(year != 2022, soil_CO2_from_2022, soil_0_20_P_CO2)
  ) %>%
  # Die Hilfsspalten werden nicht mehr benötigt und können entfernt werden.
  select(-soil_AAE10_from_2022, -soil_CO2_from_2022)


nlme.coef$kPS <- nlme.coef$k * nlme.coef$PS

D <- merge(nlme.coef, Stycs_data_2022,by = "uid")
D$uuid <- interaction(D$Site, D$treatment_ID, D$rep, D$year, D$crop)

# rm FodderCrop
D <- D[D$crop != "FodderCrop",]

# set wrong 0's to NA
D[D$soil_0_20_P_CO2 %in% c(NA,0), "soil_0_20_P_CO2"] <- NA
D[D$soil_0_20_P_AAE10 %in% c(NA,0), "soil_0_20_P_AAE10"] <- NA
D[D$annual_yield_mp_DM %in% c(NA,0), "annual_yield_mp_DM"] <- NA
D[D$Ymain_rel %in% c(NA,0), "Ymain_rel"] <- NA
D[D$annual_P_uptake %in% c(NA,0), "annual_P_uptake"] <- NA
D[D$annual_P_balance %in% c(NA,0), "annual_P_balance"] <- NA

# add log-transformed versions
# D$k_logPS <- D$k * log(D$PS)
D$kPS_log <- log(D$kPS)
D$PS_log <- log(D$PS)
D$P_AAE10_log <- log(D$soil_0_20_P_AAE10)
D$P_CO2_log <- log(D$soil_0_20_P_CO2)
D$P_AAE10_CO2_loglog <- log(D$soil_0_20_P_AAE10) * log(D$soil_0_20_P_CO2)




D$Site <- as.factor(D$Site)


## Compute Ymain_norm as the median yield of P166 treatment per site, year and crop
# set 0 values of annual_yield_mp_DM, Ymain_rel, annual_P_uptake, annual_P_balance to NA
dd <- aggregate(annual_yield_mp_DM ~ Site+Treatment+year+crop, data=D, median, na.rm = TRUE, na.action = na.pass)
# only keep P166
dd <- dd[dd$Treatment == "P166",]
dd$Treatment <- NULL
tmp <- merge(D, dd, by = c("Site", "year", "crop"), suffixes = c("", ".norm"), sort=FALSE)
nrow(tmp)
tmp$Ymain_norm <- tmp$annual_yield_mp_DM / tmp$annual_yield_mp_DM.norm
# order tmp by uid s.t. it matches D
tmp <- tmp[match(D$uuid, tmp$uuid),]
# check: cbind(tmp$uid, D$uid)
D$Ymain_norm <- tmp$Ymain_norm

oxides <- read.csv(file = "~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/oxides.csv", na.strings = c("", "n.a."))

D <- left_join(D,oxides,by="ui")

IEK <- read.csv(file = "~/Documents/Master Thesis/Master-Thesis-P-kinetics/data/IEK.csv", na.strings = c("", "n.a."))

IEK_wide <- IEK %>%
  pivot_wider(
    id_cols = uid,  # The column that uniquely identifies each sample
    names_from = IEK_time, # The column containing the time points (e.g., 1, 10, 30)
    values_from = c(n,E_mod, rt_R, E_exp), # The columns whose values you want to spread out
    names_sep = "_" # Separator for the new column names (e.g., E_exp_1, E_exp_10)
  )

# Step 3: Inspect the new wide-format dataframe
# It should now have one row per 'uid'
print("First few rows of the new wide-format IEK data:")
head(IEK_wide)

print("Column names of the new wide-format IEK data:")
names(IEK_wide)
# The new names will look like: "E_exp_1", "E_exp_10", "rt_R_1", "rt_R_10", etc.

# Step 4: Merge the wide IEK data with your main dataframe D
# Now that both have one row per uid, we can perform a clean join.
D <- D %>%
  left_join(IEK_wide, by = "uid")

###################################################
# Original Comparison 1 year k,PS vs. 6 years CO2 AAE10
###################################################

D2 <- merge(nlme.coef, allP[allP$year>=2017,],by = "uid")
D2$uuid <- interaction(D2$Site, D2$treatment_ID, D2$rep, D2$year, D2$crop)

# rm FodderCrop
D2 <- D2[D2$crop != "FodderCrop",]

# set wrong 0's to NA
D2[D2$soil_0_20_P_CO2 %in% c(NA,0), "soil_0_20_P_CO2"] <- NA
D2[D2$soil_0_20_P_AAE10 %in% c(NA,0), "soil_0_20_P_AAE10"] <- NA
D2[D2$annual_yield_mp_DM %in% c(NA,0), "annual_yield_mp_DM"] <- NA
D2[D2$Ymain_rel %in% c(NA,0), "Ymain_rel"] <- NA
D2[D2$annual_P_uptake %in% c(NA,0), "annual_P_uptake"] <- NA
D2[D2$annual_P_balance %in% c(NA,0), "annual_P_balance"] <- NA

# add log-transformed versions
# D2$k_logPS <- D2$k * log(D2$PS)
D2$kPS_log <- log(D2$kPS)
D2$PS_log <- log(D2$PS)
D2$P_AAE10_log <- log(D2$soil_0_20_P_AAE10)
D2$P_CO2_log <- log(D2$soil_0_20_P_CO2)
D2$P_AAE10_CO2_loglog <- log(D2$soil_0_20_P_AAE10) * log(D2$soil_0_20_P_CO2)




D2$Site <- as.factor(D2$Site)


## Compute Ymain_norm as the median yield of P166 treatment per site, year and crop
# set 0 values of annual_yield_mp_DM, Ymain_rel, annual_P_uptake, annual_P_balance to NA
dd <- aggregate(annual_yield_mp_DM ~ Site+Treatment+year+crop, data=D2, median, na.rm = TRUE, na.action = na.pass)
# only keep P166
dd <- dd[dd$Treatment == "P166",]
dd$Treatment <- NULL
tmp <- merge(D2, dd, by = c("Site", "year", "crop"), suffixes = c("", ".norm"), sort=FALSE)
nrow(tmp)
tmp$Ymain_norm <- tmp$annual_yield_mp_DM / tmp$annual_yield_mp_DM.norm
# order tmp by uid s.t. it matches D2
tmp <- tmp[match(D2$uuid, tmp$uuid),]
# check: cbind(tmp$uid, D2$uid)
D2$Ymain_norm <- tmp$Ymain_norm





RES$D <- D
RES$D2 <- D2
RES$nlme.coef.avg <- nlme.coef.avg
RES$data <- d
saveRDS(RES, file = "./data/RES.rds")
```


```{r}

```

