---
title: "Analyzing Ultramarathon Race Performances example"
author: "A. Ghosh, C. Agostinelli and A. Basu"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{mvogammaDPD examples}
  %\VignetteEngine{knitr::knitr}
  %\VignetteEncoding{UTF-8}
---

```{r, echo = FALSE}
knitr::opts_chunk$set(
 fig.width = 8 ,
 fig.height = 12,
 fig.align ='center'
)
```

# Introduction

This file contains an example of the function ```mvogammaDPD```. 
It reproduces all the figures and tables in section "Application: Analyzing Ultramarathon Race Performances" of A. Ghosh, C. Agostinelli and A. Basu (2026) A Composite Divergence Approach to Robust Multivariate Estimation under Cellwise and Casewise Contamination, arXiv:2608.18914, https://arxiv.org/abs/2608.18914.

The dataset was obtained from the publicly available data on the website [www.kaggle.com/datasets/fatihyavuzz/two-centuries-of-um-races](https://www.kaggle.com/datasets/fatihyavuzz/two-centuries-of-um-races) which contains ace results from ultramarathon events held worldwide over approximately two centuries.

We restrict to mens' 50-mile races. For each race event (identified by the event name and date), finishing times were converted to hours and the ten fastest finishers were retained, producing one ordered $10$-dimensional observation per race. Events with fewer than ten recorded finishers or with tied finishing times among the top ten were discarded, yielding a collection of complete ordered observations of 4683 races.
  
```{r}
library("mvdpd")
library("ggplot2")
library("reshape2")
library(tidyr)
library(knitr)  
```
  
## Load data and prepare datasets

```{r}  
data("twocenturies50")
ordered_data <- twocenturies50[order(twocenturies50[,1]),]
n <- nrow(twocenturies50)
p <- ncol(twocenturies50)
```
### Best 200 races in term of the first finisher (clean data?)

```{r}  
best0 <- ordered_data[1:200,]
best0 <- best0/3600
```
### Casewise contaminated data

```{r}    
best1 <- ordered_data[c(1:180, (n-19):n),]
best1 <- best1/3600
```

### Cellwise contaminated data

```{r}    
best2 <- ordered_data[1:200,]
set.seed(1234)
for (i in 1:200) {
  nn <- rbinom(n=1, size=10, prob=0.1)
  if (nn){
    pos <- sample(1:p, size=nn)
    best2[i,pos] <- runif(1, 0, 3)*best0[i,pos]
    best2[i,] <- sort(best2[i,]) 
  }
}
best2 <- best2/3600
```

## PLot Heatmaps of data in common scale

```{r}  
min_val <- min(best0, best1, best2) 
max_val <- max(best0, best1, best2) 

heatplot <- function(mat) {
  df <- melt(mat)
  colnames(df) <- c("Event", "FinishingTime", "Hours")  
  df$FinishingTime <- as.numeric(as.factor(df$FinishingTime))
  
  ggplot(df, aes(FinishingTime, Event, fill = Hours)) +
    geom_tile() +
    scale_x_continuous(
      breaks = 1:10,
      labels = 1:10,
      expand = expansion(mult = c(0.01, 0.01))
    ) +
    scale_y_discrete(
      expand = expansion(mult = c(0.03, 0.03))
    ) +
    scale_fill_gradient2(
      low = "blue",
      mid = "white",
      high = "red",
      midpoint = (min_val + max_val) / 2,
      limits = c(min_val, max_val),
      breaks = round(seq(min_val, max_val, length.out = 5),0)
    ) +
    labs(
      x = "Top 10 Finishing Times",
      y = "Events",
      fill = "Hours"
    ) +
    theme_classic(base_size = 14)
}
```

### Figure 5
  
```{r, fig=TRUE}    
heatplot(best0)
```

```{r, fig=TRUE}      
heatplot(best1)
```

```{r, fig=TRUE}      
heatplot(best2)
```
  
## Estimation of parametrs
   
```{r}  
perform.analysis <- function(x) { 
  # Maximum likelihood
  ml <- mvogammaML(x=x, method="multivariate",
    initial=list(delta=c(18, rep(0.5, 9)), lambda=5))
  # Maximum composite likelihood (minimum CDPD with beta=0)
  cml <- mvogammaML(x=x, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum composite DPD
  cdpd1 <- mvogammaDPD(x=x, beta=0.1, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd2 <- mvogammaDPD(x=x, beta=0.2, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd3 <- mvogammaDPD(x=x, beta=0.3, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd5 <- mvogammaDPD(x=x, beta=0.5, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum (multivariate) DPD
  mdpd1 <- mvogammaDPD(x=x, beta=0.1, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd2 <- mvogammaDPD(x=x, beta=0.2, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd3 <- mvogammaDPD(x=x, beta=0.3, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd5 <- mvogammaDPD(x=x, beta=0.5, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))

  results.multivariate <- cbind(
    c(ml$delta, ml$lambda),
    c(mdpd1$delta, mdpd1$lambda),
    c(mdpd2$delta, mdpd2$lambda),
    c(mdpd3$delta, mdpd3$lambda),
    c(mdpd5$delta, mdpd5$lambda)
  )

  results.composite <- cbind(
    c(cml$delta, cml$lambda),
    c(cdpd1$delta, cdpd1$lambda),
    c(cdpd2$delta, cdpd2$lambda),
    c(cdpd3$delta, cdpd3$lambda),
    c(cdpd5$delta, cdpd5$lambda)
  )

  colnames(results.multivariate) <- c("ML", "DPD(0.1)", "DPD(0.2)",
                                     "DPD(0.3)", "DPD(0.5)")
  colnames(results.composite) <- c("CML", "CDPD(0.1)", "CDPD(0.2)",
                                  "CDPD(0.3)", "CDPD(0.5)")
  rownames(results.multivariate) <- c(paste0("delta", 1:10), "lambda")
  rownames(results.composite) <- rownames(results.multivariate) 
  res <- list(MDPDE=results.multivariate, MCDPDE=results.composite)
  return(res)
}
```
  
```{r}
results0 <- perform.analysis(best0)
results1 <- perform.analysis(best1)
results2 <- perform.analysis(best2)
```

### Figure 5

### Compute absolute relative changes (in %)
    
```{r}  
casewiseCDPD <- 100*abs(results1$MCDPDE-results0$MCDPDE)/results0$MCDPDE
cellwiseCDPD <- 100*abs(results2$MCDPDE-results0$MCDPDE)/results0$MCDPDE
```

### Grouped Bar Plot for relative changes 

```{r}  
# Define the methods (X-axis categories) in the exact order from the image
methods <- c("ML", "CML", "DPD(0.1)", "DPD(0.2)", "DPD(0.3)", 
             "DPD(0.5)", "CDPD(0.1)", "CDPD(0.2)", "CDPD(0.3)", "CDPD(0.5)")

mybarplot <- function(X1, X2) {
  # Create a data frame with the extracted percentages
df <- data.frame(Method = factor(methods, levels = methods),
                 Row_1 = X1, Row_2 = X2)

# Convert data to long format for ggplot2
df_long <- pivot_longer(df, cols = c(Row_1, Row_2), 
                        names_to = "Parameter", values_to = "Percentage")

# Create the Plot
p <- ggplot(df_long, aes(x = Method, y = Percentage, fill = Parameter)) +
  
  # Add grouped bars with black outlines
  # Matching position_dodge width and bar width removes the inner gap
  geom_bar(
    stat = "identity",
    position = position_dodge(width = 0.8), # Changed from 0.85 
    color = "black",       
    width = 0.8,                            # Changed from 0.75 to match dodge width
    linewidth = 0.6
  ) +
  
  # Add data labels on top of the bars
  geom_text(
    aes(label = sprintf("%.0f", Percentage)),
    position = position_dodge(width = 0.8), # Must match the dodge width in geom_bar
    vjust = -0.5,
    size = 4,           
    color = "black"
  ) +
  
  # Map colors and use mathematical expressions for the legend labels (lambda and delta)
  scale_fill_manual(
    values = c("Row_1" = "#4C72B0", "Row_2" = "#DD8452"),
    labels = c(expression(lambda), expression(delta[1]))
  ) +
  
  # Configure Y-axis: force origin at 0 and set max limit slightly above data
  scale_y_continuous(
    limits = c(0, max(df_long$Percentage) * 1.15), 
    breaks = seq(0, 100, 20),
    expand = expansion(mult = c(0, 0)) 
  ) +
  
  # Axis and Legend Labels
  labs(
    x = "Estimation Method",
    y = "Absolute relative changes (in %)",
    fill = "Parameter"
  ) +
  
  # Theme Modifications
  theme_classic(base_size = 14) +
  
  theme(
    # Axis tick values bold and larger
    axis.text.x = element_text(angle = 45, hjust = 1, color = "black", size = 14, face = "bold"),
    axis.text.y = element_text(color = "black", size = 14, face = "bold"),
    
    # Axis titles
    axis.title = element_text(face = "bold", size = 16),
    
    axis.line = element_line(color = "black", linewidth = 0.7),
    axis.ticks = element_line(color = "black", linewidth = 0.7),
    
    # Legend formatting
    legend.position = c(0.88, 0.88),
    legend.background = element_rect(fill = "white", color = NA),
    legend.box.background = element_rect(color = "black", linewidth = 0.5),
    legend.title = element_text(face = "bold", size = 14),
    legend.text = element_text(size = 15), 
    
    # Plot margins
    plot.margin = margin(t = 15, r = 15, b = 10, l = 10)
  )

  # Display the plot
  print(p)
}
```

### Casewise contamination 

```{r, fig=TRUE}      
mybarplot(casewiseCDPD[11,], casewiseCDPD[1,])
```


### Cellwise contamination  
  

```{r, fig=TRUE}      
mybarplot(cellwiseCDPD[11,], cellwiseCDPD[1,])
```

### Table S3  

### Clean data

```{r}
kable(cbind(results0$MDPDE, results0$MCDPDE))
```
  
### Casewise contamination 
  
```{r}
kable(cbind(results1$MDPDE, results1$MCDPDE))
```

### Cellwise contamination  
  
```{r}
kable(cbind(results2$MDPDE, results2$MCDPDE))
```