Experiment 2: Relationship between PC gates and output gates

Introduction

We need to decide whether our “gates” methodology still makes sense when applied in the principal components (PC) coordinate system. In this experiment we compare runs filtered by output gates to runs filtered by PC gates to see how well they correspond.

Setup

library('dplyr', warn.conflicts = FALSE)
library('ggplot2')
library('ggthemes')
library('foreach')
library('hectorcal')

Evaluation of models that pass the gates in the PC space

Concentration-driven runs

Comparison of ESM and Hector PC ranges

For the concentration-driven runs most of the variation is carried in the first two PCs. The box marked out by the ESMs in this two-dimensional space gives a first approximation to our calibration region. How do models in this region perform relative to our old-style gates defined in the model output space?

Although PC1 and PC2 account for most of the variation in the hector results, there is still some variation in the other components, so we’ll take a look at different levels for PC 3 and 4 as well.

## We need the PC projections for the ensemble
hector_proj <- compute_pc(hector_conc_ensemble, c('historical','rcp26','rcp45','rcp60','rcp85'), 
                          'tas', 2006:2100, 1861:2005, retx=TRUE)$x
## ESM projections are stored as package data: cmip_conc_pcproj

## Tables of min, max, mean values by PC, for both sets
maxpc <- 10
esm_pc_stats <- group_by(cmip_conc_pcproj, PC) %>% 
    summarise(min=min(value), max=max(value), mean=mean(value), sd=sd(value)) %>%
    filter(PC <= maxpc)
stats_matrix <- apply(hector_proj[ , 1:maxpc], 2, function(x) {c(min(x), max(x), mean(x), sd(x))})
hector_pc_stats <- data.frame(PC=seq_along(stats_matrix[1,]), 
                              min=stats_matrix[1,],
                              max=stats_matrix[2,],
                              mean=stats_matrix[3,],
                              sd=stats_matrix[4,])
bind_cols(select(esm_pc_stats, PC, esm.min=min, esm.max=max, esm.mean=mean, esm.sd=sd),
          select(hector_pc_stats, hector.min=min, hector.max=max, hector.mean=mean, hector.sd=sd))

The first thing we can see from these statistics is that though the ESM statistics exclude a lot of the range covered by the Hector ensemble, the ensemble average passes all of the gates for the first 10 PCs, except for PC7, which marginally excludes zero. For PCs \(\geq\) 9, the coefficients of all of the hector ensemble members are all very close to zero, indicating that for this ensemble there is definitely no useful calibration to be done using those components because for any such test either the enitre parameter space will pass, or the entire space will fail.

We are now in a position to evaluate how the joint PC1-PC2 gate deprojects into the output plane. Here is a graphical depiction of the gate, with the positions of the ESMs labeled.

esm_pcplot(cmip_conc_pcproj) + geom_abline(slope=-0.6, intercept=-3, color='lightgrey', size=1.5, alpha=0.5)

One noticeable trend is that although these dimensions are uncorrelated in the Hector ensemble, the ESMs appear to cluster around a line with a slope of about \(-0.6\). The sample here is very small, so it’s not clear how seriously we should take this possible relationship. Is there no model around \((-7.5, -4)\) because that model is excluded by the ESMs’ calibration process, or because nobody happens to have written the model with that behavior? We are not in a position to answer that question right now, but it’s worth keeping in mind.

Evaluating the relationship between PCs and outputs

To get a sense of where the points in this plane fall in output space, we will define low, medium, and high levels for PC1 and PC2, defined by the min, mean, and max values, respectively, for the ESM projections. PC1 levels will be depicted by color in the plots, and PC2 levels will be depicted by line type. We will also define low, medium, and high levels for PC3 and PC4 using this time the min, mean, and max for the Hector ensemble, reflecting the fact that for now we are not planning to filter on those PCs. The PC3 and PC4 levels will be presented in facet plots. Finally, the results will be broken out by experiment, with each experiment presented as a separate family of plots.

idcols <- c('min','mean','max')
levels <- c('low','med','zhigh')
pc1lvls <- as.data.frame(esm_pc_stats)[1, idcols]
pc2lvls <- as.data.frame(esm_pc_stats)[2, idcols]
pc3lvls <- hector_pc_stats[3, idcols]
pc4lvls <- hector_pc_stats[4, idcols]
names(pc1lvls) <- names(pc2lvls) <- names(pc3lvls) <- names(pc4lvls) <- levels
### Create the tables for each combination of levels
climout <- 
    foreach(lpc1=levels, .combine=bind_rows) %do% {
        foreach(lpc2=levels, .combine=bind_rows) %do% {
            foreach(lpc3=levels, .combine=bind_rows) %do% {
                foreach(lpc4=levels, .combine=bind_rows) %do% {
                    coef <- c(pc1lvls[[lpc1]], pc2lvls[[lpc2]], pc3lvls[[lpc3]], pc4lvls[[lpc4]])
                    reconstruct_climate(coef, pc_conc, 4) %>%
                        mutate(PC1=lpc1, PC2=lpc2, PC3=lpc3, PC4=lpc4)
                }
            }
        }
    }
## Order the levels so that they are in order in the plots
#climout <- mutate(climout, 
#                  PC1 = ordered(PC1, labels=levels),
#                  PC2 = ordered(PC2, labels=levels),
#                  PC3 = ordered(PC3, labels=levels),
#                  PC4 = ordered(PC4, labels=levels))
## Get a data frame of gates for an experiment
get_gates <- function(expt, var, minyear=1861, maxyear=2100) {
    if(grepl('[Hh]istorical', expt)) {
        maxyear <- 2005
    }
    else {
        minyear <- 2006
    }
    
    exptdata <- filter(esm_comparison, experiment==expt, variable==var, 
                       year <= maxyear, year >= minyear)
    yrs <- c(min(exptdata$year),
             max(exptdata$year[exptdata$year <= median(exptdata$year)]),
             max(exptdata$year))
    filter(exptdata, year %in% yrs)
}
for(expt in unique(climout$experiment)) {
    print(
        ggplot(data=filter(climout, experiment==expt),
               aes(x=year, y=value, color=PC1, linetype=PC2)) +
            geom_line(size=0.75) +
            geom_errorbar(data=get_gates(expt, 'tas'), aes(x=year, ymin=mina, ymax=maxb), 
                          linetype=2, color='lightgrey', inherit.aes = FALSE, width=1) +
            facet_grid(PC3~PC4) + 
            ggtitle(paste('Experiment:', expt)) + ylab('Temperature') +
            theme_solarized_2(light=FALSE, base_size = 18) + 
            scale_color_solarized()
    )
}

Looking at the data this way makes it a lot more clear what each of the components is doing. For the future runs (the behavior for the historical runs is a bit different), PC1 shifts the temperature nearly uniformly up and down (look, for example, at the series of solid lines in the plots). Meanwhile, PC2 appears to control the average slope of the temperature curve over time (compare the sequence of lines in any single color).

PC3 (horizontal facets in the grid) appears to control the tendency of the temperature curve to bend downward in the late 21st century. Consider, for example, the RCP 2.6 experiment. With PC3 at the high level (right column), none of these scenarios show the peak and decline that we normally associate with RCP 2.6. PC3 also has a clear effect in the historical experiment, where it is clearly related to the strength of the aerosol effect.

The effect of PC4 (vertical facets) is a little hard to discern in the future runs. It looks to be a little duplicative of PC3, inducing some downward curvature late in the century. However, in the historical period PC4 is clearly influencing the behavior in the early part of the run. With PC4 low, scenarios tend to be below or at the low end of the 1933 gate. With PC4 medium, all scenarios go through the 1933 gate, irrespective of the other PC values, and with PC4 high, the scenarios tend to be high at the 1933 gate.

The conclusion from all of this is that notwithstanding the small fraction of the total variance accounted for by PC3 and PC4, if we want our PC methodology to correspond, at least roughly, to the output gate methodology, then we are going to have to constrain PC3 and PC4 using the ESM results.

If we constrain PC3 and PC4 with the ESM values and rerun the experiment, then we get the following.

pc1lvls_pc34 <- as.data.frame(esm_pc_stats)[1, idcols]
pc2lvls_pc34 <- as.data.frame(esm_pc_stats)[2, idcols]
pc3lvls_pc34 <- as.data.frame(esm_pc_stats)[3, idcols]
pc4lvls_pc34 <- as.data.frame(esm_pc_stats)[4, idcols]
names(pc1lvls_pc34) <- names(pc2lvls_pc34) <- 
    names(pc3lvls_pc34) <- names(pc4lvls_pc34) <- levels
### Create the tables for each combination of levels
climout_pc34 <- 
    foreach(lpc1=levels, .combine=bind_rows) %do% {
        foreach(lpc2=levels, .combine=bind_rows) %do% {
            foreach(lpc3=levels, .combine=bind_rows) %do% {
                foreach(lpc4=levels, .combine=bind_rows) %do% {
                    coef <- c(pc1lvls_pc34[[lpc1]], pc2lvls_pc34[[lpc2]], 
                              pc3lvls_pc34[[lpc3]], pc4lvls_pc34[[lpc4]])
                    reconstruct_climate(coef, pc_conc, 4) %>%
                        mutate(PC1=lpc1, PC2=lpc2, PC3=lpc3, PC4=lpc4)
                }
            }
        }
    }
for(expt in unique(climout_pc34$experiment)) {
    print(
        ggplot(data=filter(climout_pc34, experiment==expt),
               aes(x=year, y=value, color=PC1, linetype=PC2)) +
            geom_line(size=0.75) +
            geom_errorbar(data=get_gates(expt, 'tas'), aes(x=year, ymin=mina, ymax=maxb), 
                          linetype=2, color='lightgrey', inherit.aes = FALSE, width=1) +
            facet_grid(PC3~PC4) + 
            ggtitle(paste('Experiment:', expt)) + ylab('Temperature') +
            theme_solarized_2(light=FALSE, base_size = 18) + 
            scale_color_solarized()
    )
}

Not surprisingly, these constraints put most of the scenarios a lot closer to the output gates. We’re still missing sometimes, but at this point one starts to wonder, how firm are these gates. We don’t have that many ESMs, and so it’s not too hard to believe that scenarios a little outside of the gates are still reasonable. The only cases that are a little concerning are the RCP 8.5 scenarios for the (low/med, high, any, low) combinations. Some of those scenarios are \(1.5^{\circ}C\) below the lower edge of the 2100 gate, which seems a little much to me.

Note that for all of these scenarios, the projection onto PC5 is pegged at zero, which is within the ESM constraints, so constraining that parameter won’t affect anything in these plots. However, we might well wonder what happens to these results if PC5 is near the high or low end of its range. At the risk of inducing plot overload, here is the same set with PC5 low.

### -2.76 is the low value for PC5 (from the table above)
climout_pc5lo <- 
    foreach(lpc1=levels, .combine=bind_rows) %do% {
        foreach(lpc2=levels, .combine=bind_rows) %do% {
            foreach(lpc3=levels, .combine=bind_rows) %do% {
                foreach(lpc4=levels, .combine=bind_rows) %do% {
                    coef <- c(pc1lvls_pc34[[lpc1]], pc2lvls_pc34[[lpc2]], 
                              pc3lvls_pc34[[lpc3]], pc4lvls_pc34[[lpc4]], -2.76)
                    reconstruct_climate(coef, pc_conc, 5) %>%
                        mutate(PC1=lpc1, PC2=lpc2, PC3=lpc3, PC4=lpc4)
                }
            }
        }
    }
for(expt in unique(climout_pc5lo$experiment)) {
    print(
        ggplot(data=filter(climout_pc5lo, experiment==expt),
               aes(x=year, y=value, color=PC1, linetype=PC2)) +
            geom_line(size=0.75) +
            geom_errorbar(data=get_gates(expt, 'tas'), aes(x=year, ymin=mina, ymax=maxb), 
                          linetype=2, color='lightgrey', inherit.aes = FALSE, width=1) +
            facet_grid(PC3~PC4) + 
            ggtitle(paste('Experiment:', expt)) + ylab('Temperature') +
            theme_solarized_2(light=FALSE, base_size = 18) + 
            scale_color_solarized()
    )
}

Some of these are way outside the gates; the RCP 2.6 scenarios, especially, are very high at the 2100 gate.

### 0.382 is the high value for PC5 (from the table above)
climout_pc5hi <- 
    foreach(lpc1=levels, .combine=bind_rows) %do% {
        foreach(lpc2=levels, .combine=bind_rows) %do% {
            foreach(lpc3=levels, .combine=bind_rows) %do% {
                foreach(lpc4=levels, .combine=bind_rows) %do% {
                    coef <- c(pc1lvls_pc34[[lpc1]], pc2lvls_pc34[[lpc2]], 
                              pc3lvls_pc34[[lpc3]], pc4lvls_pc34[[lpc4]], 0.382)
                    reconstruct_climate(coef, pc_conc, 5) %>%
                        mutate(PC1=lpc1, PC2=lpc2, PC3=lpc3, PC4=lpc4)
                }
            }
        }
    }
for(expt in unique(climout_pc5hi$experiment)) {
    print(
        ggplot(data=filter(climout_pc5hi, experiment==expt),
               aes(x=year, y=value, color=PC1, linetype=PC2)) +
            geom_line(size=0.75) +
            geom_errorbar(data=get_gates(expt, 'tas'), aes(x=year, ymin=mina, ymax=maxb), 
                          linetype=2, color='lightgrey', inherit.aes = FALSE, width=1) +
            facet_grid(PC3~PC4) + 
            ggtitle(paste('Experiment:', expt)) + ylab('Temperature') +
            theme_solarized_2(light=FALSE, base_size = 18) + 
            scale_color_solarized()
    )
}

These actually don’t look too bad. However, given the poor performance at the low end, it’s pretty clear that we will need to constrain PC5.

We could continue this line of analysis, but if we look at the range for PC6, the ESM maximum doesn’t constrain the ensemble at all because the ensemble maximum is actually lower than the ESM maximum. There is some constraint at the low end, but it’s a much weaker constraint than we see for PC5, so I’m pretty confident that constraining PC6 won’t make much difference. Note also that for PCs 8-10 the ESM range includes the entire ensemble range, so constraining these components definitely won’t have any effect. Therefore, we will probably end up running MCMC runs with constraints on PC1-5, 1-6, and 1-7, and we’ll compare the parameter distributions that we get from these to see if there is any significant difference.

Emissions-driven runs

Reliability of the small concentration-driven ESM sample

Here we go again. My main concern here is that with only three valid models to use in the calibration, the statistics for our calibration data may be not be very good. In particular we could be seeing gates that are way too narrow if the models that we happen to have available to us aren’t near the edges of the true distribution of model outcomes. We can get some idea of the extent to which this happens by comparing the temperature output gates in the historical and RCP 8.5 experiments for the concentration and emissions driven ESM runs.

bind_rows(
    get_gates('historical', 'tas'),
    get_gates('esmHistorical', 'tas'))

bind_rows(
    get_gates('rcp85', 'tas'),
    get_gates('esmrcp85', 'tas'))

It looks as if the gates for the emissions-driven runs are at their most restrictive when we’re close to the historical-future change over. The gates for the 1861 and 2100 don’t look much smaller than their concentration-driven counterparts. However, we would expect the emissions-driven runs to show a greater dispersion in outcomes than the concentration- driven runs, since they have more dimensions along which they can be different from one another. That we don’t see that additional dispersion indicates that our gates in the emission-driven runs probably are narrower than they should be.

I am thinking that when we go to do the calibration we will need to correct for this by widening the likelihood function. We can do this either by increasing the width of the gates, or by increasing the smoothing parameter. I haven’t yet decided which way is better, but I am leaning toward using the smoothing parameter, since it better expresses the idea that we have softened the edge of the gate because we don’t really know where that edge is. How do we know how much to soften the gates? Here we can take advantage of the better sampling of the concentration-driven model space. If we take the marginal variance of the temperature parameters as being relatively reliably determined by the concentration- driven MCMC, then we can adjust the softening in the emissions-driven MCMC until the marginal variances of the temperature parameters are reasonably close to the values from the concentration-driven MCMC.

Comparison of Hector and ESM PC ranges

## We need the PC projections for the ensemble
hector_emiss_proj <- compute_pc(hector_emiss_ensemble, c('esmHistorical','esmrcp85'), 
                                c('tas', 'co2'), 2006:2100, 1861:2005, retx=TRUE)$x

## Tables of min, max, mean values by PC, for both sets
maxpc <- 20
esm_emiss_pc_stats <- group_by(cmip_emiss_pcproj, PC) %>% 
    summarise(min=min(value), max=max(value), mean=mean(value), sd=sd(value)) %>%
    filter(PC <= maxpc)
stats_matrix <- apply(hector_emiss_proj[ , 1:maxpc], 2, function(x) {c(min(x), max(x), mean(x), sd(x))})
hector_emiss_pc_stats <- data.frame(PC=seq_along(stats_matrix[1,]), 
                                    min=stats_matrix[1,],
                                    max=stats_matrix[2,],
                                    mean=stats_matrix[3,],
                                    sd=stats_matrix[4,])
bind_cols(select(esm_emiss_pc_stats, PC, esm.min=min, esm.max=max, esm.mean=mean, esm.sd=sd),
          select(hector_emiss_pc_stats, hector.min=min, hector.max=max, hector.mean=mean, hector.sd=sd))

Between the variance fraction analysis, the results from the concentration-driven runs, and the table above it seems likely that we are going to need to constrain at least the first 4 PCs, so I’ve dispensed with the version of the plot where PC3 and PC4 are unconstrained.

pc1lvls_emiss <- as.data.frame(esm_emiss_pc_stats)[1, idcols]
pc2lvls_emiss <- as.data.frame(esm_emiss_pc_stats)[2, idcols]
pc3lvls_emiss <- as.data.frame(esm_emiss_pc_stats)[3, idcols]
pc4lvls_emiss <- as.data.frame(esm_emiss_pc_stats)[4, idcols]
names(pc1lvls_emiss) <- names(pc2lvls_emiss) <- 
    names(pc3lvls_emiss) <- names(pc4lvls_emiss) <- levels
### Create the tables for each combination of levels
climout_emiss <- 
    foreach(lpc1=levels, .combine=bind_rows) %do% {
        foreach(lpc2=levels, .combine=bind_rows) %do% {
            foreach(lpc3=levels, .combine=bind_rows) %do% {
                foreach(lpc4=levels, .combine=bind_rows) %do% {
                    coef <- c(pc1lvls_emiss[[lpc1]], pc2lvls_emiss[[lpc2]], 
                              pc3lvls_emiss[[lpc3]], pc4lvls_emiss[[lpc4]])
                    reconstruct_climate(coef, pc_emiss, 4) %>%
                        mutate(PC1=lpc1, PC2=lpc2, PC3=lpc3, PC4=lpc4)
                }
            }
        }
    }
for(expt in unique(climout_emiss$experiment)) {
    for(var in unique(climout_emiss$variable)) {
        if(var=='tas') {
            yl <- 'Temperature'
        }
        else {
            yl <- 'CO2 concentration'
        }
        print(
            ggplot(data=filter(climout_emiss, experiment==expt, variable==var),
                   aes(x=year, y=value, color=PC1, linetype=PC2)) +
                geom_line(size=0.75) +
                geom_errorbar(data=get_gates(expt, var), aes(x=year, ymin=mina, ymax=maxb), 
                              linetype=2, color='lightgrey', inherit.aes = FALSE, width=1) +
                facet_grid(PC3~PC4) + 
                ggtitle(paste('Experiment:', expt)) + ylab(yl) +
                theme_solarized_2(light=FALSE, base_size = 18) + 
                scale_color_solarized()
        )
    }
}

The effects of the various components are a little harder to tease out here, partly because we have two variables to look at now, and partly because we don’t have the lower emissions scenarios to work with. It does look like PC1 and PC2 are felt mainly in the temperature, while PC3 and PC4 seem to affect both variables.

With just these four components constrained, we are actually getting pretty good agreement with the output gates. Looking at the table of principal component stats suggests that the ESMs will provide substantial constraints up through PC7, but for the reasons discussed above, it’s not clear how seriously we should take this with only three ESMs in the sample. Notice that for PC8 the ESMs span a very narrow range of coefficients. This seems almost certain to be an artifact of the small sample size. Therefore, I’m thinking our experimental range for the number of PCs to constrain in the emissions-driven runs should be 5-7.

It would be nice if we could come up with an objective formula for deciding the number of PCs to constrain in calibration. The cumulative variance fraction for the PCA doesn’t seem like a very good guide. I was thinking of something related to the variance in the ESM projections, but that doesn’t look very promising either. For now, the best answer seems to be to look at the output performance relative to the ESMs (as we’ve done above), but that analysis gets murky as we go to higher dimensions.

Evaluation of models that pass the gates in output space

What we want to determine here is, does the analysis in the PCA space miss any models that we would have deemed good when we looked at them in the output space? Stay tuned.