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 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.

library("mvdpd")
library("ggplot2")
library("reshape2")
library(tidyr)

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:reshape2':
## 
##     smiths

## The following object is masked from 'package:cellWise':
## 
##     unpack

library(knitr)  

Load data and prepare datasets

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?)

best0 <- ordered_data[1:200,]
best0 <- best0/3600

Casewise contaminated data

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

Cellwise contaminated data

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

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

heatplot(best0)

heatplot(best1)

heatplot(best2)

Estimation of parametrs

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)
}

results0 <- perform.analysis(best0)
results1 <- perform.analysis(best1)
results2 <- perform.analysis(best2)

Figure 5

Compute absolute relative changes (in %)

casewiseCDPD <- 100*abs(results1$MCDPDE-results0$MCDPDE)/results0$MCDPDE
cellwiseCDPD <- 100*abs(results2$MCDPDE-results0$MCDPDE)/results0$MCDPDE

Grouped Bar Plot for relative changes

# 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

mybarplot(casewiseCDPD[11,], casewiseCDPD[1,])

## Warning in data.frame(Method = factor(methods, levels = methods), Row_1 = X1, :
## row names were found from a short variable and have been discarded

## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Cellwise contamination

mybarplot(cellwiseCDPD[11,], cellwiseCDPD[1,])

## Warning in data.frame(Method = factor(methods, levels = methods), Row_1 = X1, :
## row names were found from a short variable and have been discarded

Table S3

Clean data

kable(cbind(results0$MDPDE, results0$MCDPDE))

	ML	DPD(0.1)	DPD(0.2)	DPD(0.3)	DPD(0.5)	CML	CDPD(0.1)	CDPD(0.2)	CDPD(0.3)	CDPD(0.5)
delta1	35.5937929	47.5632934	55.4770160	57.1681311	85.2023908	37.2687001	43.7593157	50.3258020	55.8767883	63.1306796
delta2	1.3727117	1.3113695	1.1431913	0.9323130	0.5908685	1.5180410	1.6523939	1.7781186	1.8720803	1.9621139
delta3	1.0620231	1.1174047	1.0635493	0.9350123	0.5890869	1.1137143	1.2535921	1.3815005	1.4783155	1.5865120
delta4	0.9595409	0.9773686	0.9669776	0.9156475	1.6739878	0.9532154	1.0460214	1.1336777	1.2025272	1.2811915
delta5	0.8816706	0.9569257	0.9840510	0.9582662	1.4554821	0.8271594	0.9164400	1.0027796	1.0717456	1.1532131
delta6	0.8861077	0.9459913	0.9675964	0.9386188	0.8378811	0.8058590	0.9025761	0.9935286	1.0654554	1.1509209
delta7	0.8175404	0.8643808	0.8881577	0.9013135	2.4527579	0.7457256	0.8304402	0.9120202	0.9781022	1.0592739
delta8	0.7580196	0.8313041	0.8502355	0.7951094	0.5875123	0.6682450	0.7542774	0.8364581	0.9026133	0.9845109
delta9	0.7672345	0.8051608	0.8127383	0.7725575	0.5888432	0.6698725	0.7467534	0.8223565	0.8835308	0.9593029
delta10	0.7567449	0.7612948	0.7606741	0.7503468	0.7394154	0.6953653	0.7578130	0.8212529	0.8710231	0.9291434
lambda	11.5272861	15.5548295	18.3712501	19.1444286	27.3300750	11.9365827	14.0126824	16.1170649	17.9008318	20.2443890

Casewise contamination

kable(cbind(results1$MDPDE, results1$MCDPDE))

	ML	DPD(0.1)	DPD(0.2)	DPD(0.3)	DPD(0.5)	CML	CDPD(0.1)	CDPD(0.2)	CDPD(0.3)	CDPD(0.5)
delta1	18.0192523	47.8143091	55.5464007	56.8379667	85.3278150	7.8325930	39.9103882	49.8146072	55.2221501	61.2740571
delta2	0.7698542	1.2564433	1.1008747	0.9075394	0.5908794	0.3688863	1.4328984	1.6443228	1.7210093	1.7773748
delta3	0.6735827	1.0949291	1.0436298	0.9198889	0.5891781	0.2989162	1.1226493	1.3257939	1.4159214	1.4963618
delta4	0.6020602	0.9580732	0.9557754	0.9096475	1.6760233	0.2578770	0.9410788	1.0947518	1.1638677	1.2340627
delta5	0.5776649	0.9550089	0.9875142	0.9602607	1.4574815	0.2368911	0.8414102	0.9821623	1.0494165	1.1211363
delta6	0.5784387	0.9618857	0.9849749	0.9451675	0.8384121	0.2289897	0.8299240	0.9827797	1.0569518	1.1361240
delta7	0.5330029	0.8526394	0.8815506	0.8982681	2.4574553	0.2203913	0.7502061	0.8797155	0.9426386	1.0118011
delta8	0.5399201	0.8465783	0.8527836	0.7894589	0.5876788	0.2243305	0.7157576	0.8479096	0.9164131	0.9922036
delta9	0.5363224	0.8104838	0.8146358	0.7690818	0.5889898	0.2329905	0.7011945	0.8220302	0.8851415	0.9577589
delta10	0.5299473	0.7691254	0.7673829	0.7531313	0.7398032	0.2604165	0.7243400	0.8285250	0.8787054	0.9317627
lambda	5.1246558	15.7137708	18.4768953	19.0852189	27.3698654	2.1824415	12.8200050	16.0382932	17.7877824	19.7588087

Cellwise contamination

kable(cbind(results2$MDPDE, results2$MCDPDE))

	ML	DPD(0.1)	DPD(0.2)	DPD(0.3)	DPD(0.5)	CML	CDPD(0.1)	CDPD(0.2)	CDPD(0.3)	CDPD(0.5)
delta1	0.2498760	0.2704943	0.2952131	0.0048965	1.7602844	0.2513854	0.2700291	0.4071210	0.1374468	0.1327150
delta2	0.2927930	0.2605365	0.2620843	0.4007590	0.1736728	0.4388801	0.3053431	0.3022416	0.1801788	0.0804210
delta3	0.3635870	0.3175717	0.2812998	0.4391833	0.1270807	0.4863957	0.6945053	59.9926168	16.6969810	0.4415987
delta4	0.3410185	0.3328679	0.2952049	0.3499743	0.1240661	0.2220891	0.2643207	1.4313818	1.0364156	45.2481402
delta5	0.3384358	0.3673984	0.3604647	0.3517751	2257.8928302	0.1662756	0.2078227	1.3855447	0.1568331	42.4343598
delta6	0.3159736	0.3598314	0.4309145	651.4184000	1.0128460	0.1371503	0.1852244	1.2097813	0.1743974	2.4030188
delta7	0.3139827	0.3512311	0.3921932	7.1171808	26.9175125	0.1291231	0.1779783	1.2286907	0.1595250	0.1343237
delta8	0.3024390	0.3451335	0.3971786	38.2496338	46.7678834	0.1261642	0.1769319	1.0896583	0.7132413	2.5822693
delta9	0.2887453	0.3228545	0.3486171	1.6106209	53.7341932	0.1293902	0.1832293	0.9866494	0.4017231	0.2275948
delta10	0.2869471	0.3255191	0.3664096	3.8686319	27.3799278	0.1482400	0.2054014	0.9353893	0.5212460	0.2278015
lambda	0.8165636	0.6756489	0.6075678	198.3140329	822.0271175	0.5221463	0.5552243	18.5775272	5.2093248	25.9419188