Analysis_manuscript.Rmd

---
title: "Analysis_manuscript"
author: "Jinshi"
date: "7/16/2021"
output: html_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(message = TRUE,
                      #results = 'hide',
                      include = TRUE, # results appear in the finished file. 
                      echo = TRUE, # show code
                      warning = FALSE, # prevents warnings that are generated by code from appearing in the finished.
                      cache = TRUE, # save run results in the cache, it saves time
                      fig.height = 4, fig.width = 8)
```


## Load packages
```{r load packages}
# install.packages("dplyr")
library(dplyr)
library(ggplot2)
theme_set(theme_bw())
library(cowplot)
library(ggmap)
library(maps)
library(mapdata)
# install.packages("tidyverse")
library(tidyverse)
library(patchwork)
source("MCMC_function.R")
```

## Create functions
```{r create functions}
predict_vwc = function(theta_r, theta_s, alpha, n, h) {
  theta_r + (theta_s-theta_r)/ (1+(alpha*h)^n)^(1-1/n)
}
```


## Run and get parameters of VG model for NCSS samples  
* This chunk takes several hours, the outputs already saved in the outputs, no need to re-run again
```{r get parameters for NCSS}
# NCSSInv <- read.csv('data/NCSS_MCMC2_reorder.csv')
# head(NCSSInv)
# sapply(NCSSInv, class)
# length(levels(NCSSInv$labsampnum))
# length(unique(NCSSInv$labsampnum))
# nrow(NCSSInv)
# 
# 
# pdf(file='outputs/MCMC_rosseta2_final.pdf', width = 8, height = 12)
# results <- van_cal_param (NCSSInv) 
# dev.off()
# 
# write.csv(results,"outputs/NCSS_MCMC2_final.csv", row.names = F)

# load data
ncssfinal <- read.csv("outputs/NCSS_MCMC2_final.csv")
rosettaData <- read.csv('data/RosettaFittedVG.csv', header = T)
lu_2008_2013_data <- read.csv("data/Lu_2008_2013_Data.csv")

USA <- map_data("world", region = c("USA"), exact = FALSE) %>% 
  filter(long < 0)
# USA <- map_data("state")
```

## Test MCMC for other data (data sources from literature)
```{r test other data, include = FALSE}
# this chunk is already done, no need to run it again
# when all data included
soil_all_param1 <- general_van_param(lu_2008_2013_data %>% filter(soil == "I"))
soil_all_param2 <- general_van_param(lu_2008_2013_data %>% filter(soil == "II"))
soil_all_param3 <- general_van_param(lu_2008_2013_data %>% filter(soil == "III"))
soil_all_param4 <- general_van_param(lu_2008_2013_data %>% filter(soil == "IV"))
soil_all_param5 <- general_van_param(lu_2008_2013_data %>% filter(soil == "V"))
soil_all_param6 <- general_van_param(lu_2008_2013_data %>% filter(soil == "VI"))
soil_all_param7 <- general_van_param(lu_2008_2013_data %>% filter(soil == "VII"))
soil_all_param8 <- general_van_param(lu_2008_2013_data %>% filter(soil == "VIII"))


write.csv(soil_all_param1, "outputs/soil_all_param1.csv", row.names = FALSE)
write.csv(soil_all_param2, "outputs/soil_all_param2.csv", row.names = FALSE)
write.csv(soil_all_param3, "outputs/soil_all_param3.csv", row.names = FALSE)
write.csv(soil_all_param4, "outputs/soil_all_param4.csv", row.names = FALSE)
write.csv(soil_all_param5, "outputs/soil_all_param5.csv", row.names = FALSE)
write.csv(soil_all_param6, "outputs/soil_all_param6.csv", row.names = FALSE)
write.csv(soil_all_param7, "outputs/soil_all_param7.csv", row.names = FALSE)
write.csv(soil_all_param8, "outputs/soil_all_param8.csv", row.names = FALSE)

# only 5 data points used
lu_2008_2013_data %>%
  filter(soil == "I", h_cm %in% c(60, 102, 142.9, 403.1, 15306)) -> lu_soil_I

lu_2008_2013_data %>%
  filter(soil == "II", h_cm %in% c(51, 102, 295.9, 510.2, 15101.9)) -> lu_soil_II

lu_2008_2013_data %>%
  filter(soil == "III", h_cm %in% c(61.2, 102, 214.3, 510.2, 15306)) -> lu_soil_III

lu_2008_2013_data %>%
  filter(soil == "IV", h_cm %in% c(71.4, 102, 204.1, 510.2, 15306)) -> lu_soil_IV

lu_2008_2013_data %>%
  filter(soil == "V", h_cm %in% c(61.2, 102, 214.3, 510.2, 15306)) -> lu_soil_V

lu_2008_2013_data %>%
  filter(soil == "VI", h_cm %in% c(51, 102, 295.9, 510.2, 15101.9)) -> lu_soil_VI

lu_2008_2013_data %>%
  filter(soil == "VII", h_cm %in% c(51, 102, 295.9, 510.2, 15101.9)) -> lu_soil_VII

lu_2008_2013_data %>%
  filter(soil == "VIII", h_cm %in% c(50, 95.2, 203.9, 407.9, 15116.4)) -> lu_soil_VIII

soil_I_param1 <- general_van_param(lu_soil_I)
soil_I_param2 <- general_van_param(lu_soil_II)
soil_I_param3 <- general_van_param(lu_soil_III)
soil_I_param4 <- general_van_param(lu_soil_IV)
soil_I_param5 <- general_van_param(lu_soil_V)
soil_I_param6 <- general_van_param(lu_soil_VI)
soil_I_param7 <- general_van_param(lu_soil_VII)
soil_I_param8 <- general_van_param(lu_soil_VIII)

write.csv(soil_I_param1, "outputs/soil_I_param1.csv", row.names = FALSE)
write.csv(soil_I_param2, "outputs/soil_I_param2.csv", row.names = FALSE)
write.csv(soil_I_param3, "outputs/soil_I_param3.csv", row.names = FALSE)
write.csv(soil_I_param4, "outputs/soil_I_param4.csv", row.names = FALSE)
write.csv(soil_I_param5, "outputs/soil_I_param5.csv", row.names = FALSE)
write.csv(soil_I_param6, "outputs/soil_I_param6.csv", row.names = FALSE)
write.csv(soil_I_param7, "outputs/soil_I_param7.csv", row.names = FALSE)
write.csv(soil_I_param8, "outputs/soil_I_param8.csv", row.names = FALSE)
```



## Prepare data
```{r load data}
# calculate swc at different head pressure
h <- c(60,100,330,15000)
ncssfinal %>% 
  mutate(pred_vmc60 = predict_vwc(theta_r_mcmc, theta_s_mcmc, alpha_mcmc, n_mcmc, h[1]),
         pred_vmc100 = predict_vwc(theta_r_mcmc, theta_s_mcmc, alpha_mcmc, n_mcmc, h[2]),
         pred_vmc330 = predict_vwc(theta_r_mcmc, theta_s_mcmc, alpha_mcmc, n_mcmc, h[3]),
         pred_vmc15000 = predict_vwc(theta_r_mcmc, theta_s_mcmc, alpha_mcmc, n_mcmc, h[4]),) ->
  ncssfinal

# prepare data for scatter plot
bind_rows(
  ncssfinal %>% 
    dplyr::select(w6clodVMW, pred_vmc60) %>% 
    rename(mvwc = w6clodVMW, pvwc = pred_vmc60) %>%
    mutate(group = "(a)"),
  
  ncssfinal %>% 
    dplyr::select(w10cldVMW, pred_vmc100) %>%
    rename(mvwc = w10cldVMW, pvwc = pred_vmc100) %>%
    mutate(group = "(b)"),
  
  ncssfinal %>% 
    dplyr::select(w3cldVMW, pred_vmc330) %>%
    rename(mvwc = w3cldVMW, pvwc = pred_vmc330) %>%
    mutate(group = "(c)"),
  
  ncssfinal %>% 
    dplyr::select(w15l2VMW, pred_vmc15000) %>% 
    rename(mvwc = w15l2VMW, pvwc = pred_vmc15000) %>%
    mutate(group = "(d)")) -> 
  ncss_cbind

soil_I_param1 <- read.csv("outputs/soil_I_param1.csv")
soil_I_param2 <- read.csv("outputs/soil_I_param2.csv")
soil_I_param3 <- read.csv("outputs/soil_I_param3.csv")
soil_I_param4 <- read.csv("outputs/soil_I_param4.csv")
soil_I_param5 <- read.csv("outputs/soil_I_param5.csv")
soil_I_param6 <- read.csv("outputs/soil_I_param6.csv")
soil_I_param7 <- read.csv("outputs/soil_I_param7.csv")
soil_I_param8 <- read.csv("outputs/soil_I_param8.csv")

soil_all_param1 <- read.csv("outputs/soil_all_param1.csv")
soil_all_param2 <- read.csv("outputs/soil_all_param2.csv")
soil_all_param3 <- read.csv("outputs/soil_all_param3.csv")
soil_all_param4 <- read.csv("outputs/soil_all_param4.csv")
soil_all_param5 <- read.csv("outputs/soil_all_param5.csv")
soil_all_param6 <- read.csv("outputs/soil_all_param6.csv")
soil_all_param7 <- read.csv("outputs/soil_all_param7.csv")
soil_all_param8 <- read.csv("outputs/soil_all_param8.csv")

```


# get VG model parameters for lu_2008_2013
```{r get VG model parameters}
lu_2008_2013_data %>% 
  mutate(n = case_when(soil == "I" ~ mean(soil_I_param1$n),
                       soil == "II" ~ mean(soil_I_param2$n),
                       soil == "III" ~ mean(soil_I_param3$n),
                       soil == "IV" ~ mean(soil_I_param4$n),
                       soil == "V" ~ mean(soil_I_param5$n),
                       soil == "VI" ~ mean(soil_I_param6$n),
                       soil == "VII" ~ mean(soil_I_param7$n),
                       soil == "VIII" ~ mean(soil_I_param8$n)),
         alpha = case_when(soil == "I" ~ mean(soil_I_param1$alpha),
                           soil == "II" ~ mean(soil_I_param2$alpha),
                           soil == "III" ~ mean(soil_I_param3$alpha),
                           soil == "IV" ~ mean(soil_I_param4$alpha),
                           soil == "V" ~ mean(soil_I_param5$alpha),
                           soil == "VI" ~ mean(soil_I_param6$alpha),
                           soil == "VII" ~ mean(soil_I_param7$alpha),
                           soil == "VIII" ~ mean(soil_I_param8$alpha)),
         theta_s = case_when(soil == "I" ~ mean(soil_I_param1$theta_s),
                             soil == "II" ~ mean(soil_I_param2$theta_s),
                             soil == "III" ~ mean(soil_I_param3$theta_s),
                             soil == "IV" ~ mean(soil_I_param4$theta_s),
                             soil == "V" ~ mean(soil_I_param5$theta_s),
                             soil == "VI" ~ mean(soil_I_param6$theta_s),
                             soil == "VII" ~ mean(soil_I_param7$theta_s),
                             soil == "VIII" ~ mean(soil_I_param8$theta_s)),
         theta_r = case_when(soil == "I" ~ mean(soil_I_param1$theta_r),
                             soil == "II" ~ mean(soil_I_param2$theta_r),
                             soil == "III" ~ mean(soil_I_param3$theta_r),
                             soil == "IV" ~ mean(soil_I_param4$theta_r),
                             soil == "V" ~ mean(soil_I_param5$theta_r),
                             soil == "VI" ~ mean(soil_I_param6$theta_r),
                             soil == "VII" ~ mean(soil_I_param7$theta_r),
                             soil == "VIII" ~ mean(soil_I_param8$theta_r)),
         
         n_all = case_when(soil == "I" ~ mean(soil_all_param1$n),
                           soil == "II" ~ mean(soil_all_param2$n),
                           soil == "III" ~ mean(soil_all_param3$n),
                           soil == "IV" ~ mean(soil_all_param4$n),
                           soil == "V" ~ mean(soil_all_param5$n),
                           soil == "VI" ~ mean(soil_all_param6$n),
                           soil == "VII" ~ mean(soil_all_param7$n),
                           soil == "VIII" ~ mean(soil_all_param8$n)),
         alpha_all = case_when(soil == "I" ~ mean(soil_all_param1$alpha),
                               soil == "II" ~ mean(soil_all_param2$alpha),
                               soil == "III" ~ mean(soil_all_param3$alpha),
                               soil == "IV" ~ mean(soil_all_param4$alpha),
                               soil == "V" ~ mean(soil_all_param5$alpha),
                               soil == "VI" ~ mean(soil_all_param6$alpha),
                               soil == "VII" ~ mean(soil_all_param7$alpha),
                               soil == "VIII" ~ mean(soil_all_param8$alpha)),
         theta_s_all = case_when(soil == "I" ~ mean(soil_all_param1$theta_s),
                                 soil == "II" ~ mean(soil_all_param2$theta_s),
                                 soil == "III" ~ mean(soil_all_param3$theta_s),
                                 soil == "IV" ~ mean(soil_all_param4$theta_s),
                                 soil == "V" ~ mean(soil_all_param5$theta_s),
                                 soil == "VI" ~ mean(soil_all_param6$theta_s),
                                 soil == "VII" ~ mean(soil_all_param7$theta_s),
                                 soil == "VIII" ~ mean(soil_all_param8$theta_s)),
         theta_r_all = case_when(soil == "I" ~ mean(soil_all_param1$theta_r),
                                 soil == "II" ~ mean(soil_all_param2$theta_r),
                                 soil == "III" ~ mean(soil_all_param3$theta_r),
                                 soil == "IV" ~ mean(soil_all_param4$theta_r),
                                 soil == "V" ~ mean(soil_all_param5$theta_r),
                                 soil == "VI" ~ mean(soil_all_param6$theta_r),
                                 soil == "VII" ~ mean(soil_all_param7$theta_r),
                                 soil == "VIII" ~ mean(soil_all_param8$theta_r))   ) -> lu_2008_2013_data_VG_param

# precicted vwc
lu_2008_2013_data_VG_param %>% 
  mutate(pvwc_all = theta_r_all + (theta_s_all-theta_r_all)/ (1+(alpha_all*h_cm)^n_all)^(1-1/n_all),
         pvwc_5 = theta_r + (theta_s-theta_r)/ (1+(alpha*h_cm)^n)^(1-1/n)) -> lu_2008_2013_data_VG_param
  
write.csv(lu_2008_2013_data_VG_param, "outputs/lu_2008_2013_data_VG_param.csv", row.names = FALSE)
  
```


```{r}
lu_2008_2013_data_VG_param <- read.csv("outputs/lu_2008_2013_data_VG_param.csv")
```

## load RecT model results
```{r RecT data}
lu_2008_rect <- read.csv("outputs/lu_2008_2013_data_VG_param_rect.csv")

# precicted vwc by RecT
lu_2008_rect %>% 
  mutate(pvwc_all_rect = theta_r_rect_all + (theta_s_rect_all-theta_r_rect_all)/ (1+(alpha_rect_all*h_cm)^n_rect_all)^(1-1/n_rect_all),
         pvwc_5_rect = theta_r_rect5 + (theta_s_rect5-theta_r_rect5)/ (1+(alpha_rect5*h_cm)^n_rect5)^(1-1/n_rect5)) -> lu_2008_rect
```


```{r evaluation-scatter, fig.height=4, fig.width=8}
lu_2008_2013_data_VG_param %>% 
  ggplot(aes(pvwc_all, VWC)) + 
  geom_point() +
  geom_abline(slope = 1, intercept = 0) +
  geom_smooth(method = "lm") +
  facet_wrap(.~soil, nrow = 2) +
  labs(y=expression(Measured~SWC~('%')),
       x=expression(Predicted~SWC~('%')))

# ggsave("outputs/Figure_S5.jpg", width = 8, height = 4, units = "in")


lu_2008_2013_data_VG_param %>% 
  ggplot(aes(pvwc_5, VWC)) + 
  geom_point() +
  geom_abline(slope = 1, intercept = 0) +
  geom_smooth(method = "lm") +
  facet_wrap(.~soil, nrow = 2) +
  labs(y=expression(Measured~SWC~('%')),
       x=expression(Predicted~SWC~('%')))

# ggsave("outputs/Figure_S6.jpg", width = 8, height = 4, units = "in")
```


```{r Table2:model matrix}
model_matrix <- function(sdata, soil){
  lm_model <- lm(sdata$VWC ~ sdata$pvwc_all)
  adj_r2 = summary(lm_model)$adj.r.squared
  rmse = sqrt(sum(lm_model$residuals^2) / lm_model$df)
  me = sum(sdata$pvwc_all - sdata$VWC) / nrow(sdata)
  out = tibble(adj_r2 = adj_r2,
               rmse = rmse,
               me = me,
               obs = nrow(sdata),
               soil = soil)
  return(out)
}

model_matrix2 <- function(sdata, soil){
  lm_model <- lm(sdata$VWC ~ sdata$pvwc_5)
  adj_r2 = summary(lm_model)$adj.r.squared
  rmse = sqrt(sum(lm_model$residuals^2) / lm_model$df)
  me = sum(sdata$pvwc_5 - sdata$VWC) / nrow(sdata)
  out = tibble(adj_r2 = adj_r2,
               rmse = rmse,
               me = me,
               obs = nrow(sdata),
               soil = soil)
  return(out)
}

model_matrix3 <- function(sdata, soil){
  lm_model <- lm(sdata$VWC ~ sdata$pvwc_all_rect)
  adj_r2 = summary(lm_model)$adj.r.squared
  rmse = sqrt(sum(lm_model$residuals^2) / lm_model$df)
  me = sum(sdata$pvwc_all_rect - sdata$VWC) / nrow(sdata)
  out = tibble(adj_r2 = adj_r2,
               rmse = rmse,
               me = me,
               obs = nrow(sdata),
               soil = soil)
  return(out)
}

model_matrix4 <- function(sdata, soil){
  lm_model <- lm(sdata$VWC ~ sdata$pvwc_5_rect)
  adj_r2 = summary(lm_model)$adj.r.squared
  rmse = sqrt(sum(lm_model$residuals^2) / lm_model$df)
  me = sum(sdata$pvwc_5_rect - sdata$VWC) / nrow(sdata)
  out = tibble(adj_r2 = adj_r2,
               rmse = rmse,
               me = me,
               obs = nrow(sdata),
               soil = soil)
  return(out)
}

bind_rows(
  model_matrix(lu_2008_rect %>% filter(soil == "I"), "all Soil 1"),
  model_matrix(lu_2008_rect %>% filter(soil == "II"), "all Soil 2"),
  model_matrix(lu_2008_rect %>% filter(soil == "III"), "all Soil 3"),
  model_matrix(lu_2008_rect %>% filter(soil == "IV"), "all Soil 4"),
  model_matrix(lu_2008_rect %>% filter(soil == "V"), "all Soil 5"),
  model_matrix(lu_2008_rect %>% filter(soil == "VI"), "all Soil 6"),
  model_matrix(lu_2008_rect %>% filter(soil == "VII"), "all Soil 7"),
  model_matrix(lu_2008_rect %>% filter(soil == "VIII"), "all Soil 8"),
  
  # 5 points used 
  model_matrix2(lu_2008_rect %>% filter(soil == "I"), "5 Soil 1"),
  model_matrix2(lu_2008_rect %>% filter(soil == "II"), "5 Soil 2"),
  model_matrix2(lu_2008_rect %>% filter(soil == "III"), "5 Soil 3"),
  model_matrix2(lu_2008_rect %>% filter(soil == "IV"), "5 Soil 4"),
  model_matrix2(lu_2008_rect %>% filter(soil == "V"), "5 Soil 5"),
  model_matrix2(lu_2008_rect %>% filter(soil == "VI"), "5 Soil 6"),
  model_matrix2(lu_2008_rect %>% filter(soil == "VII"), "5 Soil 7"),
  model_matrix2(lu_2008_rect %>% filter(soil == "VIII"), "5 Soil 8"),
  
   # rect
  model_matrix3(lu_2008_rect %>% filter(soil == "I"), "all Soil 1"),
  model_matrix3(lu_2008_rect %>% filter(soil == "II"), "all Soil 2"),
  model_matrix3(lu_2008_rect %>% filter(soil == "III"), "all Soil 3"),
  model_matrix3(lu_2008_rect %>% filter(soil == "IV"), "all Soil 4"),
  model_matrix3(lu_2008_rect %>% filter(soil == "V"), "all Soil 5"),
  model_matrix3(lu_2008_rect %>% filter(soil == "VI"), "all Soil 6"),
  model_matrix3(lu_2008_rect %>% filter(soil == "VII"), "all Soil 7"),
  model_matrix3(lu_2008_rect %>% filter(soil == "VIII"), "all Soil 8"),
  
  model_matrix4(lu_2008_rect %>% filter(soil == "I"), "5 Soil 1"),
  model_matrix4(lu_2008_rect %>% filter(soil == "II"), "5 Soil 2"),
  model_matrix4(lu_2008_rect %>% filter(soil == "III"), "5 Soil 3"),
  model_matrix4(lu_2008_rect %>% filter(soil == "IV"), "5 Soil 4"),
  model_matrix4(lu_2008_rect %>% filter(soil == "V"), "5 Soil 5"),
  model_matrix4(lu_2008_rect %>% filter(soil == "VI"), "5 Soil 6"),
  model_matrix4(lu_2008_rect %>% filter(soil == "VII"), "5 Soil 7"),
  model_matrix4(lu_2008_rect %>% filter(soil == "VIII"), "5 Soil 8")
)


lm(lu_2008_2013_data_VG_param$VWC ~ lu_2008_2013_data_VG_param$pvwc_all) -> lm_model
lm_model$residuals
summary(lm_model)

lu_2008_2013_data_VG_param %>% filter(soil == "I") %>% nrow()
lu_2008_2013_data_VG_param %>% filter(soil == "II") %>% nrow()

```


## Plot map
```{r site map}
ggplot(data = USA) + 
  geom_polygon(aes(x = long, y = lat, group=group), fill = "gray", color = "white") +
  geom_point(aes(x = longitude_decimal_degrees, y = latitude_decimal_degrees),
             shape = 3, color = "blue", size = 2,
             data = ncssfinal) +
  labs(x = "Longitude", y = "Latitude") +
  xlim(c(-170, -60))

# ggsave("outputs/Figure1a_Sites.jpg", width = 8, height = 5, units = "in")
```

## Prior distribution of parameters
```{r Prior distribution}
ncssfinal %>% 
  dplyr::select(ten_alpha, ten_npar, theta_r, theta_s) %>% 
  gather() %>% 
  mutate(Groups = case_when(key == "ten_alpha" ~ "(d)",
                            key == "ten_npar" ~ "(c)",
                            key == "theta_r" ~ "(b)",
                            TRUE ~ "(a)")) %>% 
  ggplot(aes(value)) +
  geom_density(fill = "gray") +
  facet_wrap(.~Groups, scales = "free") +
  labs(x = expression(Value),
       y = expression(Density))

# add verticle line for each panel
vline_low <- tibble(X = c("(a)", "(b)", "(c)", "(d)"), Z = c(0.3, 0.0003, 0.0001, 1)) 
  

rosettaData %>% 
  dplyr::select(Alpha, Npar, Theta_r, Theta_s) %>% 
  gather() %>% 
  mutate(Groups = case_when(key == "Theta_s" ~ "(a)",
                            key == "Theta_r" ~ "(b)",
                            key == "Npar" ~ "(c)",
                            TRUE ~ "(d)")) %>% 
  mutate(v_low = case_when (key == "Theta_s" ~ '0.3',
                            key == "Theta_r" ~ '0.0003',
                            key == "Npar" ~ '1',
                            TRUE ~ '0.0001')) %>% 
  mutate(v_high = case_when (key == "Theta_s" ~ '0.8',
                            key == "Theta_r" ~ '0.3',
                            key == "Npar" ~ '10',
                            TRUE ~ '1')) %>% 
  ggplot(aes(value)) +
  geom_density(fill = "gray") +
  geom_rug(alpha = 0.5) +
  facet_wrap(.~Groups, scales = "free") +
  geom_vline(aes(xintercept = as.numeric(v_low)), col = "red", size = 1, linetype = "dotted") +
  geom_vline(aes(xintercept = as.numeric(v_high)), col = "red", size = 1, linetype = "dotted") +
  labs(x = expression(VG~model~parameters),
       y = expression(Density))

ggsave("outputs/Figure2.jpg", width = 8, height = 6, units = "in")
  
```

# Lu 2008 data analysis
## all data from Lu_2008 used
```{r}
# all data used
plot_vg_all <- function(sdata, sdata2) {
  h <- seq(0,1700000,5)
  theta_mcmc <- sdata$theta_r_all[1] + (sdata$theta_s_all[1]-sdata$theta_r_all[1])/ (1+(sdata$alpha_all[1]*h)^sdata$n_all[1])^(1-1/sdata$n_all[1])
  
  theta_rect <- sdata$theta_r_rect_all[1] + (sdata$theta_s_rect_all[1]-sdata$theta_r_rect_all[1])/ (1+(sdata$alpha_rect_all[1]*h)^sdata$n_rect_all[1])^(1-1/sdata$n_rect_all[1])
  
  sdata %>% 
    ggplot(aes(log(h_cm), VWC)) + 
    geom_point(shape = 1, size = 3.5) +
    # geom_point(aes(log(h_cm), VWC),
    #            data = sdata2, col = "red",
    #            shape = 3,
    #            size = 5) +
    geom_line(aes(log(h), theta_mcmc), data = tibble(h, theta_mcmc),
              col = "blue", size = 1.2) +
    geom_line(aes(log(h), theta_rect), data = tibble(h, theta_rect),
              col = "red", size = 1.2, alpha = 0.5) +
    theme(axis.title = element_blank()) +
    scale_x_continuous(breaks = seq(0,15,2.5))}

plot_vg_all(lu_2008_rect %>% filter(soil == "I"), lu_soil_I) -> plot_vg_all1
plot_vg_all(lu_2008_rect %>% filter(soil == "II"), lu_soil_II) -> plot_vg_all2
plot_vg_all(lu_2008_rect %>% filter(soil == "III"), lu_soil_III) -> plot_vg_all3
plot_vg_all(lu_2008_rect %>% filter(soil == "IV"), lu_soil_IV) -> plot_vg_all4
plot_vg_all(lu_2008_rect %>% filter(soil == "V"), lu_soil_V) -> plot_vg_all5
plot_vg_all(lu_2008_rect %>% filter(soil == "VI"), lu_soil_VI) -> plot_vg_all6
plot_vg_all(lu_2008_rect %>% filter(soil == "VII"), lu_soil_VII) -> plot_vg_all7
plot_vg_all(lu_2008_rect %>% filter(soil == "VIII"), lu_soil_VIII) -> plot_vg_all8

plot_vg_all1/plot_vg_all2/plot_vg_all3/plot_vg_all4/plot_vg_all5/plot_vg_all6/plot_vg_all7/plot_vg_all8 -> plot_vg_all
```


```{r plot trace of parameters}
# create function
plot_trace <- function(sdata, var_col) {
  dat <- tibble(x = (2000:10000),
                y = sdata[, var_col])
  dat %>% 
    ggplot(aes(x, y)) +
    geom_line() +
    scale_x_continuous(breaks = seq(2000,10000,3000))+
    theme(axis.title = element_blank())}

# trace of n
plot_trace(soil_all_param1, 1) -> trace_n_all1
plot_trace(soil_all_param2, 1) -> trace_n_all2
plot_trace(soil_all_param3, 1) -> trace_n_all3
plot_trace(soil_all_param4, 1) -> trace_n_all4
plot_trace(soil_all_param5, 1) -> trace_n_all5
plot_trace(soil_all_param6, 1) -> trace_n_all6
plot_trace(soil_all_param7, 1) -> trace_n_all7
plot_trace(soil_all_param8, 1) -> trace_n_all8
trace_n_all1/trace_n_all2/trace_n_all3/trace_n_all4/trace_n_all5/trace_n_all6/trace_n_all7/trace_n_all8 -> trace_n_all

# trace of alpha
plot_trace(soil_all_param1, 4) -> trace_a_all1
plot_trace(soil_all_param2, 4) -> trace_a_all2
plot_trace(soil_all_param3, 4) -> trace_a_all3
plot_trace(soil_all_param4, 4) -> trace_a_all4
plot_trace(soil_all_param5, 4) -> trace_a_all5
plot_trace(soil_all_param6, 4) -> trace_a_all6
plot_trace(soil_all_param7, 4) -> trace_a_all7
plot_trace(soil_all_param8, 4) -> trace_a_all8
trace_a_all1/trace_a_all2/trace_a_all3/trace_a_all4/trace_a_all5/trace_a_all6/trace_a_all7/trace_a_all8 -> trace_a_all

# trace of theta_s
plot_trace(soil_all_param1, 2) -> trace_s_all1
plot_trace(soil_all_param2, 2) -> trace_s_all2
plot_trace(soil_all_param3, 2) -> trace_s_all3
plot_trace(soil_all_param4, 2) -> trace_s_all4
plot_trace(soil_all_param5, 2) -> trace_s_all5
plot_trace(soil_all_param6, 2) -> trace_s_all6
plot_trace(soil_all_param7, 2) -> trace_s_all7
plot_trace(soil_all_param8, 2) -> trace_s_all8
trace_s_all1/trace_s_all2/trace_s_all3/trace_s_all4/trace_s_all5/trace_s_all6/trace_s_all7/trace_s_all8 -> trace_s_all

# trace of theta_r
plot_trace(soil_all_param1, 3) -> trace_r_all1
plot_trace(soil_all_param2, 3) -> trace_r_all2
plot_trace(soil_all_param3, 3) -> trace_r_all3
plot_trace(soil_all_param4, 3) -> trace_r_all4
plot_trace(soil_all_param5, 3) -> trace_r_all5
plot_trace(soil_all_param6, 3) -> trace_r_all6
plot_trace(soil_all_param7, 3) -> trace_r_all7
plot_trace(soil_all_param8, 3) -> trace_r_all8
trace_r_all1/trace_r_all2/trace_r_all3/trace_r_all4/trace_r_all5/trace_r_all6/trace_r_all7/trace_r_all8 -> trace_r_all
```


```{r comparisom plot, fig.height=12, fig.width=10}
plot_vg_all | trace_s_all | trace_r_all | trace_n_all | trace_a_all 
# ggsave("outputs/Figure3_lu2008_evaluation.jpg", width = 10, height = 12, units = "in")
```

## Distribution and ACF plot
```{r Histgram of parameter, fig.height=12, fig.width=10}
#plot hist of parameter
tiff(paste("outputs/FigureS2.HistGram.tiff"), width = 8, height = 9, pointsize = 1/300, units = 'in', res = 300)

par( mar=c(2, 0.5, 0.5, 0.5)
     , mai=c(0.15, 0.3, 0.1, 0.1)  # by inches, inner margin
     , omi = c(0.3, 0.3, 0.1, 0.1)  # by inches, outer margin 
     , mgp = c(0, 0.5, 0) # set distance of axis
     , tcl = 0.2
     , cex.axis = 1
     , las = 1
     , mfrow=c(8,4))

# theta_s
hist_theta_s <- function(sdata){
  truehist(sdata$theta_s, col='gray', xlab='', ylab='', main='', las=1)
  lines(density(sdata$theta_s, col='blue', lwd=2))
  box()}
# theta_r
hist_theta_r <- function(sdata){
  truehist(sdata$theta_r, col='gray', xlab='', ylab='', main='', las=1)
  lines(density(sdata$theta_r, col='blue', lwd=2))
  box()}
# n
hist_n <- function(sdata){
  truehist(sdata$n, col='gray', xlab='', ylab='', main='', las=1)
  lines(density(sdata$n, col='blue', lwd=2))
  box()}
# alpha
hist_alpha <- function(sdata){
  truehist(sdata$alpha, col='gray', xlab='', ylab='', main='', las=1)
  lines(density(sdata$alpha, col='blue', lwd=2))
  box()}

# soil 1
hist_theta_s(soil_all_param1)
hist_theta_r(soil_all_param1)
hist_n(soil_all_param1)
hist_alpha(soil_all_param1)

# soil 2
hist_theta_s(soil_all_param2)
hist_theta_r(soil_all_param2)
hist_n(soil_all_param2)
hist_alpha(soil_all_param2)

# soil 3
hist_theta_s(soil_all_param3)
hist_theta_r(soil_all_param3)
hist_n(soil_all_param3)
hist_alpha(soil_all_param3)

# soil 4
hist_theta_s(soil_all_param4)
hist_theta_r(soil_all_param4)
hist_n(soil_all_param4)
hist_alpha(soil_all_param4)

# soil 5
hist_theta_s(soil_all_param5)
hist_theta_r(soil_all_param5)
hist_n(soil_all_param5)
hist_alpha(soil_all_param5)

# soil 6
hist_theta_s(soil_all_param6)
hist_theta_r(soil_all_param6)
hist_n(soil_all_param6)
hist_alpha(soil_all_param6)

# soil 7
hist_theta_s(soil_all_param7)
hist_theta_r(soil_all_param7)
hist_n(soil_all_param7)
hist_alpha(soil_all_param7)

# soil 8
hist_theta_s(soil_all_param8)
mtext(side = 1, text = expression("Histgram of "* theta[s]), line = 2.0, cex=1, outer = F)
hist_theta_r(soil_all_param8)
mtext(side = 1, text = expression("Histgram of "* theta[r]), line = 2.0, cex=1, outer = F)
hist_n(soil_all_param8)
mtext(side = 1, text = expression("Histgram of n"), line = 2.0, cex=1, outer = F)
hist_alpha(soil_all_param8)
mtext(side = 1, text = expression("Histgram of "* alpha), line = 2.0, cex=1, outer = F)
mtext(side = 2, text = expression("Frequency (n)"), line = 0, cex=1, outer = T, las = 0)
dev.off()
```

```{r ACF, fig.height=9, fig.width=8}

tiff(paste("outputs/FigureS3.ACF.tiff"), width = 8, height = 9, pointsize = 1/300, units = 'in', res = 300)

par( mar=c(2, 0.5, 0.5, 0.5)
     , mai=c(0.15, 0.3, 0.1, 0.1)  # by inches, inner margin
     , omi = c(0.3, 0.3, 0.1, 0.1)  # by inches, outer margin 
     , mgp = c(0, 0.5, 0) # set distance of axis
     , tcl = 0.2
     , cex.axis = 1
     , las = 1
     , mfrow=c(8,4))

# soil 1
acf(soil_all_param1$theta_s, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 2, text = expression("ACF"), line = 2, cex=1, outer = F, las=0)
acf(soil_all_param1$theta_r, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param1$n, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param1$alpha, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")

# soil 2
acf(soil_all_param2$theta_s, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 2, text = expression("ACF"), line = 2, cex=1, outer = F, las=0)
acf(soil_all_param2$theta_r, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param2$n, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param2$alpha, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")

# soil 3
acf(soil_all_param3$theta_s, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 2, text = expression("ACF"), line = 2, cex=1, outer = F, las=0)
acf(soil_all_param3$theta_r, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param3$n, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param3$alpha, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")

# soil 4
acf(soil_all_param4$theta_s, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 2, text = expression("ACF"), line = 2, cex=1, outer = F, las=0)
acf(soil_all_param4$theta_r, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param4$n, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param4$alpha, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")

# soil 5
acf(soil_all_param5$theta_s, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 2, text = expression("ACF"), line = 2, cex=1, outer = F, las=0)
acf(soil_all_param5$theta_r, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param5$n, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param5$alpha, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")

# soil 6
acf(soil_all_param6$theta_s, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 2, text = expression("ACF"), line = 2, cex=1, outer = F, las=0)
acf(soil_all_param6$theta_r, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param6$n, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param6$alpha, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")

# soil 7
acf(soil_all_param7$theta_s, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 2, text = expression("ACF"), line = 2, cex=1, outer = F, las=0)
acf(soil_all_param6$theta_r, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param6$n, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
acf(soil_all_param6$alpha, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")

# soil 8
acf(soil_all_param8$theta_s, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 2, text = expression("ACF"), line = 2, cex=1, outer = F, las=0)
mtext(side = 1, text = expression("Lag of "*theta[s]), line = 2.0, cex=1, outer = F)
acf(soil_all_param6$theta_r, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 1, text = expression("Lag of "*theta[r]), line = 2.0, cex=1, outer = F)
acf(soil_all_param6$n, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 1, text = expression("Lag of n"), line = 2.0, cex=1, outer = F)
acf(soil_all_param8$alpha, lag.max = 300, plot = TRUE, main = "",xlab = "",ylab = "")
mtext(side = 1, text = expression("Lag of "*alpha), line = 2.0, cex=1, outer = F)

dev.off()
```



# NCSS data analysis
## Evaluation: scatter plot
```{r Evaluation-scatter plot, fig.height=7}
ncss_cbind %>% 
  mutate(residuals = abs(mvwc - pvwc)) %>% 
  ggplot(aes(pvwc, mvwc)) +
  geom_point(aes(size = residuals), alpha = 0.5) +  
  geom_smooth(method = "lm") +
  facet_wrap(.~group, ncol = 2) +
  geom_abline(intercept = 0, slope = 1, col = "red", linetype = "dashed", size = 1) +
  labs(x = expression(Predicted~soil~water~contend~('%')),
       y = expression(Measured~soil~water~contend~('%'))) +
  guides(size=guide_legend(title="Abs (residuals)"))

# ggsave("outputs/Figure4_scatter.jpg", width = 8, height = 7, units = "in")

# [1] are measured swc and [2] are predicted swc
model_matrix_ncss <- function(sdata){
  colnames(sdata) <- c("swc", "pswc")
  lm_model <- lm(sdata$swc ~ sdata$pswc)
  adj_r2 = summary(lm_model)$adj.r.squared
  rmse = sqrt(sum(lm_model$residuals^2) / lm_model$df)
  me = sum(sdata$pswc - sdata$swc) / nrow(sdata)
  out = tibble(adj_r2 = adj_r2,
               rmse = rmse,
               me = me,
               obs = nrow(sdata))
  return(out)
}

model_matrix_ncss(ncssfinal %>% dplyr::select(w6clodVMW, pred_vmc60))
model_matrix_ncss(ncssfinal %>% dplyr::select(w10cldVMW, pred_vmc100))
model_matrix_ncss(ncssfinal %>% dplyr::select(w3cldVMW, pred_vmc330))
model_matrix_ncss(ncssfinal %>% dplyr::select(w15l2VMW, pred_vmc15000))

```

## Evaluation: compare with benchmark
```{r Evaluation-benckmark, fig.height=8}
# theta s
ncssfinal %>% 
  mutate(MCMC = theta_s_mcmc, Measured = w0) %>% 
  dplyr::select(MCMC, Measured) %>% 
  gather() %>% 
  ggplot(aes(value, fill = key))+
  geom_density(alpha = 0.5) +
  theme(legend.position = c(0.75, 0.5),
        legend.title = element_blank()) ->
  panel_s

# theta r
ncssfinal %>% 
  mutate(MCMC = theta_r_mcmc, Measured = w15l2VMW) %>%
  dplyr::select(MCMC, Measured) %>% 
  gather() %>% 
  ggplot(aes(value, fill = key))+
  geom_density(alpha = 0.5)+
  theme(legend.position = c(0.75, 0.5),
        legend.title = element_blank()) ->
  panel_r

# theta_r compared with theta_r from UNSODA
bind_rows(
  ncssfinal %>% 
    mutate(key = 'MCMC', value = theta_r_mcmc) %>% 
    dplyr::select(value, key), 
    
  rosettaData %>% 
    mutate(key = 'UNSODA', value = Theta_r) %>%
    dplyr::select(value, key)) %>% 
  ggplot(aes(value, fill = key))+
  geom_density(alpha = 0.5)+
  theme(legend.position = c(0.75, 0.5),
        legend.title = element_blank())

# compare with ROSETTA theta_r of NCSS
ncssfinal %>% 
  mutate(MCMC = theta_r_mcmc, Measured = theta_r) %>%
  dplyr::select(MCMC, Measured) %>% 
  gather() %>% 
  ggplot(aes(value, fill = key))+
  geom_density(alpha = 0.5)+
  theme(legend.position = c(0.75, 0.5),
        legend.title = element_blank())

# n
bind_rows(
  ncssfinal %>% 
    mutate(key = 'MCMC', value = n_mcmc) %>% 
    dplyr::select(value, key), 
    
  rosettaData %>% 
    mutate(key = 'UNSODA', value = Npar) %>%
    dplyr::select(value, key)) %>% 
  ggplot(aes(value, fill = key))+
  geom_density(alpha = 0.5)+
  theme(legend.position = c(0.75, 0.5),
        legend.title = element_blank()) ->
  panel_n

# alpha
bind_rows(
  ncssfinal %>% 
    mutate(key = 'MCMC', value = alpha_mcmc) %>% 
    dplyr::select(value, key), 
    
  rosettaData %>% 
    mutate(key = 'UNSODA', value = Alpha) %>%
    dplyr::select(value, key)) %>% 
  ggplot(aes(value, fill = key))+
  geom_density(alpha = 0.5) +
  xlim(c(0, 0.4))+
  theme(legend.position = c(0.75, 0.5),
        legend.title = element_blank()) ->
  panel_alpha

(panel_s | panel_r) / (panel_n | panel_alpha) +
  plot_annotation(tag_levels = "a",
                  tag_prefix = "(",
                  tag_suffix = ")")
  

# ggsave("outputs/Figure5_evaluation.jpg", width = 8, height = 6, units = "in")
```


## lu2008 data: only 5 points used
```{r lu subset-test other data-plot}
# plot Measured points and fit line
# only 5 data point used
plot_vg <- function(sdata, sdata2) {
  h <- seq(0,1700000,5)
  theta_mcmc <- sdata$theta_r[1] + (sdata$theta_s[1]-sdata$theta_r[1])/ (1+(sdata$alpha[1]*h)^sdata$n[1])^(1-1/sdata$n[1])
  
  theta_rect <- sdata$theta_r_rect5[1] + (sdata$theta_s_rect5[1]-sdata$theta_r_rect5[1])/ (1+(sdata$alpha_rect5[1]*h)^sdata$n_rect5[1])^(1-1/sdata$n_rect5[1])
  
  sdata %>% 
    ggplot(aes(log(h_cm), VWC)) + 
    geom_point(shape = 1, size = 3.5) +
    geom_point(aes(log(h_cm), VWC),
               data = sdata2, col = "red",
               shape = 3,
               size = 5,
               stroke = 2) +
    geom_line(aes(log(h), theta_mcmc), data = tibble(h, theta_mcmc),
              col = "blue", size = 1.2) +
    geom_line(aes(log(h), theta_rect), data = tibble(h, theta_rect),
              col = "red", size = 1.2, alpha = 0.5) +
    theme(axis.title = element_blank()) +
    scale_x_continuous(breaks = seq(0,15,2.5))}

plot_vg(lu_2008_rect %>% filter(soil == "I"), lu_soil_I) -> plot_vg_soil1
plot_vg(lu_2008_rect %>% filter(soil == "II"), lu_soil_II) -> plot_vg_soil2
plot_vg(lu_2008_rect %>% filter(soil == "III"), lu_soil_III) -> plot_vg_soil3
plot_vg(lu_2008_rect %>% filter(soil == "IV"), lu_soil_IV) -> plot_vg_soil4
plot_vg(lu_2008_rect %>% filter(soil == "V"), lu_soil_V) -> plot_vg_soil5
plot_vg(lu_2008_rect %>% filter(soil == "VI"), lu_soil_VI) -> plot_vg_soil6
plot_vg(lu_2008_rect %>% filter(soil == "VII"), lu_soil_VII) -> plot_vg_soil7
plot_vg(lu_2008_rect %>% filter(soil == "VIII"), lu_soil_VIII) -> plot_vg_soil8

plot_vg_soil1/plot_vg_soil2/plot_vg_soil3/plot_vg_soil4/plot_vg_soil5/plot_vg_soil6/plot_vg_soil7/plot_vg_soil8 -> plot_vg_soil

```



```{r plot trace of parameters}
# create function
# 5 data point used
# trace of n
plot_trace(soil_I_param1, 1) -> trace_n_soil1
plot_trace(soil_I_param2, 1) -> trace_n_soil2
plot_trace(soil_I_param3, 1) -> trace_n_soil3
plot_trace(soil_I_param4, 1) -> trace_n_soil4
plot_trace(soil_I_param5, 1) -> trace_n_soil5
plot_trace(soil_I_param6, 1) -> trace_n_soil6
plot_trace(soil_I_param7, 1) -> trace_n_soil7
plot_trace(soil_I_param8, 1) -> trace_n_soil8
trace_n_soil1/trace_n_soil2/trace_n_soil3/trace_n_soil4/trace_n_soil5/trace_n_soil6/trace_n_soil7/trace_n_soil8 -> trace_n_soil

# trace of alpha
plot_trace(soil_I_param1, 4) -> trace_a_soil1
plot_trace(soil_I_param2, 4) -> trace_a_soil2
plot_trace(soil_I_param3, 4) -> trace_a_soil3
plot_trace(soil_I_param4, 4) -> trace_a_soil4
plot_trace(soil_I_param5, 4) -> trace_a_soil5
plot_trace(soil_I_param6, 4) -> trace_a_soil6
plot_trace(soil_I_param7, 4) -> trace_a_soil7
plot_trace(soil_I_param8, 4) -> trace_a_soil8
trace_a_soil1/trace_a_soil2/trace_a_soil3/trace_a_soil4/trace_a_soil5/trace_a_soil6/trace_a_soil7/trace_a_soil8 -> trace_a_soil

# trace of theta_s
plot_trace(soil_I_param1, 2) -> trace_s_soil1
plot_trace(soil_I_param2, 2) -> trace_s_soil2
plot_trace(soil_I_param3, 2) -> trace_s_soil3
plot_trace(soil_I_param4, 2) -> trace_s_soil4
plot_trace(soil_I_param5, 2) -> trace_s_soil5
plot_trace(soil_I_param6, 2) -> trace_s_soil6
plot_trace(soil_I_param7, 2) -> trace_s_soil7
plot_trace(soil_I_param8, 2) -> trace_s_soil8
trace_s_soil1/trace_s_soil2/trace_s_soil3/trace_s_soil4/trace_s_soil5/trace_s_soil6/trace_s_soil7/trace_s_soil8 -> trace_s_soil

# trace of theta_r
plot_trace(soil_I_param1, 3) -> trace_r_soil1
plot_trace(soil_I_param2, 3) -> trace_r_soil2
plot_trace(soil_I_param3, 3) -> trace_r_soil3
plot_trace(soil_I_param4, 3) -> trace_r_soil4
plot_trace(soil_I_param5, 3) -> trace_r_soil5
plot_trace(soil_I_param6, 3) -> trace_r_soil6
plot_trace(soil_I_param7, 3) -> trace_r_soil7
plot_trace(soil_I_param8, 3) -> trace_r_soil8
trace_r_soil1/trace_r_soil2/trace_r_soil3/trace_r_soil4/trace_r_soil5/trace_r_soil6/trace_r_soil7/trace_r_soil8 -> trace_r_soil
```


```{r comparisom plot, fig.height=12, fig.width=10}
plot_vg_soil | trace_s_soil | trace_r_soil | trace_n_soil | trace_a_soil
ggsave("outputs/FigureSX_lu2008_evaluation.jpg", width = 10, height = 12, units = "in")
```

## Changing n and alpha and their effects on the curve
```{r Testing n and alpha, fig.height=8}

# tiff("outputs/Figure_S1.tiff", width = 8, height = 8, pointsize = 1/300, units = 'in', res = 300)
par( mar=c(2, 0.2, 0.2, 0.2)
     , mai=c(0.2, 0.3, 0.1, 0.1)  # by inches, inner margin
     , omi = c(0.4, 0.5, 0., 0.1)  # by inches, outer margin
     , mgp = c(0, 0.3, 0) # set distance of axis
     , tcl = 0.2
     , cex.axis = 1.25
     , las = 1
     , mfrow=c(2,1))

n <- mean(ncssfinal$ten_npar, na.rm = T)
alpha <- mean(ncssfinal$ten_alpha, na.rm = T)
theta_s <- mean(ncssfinal$theta_s, na.rm = T)
theta_r <- mean(ncssfinal$theta_r, na.rm = T)


h <- seq(0,20000,1)

SWC <- theta_r + (theta_s-theta_r)/ (1+(alpha*h)^n)^(1-1/n)

plot(log(h),SWC, type="l", col="green", lwd=2
     , xlab = ""
     , ylab = ""
     ,ylim=c(theta_r, theta_s))

#****************************************************************************************************
# test n value
var_n <- seq(from=1.1, to=3.5, by=0.3)

for(i in 1:length(var_n)){
  var_SWC <- theta_r + (theta_s-theta_r)/ (1+(alpha*h)^var_n[i])^(1-1/var_n[i])
  points(log(h),var_SWC, type="l", col=(i+1), lwd=2, lty=i+1)
}

text(9,0.4, expression('( a )'), cex = 1.5)

#****************************************************************************************************
# test alpha
h <- seq(0,20000,1)
SWC <- theta_r + (theta_s-theta_r)/ (1+(alpha*h)^n)^(1-1/n)

plot(log(h),SWC, type="l", col="green", lwd=2
     , xlab = ""
     , ylab = ""
     ,ylim=c(theta_r, theta_s))

# points(log(h)
#        , theta_r + (theta_s-theta_r)/ (1+(0.005*h)^n)^(1-1/n)
#        , type="l", col=(i+1), lwd=2, lty=i+1)

# points(log(h)
#        , theta_r + (theta_s-theta_r)/ (1+(0.05*h)^n)^(1-1/n)
#        , type="l", col=(i+1), lwd=2, lty=i+1)


# change n value
var_alpha <- seq(from=0.005, to=0.75, by=0.05)

for(i in 1:length(var_alpha)){
  var_SWC <- theta_r + (theta_s-theta_r)/ (1+(var_alpha[i]*h)^n)^(1-1/n)
  points(log(h),var_SWC, type="l", col=(i+1), lwd=2, lty=i+1)
}

text(9,0.4, expression('( b )'), cex = 1.5)

mtext(side = 1, text = expression("log (h)"), line = 0.75, cex=1.5, outer = T, las = 0)
mtext(side = 2, text = expression("SWC (%)"), line = 0.5, cex=1.5, outer = T, las = 0)

# dev.off()

```

```{r Table A1}
lu_2008_rect %>% 
  dplyr::select(soil, theta_r_all, theta_s_all, alpha_all, n_all,
         theta_r_rect_all, theta_s_rect_all, alpha_rect_all, n_rect_all) %>% 
  distinct()

lu_2008_rect %>% 
  dplyr::select(soil, theta_r, theta_s, alpha, n,
         theta_r_rect5, theta_s_rect5, alpha_rect5, n_rect5) %>% 
  distinct()
```


```{r diagnose convergence}
library(coda)
post=as.matrix(soil_all_param1)
postmcmc=as.mcmc(post)
raftery.diag(postmcmc,q=0.025)
geweke.diag(post)

post=as.matrix(soil_all_param2)
postmcmc=as.mcmc(post)
raftery.diag(postmcmc,q=0.025)
geweke.diag(post)
```

```{r sdv testing}
mov_mn_sdv <- function(sdata, row_n) {
  MN <- rep(NA, nrow(sdata)-5)
  SDV <- rep(NA, nrow(sdata)-5)
  
  for (i in 5:nrow(sdata)) {
    subdata = sdata[1:i,row_n]
    mn = mean(subdata)
    sdv = sd(subdata)
    
    MN[i-5] = mn
    SDV[i-5] = sdv
    
    print(i)
  }
  mov_m_sd <- as.data.frame(cbind(MN, SDV))
  return (mov_m_sd)
}

# mov_m_sd <- mov_mn_sdv(soil_all_param1_test, 1)
# mov_m_sd_s <- mov_mn_sdv(soil_all_param1_test, 2)
# mov_m_sd_r <- mov_mn_sdv(soil_all_param1_test, 3)
# mov_m_sd_a <- mov_mn_sdv(soil_all_param1_test, 4)

```




```{r}
mov_m_sd %>% 
  ggplot(aes(x = 1:nrow(mov_m_sd), y = MN)) +
  geom_line() +
  geom_line(aes(x = 1:nrow(mov_m_sd), y = SDV), col = "red") +
  labs(x = expression(Sampling~time),
       y = expression(n)) -> plot_n_mov

mov_m_sd_s %>% 
  ggplot(aes(x = 1:nrow(mov_m_sd), y = MN)) +
  geom_line() +
  geom_line(aes(x = 1:nrow(mov_m_sd), y = SDV), col = "red") +
  labs(x = expression(Sampling~time),
       y = expression(theta_s)) -> plot_s_mov

mov_m_sd_r %>% 
  ggplot(aes(x = 1:nrow(mov_m_sd), y = MN)) +
  geom_line() +
  geom_line(aes(x = 1:nrow(mov_m_sd), y = SDV), col = "red") +
  labs(x = expression(Sampling~time),
       y = expression(theta_r)) -> plot_r_mov

mov_m_sd_a %>% 
  ggplot(aes(x = 1:nrow(mov_m_sd), y = MN)) +
  geom_line() +
  geom_line(aes(x = 1:nrow(mov_m_sd), y = SDV), col = "red") +
  labs(x = expression(Sampling~time),
       y = expression(alpha)) -> plot_a_mov

(plot_n_mov | plot_s_mov) / (plot_r_mov | plot_a_mov)
```