R Notebook

Objective

Now that we have heat flux data we are ready to start fitting Hector runs to the ESM data! I started with the 24 models that had both temperature and heat flux data for each concentration driven experiment / ensemble member. I think that some of the results look promising but other fits are a bit um wonky.

library(gridExtra)
library(dplyr)
library(tidyr)
library(ggplot2)
library(rlist)
BASE_DIR <-  "/Users/dorh012/Documents/2019/hectorcal"
INPUT_DIR <- file.path(BASE_DIR, 'analysis', 'best_fit', 'temp-heat-1')
files <- list.files(INPUT_DIR, '.rds', full.names = TRUE) 
model_names <- gsub(pattern = '.rds', replacement = '', x = basename(files))    
    
comparison_plots <- lapply(files, function(input){
    
    x <- readRDS(input)
    
    if('message' %in% names(x)){
        
        NULL
    } else {
        
        x[[1]]
    }
    
    })
names(comparison_plots) <- model_names
did_not_converge <- model_names[which(unlist(lapply(comparison_plots, is.null)))]
comparison_plots <- list.clean(comparison_plots, fun = is.null)
residual_plots <- lapply(files, function(input){    
    x <- readRDS(input)
    
    if('message' %in% names(x)){
        
        NULL
    } else {
        
        x[[2]] 
    }})
names(residual_plots) <- model_names
residual_plots <- list.clean(residual_plots, fun = is.null)
MSE <- bind_rows(mapply(function(x, y){
    
    input <- readRDS(x)
    
    if('message' %in% names(input)){
        
        NULL
    } else {
        
        input[[3]] %>%  
            mutate(name = y)
    }}, x = files, y = model_names))
fit_params <- bind_rows(mapply(function(x, y){
    
    input <- readRDS(x)
   if('message' %in% names(input)){
        
        NULL
    } else {
        
       m <- matrix(input[[4]]$par, nrow = 1)  
       colnames(m) <- names(input[[4]]$par)
       tibble::as_tibble(m) %>% 
           bind_cols(value = input[[4]]$value, count = input[[4]]$counts[['function']]) %>% 
           mutate(name = y)
    }}, x = files, y = model_names))

Plot the residuals between the Hector and ESM data for heat flux for each model.

# Format the residual pltos to 
residual_plots %>% 
    lapply(function(input){
        
        input + 
            theme_bw(base_size = 10) + 
            theme(legend.position = 'none') + 
            labs(y = NULL, x = NULL) 
    }) -> 
    to_plot
grid.arrange(grobs = to_plot, ncol = 2)

Hmmm… so there are larger patterns in the residuals than I was expecting. Previously the only pattern we had seen in the residuals was that there was a greater spread in the residuals during the historical experiment that narrowed down in the future scenarios but was centered around 0. We still see that in some models (CCSM4) but then there is also all sorts of wonky behavior going on in the residuals for the other fits (FGOALS).

Comparison of the fits

# Format the residual pltos to 
comparison_plots %>% 
    lapply(function(input){
        
        other_model <- unique(pull(filter(input$data, model != 'hector'), model))
        
        colors <- c('grey', 'blue')
        names(colors) <- c(other_model, 'hector')
        
        alpha <- c(0, 1)
        names(alpha) <- c(other_model, 'hector')
        
        input + 
            theme_bw(base_size = 10) + 
            theme(legend.position = 'none') + 
            labs(y = NULL, x = NULL) + 
            scale_color_manual(values = colors) 
    }) -> 
    to_plot
grid.arrange(grobs = to_plot, ncol = 2)

Well it looks like we have some not so great fits. Like what is going on with incm4, I have no idea.

1. Could Hector be struggling to capture the variability in the Heat Flux data that we don’t really expect for it to be able to capture?
2. We had to reduce the number of experiments / ensemble members being used as comparison data so that there was parity for the heat flux and temperature data.

Should we include more temp data than heat flux data?

What do the paramters look like?

fit_params %>%
    gather(param, param_value, S, diff, alpha, volscl) -> 
    params_long
params_long %>% 
    dplyr::group_by(param) %>%  
    dplyr::summarise(min = min(param_value), 
              mean = mean(param_value), 
              max = max(param_value), 
              sd = sd(param_value)) %>% 
    knitr::kable(digits = 3) %>%  
    kableExtra::kable_styling(bootstrap_options = c("striped"))

param	min	mean	max	sd
alpha	-21.462	-0.892	1.722	4.666
diff	0.996	104.819	1553.504	325.005
S	0.232	8.614	139.990	28.676
volscl	-6.465	2.976	60.251	12.669

So it looks like we are still getting some unrealistic parameter values. I wonder if they are for the models that only contain the historical runs? Then I suspect that the heat flux data is unable to help resolve the issues of identifying S and diff.

# Determine which models only have historical data. 
MSE %>%
    separate(experiment, c("cmip_experiment", "ensemble"), sep = '_') %>%  
    mutate(model = gsub('conc_', '', name)) %>% 
    select(cmip_experiment, model) %>% 
    distinct %>% 
    mutate(historical = if_else(cmip_experiment == 'historical', 1, 0)) %>%  
    group_by(model) %>%  
    summarise(exp_count = n(), 
              historical = sum(historical)) %>%  
    filter(exp_count == 1 & historical == 1) %>% 
    select(model) %>% 
    mutate(only_history = TRUE) -> 
    history_id_tibble
params_long %>%  
        mutate(model = gsub('conc_', '', name)) %>% 
    left_join(history_id_tibble) -> 
    to_plot
ggplot(data = to_plot) + 
    geom_dotplot(aes(param_value, fill = only_history), stackgroups = TRUE, binpositions = "all") + 
    facet_wrap(~param, scales = 'free')

Hmmm that was not the pattern I was expecting.

Take Aways

Well we have both not so great parameter values and Hector fits. Do we want to increase the amount of temperature data being used in the fits? Do we want to limit heat flux comparison data to year 2100? Or do we want to smooth out the heat flux data so Hector is not being contorted to fit the peaks (I think that is where some of the off volscl and alpha values are coming from).