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

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.4770158	57.1681314	85.2730287	37.2687001	43.7593157	50.3258019	55.8767889	63.1306796
delta2	1.3727117	1.3113695	1.1431913	0.9323130	0.5909875	1.5180410	1.6523939	1.7781186	1.8720803	1.9621139
delta3	1.0620231	1.1174047	1.0635492	0.9350123	0.5892075	1.1137143	1.2535921	1.3815005	1.4783155	1.5865120
delta4	0.9595409	0.9773686	0.9669776	0.9156475	1.6750364	0.9532154	1.0460214	1.1336777	1.2025272	1.2811915
delta5	0.8816706	0.9569257	0.9840510	0.9582662	1.4560956	0.8271594	0.9164400	1.0027796	1.0717456	1.1532131
delta6	0.8861077	0.9459913	0.9675964	0.9386188	0.8380347	0.8058590	0.9025761	0.9935286	1.0654555	1.1509209
delta7	0.8175404	0.8643808	0.8881577	0.9013135	2.4540202	0.7457256	0.8304402	0.9120202	0.9781022	1.0592739
delta8	0.7580196	0.8313041	0.8502355	0.7951094	0.5876140	0.6682450	0.7542774	0.8364581	0.9026133	0.9845109
delta9	0.7672345	0.8051608	0.8127383	0.7725575	0.5889167	0.6698725	0.7467534	0.8223565	0.8835308	0.9593029
delta10	0.7567449	0.7612948	0.7606741	0.7503468	0.7395887	0.6953653	0.7578130	0.8212529	0.8710231	0.9291434
lambda	11.5272861	15.5548295	18.3712500	19.1444287	27.3492173	11.9365827	14.0126824	16.1170649	17.9008320	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.8379668	85.3310842	7.8325930	39.9103882	49.8146071	55.2221500	61.2740569
delta2	0.7698542	1.2564433	1.1008747	0.9075394	0.5909242	0.3688863	1.4328984	1.6443228	1.7210093	1.7773748
delta3	0.6735827	1.0949291	1.0436298	0.9198889	0.5892038	0.2989162	1.1226493	1.3257939	1.4159213	1.4963618
delta4	0.6020602	0.9580732	0.9557754	0.9096475	1.6758740	0.2578770	0.9410788	1.0947518	1.1638678	1.2340627
delta5	0.5776649	0.9550089	0.9875142	0.9602607	1.4575121	0.2368911	0.8414102	0.9821623	1.0494165	1.1211363
delta6	0.5784387	0.9618857	0.9849749	0.9451675	0.8385376	0.2289897	0.8299240	0.9827797	1.0569518	1.1361240
delta7	0.5330029	0.8526394	0.8815506	0.8982681	2.4571987	0.2203913	0.7502061	0.8797155	0.9426386	1.0118011
delta8	0.5399201	0.8465783	0.8527836	0.7894589	0.5876756	0.2243305	0.7157576	0.8479096	0.9164131	0.9922036
delta9	0.5363224	0.8104838	0.8146358	0.7690818	0.5889927	0.2329905	0.7011945	0.8220302	0.8851415	0.9577589
delta10	0.5299473	0.7691254	0.7673829	0.7531313	0.7398518	0.2604165	0.7243400	0.8285250	0.8787054	0.9317627
lambda	5.1246558	15.7137708	18.4768953	19.0852188	27.3724585	2.1824415	12.8200050	16.0382931	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.0063546	1.3461534	0.2513854	0.2700291	0.4071065	0.1371442	0.1329020
delta2	0.2927930	0.2605365	0.2620843	0.3691824	0.2104163	0.4388801	0.3053431	0.3022508	0.1822355	0.0803819
delta3	0.3635870	0.3175717	0.2812998	0.4634006	0.1228359	0.4863957	0.6945053	59.9700493	19.3326694	0.4417919
delta4	0.3410185	0.3328679	0.2952049	0.3760686	0.1394769	0.2220891	0.2643207	1.4287564	1.1649873	45.2624485
delta5	0.3384358	0.3673984	0.3604647	0.3532245	1743.1405670	0.1662756	0.2078227	1.3855875	0.1661685	42.4477761
delta6	0.3159736	0.3598314	0.4309145	588.9708158	1.2864590	0.1371503	0.1852244	1.2091556	0.1731099	2.4047794
delta7	0.3139827	0.3512311	0.3921932	4.9440803	18.6369872	0.1291231	0.1779783	1.2285557	0.1586679	0.1338800
delta8	0.3024390	0.3451335	0.3971786	34.3352194	36.7160706	0.1261642	0.1769319	1.0891819	0.6468224	2.5825714
delta9	0.2887453	0.3228545	0.3486171	1.1761314	43.6180167	0.1293902	0.1832293	0.9860158	0.5666844	0.2266257
delta10	0.2869471	0.3255191	0.3664096	4.4787093	20.7763514	0.1482400	0.2054014	0.9359451	0.4861150	0.2279680
lambda	0.8165636	0.6756489	0.6075678	181.2831030	634.4833983	0.5221463	0.5552243	18.5703908	5.8164414	25.9485088

- Introduction

Analyzing Ultramarathon Race Performances example

Introduction

Load data and prepare datasets

Best 200 races in term of the first finisher (clean data?)

Casewise contaminated data

Cellwise contaminated data

PLot Heatmaps of data in common scale

Figure 5

Estimation of parametrs

Figure 5

Compute absolute relative changes (in %)

Grouped Bar Plot for relative changes

Casewise contamination

Cellwise contamination

Table S3

Clean data

Casewise contamination

Cellwise contamination