Glacier weather data – Sonnblick observatory
Weather data from Austria’s highest weather station, situated in the Austrian Central Alps on the glaciated mountain “Hoher Sonnblick”, standing at an elevation of 3106 meters above sea level. An array of dimension \((p,q,n)\), comprising \(n = 136\) observations, each represented by a \(p = 5\) times \(q = 12\) dimensional matrix. Observed parameters are monthly averages of
- air pressure (AP)
- precipitation (P)
- sunshine hours (SH)
- temperature (T)
- proportion of solid precipitation (SP)
from 1891 to 2022.
data(weather)
X <- weather[,,!apply(weather,3, function(x) any(is.na(x)))]
n<- dim(X)[3]
p<- dim(X)[1]
q<- dim(X)[2]Parameter estimation of mean, rowwise and columnwise covariance matrix using maximum likelihood estimation and robust parameter estimation based on MMCD.
par_MMLE <- mmle(X, lambda = 0)
set.seed(1)
par_MMCD <- mmcd(X, alpha = 0.5 ,lambda = 0, nsamp = 5000, nthreads = 1)Computation of (robust) Mahalanbosi distances.
MD <- as.numeric(mmd(X = X,
mu = par_MMLE$mu,
cov_row = par_MMLE$cov_row_inv,
cov_col = par_MMLE$cov_col_inv,
inverted = TRUE))
MD_rob <- as.numeric(mmd(X = X,
mu = par_MMCD$mu,
cov_row = par_MMCD$cov_row_inv,
cov_col = par_MMCD$cov_col_inv,
inverted = TRUE))
out_quant <- qchisq(0.95, p*q)
outliers <- which(MD_rob > out_quant)Distance-distance of yearly outlyingness scores, using the robust procedure, we can unmask several outliers.
names(MD_rob) <- dimnames(X)[[3]]
subs <- MD_rob < out_quant
plt_dd <- ggplot(data = NULL, aes(x = sqrt(MD), y = sqrt(MD_rob))) +
geom_point() +
geom_hline(yintercept = sqrt(out_quant)) +
geom_vline(xintercept = sqrt(out_quant)) +
geom_label_repel(aes(x = sqrt(MD[!subs]),
y = sqrt(MD_rob[!subs]),
label = names(MD_rob[!subs])),
size = 2.5, max.overlaps = 200, nudge_y = 0.5) +
theme_classic() +
xlim(c(min(sqrt(MD)),max(sqrt(MD))+1)) +
labs(x = "MMD", y = "Robust MMD")
plt_dd
Computation of cellwise Shapley values.
shv_cell <- array(dim = dim(X))
shv_cell[] <- matrixShapley(X,
mu = par_MMCD$mu,
cov_row = par_MMCD$cov_row_inv,
cov_col = par_MMCD$cov_col_inv,
inverted = TRUE, type = "cell")
dimnames(shv_cell) <- dimnames(X)Preparing data for plots.
long_weather <- function(x, val_to, var_levels = NULL) {
y <- x %>%
rownames_to_column(var = "var") %>%
pivot_longer(cols = -var, names_to = "parameter", values_to = val_to)
if(!is.null(var_levels)){
y <- y %>%
mutate(var = factor(var, levels = var_levels))
}
y
}
var_names <- c("T [°C]" = "T",
"SH [h]" = "SH",
"P [mm/m²]" = "P",
"AP [hPa]" = "AP",
"SP [%]" = "SP")
set <- cbind("var" = 1891:2022,
"Temperature [°C]" = 0,
"Preciptiation [mm/m²]" = 0,
"Solid Preciptiation [%]" = 0,
"Sunshine hours [h]" = 0,
"Air pressure [hPa]" = 0)
set[!(set[,1] %in% dimnames(X)[[3]]),-1] <- 1
set <- data.frame(set)
colnames(set) <- c("var",
"T [°C]",
"P [mm/m²]",
"SP [%]",
"SH [h]",
"AP [hPa]")
set_long <- set %>% pivot_longer(-c(var))shv_year <- apply(shv_cell, c(1,3),sum)
X_center <- array(dim = dim(X), dimnames = dimnames(X))
X_center[] <- (apply(X,3,function(x) x - par_MMCD$mu))
X_center_sign <- array(dim = dim(X), dimnames = dimnames(X))
X_center_sign[] <- sign(X_center)
sign_mat <- apply(X_center,c(1,3),function(x) sign(mean(x)))
sign_mat[,-outliers] <- 0
shv_rel <- apply(shv_year, 2, function(x) x/sum(x))
shv_rel_outliers <- shv_rel
shv_rel_outliers[,-outliers] <- 0
shv_rel_outliers[shv_rel_outliers < 0] <- 0 # ONLY PLOT POSITIVE SHAPLEY VALUES
shv_years_long <- long_weather(data.frame(t(shv_rel_outliers)), val_to = "shv_rel")
sign_years_long <- long_weather(data.frame(t(sign_mat)), val_to = "sign")
plt_data_years <- inner_join(shv_years_long,sign_years_long) %>%
mutate(val = sign*abs(shv_rel),
sign = as.character(sign),
var = as.numeric(var),
parameter = fct_recode(factor(parameter),!!! var_names))shv_long_1895 <- long_weather(data.frame(t(shv_cell[,,outliers[1]])), val_to = "shv", var_levels = dimnames(X)[[2]])
X_long_1895 <- long_weather(data.frame(t(X[,,outliers[1]])) %>% mutate(SP = SP*100), val_to = "x", var_levels = dimnames(X)[[2]])
sign_long_1895 <- long_weather(data.frame(t(sign(X[,,outliers[1]] - par_MMCD$mu))), val_to = "sign", var_levels = dimnames(X)[[2]])
plt_data_1895 <- inner_join(shv_long_1895,sign_long_1895) %>%
mutate(val = sign*abs(shv),
sign = as.character(sign)) %>%
inner_join(X_long_1895) %>%
mutate(parameter = fct_recode(factor(parameter),!!! var_names))
shv_long_2022 <- long_weather(data.frame(t(shv_cell[,,n])), val_to = "shv", var_levels = dimnames(X)[[2]])
X_long_2022 <- long_weather(data.frame(t(X[,,n])) %>% mutate(SP = SP*100), val_to = "x", var_levels = dimnames(X)[[2]])
sign_long_2022 <- long_weather(data.frame(t(sign(X[,,n] - par_MMCD$mu))), val_to = "sign", var_levels = dimnames(X)[[2]])
plt_data_2022 <- inner_join(shv_long_2022,sign_long_2022) %>%
mutate(val = sign*abs(shv),
sign = as.character(sign)) %>%
inner_join(X_long_2022) %>%
mutate(parameter = fct_recode(factor(parameter),!!! var_names))Yearly outlyingess contributions based on rowwise Shapley values. In the plot, non-outlying years are colored white, while years with missing data are gray. For the outliers, blue indicates “above average”, red indicates “below average”, and color intensity is proportional to the rowwise Shapley value.
ggplot(plt_data_years, aes(x = var, y = parameter, fill = val)) +
geom_tile(color = "lightgray") +
theme_minimal() +
scale_fill_gradient2(low = "#1F78B4", mid = "white", high = "#E31A1C",
midpoint = 0, guide = "none") +
scale_x_continuous(breaks = seq(from = 1880, to = 2020, by = 5), expand = c(0,0)) +
theme(axis.text.x = element_text(angle = 45, vjust = 1.2, hjust=1),
axis.title = element_blank(), title = element_blank()) +
geom_tile(data = set_long, aes(x = var, y = name, alpha = factor(value)), fill = "gray") +
scale_alpha_manual(values = c(0,1)) +
guides(alpha = FALSE, fill = FALSE)
Outlyingess contributions based on cellwise Shapley values are computed for 1895 and 2022 using the same color scheme as in the previous plot.
gridExtra::grid.arrange(ggplot(plt_data_1895, aes(x = var, y = parameter, fill = val, label = round(x))) +
geom_tile(color = "lightgray") +
theme_minimal() +
geom_text(size = 3) +
theme(axis.text.x = element_text(angle = 45, vjust = 1.4, hjust=1.4),
axis.title = element_blank())+
labs(title = "1895") +
scale_fill_gradient2(low = "#1F78B4", mid = "white", high = "#E31A1C",
midpoint = 0, guide = "none"),
ggplot(plt_data_2022, aes(x = var, y = parameter, fill = val, label = round(x))) +
geom_tile(color = "lightgray") +
theme_minimal() +
geom_text(size = 3) +
theme(axis.text.x = element_text(angle = 45, vjust = 1.4, hjust=1.4),
axis.title = element_blank())+
labs(title = "2022") +
scale_fill_gradient2(low = "#1F78B4", mid = "white", high = "#E31A1C",
midpoint = 0, guide = "none"),
nrow = 1)
Computing Shapley values and average proportional contributions
of the variables to the MMDs for the H and P
groups.