Settings

library(tidyverse)
library(multidplyr)
library(zoo)
library(slider)
library(pROC)

plot_theme <- theme_bw() +
    theme(text = element_text(size=16))

cores <- parallel::detectCores() - 4
cluster <- new_cluster(cores)

window <- 14
vars <- c("RD", "VPD", "T", "SOIL", "AWSSM", "RAIN", "AET", "D", "GDD")

region_name <- "Southern Rockies"

settings <- tribble(
    ~region, ~params,
    "Southern Rockies", list(data_dir =  "/home/steve/OneDrive/nothern-rockies-dryness",
                           out_dir =  "/media/steve/storage/waterbalance/fire/out/ap23/",
                           master_csv =  "SR_fire_points_extracted.csv"
                           ),
    "Middle Rockies", list(data_dir =  "/home/steve/OneDrive/nothern-rockies-dryness",
                           out_dir =  "/media/steve/storage/waterbalance/fire/out/ap23/",
                           master_csv =  "MR_fire_points_extracted.csv"
                           ),
)


settings <- settings %>% filter(region == region_name) 

## Directory containing metadata file
data_dir <- settings$params[[1]]$data_dir

## Directory with csv files from merge-csv.R
out_dir <- settings$params[[1]]$out_dir

master_csv <- settings$params[[1]]$master_csv

Southern Rockies

Preprocessing

Get MTBS database

Extract fire polys by location of interest (Ecoregion)
Create centroids for each poly (Pole of inaccessability in QGIS)

Get LANDFIRE EVT (Existing Vegetation Type data)

Use statistics tool in QGIS to get mean, median, and mode EVT by pixel count. Only mode is used going further.

Create master csv file with points of interest. Needs columns ‘Event_ID’, ‘Incid_Type’, ‘Ig_Date’, ‘x’ (=Longitude), ‘y’ (=Latitiude), and ‘LF20_EVT_5’ (Landfire 2020 EVT type). Further analysis can look at finer resolution LANDFIRE vegetation types offered.
Get Historical WB data for each centroid from Mike
Get Historical gridmet climate data for Each centroid

gridmet-historical.R - download gridmet data for each point
- pr, tmmx, tmmn, rmax, rmin, vpd
- creates one dir ‘gridmet-$var’ for each variable
merge-csv.R
- join Mike’s data and historical gridmet data into one csv for each Event_ID
- use these csvs for input to this script

Master dataframe

Read in file containing meta data for fires of interest.

setwd(data_dir)

master <- read_csv(master_csv) %>%
    filter(Incid_Type == "Wildfire") %>%
##    filter(year(Ig_Date) <= 2014) %>%
    rename("maj_veg_cl" = LF20_EVT_5,
           "Acres" = BurnBndAc) %>%
    ##  maj_class is the majority
    ##    of pixels classification from Landfire 30m pixels in burned
    ##    areas, this data was extracted in QGIS.
    mutate(maj_veg_cl = case_match(maj_veg_cl,
                                   c("Herb", "Shrub", "Sparse") ~ "non-forest",
                                   "Tree" ~ "forest",
                                   c("Agriculture", "Water") ~ "other"),
           filename = paste0(Event_ID, ".csv"),
           doy = yday(Ig_Date))

## Rows: 290 Columns: 46
## ── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (18): Event_ID, irwinID, Incid_Name, Incid_Type, Map_Prog, Asmnt_Type, Pre_ID, Post_ID, Comment, LF20_EVT_2, LF20_EVT_4, LF20_EVT_...
## dbl  (27): Map_ID, BurnBndAc, BurnBndLat, BurnBndLon, Perim_ID, dNBR_offst, dNBR_stdDv, NoData_T, IncGreen_T, Low_T, Mod_T, High_T, _EV...
## date  (1): Ig_Date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Map

library(sp)
library(sf)

coords <- data.frame(
    x = master$x,
    y = master$y
)
spdf <- SpatialPointsDataFrame(coords = coords, data = master)

ecoregions <- st_transform(read_sf("/home/steve/OneDrive/nothern-rockies-dryness/data/us_eco_l3/us_eco_l3.shp"), crs = 4326)
ecoregion <- ecoregions %>% filter(US_L3NAME == region_name)

ggplot(ecoregion) +
    geom_sf() +
    geom_point(data = master, aes(x, y, color = maj_veg_cl))

summary_tbl <- master %>% group_by(maj_veg_cl) %>% summarise(n = n())
knitr::kable(summary_tbl)

maj_veg_cl	n
forest	153
non-forest	86

Determine Fire Season

forest<-subset(master, maj_veg_cl=="forest")
forest_q<-quantile(forest$doy,c(0,1))#90% of data fall within this range 0.05-0.95# 0.27 to 0.75 or 0,1 to determine min max

nonforest<-subset(master, maj_veg_cl=="non-forest")
nonforest_q<-quantile(nonforest$doy,c(0,1))#90% of data fall within this range 0.05-0.95# 0.27 to 0.75 or 0,1 to determine min

#start and end dates of fire season by vegetation type derived from day of year in min_max table
doy_min_max<-as.data.frame(rbind(forest_q, nonforest_q))
forest_start<-doy_min_max[1,1]
nonforest_start<-doy_min_max[2,1]
forest_end<-doy_min_max[1,2]
nonforest_end<-doy_min_max[2,2]

knitr::kable(doy_min_max)

	0%	100%
forest_q	79	303
nonforest_q	7	301

Compute Percentiles

cluster_copy(cluster, "master")
cluster_copy(cluster, "window")
cluster_copy(cluster, "vars")
cluster_copy(cluster, "forest_start")
cluster_copy(cluster, "nonforest_start")
cluster_copy(cluster, "forest_end")
cluster_copy(cluster, "nonforest_end")
cluster_library(cluster, 'tidyverse')
cluster_library(cluster, 'slider')
cluster_library(cluster, 'zoo')


files <- paste0(out_dir, master$filename)
cluster_assign_partition(cluster, files = files)
cluster_send(cluster, wbdata <- vroom::vroom(files, id = "PATH"))
wbdata <- party_df(cluster, "wbdata")
master <- party_df(cluster, "master")

wbdata <- wbdata %>%
    mutate(Event_ID = str_extract(PATH, '(?<=/)[^/]*(?=\\.[^./]+$)')) %>% ## Extract file name without extension, use as Event_ID
    group_by(Event_ID) %>%
    mutate(RD = 100 - RH,
           WHC = max(SOIL), ## Guestimate soil WHC from max soil moisture observed
           AWSSM = WHC - SOIL) %>%
    left_join(select(master,
                     Event_ID, Ig_Date, maj_veg_cl, Acres),
              by = "Event_ID") %>%
    mutate(fire = if_else(Date == Ig_Date, 1, 0)) 

start <- proc.time() # Start clock
wbdata_smoothed <- select(wbdata, c(Event_ID, Ig_Date, maj_veg_cl, Acres, fire, Date, vars)) %>%
    ## mutate_at(vars[1:5],
    ##           ~ slide_dbl(.x, mean, .before = window, .after = 0, .complete = TRUE)) %>%
    ## mutate_at(vars[6:9],
    ##           ~ slide_dbl(.x, sum, .before = window, .after = 0, .complete = TRUE)) %>%
    mutate_at(vars[1:5],
              ~ rollapply(.x, window, mean, by = 1, partial = FALSE, fill = NA, align = "right")) %>%
    mutate_at(vars[6:9],
              ~ rollapply(.x, window, sum, by = 1, partial = FALSE, fill = NA, align = "right")) 
time_elapsed_parallel <- proc.time() - start # End clock

wbdata_percentiles_forest <- select(wbdata_smoothed, c(Event_ID, Ig_Date, maj_veg_cl, Acres, fire, Date, vars)) %>%
    filter(maj_veg_cl == "forest") %>%
    filter(Date <= Ig_Date,
           yday(Date) >= forest_start, ## Filter to fire season
           yday(Date) <= forest_end) %>%  ## To get this to match results from original dryness_script2, we only filter beginning of season.  Original script also didn't use doy, instead it filtered by absolute date, so would only chop off first pre-fireseason from first year in record.
##    slice(window + 1):n()) %>%
    collect() %>%
    drop_na(vars) %>% ### Drop NAs from periods not included in rolling calcs
    mutate(across(vars,
                  percent_rank))

wbdata_percentiles_nonforest <- select(wbdata_smoothed, c(Event_ID, Ig_Date, maj_veg_cl, Acres, fire, Date, vars)) %>%
    filter(maj_veg_cl == "non-forest") %>%
    filter(Date <= Ig_Date,
           yday(Date) >= nonforest_start,  ## Filter to fire season
           yday(Date) <= nonforest_end) %>%
    ##        slice((window + 1):n()) %>%
    collect() %>%
    drop_na(vars) %>% ### Drop NAs from periods not included in rolling calcs
    mutate(across(vars,
                  percent_rank)) 

ignitions_forest <- wbdata_percentiles_forest %>% filter(fire == 1)
nonignitions_forest <- wbdata_percentiles_forest %>% filter(fire == 0)
ignitions_forest_long <- ignitions_forest %>% pivot_longer(cols = vars)
ignitions_nonforest <- wbdata_percentiles_nonforest %>% filter(fire == 1)
nonignitions_nonforest <- wbdata_percentiles_nonforest %>% filter(fire == 0)
ignitions_nonforest_long <- ignitions_nonforest %>% pivot_longer(cols = vars)
ignitions_long <- rbind(ignitions_forest_long, ignitions_nonforest_long)

wbdata_smoothed_c <- wbdata_smoothed %>% collect()
master_c <- master %>% collect()

Data Visualizations

wbdata_percentiles_forest %>%
    filter(Event_ID %in% master_c$Event_ID[1:25]) %>%
    ggplot() +
    geom_point(aes(x = Date, y = Event_ID)) +
    geom_point(aes(x = Ig_Date, y = Event_ID), col = "red") +
    ggtitle("Dates included in analysis") +
    plot_theme

wbdata_smoothed_c %>%
    filter(fire == 1) %>%
    pivot_longer(vars) %>%
    ggplot(mapping = aes(x = Ig_Date, y = value, color = maj_veg_cl)) +
    geom_smooth(method = "lm") +
    geom_point() +
    facet_wrap(~name, scales = "free") +
    ggtitle(paste(region_name, "water balance on ignition date"))

## `geom_smooth()` using formula = 'y ~ x'

Different Rolling Windows

### TODO
### Refactor this so its not copy pasted from above

windows <- c(1, 3, 5, 7, 9, 11, 14, 17, 21, 31)

forest_roc_table_full <- tribble(~window, ~var, ~auc, ~auc20, ~auc10)
nonforest_roc_table_full <- tribble(~window, ~var, ~auc, ~auc20, ~auc10)

for (k in windows) {
    cluster_copy(cluster, "k")
    start <- proc.time() # Start clock
    wbdata_smoothed2 <- select(wbdata, c(Event_ID, Ig_Date, maj_veg_cl, Acres, fire, Date, vars)) %>%
        ## mutate_at(vars[1:5],
        ##           ~ slide_dbl(.x, mean, .before = window, .after = 0, .complete = TRUE)) %>%
        ## mutate_at(vars[6:9],
        ##           ~ slide_dbl(.x, sum, .before = window, .after = 0, .complete = TRUE)) %>%
        mutate_at(vars[1:5],
                  ~ rollapply(.x, k, mean, by = 1, partial = FALSE, fill = NA, align = "right")) %>%
        mutate_at(vars[6:9],
                  ~ rollapply(.x, k, sum, by = 1, partial = FALSE, fill = NA, align = "right")) 
    time_elapsed_parallel <- proc.time() - start # End clock

    wbdata_percentiles_forest2 <- select(wbdata_smoothed2, c(Event_ID, Ig_Date, maj_veg_cl, Acres, fire, Date, vars)) %>%
        filter(maj_veg_cl == "forest") %>%
        filter(Date <= Ig_Date,
               yday(Date) >= forest_start, ## Filter to fire season
               yday(Date) <= forest_end) %>%  ## To get this to match results from original dryness_script2, we only filter beginning of season.  Original script also didn't use doy, instead it filtered by absolute date, so would only chop off first pre-fireseason from first year in record.
        ##    slice(window + 1):n()) %>%
        collect() %>%
        drop_na(vars) %>% ### Drop NAs from periods not included in rolling calcs
        mutate(across(vars,
                      percent_rank))

    wbdata_percentiles_nonforest2 <- select(wbdata_smoothed2, c(Event_ID, Ig_Date, maj_veg_cl, Acres, fire, Date, vars)) %>%
        filter(maj_veg_cl == "non-forest") %>%
        filter(Date <= Ig_Date,
               yday(Date) >= nonforest_start,  ## Filter to fire season
               yday(Date) <= nonforest_end) %>%
        ##        slice((window + 1):n()) %>%
        collect() %>%
        drop_na(vars) %>% ### Drop NAs from periods not included in rolling calcs
        mutate(across(vars,
                      percent_rank))

    forest_rd_roc<-roc(wbdata_percentiles_forest2$fire,wbdata_percentiles_forest2$RD, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    forest_vpd_roc<-roc(wbdata_percentiles_forest2$fire,wbdata_percentiles_forest2$VPD, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    forest_t_roc<-roc(wbdata_percentiles_forest2$fire,wbdata_percentiles_forest2$T, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    forest_soil_roc<-roc(wbdata_percentiles_forest2$fire,wbdata_percentiles_forest2$SOIL, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    forest_gdd_roc<-roc(wbdata_percentiles_forest2$fire,wbdata_percentiles_forest2$GDD, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    forest_awssm_roc<-roc(wbdata_percentiles_forest2$fire,wbdata_percentiles_forest2$AWSSM, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    forest_rain_roc<-roc(wbdata_percentiles_forest2$fire,wbdata_percentiles_forest2$RAIN, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    forest_aet_roc<-roc(wbdata_percentiles_forest2$fire,wbdata_percentiles_forest2$AET, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    forest_d_roc<-roc(wbdata_percentiles_forest2$fire,wbdata_percentiles_forest2$D, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)

    png(filename = paste0("img/", region_name, "-forest-", k, "-days.png"), width = 800, height = 800)
    par(pty = "s") ## pty sets the aspect ratio of the plot region. Two options:
    ##                "s" - creates a square plotting region
    ##                "m" - (the default) creates a maximal plotting region
    par(bg="white")
    plot.roc(col="pink",forest_rd_roc, lwd=4, print.auc=FALSE, main = paste(region_name, "Forest ", k, " days"))
    plot.roc(add=TRUE,col="dark blue",forest_vpd_roc, main="vpd", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="green",forest_t_roc, main="temp", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="red",forest_soil_roc, main="soil", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="brown",forest_gdd_roc, main="gdd", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="magenta",forest_awssm_roc, main="awssm", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="purple",forest_rain_roc, main="rain", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="dark green",forest_aet_roc, main="aet", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="orange",forest_d_roc, main="d", lwd=4, print.auc=FALSE)
                                        # Add a legend
    legend(35, 75, legend=c(paste("RD ",round(auc(forest_rd_roc),1)),
                                        #paste("SVP ",round(auc(forest_svp_roc),1)),
                            paste("VPD ",round(auc(forest_vpd_roc),1)),
                            paste("TEMP ",round(auc(forest_t_roc),1)),
                            paste("SOIL ",round(auc(forest_soil_roc),1)),
                            paste("GDD ",round(auc(forest_gdd_roc),1)),
                            paste("SWD",round(auc(forest_awssm_roc),1)),
                            paste("RAIN ",round(auc(forest_rain_roc),1)),
                            paste("AET ",round(auc(forest_aet_roc),1)),
                            paste("CWD ",round(auc(forest_d_roc),1))),
           bty="n", col=c("pink", "dark blue","green","red","brown","magenta","purple","dark green","orange"), #col=c("black","blue", "dark blue","green","red","brown","magenta","purple","dark green","orange")
           lty=1, lwd=4, text.font=2)#cex=0.9,
    text(1,75,"AUC")
    dev.off()
    
    auc0f <- function(var) {
        roc <- paste0("forest_", var, "_roc")
        return(as.numeric(auc(get(roc))))
    }
    auc10f <- function(var) {
        roc <- paste0("forest_", var, "_roc")
        return(as.numeric(auc(get(roc),  partial.auc=c(100, 90))))
    }
    auc20f <- function(var) {
        roc <- paste0("forest_", var, "_roc")
        return(as.numeric(auc(get(roc),  partial.auc=c(100, 80))))
    }

    for (var in vars) {
        forest_roc_table_full <- forest_roc_table_full %>% add_row(window = k, var = var, auc = auc0f(tolower(var)), auc20 = auc20f(tolower(var)), auc10 = auc10f(tolower(var)))
    }

    nonforest_rd_roc<-roc(wbdata_percentiles_nonforest2$fire,wbdata_percentiles_nonforest2$RD, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    nonforest_vpd_roc<-roc(wbdata_percentiles_nonforest2$fire,wbdata_percentiles_nonforest2$VPD, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    nonforest_t_roc<-roc(wbdata_percentiles_nonforest2$fire,wbdata_percentiles_nonforest2$T, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    nonforest_soil_roc<-roc(wbdata_percentiles_nonforest2$fire,wbdata_percentiles_nonforest2$SOIL, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    nonforest_gdd_roc<-roc(wbdata_percentiles_nonforest2$fire,wbdata_percentiles_nonforest2$GDD, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    nonforest_awssm_roc<-roc(wbdata_percentiles_nonforest2$fire,wbdata_percentiles_nonforest2$AWSSM, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    nonforest_rain_roc<-roc(wbdata_percentiles_nonforest2$fire,wbdata_percentiles_nonforest2$RAIN, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    nonforest_aet_roc<-roc(wbdata_percentiles_nonforest2$fire,wbdata_percentiles_nonforest2$AET, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)
    nonforest_d_roc<-roc(wbdata_percentiles_nonforest2$fire,wbdata_percentiles_nonforest2$D, plot=FALSE, legacy.axes=TRUE, percent=TRUE, col="#377eb8", lwd=4, print.auc=TRUE)

    png(filename = paste0("img/", region_name, "-nonforest-", k, "-days.png"), width = 800, height = 800)
    par(pty = "s") ## pty sets the aspect ratio of the plot region. Two options:
    ##                "s" - creates a square plotting region
    ##                "m" - (the default) creates a maximal plotting region
    par(bg="white")
    plot.roc(col="pink",nonforest_rd_roc, lwd=4, print.auc=FALSE, main = paste(region_name, "Nonforest ", k, " days"))
    plot.roc(add=TRUE,col="dark blue",nonforest_vpd_roc, main="vpd", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="green",nonforest_t_roc, main="temp", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="red",nonforest_soil_roc, main="soil", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="brown",nonforest_gdd_roc, main="gdd", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="magenta",nonforest_awssm_roc, main="awssm", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="purple",nonforest_rain_roc, main="rain", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="dark green",nonforest_aet_roc, main="aet", lwd=4, print.auc=FALSE)
    plot.roc(add=TRUE,col="orange",nonforest_d_roc, main="d", lwd=4, print.auc=FALSE)
                                        # Add a legend
    legend(35, 75, legend=c(paste("RD ",round(auc(nonforest_rd_roc),1)),
                                        #paste("SVP ",round(auc(nonforest_svp_roc),1)),
                            paste("VPD ",round(auc(nonforest_vpd_roc),1)),
                            paste("TEMP ",round(auc(nonforest_t_roc),1)),
                            paste("SOIL ",round(auc(nonforest_soil_roc),1)),
                            paste("GDD ",round(auc(nonforest_gdd_roc),1)),
                            paste("SWD",round(auc(nonforest_awssm_roc),1)),
                            paste("RAIN ",round(auc(nonforest_rain_roc),1)),
                            paste("AET ",round(auc(nonforest_aet_roc),1)),
                            paste("CWD ",round(auc(nonforest_d_roc),1))),
           bty="n", col=c("pink", "dark blue","green","red","brown","magenta","purple","dark green","orange"), #col=c("black","blue", "dark blue","green","red","brown","magenta","purple","dark green","orange")
           lty=1, lwd=4, text.font=2)#cex=0.9,
    text(1,75,"AUC")
    dev.off()
    
    auc0nf <- function(var) {
        roc <- paste0("nonforest_", var, "_roc")
        return(as.numeric(auc(get(roc))))
    }
    auc10nf <- function(var) {
        roc <- paste0("nonforest_", var, "_roc")
        return(as.numeric(auc(get(roc),  partial.auc=c(100, 90))))
    }
    auc20nf <- function(var) {
        roc <- paste0("nonforest_", var, "_roc")
        return(as.numeric(auc(get(roc),  partial.auc=c(100, 80))))
    }

    for (var in vars) {
        nonforest_roc_table_full <- nonforest_roc_table_full %>% add_row(window = k, var = var, auc = auc0nf(tolower(var)), auc20 = auc20nf(tolower(var)), auc10 = auc10nf(tolower(var)))
    }
}

knitr::kable(forest_roc_table_full)
knitr::kable(nonforest_roc_table_full)

forest_roc_table_full %>% pivot_longer(c(auc, auc20, auc10)) %>%
    ggplot() +
    geom_line(mapping = aes(x = window, y = value, color = var), lwd=1) +
    facet_wrap(~name, ncol = 1) +
    ggtitle(paste(region_name, "Forest AUC and pAUC")) +
    plot_theme + theme(legend.position="bottom")

nonforest_roc_table_full %>% pivot_longer(c(auc, auc20, auc10)) %>%
    ggplot() +
    geom_line(mapping = aes(x = window, y = value, color = var), lwd=1) +
    facet_wrap(~name, ncol = 1) +
    ggtitle(paste(region_name, "Nonforest AUC and pAUC")) +
    plot_theme + theme(legend.position="bottom")

Corrplot

library(GGally)

ggcorr(wbdata_percentiles_forest, method = c("everything", "pearson"))

## Warning in ggcorr(wbdata_percentiles_forest, method = c("everything", "pearson")): data in column(s) 'Event_ID', 'Ig_Date', 'maj_veg_cl',
## 'Date' are not numeric and were ignored

ROC Analysis

We use the AUC (area under the curve) of the ROC (receiver operating characteristic) as a measure of goodness-of-fit to evaluate which predictors to use in the fire danger model.

window	var	auc	auc20	auc10
14	RD	83.51257	9.7954646	3.3832363
14	VPD	81.19981	8.6227823	2.9710010
14	T	71.19601	5.6287156	1.9312928
14	SOIL	74.92301	6.7278292	2.3106963
14	AWSSM	75.06209	6.7290251	2.3117579
14	RAIN	70.63633	4.2602364	1.1608188
14	AET	51.08774	0.8296199	0.1720074
14	D	80.84373	8.9306237	3.3207360
14	GDD	71.27445	5.6310038	1.9308678

window	var	auc	auc20	auc10
14	RD	86.99424	11.7117402	4.4909256
14	VPD	80.36174	10.0624360	3.9005013
14	T	73.21346	6.6842271	2.2391989
14	SOIL	79.60723	7.7959843	2.4582278
14	AWSSM	79.62269	7.7289973	2.3901110
14	RAIN	61.32740	1.0611233	0.2096575
14	AET	50.46663	0.5361996	0.0792144
14	D	81.61933	9.7841151	3.5618968
14	GDD	73.78741	6.6825864	2.2388618

demonstrate histogram of conditions associated with fire or no fire to see how much overlap there is. This helps interpret the ROC

Determine contribution of various parameters to burned area

## $gfs
##   [1] 0.000000000 0.004826023 0.002150181 0.005174935 0.002223450 0.018179481 0.015523795 0.029612524 0.020682854 0.004826135 0.005927640
##  [12] 0.004826762 0.018189864 0.024800696 0.029687974 0.065438083 0.005183065 0.002714121 0.019960089 0.025017411 0.029641161 0.025679772
##  [23] 0.005177632 0.018179491 0.023243572 0.029799860 0.058686848 0.019992445 0.024752822 0.029633130 0.025506645 0.024247270 0.048822544
##  [34] 0.026098924 0.031224772 0.022057549 0.033239298 0.005937369 0.005502943 0.023195769 0.040830947 0.029884469 0.067825169 0.005931279
##  [45] 0.019235918 0.026237603 0.033158166 0.067072621 0.023209017 0.040675105 0.029855083 0.067788006 0.048330110 0.056635528 0.082452703
##  [56] 0.034907392 0.067074928 0.065693613 0.005868741 0.022704229 0.037834780 0.030133448 0.060891226 0.020102631 0.025419005 0.030152186
##  [67] 0.026297163 0.027270237 0.048939654 0.027732882 0.035581673 0.026231703 0.036636553 0.022721960 0.037664766 0.030093986 0.060847278
##  [78] 0.045073009 0.055401116 0.074227759 0.033641124 0.060364822 0.058763196 0.027135148 0.048931720 0.027646935 0.035397614 0.026006938
##  [89] 0.036501292 0.049676260 0.026426300 0.049637814 0.035017618 0.006625369 0.023831552 0.044014235 0.033555721 0.069855687 0.023209028
## [100] 0.040868839 0.030429618 0.067874697 0.051560581 0.063697545 0.082972148 0.044032831 0.067857144 0.068194403 0.023847843 0.043850102
## [111] 0.033512590 0.069814269 0.049773220 0.060867465 0.084028619 0.037944878 0.068651945 0.067143860 0.051620556 0.063596410 0.082939036
## [122] 0.043874233 0.067823999 0.068160469 0.065287924 0.082480181 0.083155762 0.067276903 0.022722429 0.037897103 0.030692998 0.060971390
## [133] 0.047766541 0.061231164 0.074720486 0.041545768 0.060962667 0.061013524 0.027513177 0.049045079 0.028111485 0.036011166 0.026758151
## [144] 0.037194019 0.051065421 0.028013113 0.049975114 0.036645326 0.047805426 0.061155147 0.074694773 0.041378960 0.060926796 0.060971438
## [155] 0.061751804 0.074284973 0.075304467 0.060414393 0.051100317 0.027893506 0.049987427 0.036504398 0.049697412 0.023848047 0.044047944
## [166] 0.034111248 0.069907073 0.053860502 0.067000291 0.084371928 0.048086374 0.069855759 0.069977590 0.051650785 0.063782074 0.083144102
## [177] 0.044120947 0.067912968 0.068231799 0.065637356 0.084072422 0.084158056 0.068203174 0.053923495 0.066910044 0.084346126 0.047915525
## [188] 0.069814275 0.069938419 0.068951114 0.084049479 0.085261275 0.068698123 0.065600461 0.083987452 0.084106510 0.068170969 0.083172104
## [199] 0.047817706 0.061285542 0.074825574 0.041675973 0.061054468 0.061084310 0.062291497 0.076046293 0.076346615 0.061065044 0.051116030
## [210] 0.028350821 0.050016311 0.037194993 0.051100509 0.062254174 0.075974064 0.076303962 0.061029271 0.075340011 0.051140961 0.053944975
## [221] 0.067077837 0.084531120 0.048175237 0.069907140 0.070021507 0.069132616 0.085151008 0.086107485 0.069980182 0.065764840 0.084212029
## [232] 0.084349026 0.068243826 0.086059481 0.069107301 0.085081941 0.086061711 0.069940141 0.085269245 0.085924865 0.062368746 0.076126599
## [243] 0.076467185 0.061146059 0.078600961 0.051161326 0.078483355 0.069253475 0.085284853 0.086286639 0.070022679 0.087631727 0.086210453
## [254] 0.087512964 0.078690391 0.087776859
## 
## $IJ
##                 I             J       Total
## T     0.012705213 -0.0078791895 0.004826023
## SOIL  0.001540157  0.0006100237 0.002150181
## GDD   0.008486843 -0.0033119075 0.005174935
## AWSSM 0.001487961  0.0007354890 0.002223450
## AET   0.016781218  0.0013982638 0.018179481
## D     0.008168783  0.0073550117 0.015523795
## RD    0.014136702  0.0154758215 0.029612524
## VPD   0.024469982 -0.0037871281 0.020682854
## 
## $I.perc
##       ind.exp.var
## T       14.474444
## SOIL     1.754628
## GDD      9.668656
## AWSSM    1.695163
## AET     19.118043
## D        9.306306
## RD      16.105272
## VPD     27.877487
## 
## $params
## $params$full.model
## [1] "y ~ T + SOIL + GDD + AWSSM + AET + D + RD + VPD"
## 
## $params$family
## [1] "gaussian"
## 
## $params$link
## [1] "default"
## 
## $params$gof
## [1] "Rsqu"

## $gfs
##   [1] 0.0000000000 0.0007395474 0.0370956457 0.0003842961 0.0370509497 0.0438900893 0.0150547915 0.0304497778 0.0049727784 0.0828866343
##  [11] 0.0483226904 0.0828937246 0.0445356005 0.0620752157 0.0308598941 0.0268806613 0.0792798311 0.0373752377 0.0628323727 0.0449581372
##  [21] 0.0592244324 0.0435000881 0.0792877735 0.0449674416 0.0568588104 0.0306238999 0.0234500599 0.0628814919 0.0449424358 0.0591071503
##  [31] 0.0434929283 0.0554361556 0.0568299010 0.0513530749 0.0358904832 0.0227243020 0.0304791426 0.1147498671 0.0828937296 0.0839669661
##  [41] 0.0857536695 0.0951819144 0.0899925362 0.1146721807 0.0876987802 0.1028832250 0.0722763286 0.0653846120 0.0840159673 0.0857302238
##  [51] 0.0950582471 0.0899275469 0.0724569806 0.0572550752 0.0581465886 0.0636525580 0.0621941478 0.0321471654 0.0792878527 0.0808446216
##  [61] 0.0813962874 0.0920602148 0.0853311533 0.0635592208 0.0449633497 0.0619362005 0.0435014876 0.0638018127 0.0741723915 0.0632421290
##  [71] 0.0824041803 0.0450563006 0.0872064164 0.0809000898 0.0813780975 0.0919370291 0.0852744446 0.0689367862 0.0575992047 0.0562479517
##  [81] 0.0591483177 0.0570573835 0.0309659897 0.0638777084 0.0741567388 0.0633033841 0.0821722493 0.0450459978 0.0870521306 0.0628510765
##  [91] 0.0560827872 0.0581258431 0.0598729383 0.1149869668 0.1165865495 0.1188900332 0.1248943738 0.1188923679 0.0841110535 0.0857758245
## [101] 0.0958556416 0.0901558346 0.0867591322 0.0952396301 0.0907386038 0.0954919640 0.0901521633 0.0952322814 0.1165668057 0.1188031082
## [111] 0.1247180680 0.1187724122 0.1124195637 0.0971292125 0.0954753984 0.1036260038 0.1044830005 0.0723589832 0.0867707763 0.0951316830
## [121] 0.0907135567 0.0953561615 0.0900831210 0.0951069212 0.0735973943 0.0724875823 0.0592036478 0.0666239404 0.0810379347 0.0814073694
## [131] 0.0927578698 0.0854544692 0.0828893271 0.0922480371 0.0864898303 0.0928404301 0.0856134384 0.0925010413 0.0646580733 0.0742521534
## [141] 0.0640962297 0.0838966227 0.0450578090 0.0884298746 0.0847344704 0.0638818496 0.0878684715 0.0879355201 0.0829121437 0.0921515139
## [151] 0.0864801748 0.0926972720 0.0855518604 0.0923722092 0.0706185993 0.0689470123 0.0583372518 0.0637018064 0.0846282502 0.0639586579
## [161] 0.0877659925 0.0877592534 0.0700033682 0.1166034985 0.1195106773 0.1266168465 0.1196112588 0.1206265490 0.1252890891 0.1203536657
## [171] 0.1249040664 0.1194838551 0.1256083852 0.0867768953 0.0958580091 0.0907649319 0.0961848912 0.0903136608 0.0959191636 0.0955165883
## [181] 0.0908717785 0.0952720642 0.0954982506 0.1205874032 0.1251523427 0.1202909691 0.1247310171 0.1193786682 0.1254238485 0.1128883197
## [191] 0.1130044954 0.0974409523 0.1105200516 0.0953930175 0.0908445351 0.0951606710 0.0953621119 0.0743418914 0.0829451900 0.0927816554
## [201] 0.0864919653 0.0935741708 0.0857351664 0.0932464287 0.0929306191 0.0867193880 0.0925973634 0.0928931170 0.0852374603 0.0647022705
## [211] 0.0886162037 0.0891430397 0.0884917598 0.0928099896 0.0867078601 0.0924893787 0.0927508162 0.0721080882 0.0883704284 0.1208253640
## [221] 0.1266522123 0.1206302923 0.1266241417 0.1202345392 0.1274333424 0.1253305480 0.1210548945 0.1258329091 0.1260533769 0.0962070535
## [231] 0.0908988588 0.0959308885 0.0961881537 0.0955271142 0.1252012978 0.1210030883 0.1256810020 0.1258900349 0.1162097632 0.0954041829
## [241] 0.0935741786 0.0867219953 0.0932464301 0.0936410155 0.0929671313 0.0892732291 0.0928452771 0.1266667733 0.1213389841 0.1274335139
## [251] 0.1279092157 0.1263257243 0.0962081715 0.1261960144 0.0936429182 0.1279103950
## 
## $IJ
##                 I             J        Total
## T     0.026431083 -0.0256915359 0.0007395474
## SOIL  0.018178269  0.0189173772 0.0370956457
## GDD   0.023822930 -0.0234386343 0.0003842961
## AWSSM 0.018101074  0.0189498753 0.0370509497
## AET   0.013913380  0.0299767095 0.0438900893
## D     0.008815862  0.0062389295 0.0150547915
## RD    0.014426159  0.0160236190 0.0304497778
## VPD   0.004221638  0.0007511408 0.0049727784
## 
## $I.perc
##       ind.exp.var
## T       20.663749
## SOIL    14.211721
## GDD     18.624702
## AWSSM   14.151371
## AET     10.877443
## D        6.892217
## RD      11.278332
## VPD      3.300465
## 
## $params
## $params$full.model
## [1] "y ~ T + SOIL + GDD + AWSSM + AET + D + RD + VPD"
## 
## $params$family
## [1] "gaussian"
## 
## $params$link
## [1] "default"
## 
## $params$gof
## [1] "Rsqu"

Plots of wb percentile vs. area burned show weak relationships

ECDF

Fire Danger Model

The fire danger model uses ECDF curves of water balance percentiles on day of ignition. Water balance variables are selected based on AUC values in the ROC analysis above. We can use the logistic regression models evaluated in the ROC analysis to estimate probability of ignition for given water balance percentiles for a given day, but due to the very large number of non-ignition days in the dataset, the probability of ignition for any given day will be extremely low. Instead, we will look only at days that had an ignition, and compare the water balance percentiles for those days. This allows us to, for example, take a percentile of water deficit on a given day and compare that to the proportion of historic fires that burned at or below this value.

Low fire danger is taken as the water balance percentile at or below which 10% of historical fires ignited.

High fire danger is taken as the water balance percentile at or below which 40% of historical fires

A regression line is fit to the ECDF curves to develop an equation that can be used to estimate proportion of historical fires for a given cover type burning at or below an input water balance variable percentile.

selected_vars <- c("RD", "VPD", "D")

regression_formulas <- tribble(
    ~cover, ~var, ~order, ~formula, ~r2
)

precision <- 7

regression1_formula <- function(a, b) {
    ## Regression formula in string format
    a <- round(a, precision)
    b <- round(b, precision)
    return (paste0("y = ", a, "e^(", b, "x)"))
}

regression2_formula <- function(b0, b1, b2) {
    ## Regression formula in string format
    b0 <- round(b0, precision)
    b1 <- round(b1, precision)
    b2 <- round(b2, precision)
    return (paste0("y = e^(", b0, " + ", b1, "x + ", b2, "x^2)"))
}

regression3_formula <- function(b0, b1, b2, b3) {
    ## Regression formula in string format
    b0 <- round(b0, precision)
    b1 <- round(b1, precision)
    b2 <- round(b2, precision)
    b3 <- round(b3, precision)
    return (paste0("y = e^(", b0, " + ", b1, "x + ", b2, "x^2 + ", b3, "x^3)"))
}

First Order

ecdf_regression1 <- function(ecdf_fn) {
    estimated <- data.frame(
        dec = seq(from = 0.01, to = 0.99, by = 0.01)
    )

    estimated$quants <- quantile(ecdf_fn, estimated$dec)
    estimated$proportions <- ecdf_fn(estimated$quants)

    lm <- lm(log(proportions) ~ quants, data = estimated)
    estimated$predicted <- exp(predict(lm, data = estimated$quants))
    
    ##plot(lm)

    a <- exp(lm$coefficients[[1]])
    b <- lm$coefficients[[2]]

    return (list(lm = lm,
                 estimated = estimated,
                 a = a,
                 b = b))
}

ecdf_regression1_plot <- function(ecdf_fn, cover, var, estimated, a, b) {
    high_thresh <- .4
    low_thresh <- .1
    
    regression_fn <- function(x) {
        a <- a
        b <- b
        return (a*exp(b*x))
    }

    plot(ecdf_fn, main = paste(region_name, cover, var))
    points(y = estimated$proportions, x = estimated$quants, col = "blue")
    points(y = estimated$predicted, x = estimated$quants, col = "red")
    curve(regression_fn, add = TRUE, col = "red")
    x_low <- round(log(low_thresh/a)/b,2)
    segments(0, low_thresh, x1 = x_low, y1 = low_thresh) ## Low Fire Danger horiz
    text(x = .3, y = low_thresh + 0.02, label = "Low Fire Danger")
    segments(x_low, 0, x1 = x_low, y1 = low_thresh) ## Low Fire Danger vert
    text(x = x_low, y = -.02, label = x_low)
    x_high <- round(log(high_thresh/a)/b,2)
    segments(0, high_thresh, x1 = x_high, y1 = high_thresh) ## High Fire Danger horiz
    text(x = .3, y = high_thresh + 0.02, label = "High Fire Danger")
    segments(x_high, 0, x1 = x_high, y1 = high_thresh) ## High Fire Danger vert
    text(x = x_high, y = -.02, label = x_high)
}

for (var in selected_vars) {
    for (cover in c("forest", "nonforest")) {
        print(paste("First order regression for", region_name, cover, var))
        
        ignitions <- get(paste0("ignitions_", cover))
        ecdf_fn <- ecdf(ignitions[[var]])
        print(summary(ecdf_fn))

        regression_list <- ecdf_regression1(ecdf_fn)
        estimated <- regression_list$estimated
        lm <- regression_list$lm
        r2 <- summary(lm)$adj.r.squared
        a <- regression_list$a
        b <- regression_list$b

        print(summary(lm))

        ecdf_regression1_plot(ecdf_fn, cover, var, estimated, a, b)

        regression_formulas <- regression_formulas %>%
            add_row(var = var, order = "1", cover = cover, formula = regression1_formula(a, b), r2 = r2)
    }
}

## [1] "First order regression for Southern Rockies forest RD"
## Empirical CDF:     151 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01608 0.77262 0.91113 0.83293 0.96883 1.00000 
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.21425 -0.08597 -0.04050  0.09519  0.54532 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.87476    0.06189  -78.77   <2e-16 ***
## quants       4.65165    0.07184   64.75   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1353 on 97 degrees of freedom
## Multiple R-squared:  0.9774, Adjusted R-squared:  0.9772 
## F-statistic:  4192 on 1 and 97 DF,  p-value: < 2.2e-16

## [1] "First order regression for Southern Rockies nonforest RD"
## Empirical CDF:     85 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1736  0.8407  0.9507  0.8684  0.9871  1.0000 
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67227 -0.12967 -0.00609  0.14367  0.37013 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5.8348     0.1013  -57.57   <2e-16 ***
## quants        5.5608     0.1139   48.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.178 on 97 degrees of freedom
## Multiple R-squared:  0.9609, Adjusted R-squared:  0.9605 
## F-statistic:  2385 on 1 and 97 DF,  p-value: < 2.2e-16

## [1] "First order regression for Southern Rockies forest VPD"
## Empirical CDF:     151 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.02609 0.74963 0.87862 0.80937 0.96599 1.00000 
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.58243 -0.04267 -0.02428  0.04716  0.30340 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.72124    0.04326 -109.13   <2e-16 ***
## quants       4.59924    0.05157   89.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09879 on 97 degrees of freedom
## Multiple R-squared:  0.988,  Adjusted R-squared:  0.9878 
## F-statistic:  7955 on 1 and 97 DF,  p-value: < 2.2e-16

## [1] "First order regression for Southern Rockies nonforest VPD"
## Empirical CDF:     86 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.08973 0.67852 0.91600 0.80352 0.98460 1.00000 
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.94618 -0.13086 -0.00746  0.13293  0.36844 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.04158    0.06443  -62.73   <2e-16 ***
## quants       3.80344    0.07659   49.66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1752 on 97 degrees of freedom
## Multiple R-squared:  0.9622, Adjusted R-squared:  0.9618 
## F-statistic:  2466 on 1 and 97 DF,  p-value: < 2.2e-16

## [1] "First order regression for Southern Rockies forest D"
## Empirical CDF:     152 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.7097  0.9036  0.8059  0.9722  1.0000 
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.81463 -0.04907  0.00335  0.04295  0.34463 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.43992    0.04980  -89.15   <2e-16 ***
## quants       4.27929    0.05946   71.97   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.122 on 97 degrees of freedom
## Multiple R-squared:  0.9816, Adjusted R-squared:  0.9814 
## F-statistic:  5180 on 1 and 97 DF,  p-value: < 2.2e-16

## [1] "First order regression for Southern Rockies nonforest D"
## Empirical CDF:     86 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.7359  0.9078  0.8111  0.9743  1.0000 
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62200 -0.08615 -0.03224  0.10049  0.32293 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.22742    0.05296  -79.83   <2e-16 ***
## quants       3.99040    0.06251   63.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1373 on 97 degrees of freedom
## Multiple R-squared:  0.9767, Adjusted R-squared:  0.9765 
## F-statistic:  4075 on 1 and 97 DF,  p-value: < 2.2e-16

Second Order

Fits to ECDF curves may be improved by including a quadtratic term in regression

CWD Forest

ecdf_fn <- ecdf(ignitions_forest$D)
summary(ecdf_fn)

## Empirical CDF:     152 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.7097  0.9036  0.8059  0.9722  1.0000

estimated <- data.frame(
    dec = seq(from = 0.01, to = 0.99, by = 0.01)
)

estimated$quants <- quantile(ecdf_fn, estimated$dec)
estimated$proportions <- ecdf_fn(estimated$quants)

lm <- lm(log(proportions) ~ poly(quants, 2, raw = TRUE), data = estimated)
r2 <- summary(lm)$adj.r.squared
estimated$predicted <- exp(predict(lm, data = estimated$quants))

##plot(lm)

summary(lm)

## 
## Call:
## lm(formula = log(proportions) ~ poly(quants, 2, raw = TRUE), 
##     data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.76683 -0.05825 -0.00351  0.03662  0.35306 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -4.5638     0.1173 -38.920   <2e-16 ***
## poly(quants, 2, raw = TRUE)1   4.7012     0.3664  12.830   <2e-16 ***
## poly(quants, 2, raw = TRUE)2  -0.3116     0.2670  -1.167    0.246    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1218 on 96 degrees of freedom
## Multiple R-squared:  0.9819, Adjusted R-squared:  0.9815 
## F-statistic:  2600 on 2 and 96 DF,  p-value: < 2.2e-16

b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]

plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2)), add = TRUE, col = "red")

##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")

regression_formulas <- regression_formulas %>%
    add_row(var = "D", order = "2", cover = "forest", formula = regression2_formula(b0, b1, b2), r2 = r2)

Estimated function predicting proportion of historic fires burning in forests at or below given percentile of CWD on day of ignition:

$\begin{align*}f(x) = e^{-4.5638481 + 4.7012121x + -0.311572x^2}\end{align*}$

VPD Forest

ecdf_fn <- ecdf(ignitions_forest$VPD)
summary(ecdf_fn)

## Empirical CDF:     151 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.02609 0.74963 0.87862 0.80937 0.96599 1.00000

estimated <- data.frame(
    dec = seq(from = 0.01, to = 0.99, by = 0.01)
)

estimated$quants <- quantile(ecdf_fn, estimated$dec)
estimated$proportions <- ecdf_fn(estimated$quants)

lm <- lm(log(proportions) ~ poly(quants, 2, raw = TRUE), data = estimated)
r2 <- summary(lm)$adj.r.squared
estimated$predicted <- exp(predict(lm, data = estimated$quants))

##plot(lm)

summary(lm)

## 
## Call:
## lm(formula = log(proportions) ~ poly(quants, 2, raw = TRUE), 
##     data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53279 -0.04521 -0.02288  0.04875  0.31366 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -4.8406     0.1104 -43.836   <2e-16 ***
## poly(quants, 2, raw = TRUE)1   4.9907     0.3374  14.794   <2e-16 ***
## poly(quants, 2, raw = TRUE)2  -0.2846     0.2424  -1.174    0.243    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0986 on 96 degrees of freedom
## Multiple R-squared:  0.9881, Adjusted R-squared:  0.9879 
## F-statistic:  3994 on 2 and 96 DF,  p-value: < 2.2e-16

b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]

plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2)), add = TRUE, col = "red")

##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")

regression_formulas <- regression_formulas %>%
    add_row(var = "VPD", order = "2", cover = "forest", formula = regression2_formula(b0, b1, b2), r2 = r2)

Estimated function predicting proportion of historic fires burning in forests at or below given percentile of VPD on day of ignition:

$\begin{align*}f(x) = e^{-4.8405722 + 4.9907096x + -0.2846107x^2}\end{align*}$

CWD Non-forest

ecdf_fn <- ecdf(ignitions_nonforest$D)
summary(ecdf_fn)

## Empirical CDF:     86 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.7359  0.9078  0.8111  0.9743  1.0000

estimated <- data.frame(
    dec = seq(from = 0.01, to = 0.99, by = 0.01)
)

estimated$quants <- quantile(ecdf_fn, estimated$dec)
estimated$proportions <- ecdf_fn(estimated$quants)

lm <- lm(log(proportions) ~ poly(quants, 2, raw = TRUE), data = estimated)
r2 <- summary(lm)$adj.r.squared
estimated$predicted <- exp(predict(lm, data = estimated$quants))

##plot(lm)

summary(lm)

## 
## Call:
## lm(formula = log(proportions) ~ poly(quants, 2, raw = TRUE), 
##     data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72568 -0.07821 -0.02359  0.08248  0.34426 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  -4.06814    0.09911 -41.048   <2e-16 ***
## poly(quants, 2, raw = TRUE)1  3.38218    0.32740  10.330   <2e-16 ***
## poly(quants, 2, raw = TRUE)2  0.47122    0.24911   1.892   0.0616 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1355 on 96 degrees of freedom
## Multiple R-squared:  0.9776, Adjusted R-squared:  0.9771 
## F-statistic:  2093 on 2 and 96 DF,  p-value: < 2.2e-16

b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]

plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2)), add = TRUE, col = "red")

##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")

regression_formulas <- regression_formulas %>%
    add_row(var = "D", order = "2", cover = "nonforest", formula = regression2_formula(b0, b1, b2), r2 = r2)

Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of CWD on day of ignition:

$\begin{align*}f(x) = e^{-4.0681412 + 3.3821802x + 0.4712154x^2}\end{align*}$

VPD Non-forest

ecdf_fn <- ecdf(ignitions_nonforest$VPD)
summary(ecdf_fn)

## Empirical CDF:     86 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.08973 0.67852 0.91600 0.80352 0.98460 1.00000

estimated <- data.frame(
    dec = seq(from = 0.01, to = 0.99, by = 0.01)
)

estimated$quants <- quantile(ecdf_fn, estimated$dec)
estimated$proportions <- ecdf_fn(estimated$quants)

lm <- lm(log(proportions) ~ poly(quants, 2, raw = TRUE), data = estimated)
r2 <- summary(lm)$adj.r.squared
estimated$predicted <- exp(predict(lm, data = estimated$quants))

##plot(lm)

summary(lm)

## 
## Call:
## lm(formula = log(proportions) ~ poly(quants, 2, raw = TRUE), 
##     data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73402 -0.13632  0.01126  0.11806  0.38144 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -4.4289     0.1406 -31.497  < 2e-16 ***
## poly(quants, 2, raw = TRUE)1   5.1993     0.4611  11.276  < 2e-16 ***
## poly(quants, 2, raw = TRUE)2  -1.0489     0.3421  -3.066  0.00281 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.168 on 96 degrees of freedom
## Multiple R-squared:  0.9655, Adjusted R-squared:  0.9648 
## F-statistic:  1344 on 2 and 96 DF,  p-value: < 2.2e-16

b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]

plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2)), add = TRUE, col = "red")

##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")

regression_formulas <- regression_formulas %>%
    add_row(var = "VPD", order = "2", cover = "nonforest", formula = regression2_formula(b0, b1, b2), r2 = r2)

Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of VPD on day of ignition:

$\begin{align*}f(x) = e^{-4.4288633 + 5.199289x + -1.0488805x^2}\end{align*}$

Third Order

CWD Forest

ecdf_fn <- ecdf(ignitions_forest$D)
summary(ecdf_fn)

## Empirical CDF:     152 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.7097  0.9036  0.8059  0.9722  1.0000

estimated <- data.frame(
    dec = seq(from = 0.01, to = 0.99, by = 0.01)
)

estimated$quants <- quantile(ecdf_fn, estimated$dec)
estimated$proportions <- ecdf_fn(estimated$quants)

lm <- lm(log(proportions) ~ poly(quants, 3, raw = TRUE), data = estimated)
r2 <- summary(lm)$adj.r.squared
estimated$predicted <- exp(predict(lm, data = estimated$quants))

##plot(lm)

summary(lm)

## 
## Call:
## lm(formula = log(proportions) ~ poly(quants, 3, raw = TRUE), 
##     data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48901 -0.03760 -0.01246  0.05552  0.29634 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -6.4377     0.2318 -27.774  < 2e-16 ***
## poly(quants, 3, raw = TRUE)1  15.6781     1.2863  12.189  < 2e-16 ***
## poly(quants, 3, raw = TRUE)2 -18.8178     2.1281  -8.843 4.91e-14 ***
## poly(quants, 3, raw = TRUE)3   9.4812     1.0855   8.735 8.34e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09117 on 95 degrees of freedom
## Multiple R-squared:  0.9899, Adjusted R-squared:  0.9896 
## F-statistic:  3119 on 3 and 95 DF,  p-value: < 2.2e-16

b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
b3 <- lm$coefficients[[4]]

plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2) + b3*(x^3)), add = TRUE, col = "red")

##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")

regression_formulas <- regression_formulas %>%
    add_row(var = "D", order = "3", cover = "forest", formula = regression3_formula(b0, b1, b2, b3), r2 = r2)

Estimated function predicting proportion of historic fires burning in forests at or below given percentile of CWD on day of ignition:

$\begin{align*}f(x) = e^{-6.4376532 + 15.67814x + -18.8178063x^2 + 9.4811755x^3}\end{align*}$

VPD Forest

ecdf_fn <- ecdf(ignitions_forest$VPD)
summary(ecdf_fn)

## Empirical CDF:     151 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.02609 0.74963 0.87862 0.80937 0.96599 1.00000

estimated <- data.frame(
    dec = seq(from = 0.01, to = 0.99, by = 0.01)
)

estimated$quants <- quantile(ecdf_fn, estimated$dec)
estimated$proportions <- ecdf_fn(estimated$quants)

lm <- lm(log(proportions) ~ poly(quants, 3, raw = TRUE), data = estimated)
r2 <- summary(lm)$adj.r.squared
estimated$predicted <- exp(predict(lm, data = estimated$quants))

##plot(lm)

summary(lm)

## 
## Call:
## lm(formula = log(proportions) ~ poly(quants, 3, raw = TRUE), 
##     data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12681 -0.03563  0.00167  0.03105  0.20208 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -6.9555     0.1545  -45.03   <2e-16 ***
## poly(quants, 3, raw = TRUE)1  16.7469     0.8109   20.65   <2e-16 ***
## poly(quants, 3, raw = TRUE)2 -19.5730     1.3019  -15.03   <2e-16 ***
## poly(quants, 3, raw = TRUE)3   9.7347     0.6536   14.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05428 on 95 degrees of freedom
## Multiple R-squared:  0.9964, Adjusted R-squared:  0.9963 
## F-statistic:  8861 on 3 and 95 DF,  p-value: < 2.2e-16

b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
b3 <- lm$coefficients[[4]]

plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2) + b3*(x^3)), add = TRUE, col = "red")

##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")

regression_formulas <- regression_formulas %>%
    add_row(var = "VPD", order = "3", cover = "forest", formula = regression3_formula(b0, b1, b2, b3), r2 = r2)

Estimated function predicting proportion of historic fires burning in forests at or below given percentile of VPD on day of ignition:

$\begin{align*}f(x) = e^{-6.955476 + 16.7469326x + -19.5729642x^2 + 9.7346975x^3}\end{align*}$

CWD Non-forest

ecdf_fn <- ecdf(ignitions_nonforest$D)
summary(ecdf_fn)

## Empirical CDF:     86 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.7359  0.9078  0.8111  0.9743  1.0000

estimated <- data.frame(
    dec = seq(from = 0.01, to = 0.99, by = 0.01)
)

estimated$quants <- quantile(ecdf_fn, estimated$dec)
estimated$proportions <- ecdf_fn(estimated$quants)

lm <- lm(log(proportions) ~ poly(quants, 3, raw = TRUE), data = estimated)
r2 <- summary(lm)$adj.r.squared
estimated$predicted <- exp(predict(lm, data = estimated$quants))

##plot(lm)

summary(lm)

## 
## Call:
## lm(formula = log(proportions) ~ poly(quants, 3, raw = TRUE), 
##     data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35124 -0.03386 -0.00965  0.04239  0.37073 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -5.05656    0.09585  -52.75   <2e-16 ***
## poly(quants, 3, raw = TRUE)1  10.94736    0.60845   17.99   <2e-16 ***
## poly(quants, 3, raw = TRUE)2 -14.10398    1.11966  -12.60   <2e-16 ***
## poly(quants, 3, raw = TRUE)3   8.09566    0.61636   13.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08118 on 95 degrees of freedom
## Multiple R-squared:  0.992,  Adjusted R-squared:  0.9918 
## F-statistic:  3946 on 3 and 95 DF,  p-value: < 2.2e-16

b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
b3 <- lm$coefficients[[4]]

plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2) + b3*(x^3)), add = TRUE, col = "red")

##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")

regression_formulas <- regression_formulas %>%
    add_row(var = "D", order = "3", cover = "nonforest", formula = regression3_formula(b0, b1, b2, b3), r2 = r2)

Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of CWD on day of ignition:

$\begin{align*}f(x) = e^{-5.0565568 + 10.9473619x + -14.1039814x^2 + 8.0956582x^3}\end{align*}$

VPD Non-forest

ecdf_fn <- ecdf(ignitions_nonforest$VPD)
summary(ecdf_fn)

## Empirical CDF:     86 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.08973 0.67852 0.91600 0.80352 0.98460 1.00000

estimated <- data.frame(
    dec = seq(from = 0.01, to = 0.99, by = 0.01)
)

estimated$quants <- quantile(ecdf_fn, estimated$dec)
estimated$proportions <- ecdf_fn(estimated$quants)

lm <- lm(log(proportions) ~ poly(quants, 3, raw = TRUE), data = estimated)
r2 <- summary(lm)$adj.r.squared
estimated$predicted <- exp(predict(lm, data = estimated$quants))

##plot(lm)

summary(lm)

## 
## Call:
## lm(formula = log(proportions) ~ poly(quants, 3, raw = TRUE), 
##     data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27312 -0.05405 -0.02523  0.06138  0.39568 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -6.3172     0.1402  -45.07   <2e-16 ***
## poly(quants, 3, raw = TRUE)1  17.5460     0.8152   21.52   <2e-16 ***
## poly(quants, 3, raw = TRUE)2 -23.4986     1.4265  -16.47   <2e-16 ***
## poly(quants, 3, raw = TRUE)3  12.1065     0.7631   15.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08843 on 95 degrees of freedom
## Multiple R-squared:  0.9906, Adjusted R-squared:  0.9903 
## F-statistic:  3321 on 3 and 95 DF,  p-value: < 2.2e-16

b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
b3 <- lm$coefficients[[4]]

plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2) + b3*(x^3)), add = TRUE, col = "red")

##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")

regression_formulas <- regression_formulas %>%
    add_row(var = "VPD", order = "3", cover = "nonforest", formula = regression3_formula(b0, b1, b2, b3), r2 = r2)

Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of VPD on day of ignition:

$\begin{align*}f(x) = e^{-6.3171867 + 17.5459699x + -23.4985919x^2 + 12.1064895x^3}\end{align*}$

RD Non-forest

ecdf_fn <- ecdf(ignitions_nonforest$RD)
summary(ecdf_fn)

## Empirical CDF:     85 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1736  0.8407  0.9507  0.8684  0.9871  1.0000

estimated <- data.frame(
    dec = seq(from = 0.01, to = 0.99, by = 0.01)
)

estimated$quants <- quantile(ecdf_fn, estimated$dec)
estimated$proportions <- ecdf_fn(estimated$quants)

lm <- lm(log(proportions) ~ poly(quants, 3, raw = TRUE), data = estimated)
r2 <- summary(lm)$adj.r.squared
estimated$predicted <- exp(predict(lm, data = estimated$quants))

##plot(lm)

summary(lm)

## 
## Call:
## lm(formula = log(proportions) ~ poly(quants, 3, raw = TRUE), 
##     data = estimated)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.240049 -0.032203 -0.007788  0.039526  0.164608 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  -18.9616     0.5042  -37.60   <2e-16 ***
## poly(quants, 3, raw = TRUE)1  64.9844     2.2037   29.49   <2e-16 ***
## poly(quants, 3, raw = TRUE)2 -85.2061     3.0886  -27.59   <2e-16 ***
## poly(quants, 3, raw = TRUE)3  39.0603     1.3977   27.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05886 on 95 degrees of freedom
## Multiple R-squared:  0.9958, Adjusted R-squared:  0.9957 
## F-statistic:  7537 on 3 and 95 DF,  p-value: < 2.2e-16

b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
b3 <- lm$coefficients[[4]]

plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2) + b3*(x^3)), add = TRUE, col = "red")

##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")

regression_formulas <- regression_formulas %>%
    add_row(var = "RD", order = "3", cover = "nonforest", formula = regression3_formula(b0, b1, b2, b3), r2 = r2)

Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of 14 day rolling mean of RD on day of ignition:

$\begin{align*}f(x) = e^{-18.9615581 + 64.9843914x + -85.2061x^2 + 39.0603233x^3}\end{align*}$

Estimated Regression Formulas

cover	var	order	formula	r2
forest	RD	1	y = 0.007637e^(4.6516454x)	0.9771518
nonforest	RD	1	y = 0.002924e^(5.5607905x)	0.9605097
forest	VPD	1	y = 0.0089041e^(4.5992425x)	0.9878286
nonforest	VPD	1	y = 0.0175697e^(3.803439x)	0.9617630
forest	D	1	y = 0.0117969e^(4.2792856x)	0.9814296
nonforest	D	1	y = 0.01459e^(3.9903979x)	0.9765084
forest	D	2	y = e^(-4.5638481 + 4.7012121x + -0.311572x^2)	0.9814985
forest	VPD	2	y = e^(-4.8405722 + 4.9907096x + -0.2846107x^2)	0.9878759
nonforest	D	2	y = e^(-4.0681412 + 3.3821802x + 0.4712154x^2)	0.9771166
nonforest	VPD	2	y = e^(-4.4288633 + 5.199289x + -1.0488805x^2)	0.9648114
forest	D	3	y = e^(-6.4376532 + 15.67814x + -18.8178063x^2 + 9.4811755x^3)	0.9896312
forest	VPD	3	y = e^(-6.955476 + 16.7469326x + -19.5729642x^2 + 9.7346975x^3)	0.9963264
nonforest	D	3	y = e^(-5.0565568 + 10.9473619x + -14.1039814x^2 + 8.0956582x^3)	0.9917883
nonforest	VPD	3	y = e^(-6.3171867 + 17.5459699x + -23.4985919x^2 + 12.1064895x^3)	0.9902560
nonforest	RD	3	y = e^(-18.9615581 + 64.9843914x + -85.2061x^2 + 39.0603233x^3)	0.9956838

R Version

sessionInfo()

## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 11 (bullseye)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
##  [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/Boise
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] hier.part_1.0-6  GGally_2.1.2     sf_1.0-12        sp_1.6-0         pROC_1.18.0      slider_0.3.0     zoo_1.8-12      
##  [8] multidplyr_0.1.3 lubridate_1.9.2  forcats_1.0.0    stringr_1.5.0    dplyr_1.1.2      purrr_1.0.1      readr_2.1.4     
## [15] tidyr_1.3.0      tibble_3.2.1     ggplot2_3.4.2    tidyverse_2.0.0 
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.0    farver_2.1.1        warp_0.2.0          fastmap_1.1.1       reshape_0.8.9       stringfish_0.15.7  
##  [7] digest_0.6.31       timechange_0.2.0    lifecycle_1.0.3     processx_3.8.1      magrittr_2.0.3      compiler_4.3.0     
## [13] rlang_1.1.0         sass_0.4.5          tools_4.3.0         utf8_1.2.3          yaml_2.3.7          knitr_1.42         
## [19] labeling_0.4.2      bit_4.0.5           classInt_0.4-9      plyr_1.8.8          RColorBrewer_1.1-3  KernSmooth_2.23-20 
## [25] withr_2.5.0         nnet_7.3-18         stats4_4.3.0        grid_4.3.0          fansi_1.0.4         e1071_1.7-13       
## [31] colorspace_2.1-0    MASS_7.3-58.3       scales_1.2.1        gtools_3.9.4        betareg_3.1-4       cli_3.6.1          
## [37] rmarkdown_2.21      crayon_1.5.2        generics_0.1.3      RcppParallel_5.1.7  tzdb_0.3.0          DBI_1.1.3          
## [43] cachem_1.0.7        proxy_0.4-27        RApiSerialize_0.1.2 modeltools_0.2-23   splines_4.3.0       parallel_4.3.0     
## [49] vctrs_0.6.2         Matrix_1.5-3        sandwich_3.0-2      jsonlite_1.8.4      callr_3.7.3         hms_1.1.3          
## [55] bit64_4.0.5         Formula_1.2-5       jquerylib_0.1.4     units_0.8-2         glue_1.6.2          ps_1.7.5           
## [61] stringi_1.7.12      gtable_0.3.3        lmtest_0.9-40       munsell_0.5.0       pillar_1.9.0        htmltools_0.5.5    
## [67] R6_2.5.1            vroom_1.6.1         evaluate_0.20       lattice_0.20-45     highr_0.10          qs_0.25.5          
## [73] bslib_0.4.2         class_7.3-21        Rcpp_1.0.10         flexmix_2.3-19      nlme_3.1-162        mgcv_1.8-42        
## [79] xfun_0.39           pkgconfig_2.0.3

Rockies Fire Analysis

Stephen Huysman

2023-05-18

Settings

Southern Rockies

Preprocessing

Master dataframe

Map

Determine Fire Season

Compute Percentiles

Data Visualizations

Different Rolling Windows

Corrplot

ROC Analysis

ECDF

Fire Danger Model

First Order

Second Order

CWD Forest

VPD Forest

CWD Non-forest

VPD Non-forest

Third Order

CWD Forest

VPD Forest

CWD Non-forest

VPD Non-forest

RD Non-forest

Estimated Regression Formulas

R Version