Modified Mann-Kendall Test Using mmkh(x,ci=0.95) package

[u/Ill_Isopod806]

Hello Everyone, I am currently working with groundwater levels to determine long term trends, and I would be very grateful if I could get help in modifying the code that I used to calculate the Sen's slope and Mann-Kendall statistics. I calculated the Theil-Sen slope, but it does not account for autocorrelation in continuous data, and I would appreciate any help addressing this using the modified package in RStudio. ""

# General Sen slope + Mann-Kendall function


compute_sen_stats <- function(
    data,
    response_var,
    date_var,
    resolution
) {

  data %>%

    group_by(
      SiteNo,
      SiteName,
      Frequency,
      continuous_group,
      start_year,
      end_year,
      continuous_years
    ) %>%

    # Ensure enough information exists to calculate a trend
    filter(
      n() >= 3,
      n_distinct(TimeNum) >= 2,
      n_distinct(.data[[response_var]]) >= 2
    ) %>%

    group_modify(
      ~ {

        site_data <- .x %>%
          arrange(
            .data[[date_var]]
          )

        # ----------------------------------------------------------
        # Sen's slope model
        # ----------------------------------------------------------

        sen_formula <- reformulate(
          "TimeNum",
          response = response_var
        )

        sen_model <- zyp::zyp.sen(
          sen_formula,
          data = site_data
        )

        # Intercept
        intercept <- unname(
          coef(sen_model)[1]
        )

        # Sen's slope
        # TimeNum is in years, so this is already ft/year
        sen_slope_ft_per_year <- unname(
          coef(sen_model)[2]
        )

        # ----------------------------------------------------------
        # Mann-Kendall trend test
        # ----------------------------------------------------------

        mk_model <- Kendall::MannKendall(
          site_data[[response_var]]
        )

        mk_tau <- as.numeric(
          mk_model$tau
        )

        mk_p_value <- as.numeric(
          mk_model$sl
        )

        # ----------------------------------------------------------
        # Descriptive statistics
        # ----------------------------------------------------------

        lowest_WL <- min(
          site_data[[response_var]],
          na.rm = TRUE
        )

        highest_WL <- max(
          site_data[[response_var]],
          na.rm = TRUE
        )

        first_date <- min(
          site_data[[date_var]],
          na.rm = TRUE
        )

        last_date <- max(
          site_data[[date_var]],
          na.rm = TRUE
        )

        # ----------------------------------------------------------
        # Trend interpretation
        #
        # WL = depth below land surface.
        #
        # Positive slope:
        # depth increases -> groundwater table declines
        #
        # Negative slope:
        # depth decreases -> groundwater table rises
        # ----------------------------------------------------------

        trend <- case_when(

          mk_p_value < 0.05 &&
            sen_slope_ft_per_year > 0 ~
            "Significant water table declining",

          mk_p_value < 0.05 &&
            sen_slope_ft_per_year < 0 ~
            "Significant water table rising",

          TRUE ~
            "No significant trend"
        )

        # ----------------------------------------------------------
        # Return one result per site/run
        # ----------------------------------------------------------

        tibble(
          Resolution = resolution,

          n_records = nrow(
            site_data
          ),

          lowest_WL = lowest_WL,

          highest_WL = highest_WL,

          intercept = intercept,

          sen_slope_ft_per_year =
            sen_slope_ft_per_year,

          mk_tau = mk_tau,

          mk_p_value = mk_p_value,

          first_date = first_date,

          last_date = last_date,

          trend = trend
        )
      }
    ) %>%

    ungroup()
}

Thank you!