Technical Memorandum

September 5, 2014 STI-910507-6062

To: Dennis Doll and Mark Evangelista, U.S. Environmental Protection Agency

From: Dianne Miller, Jennifer DeWinter, and Stephen Reid

Re: Documentation of Data Portal and Case Study to Support Analysis of Fire Impacts 

on Ground-Level Ozone Concentrations

Introduction

The U.S. EPA’s Exceptional Events Rule (72 FR 13560) defines exceptional events as 

unusual or naturally occurring events that can affect air quality but are not reasonably 

controllable using techniques that tribal, state, or local air agencies may implement to attain and 

maintain the National Ambient Air Quality Standards (NAAQS). Ambient data can be excluded 

from use in regulatory decision-making if an air agency demonstrates to EPA’s satisfaction that 

measurements were influenced by exceptional events. Technical evidence of such influences 

must be submitted to EPA in the form of demonstration packages that address the following 

criteria:

 The event affects air quality.

 The event is not reasonably controllable or preventable.

 The event is caused by human activity that is unlikely to recur at a particular location, or 

is a natural event.

 There is a clear, causal relationship between the measured concentrations at the 

impacted monitoring site and the specified event.

 The measured values were in excess of normal historical fluctuations.

 No exceedance or violation would have occurred but for the event.

Examples of exceptional events include large wildfires, which emit pollutants such as 

nitrogen oxides (NOx), volatile organic compounds (VOC), and fine particulate matter (PM2.5), 

important contributors to enhanced ozone and PM2.5 concentrations in the atmosphere. 

Because ozone is a secondary pollutant formed by reactions involving precursor species 

emitted by fires and other sources, significant technical challenges are associated with meeting 

the “but for” criterion listed in the last bullet above. Rather than simply establish a line of 

transport from the fire to the impacted ambient monitor, “but for” demonstrations must evaluate 

ozone levels that would be typical for a particular time and place in the absence of fire 
September 5, 2014 Page 2

emissions, and must provide a quantitative assessment of ozone levels with and without the 

influence of fire. 

In 2013, under Work Assignment 4-17, we completed an assessment of analysis 

techniques available to address the “but for” criterion for exceptional events related to 

fire-influenced ozone exceedances and made recommendations for developing tools or 

methods that would assist states in preparing fire-related exceptional events demonstration 

packages. A key recommendation was the development of a data portal that would provide 

state agencies with access to EPA’s database of quality-controlled meteorological data for the 

continental United States, known as the Omnibus Meteorological Database, or MetDat. The 

MetDat database includes a wide variety of parameters that are useful for developing statistical 

regression equations for ozone forecasting, one of the principal methods for estimating the 

impact of fire emissions on ambient ozone concentrations.

Under the current work assignment (WA 5-07), we developed a web-based data 

acquisition tool for MetDat and integrated this tool into the AirNow-Tech software system. We

also completed a sample case study to demonstrate the process of using MetDat data to 

develop a statistical regression tool and apply it to an exceptional event analysis of wildfire 

smoke. In addition, we produced R scripts that can be modified by individual agencies to 

perform regression development and testing. These work steps are documented in this 

technical memorandum, which includes

1. An overview of the MetDat software.

2. Documentation of a sample case study for developing a statistical regression equation to 

predict 8-hr ozone concentrations in Los Angeles using data acquired from MetDat.

3. Documentation of the R scripts used to automate equation development and testing. 

MetDat Overview

We worked with EPA to develop MetDat, which includes a synthesis of raw and 

computed meteorological parameters from surface and upper air observations, model 

reanalyses, and trajectory calculations. EPA has used MetDat to develop statistical models 

(Camalier et al., 2007) for adjusting air quality trends to account for year-to-year meteorological 

variability (U.S. Environmental Protection Agency, 2008). The original MetDat database 

covered 1995 through 2006. STI recently updated the database to also cover 2007 through 

2011. This updated database, along with the corresponding data processing codes, is known 

as MetDat2. The MetDat2 database provides data for over 700 sites across the United States

(Craig et al., 2012).

To use the MetDat data portal, which currently includes data for 2007 through 2011, log 

in to www.airnowtech.org. Users can register for a new AirNow-Tech account on this page. 

MetDat is linked under the Data page. Once on the MetDat page (Figure 1), choose the 

surface meteorological monitoring site for which to download data. Sites are grouped 

alphabetically by state. Next, select the meteorological parameters to download for the chosen 

site. Click Select All to get the entire list of parameters, or customize the parameter list by 

checking individual parameter boxes. Hover over any parameter to view a full description. 
September 5, 2014 Page 3

Parameters marked with an asterisk are a suggested subset to develop ozone forecasting 

regression equations. The resources page 

(www.airnowtech.org/Resources/MetDatResources.cfm) contains parameter definition tables, a 

data inventory table, and more information about the suggested parameter subset. The 

parameter definition table is also available in Appendix A.

Figure 1. MetDat portal in AirNow-Tech (circled in red).

Many of these parameters are calculated from hourly surface observations. 

Occasionally, fewer than 75% of the potential hours of data are available for calculating the 

metrics. To include these incomplete data in the data set query, check the box next to “Include 

incomplete data.”

Next, select the desired range of years and months. The default years are the first and 

last years of data available for the site selected. After selecting the years, a subset of months of 

data to download may be selected. For example, if summer-season data for each year is 

preferred, change the first month to June and the last month to September. 

Click Run Query to submit the data request. File links will appear in the Results pane 

when the data set is ready to be downloaded. A link for an unpivoted data file will appear, which 

includes columns for date, parameter, value, and completeness status (Figure 2). If the 
September 5, 2014 Page 4

“Include incomplete data” box was not checked, a link to a pivoted file will also appear, which 

displays each parameter in a separate column (Figure 3).

Figure 2. Unpivoted MetDat data file (excerpt).

Figure 3. Pivoted MetDat data file (excerpt).
September 5, 2014 Page 5

Los Angeles Case Study 

We completed a case study to demonstrate how MetDat can be used to develop a 

statistical regression equation to evaluate the impact of smoke on ozone concentrations. We 

selected potential cities for the case study based on the following criteria:

 Western states 

 Urbanized

 Cover different types of terrain (e.g., simple flat; complex terrain, coastal, valley)

 Cover wide range of ozone exceedances (i.e., ozone levels slightly higher than NAAQS 

[e.g., 80 ppb] to ozone levels much higher than NAAQS [e.g., 150 ppb]). 

 Occurrence of high wildland fire emissions in the past five years

Based on these criteria, we selected six candidate cities: Albuquerque, Salt Lake City, 

Las Vegas, Seattle, Sacramento, and Los Angeles. The next step was to screen the potential 

cities with the following steps:

 Retrieve Air Quality System (AQS) data.

– Retrieve site information (AMP500) and 8-hr ozone1(AMP350NW) from AQS. 

– Identify ozone exceedance days for each site. 

– Prepare site location list for further analysis. 

 Review BlueSky smoke impact plots to identify impacts from fires. 

 Review satellite images.

– Hazard Mapping System (HMS) smoke product (only graphics from the last 6 months 

are available).

– The Smog Blog,

2 which discusses impacts on air quality from noticeable fires.

– Moderate Resolution Imaging Spectroradiometer (MODIS) visible satellite imagery.

 Run HMS Plumes/AQS Intersection tool to identify whether there is any smoke impact 

for each ozone exceedance event. The HMS Plumes/AQS Intersection tool is an

internally developed GIS-python program based on: (1) the digitalized, date-dependent 

HMS plumes in a 1-km resolution grid converted from HMS polygons from the National 

Oceanic and Atmospheric Administration (NOAA) and (2) the AQS site-specific data 

prepared above.

3

 The workflow of the program is:

– Import AQS site-specific data into ArcGIS.

– Query the AQS data, and determine whether plume rasters were present at the 

same date and location as the air quality exceedances. 

– Write results to two separate output files. If no digitized rasters exist for a query 

date, the date is recorded in a “no plumes” output file. For rasters that exist, ozone 

exceedance information, site latitude and longitude, and the value of the plume 

raster (Y for intersected with plume, N for not intersected with plume) is written to a 

 

1 Duration code: W – 8-hr run average begin hour. 

2

http://alg.umbc.edu/usaq/archives/2012_07.html.

3

Analysts can use AirNow-Tech to manually review HMS smoke plumes relative to their sites on days with high PM2.5

concentrations.
September 5, 2014 Page 6

“results” output file. The plume raster values are then used to identify smoke days or 

non-smoke days. 

After completing these steps, we determined that Los Angeles was the best candidate 

city for our case study, based largely on the number of sufficient exceedance days with and 

without smoke influence. It is important to have several non-smoke exceedance days in a data 

set when developing an equation so the equation is tuned to capture high ozone concentrations.

Background on Regression Equations

Regression is a statistical method for describing relationships among variables. For 

estimating air quality concentrations, regression equations are developed to describe the 

relationship between pollutant concentrations (referred to as the predictand, what is being 

predicted) and primarily meteorological variables (referred to as the predictors). Regression 

equations have been successfully used to predict daily or sub-daily pollutant concentrations in 

many areas of the country (Cassmassi, 1987; Hubbard and Cobourn, 1997; Ryan, 1994; Dye et 

al., 1996; U.S. Environmental Protection Agency, 2003). 

Because regression equations are developed with several years of data, they represent 

the relationship between air quality and meteorology under typical emission patterns. 

Therefore, the difference between the predictions and observations can provide a reasonable 

estimate of the air pollution caused by unusual emissions (e.g., emissions from wildland fires). 

How Regression Equations Work

If two variables are correlated, a mathematical equation can generate a line or a curve 

that depicts the relationship between those variables. With this mathematical equation, one 

variable (e.g., ozone) can be predicted from other variables (e.g., meteorology). Multi-linear 

regression is most commonly used to predict ozone (Equation 1) from multiple variables. 

However, curvilinear regression (Equation 2) is also useful in predicting ozone because it 

captures the non-linear relationships of ozone and predictor variables. 

Ozone = c1 V1 + c2 V2 ....... cn Vn + constant (1)

Ozone = c1 V1+ c2 V2

2

+ c3 V3

3 ......, cn Vnn

+ constant (2)

where:

Ozone = predictand

c = coefficients (weighting factors)

V = meteorological predictor variables
September 5, 2014 Page 7

Regression Equation Development for Los Angeles

We developed a sample regression equation to evaluate the impact of smoke on ozone 

concentrations for the Los Angeles metropolitan area using data from MetDat and EPA’s AQS. 

The first task was to build a data set as follows: 

 Obtain and process five years (2007-2011) of April through October ozone data from 

AQS for the Los Angeles metropolitan area. 

 Obtain five years (2007-2011) of surface and upper-air meteorological data from the 

MetDat data portal, with the output file in pivoted format. See Appendix A for parameter 

definitions. 

 Perform general quality control of the air quality and meteorological data.

 Calculate the daily maximum 8-hr ozone concentration across all sites in the Los 

Angeles region. Depending on an individual agency’s needs, a single site could be used 

instead of a regional maximum. 

 Identify and flag all days likely impacted by smoke (as described at the beginning of this 

section). 

 Merge the air quality data and daily smoke impact flags with the MetDat output file, 

linking records by date. 

Next, we imported the database into SYSTAT statistical software and performed some 

basic data exploration on the non-smoke days, including analysis of scatter plots and 

histograms, to assess how each potential variable relates to ozone. We also ran a factor 

analysis to identify statistically unique variables. 

Factor analysis groups related variables that are statistically unique so analysts can 

select representative variables from each group. This approach helps keep the variable list 

small and manageable. The factor analysis output for non-smoke days in Los Angeles is shown 

in Tables 1 and 2. Table 1 shows the eigenvalues for each factor, or column, shown in Table 2. 

Typically, only factors with an eigenvalue greater than 1 are used. The factor loadings in 

Table 2 show the importance of each variable in each factor. Higher loadings indicate a 

stronger relationship among the variables in that factor. Therefore, no more than two variables 

with a high loading in the factor should be selected as variables for regression. For the MetDat 

meteorological data set, 0.5 is a suggested cutoff loading value for determining which variables 

are statistically similar. 

Once the data were sorted into factors, we selected representative variables from each 

group. Considerations in variable selection included data completeness, ease of obtaining the 

parameter on a forecast basis, and the conceptual understanding of the parameter’s importance 

to ozone concentrations.
September 5, 2014 Page 8

Table 1. Eigenvalues for factor analysis.

Factor 1 2 3 4 5 6 7 8 9 10 11 12 13 14

Eigenvalue 16.286 7.316 4.928 3.508 2.649 2.347 1.995 1.523 1.380 1.305 1.150 1.091 1.060 1.014

Table 2. Factor loadings. Identified factor groupings are shown in bold italics. 

Parameter definitions can be found in Appendix A. 

Parameter 1 2 3 4 5 6 7 8 9 10 11 12 13 14

T850 0.94 0.17 0.02 -0.06 -0.03 -0.07 -0.04 -0.08 -0.01 0.02 0.12 -0.01 0.01 0.10

H700 0.92 0.06 0.00 -0.13 -0.22 0.14 0.07 0.03 0.06 -0.05 0.11 -0.08 0.10 0.07

H500 0.91 0.15 0.01 -0.12 -0.22 -0.04 0.08 -0.03 0.08 0.00 0.02 -0.07 0.18 0.09

T925 0.91 -0.11 0.08 -0.03 0.05 0.16 0.05 -0.09 -0.18 0.01 0.08 -0.10 0.02 -0.02

TAAVG 0.87 -0.30 0.14 0.24 0.09 0.00 0.07 0.10 0.11 0.06 0.08 0.06 -0.01 0.00

MRMAX 0.87 0.14 -0.05 0.23 0.16 0.08 0.14 0.15 0.08 0.18 0.00 0.02 -0.04 -0.09

DPMAX 0.87 0.17 -0.08 0.23 0.13 0.08 0.15 0.16 0.09 0.19 -0.01 0.02 -0.02 -0.10

TAMIN 0.86 0.00 0.10 0.32 0.15 -0.01 0.05 0.16 0.14 0.06 0.10 0.06 -0.03 0.01

T700 0.84 0.24 0.01 -0.08 -0.11 -0.28 0.01 -0.08 0.06 0.00 0.00 -0.01 0.21 0.10

DT925 0.82 0.04 0.03 -0.17 0.01 0.17 0.03 -0.16 -0.27 -0.04 0.03 -0.14 0.04 0.00

TAMAX 0.82 -0.48 0.14 0.15 0.08 0.03 0.05 0.07 0.10 0.06 0.07 0.02 -0.03 0.03

DT850 0.78 0.40 -0.07 -0.22 -0.10 -0.11 -0.08 -0.14 -0.05 -0.03 0.06 -0.02 0.04 0.15

TAVGAM 0.76 -0.56 0.17 0.15 0.04 0.00 0.04 0.05 0.07 0.04 0.08 0.03 -0.02 0.00

WS300 -0.76 -0.20 0.06 0.12 0.08 0.13 0.07 0.09 -0.16 0.33 0.11 0.09 0.28 0.08

TAVGPM 0.76 -0.53 0.23 0.16 0.03 -0.02 0.05 0.03 0.10 0.02 0.06 0.06 0.02 0.03

DPAVG 0.75 0.44 -0.12 0.23 0.08 0.03 0.13 0.22 0.07 0.21 -0.03 0.04 -0.05 -0.09

WS250 -0.73 -0.20 0.05 0.12 0.06 0.11 0.12 0.12 -0.17 0.38 0.17 0.11 0.32 0.12

TMAX 0.72 -0.63 0.17 0.10 0.06 0.03 0.03 0.02 0.09 0.02 0.08 0.01 -0.02 0.05

WS200 -0.71 -0.19 0.03 0.06 0.05 0.08 0.15 0.10 -0.16 0.35 0.19 0.11 0.33 0.15

T500 0.67 0.19 0.02 -0.07 -0.16 -0.25 0.11 -0.10 0.10 0.18 -0.14 -0.07 0.25 0.11

H850 0.64 -0.17 -0.05 -0.20 -0.36 0.47 0.16 0.12 0.08 -0.09 0.15 -0.15 0.07 0.01

MHMAX -0.61 -0.07 0.08 0.20 -0.05 -0.30 -0.15 0.13 0.20 0.10 -0.14 0.08 0.07 -0.14

SLPAVG -0.55 -0.07 -0.19 -0.25 -0.47 0.42 0.08 0.24 0.24 -0.11 0.05 -0.11 0.07 -0.04

STPAVG -0.54 -0.08 -0.19 -0.25 -0.47 0.42 0.08 0.24 0.24 -0.11 0.06 -0.11 0.07 -0.04

RHAVGMID 0.05 0.87 -0.34 0.05 0.07 0.06 0.06 0.13 0.03 0.16 -0.05 -0.01 -0.05 -0.08

RHAVG 0.11 0.87 -0.33 -0.01 0.05 0.08 0.08 0.14 -0.01 0.18 -0.09 -0.03 -0.04 -0.08

TDELTA 0.20 -0.85 0.11 -0.20 -0.05 0.07 0.01 -0.11 -0.03 -0.01 0.00 -0.05 0.01 0.06

RHAVGNIGHT 0.16 0.72 -0.31 -0.01 -0.02 0.09 0.10 0.17 -0.01 0.22 -0.13 -0.05 -0.10 -0.11

WDAVG -0.08 0.61 0.49 0.04 0.05 0.15 0.05 0.04 0.12 -0.15 0.16 0.06 -0.05 0.04

WDAVGAM 0.10 0.57 0.53 0.10 -0.08 0.12 0.06 -0.08 -0.18 -0.17 0.22 0.15 0.09 -0.24

DT700 0.43 0.54 -0.11 -0.29 -0.21 -0.36 -0.04 -0.15 0.01 -0.07 -0.12 -0.03 0.26 0.16

HAZEHRS 0.15 0.51 -0.35 -0.23 0.19 0.06 -0.20 -0.20 -0.03 0.08 0.10 0.00 0.00 0.14

UAVG -0.12 0.50 0.74 -0.05 -0.13 0.13 0.02 0.04 0.05 0.03 0.03 -0.03 -0.09 0.05

CCAVTAM -0.32 0.50 -0.27 0.02 0.27 0.06 -0.07 0.09 0.24 -0.14 0.24 0.08 0.00 0.37
September 5, 2014 Page 9

Parameter 1 2 3 4 5 6 7 8 9 10 11 12 13 14

UAVGPM -0.25 0.15 0.77 -0.16 -0.10 0.06 0.02 -0.03 0.17 0.27 -0.17 -0.15 -0.15 0.20

WSAVGPM -0.35 0.12 0.72 0.00 -0.15 0.04 0.04 -0.10 0.15 0.24 -0.14 -0.20 -0.15 0.25

UAVGAM 0.13 0.43 0.65 0.02 -0.18 0.06 0.02 -0.02 -0.22 -0.08 0.13 0.11 0.05 -0.28

WSAVG -0.49 0.29 0.64 0.23 -0.05 0.08 0.09 -0.13 0.15 0.04 0.03 -0.13 -0.05 0.14

WSAVGAM -0.24 0.33 0.52 0.35 -0.07 0.07 0.12 -0.26 -0.10 -0.13 0.15 0.00 0.17 -0.26

WDAVGPM 0.13 0.09 0.45 -0.59 0.27 0.04 -0.05 0.19 0.12 0.12 -0.15 0.10 0.00 -0.05

VAVG -0.04 0.08 -0.28 0.59 -0.52 -0.23 0.05 0.04 -0.11 0.05 0.08 -0.17 -0.11 0.09

VAVGPM -0.15 -0.03 -0.20 0.57 -0.52 -0.13 -0.03 -0.06 -0.14 0.04 0.11 -0.15 -0.15 0.13

VAVGAM -0.19 0.09 -0.25 0.55 -0.34 -0.19 0.06 0.08 -0.07 0.00 0.03 -0.14 -0.02 -0.02

DH500 0.15 -0.06 -0.22 -0.28 -0.59 0.20 0.14 -0.05 0.09 0.11 -0.08 0.40 -0.03 0.00

WD250 -0.25 -0.03 0.01 -0.27 0.00 -0.36 0.66 0.09 0.01 -0.15 0.11 0.00 0.00 -0.02

WD300 -0.30 -0.03 0.01 -0.30 0.01 -0.44 0.64 0.02 0.00 -0.03 0.14 -0.02 -0.04 0.03

WD200 -0.29 0.03 -0.06 -0.32 0.14 -0.36 0.58 0.05 -0.01 -0.05 0.12 -0.01 -0.20 -0.02

PRECIP -0.19 0.00 -0.27 0.33 0.15 0.28 0.35 -0.56 0.31 0.01 -0.16 0.02 0.08 -0.03

RAINHRS -0.24 0.03 -0.29 0.36 0.21 0.31 0.35 -0.50 0.29 -0.01 -0.11 0.05 0.06 -0.03

D24TD850 0.09 0.06 0.10 0.17 0.02 0.18 0.21 0.20 -0.33 -0.37 -0.47 0.20 0.15 0.36

CCAVGPM -0.09 0.39 -0.28 0.10 0.28 0.29 -0.06 0.04 0.02 -0.18 0.41 0.17 -0.13 0.30

MRMHI -0.14 0.09 -0.03 0.10 0.06 -0.21 -0.19 0.17 0.44 -0.22 -0.02 -0.08 0.46 -0.18

LOD 0.15 0.33 0.26 0.22 -0.33 -0.31 -0.18 -0.07 0.15 -0.02 0.05 0.37 0.08 0.11

D24T850 0.22 -0.09 -0.21 -0.35 -0.42 0.09 -0.02 -0.27 -0.08 0.23 -0.03 0.45 -0.16 -0.09

TD850 0.34 0.10 0.15 0.36 0.06 0.29 0.30 0.23 -0.22 -0.19 -0.35 0.06 0.07 0.07

FOGHRS 0.12 0.19 -0.16 -0.30 0.14 0.20 -0.01 -0.18 -0.26 0.14 0.04 -0.39 0.18 -0.06

After data exploration, we extracted a subset of data for independent testing. Our 

approach to creating the subset was to sort the data in descending order by daily peak ozone 

concentration and then flag every fourth point. These flagged data are then extracted and 

reserved as the independent data set for later testing of the regression equation. This 

methodology ensures that there are independent test samples representative of all 

concentrations in the data set. The nonflagged data are the development data set. We then 

used the regression algorithm in SYSTAT to generate an equation with the remaining data for 

non-smoke days. As part of this step, several input metrics can be set to ensure statistically 

sound regression equations. These metrics were as follows:

 Tolerance – limits the similarity or collinearity of input variables. For example, 950-mb 

temperature and 925-mb temperature are closely related; therefore, only one factor will 

be allowed into the equation. We generally recommend tolerance values of 0.3, but 

values as low as 0.1 can be used if too few variables are selected for the equation. For 

our Los Angeles example, we used 0.1. 

 Estimation – controls the method used to enter and remove variables from the 

equation. We used Stepwise estimation, in which variables are entered or removed from 
September 5, 2014 Page 10

the model one at a time. The alternative approach is Complete, which will include every 

variable in the initial configuration in the final equation. 

 Probability – specifies probabilities to enter and remove a variable from the model. A 

variable is entered into the model if its alpha value is less than the specified Enter value, 

and it is removed from the model if its alpha value is greater than the specified Remove 

value. We used probabilities of 0.05. An alpha of 0.05 is commonly used to prevent 

nonsignificant variables from being entered and to force out nonsignificant variables from 

the model. Basically, the probability metric assesses whether the variable improves the 

equation enough to warrant inclusion.

 MaxStep – dictates the maximum number of steps (variables) allowed. We used a 

maximum of 10 variables. For a record set this size, having more than 10 variables 

results in overfitting the data. This means that the model has been overly trained to the 

developmental data set to the point where it won’t have good forecasting ability with new 

data. Such an equation would be modeling noise or slight inaccuracies, which could 

become substantial errors when applied to different data records.

 Force – forces the first n variables listed to remain in the equation. We did not force any 

variables into the equations (i.e., n = 0).

Our first iteration used only 850-mb temperature and length of day as variables. We 

know that these variables alone do not capture all of the important atmospheric processes for 

ozone formation, so to increase the number of variables, we lowered the tolerance from 0.3 to 

0.2 and 0.1. In addition, because ozone is log-normally distributed, using the natural log of 

ozone concentrations as the dependent variable can improve performance, which proved to be 

the case for the Los Angeles data set. Once we found an equation with a sufficient number of 

variables, we excluded data points the software flagged as outliers and those with high leverage 

and reran the algorithm. This process was iterated until no data points were overly influencing 

the results. Then, we evaluated the output variables for physical sense (i.e., higher ozone 

should be predicted with lower wind speeds and higher temperatures) and evaluated the 

statistical strength of the equation. Output metrics included 

 Standard error – shows the average deviation of a sample from the expected mean; 

small error is best.

 T-test – shows whether each variable in the equation is statistically significant by 

comparing the mean of the predictions with and without the new variable. 

 P2-tail – shows the likelihood of the equation performing well with a random sample,

such as when used operationally. Values less than 0.05 are desired, indicating a 95% 

probability that the model will work with a random sample.

 Standard coefficient of a variable – shows how much the equation output would 

change by altering the variable by one standard deviation. In other words, it shows how 

heavily each variable contributes to the overall equation.

After several iterations, we achieved an equation that met the desired statistical strength 

metrics outlined above. Note, when developing a tool for forecasting or assessing smoke 
September 5, 2014 Page 11

impacts, several potential equations may be evident at this point in the process. For this sample 

case study, we are continuing with only one equation, shown below. 

 ( ) ( ) ( ) ( ) ( )

 ( )

Where: 

O3 = ozone concentration in ppb

850-mb temp = morning (1200 UTC) temperature at 850 mb (degK)

925-mb temp = morning (1200 UTC) temperature at 925 mb (degK)

SLPAvg = average sea level pressure (mb)

TAVGPM = afternoon (1-4 p.m. LST) average temperature (degK)

WSAVGPM = average afternoon (1-4 p.m. LST) wind speed (m/s)

Regression Equation Testing

Immediately after developing the equation, we evaluated its performance by comparing 

the predictions from the equations to the observations from both the data set used to develop 

the equations and the data subset reserved for testing. The results for the correlation coefficient 

(r2), bias, and absolute error values are shown in Table 3. Typically, we would expect to see an 

r

2

of at least 0.6 for this type of tool. However, since this example was prepared for 

demonstration purposes only, we did not continue iterating through the development steps to 

reach this standard. 

Table 3. Los Angeles 8-hr ozone equation test results.

Metric Result

Development data set r20.4473

Development data set bias 0.0043 ppm

Development data set absolute error 0.0112 ppm

Independent data set r20.4904

Independent data set bias 0.0013 ppm

Independent data set absolute error 0.0104 ppm

Application to Exceptional Event Analysis

After a statistically sound equation has been finalized, it can be used to estimate the 

wildland fire contribution to ozone concentrations. Using our demonstration equation for Los 

Angeles, the following steps were performed in Excel:

1. Calculate the difference between the daily ozone predictions by the regression equation 

and the daily observed ozone concentrations (error) for the non-smoke days.
September 5, 2014 Page 12

2. Plot the distribution of error to determine whether the errors are normally distributed 

about the mean (i.e., the distribution is symmetric). This indicates there is no directional 

bias and the errors are random. A histogram of the errors in our example indicates that 

the majority of days showed relatively small prediction errors, while very few days 

showed large prediction errors (Figure 4).

3. Calculate the value at which 68% and 95% of the positive daily errors, and 68% and 

95% of the negative daily errors, fall within the mean error for each prediction method. 

For a normal distribution of data, the 68th and 95th percentiles correspond to roughly one 

and two standard deviations from the mean, respectively. The use of the 95th percentile 

provides a conservative upper and lower bounds of the predictions. These error 

percentiles are summarized in Table 4. 

4. Add the upper bound (95th percentile of the positive error) to the predicted ozone value 

and subtract that result from the observed value to provide a conservative estimate of 

the wildfire smoke contribution to ozone each day. Table 5 shows the calculation steps

for the Los Angeles area. Because the sum of the prediction and error exceeded the 

observation on all days, no ozone can definitively be attributed to smoke. 

For this example, additional development work would be needed to improve the ozone 

prediction based on meteorological variables to perhaps better estimate ozone concentrations 

from wildfire. For example, next steps could include selecting some monitoring sites instead of 

all or further refining which variables to include.

Figure 4. Distribution of regression error on non-smoke days.
September 5, 2014 Page 13

Table 4. Summary of error values.

Error Metric Value

Error value (ppm) at which 95% of positive errors are within the mean +0.027

Error value (ppm) at which 68% of positive errors are within the mean +0.014

Error value (ppm) at which 68% of negative errors are within the mean -0.011

Error value (ppm) at which 95% of negative errors are within the mean -0.030

Table 5. Analysis of smoke days.

Date Observed O3

(ppm)

Predicted O3

(ppm)

95th

Percentile 

Error (ppm)

Maximum 

Expected 

Concentration 

(ppm)

Potential 

Smoke 

Contribution 

(ppm)

7/29/2007 0.086 0.084 0.027 0.111 -0.025

7/31/2007 0.08 0.089 0.027 0.116 -0.036

8/16/2007 0.088 0.094 0.027 0.121 -0.033

8/17/2007 0.082 0.091 0.027 0.118 -0.036

8/22/2007 0.085 0.090 0.027 0.117 -0.032

6/18/2008 0.1 0.096 0.027 0.123 -0.023

6/24/2008 0.095 0.081 0.027 0.108 -0.013

6/25/2008 0.094 0.074 0.027 0.101 -0.007

6/27/2008 0.086 0.075 0.027 0.102 -0.016

7/3/2008 0.108 0.099 0.027 0.126 -0.018

7/4/2008 0.112 0.097 0.027 0.124 -0.012

7/10/2008 0.091 0.094 0.027 0.121 -0.030

7/17/2008 0.097 0.079 0.027 0.106 -0.009

8/10/2009 0.084 0.077 0.027 0.104 -0.020

8/14/2009 0.086 0.081 0.027 0.108 -0.022

8/19/2009 0.085 0.088 0.027 0.115 -0.030

8/26/2009 0.095 0.085 0.027 0.112 -0.017

8/27/2009 0.089 0.090 0.027 0.117 -0.028

8/30/2009 0.118 0.105 0.027 0.132 -0.014

8/31/2009 0.108 0.095 0.027 0.122 -0.014

9/2/2009 0.093 0.079 0.027 0.106 -0.013

9/3/2009 0.102 0.090 0.027 0.117 -0.015

9/4/2009 0.087 0.092 0.027 0.119 -0.032

9/6/2009 0.085 0.077 0.027 0.104 -0.019

7/30/2010 0.085 0.080 0.027 0.107 -0.022

9/4/2010 0.099 0.103 0.027 0.130 -0.031

9/5/2011 0.085 0.093 0.027 0.120 -0.035
September 5, 2014 Page 14

R Scripts

To simplify the equation development and testing process for MetDat users, we created 

R scripts that perform the same steps described previously that were performed using Systat 

and Excel. By using R open source code, analysts will not be reliant on specific statistical 

packages. Appendix B summarizes the process used to develop and test the R scripts.

References

Camalier L., Cox W., and Dolwick P. (2007) The effects of meteorology on ozone in urban areas 

and their use in assessing ozone trends. Atmos. Environ., 41, ed., 7127-7137. doi: 

10.1016/j.atmosenv.2007.04.061. 

Cassmassi J.C. (1987): Development of an objective ozone forecast model for the South Coast 

Air Basin: Air Pollution Control Association 80th Annual Meeting, New York, NY, June 

21-26.

Craig K.J., Erdakos G.B., and Reid S.B. (2012) Report on meteorological input data for 

MetDat2. Technical memorandum prepared for the U.S. Environmental Protection 

Agency by Sonoma Technology, Inc. Petaluma, CA, 910319-5373-TM, May. 

Dye T.S., Ray S.E., Lindsey C.G., Arthur M., and Chinkin L.R. (1996) Summary of ozone 

forecasting and equation development for the air districts of Sacramento, Yolo-Solano, 

and Placer. Vol. I: ozone forecasting. Vol. II: equation development. Final report 

prepared for Sacramento Metropolitan Air Quality Management District, Sacramento, 

CA, by Sonoma Technology, Inc., Santa Rosa, CA, STI-996210-1701-FR, December. 

Hubbard M.C. and Cobourn W.G. (1997) Development of a regression model to forecast 

ground-level ozone concentration in Louisville, KY. Atmos. Environ., (submitted), ed. 

Ryan W.F. (1994) Forecasting severe ozone episodes in the Baltimore metropolitan area. 

Atmos. Environ., 29, ed., 2387-2398. 

U.S. Environmental Protection Agency (2003) Guidelines for developing an air quality (ozone 

and PM2.5) forecasting program. Prepared by Dye T.S., MacDonald C.P., Anderson C.B., 

Hafner H.R., Wheeler N.J.M., and Chan A.C. at Sonoma Technology, Inc., Petaluma, 

CA, STI-902461-2295-FR, EPA-456/R-03-002, June. Available on the Internet at 

http://www.epa.gov/airnow/aq_forecasting_guidance-1016.pdf.

U.S. Environmental Protection Agency (2008) National air quality status and trends through 

2007. Report prepared by the Office of Air Quality Planning and Standards, Air Quality 

Assessment Division, Research Triangle Park, NC, EPA-454/R-08-006, November. 

Available on the Internet at 

http://www.epa.gov/airtrends/2008/report/TrendsReportfull.pdf.
September 5, 2014 Page A-1

Appendix A

Parameter List

Name Description

CCAVGPM Afternoon (1-4 p.m. LST) average fractional cloud cover (tenths)

CCAVTAM Morning (7-10 a.m. LST) average fractional cloud cover (tenths)

D24T850 24-hr change in morning (1200 UTC) temperature at 850 mb (degK)

D24TD850 24-hr change in morning (1200 UTC) dewpoint temperature at 850 mb (degK)

DEVH500 Deviation in morning (1200 UTC) height of 500 mb surface (m) from the 10-year monthly mean (degK)

DEVH700 Deviation in morning (1200 UTC) height of 700 mb surface (m) from the 10-year monthly mean (degK)

DEVH850 Deviation in morning (1200 UTC) height of 850 mb surface (m) from the 10-year monthly mean (degK)

DEVT500 Deviation in morning (1200 UTC) temperature (K) at 500 mb from the 10-year monthly mean (degK)

DEVT700 Deviation in morning (1200 UTC) temperature (K) at 700 mb from the 10-year monthly mean (degK)

DEVT850 Deviation in morning (1200 UTC) temperature (K) at 850 mb from the 10-year monthly mean (degK)

DH500 24-hr 500 mb height change (m)

DPAVG Average dew point temperature (degK)

DPMAX Maximum dew point temperature (degK)

DT700 Morning (1200 UTC) 700 mb-surface temperature difference (degK)

DT850 Morning (1200 UTC) 850 mb-surface temperature difference (degK)

DT925 Morning (1200 UTC) 925 mb-surface temperature difference (degK)

FOGHRS Occurrence of fog (hours)

H1000 Morning (1200 UTC) height of 1000 mb surface (m)

H500 Morning (1200 UTC) height of 500 mb surface (m)

H700 Morning (1200 UTC) height of 700 mb surface (m)

H850 Morning (1200 UTC) height of 850 mb surface (m)

HAZEHRS Occurrence of haze (hours)

MHMAX Maximum mixing height (m) 4 a.m.-4 p.m. LST

MRMAX Maximum water vapor mixing ratio (g/kg)

MRMHI Maximum rate of mixing height increase (m/hr) 6 a.m.-12 p.m. LST

PRECIP Total precipitation (inches)

RAINHRS Occurrence of rain (hours)

RHAVG Average daily relative humidity (%)

RHAVGMID Midday (10 a.m.-4 p.m. LST) average relative humidity (%)

RHAVGNIGHT Nighttime (8 p.m.-7 a.m. LST) average relative humidity (%)

SLPAVG Average Sea level pressure (mb)

STPAVG Avereage Station pressure (mb)

T500 Morning (1200 UTC) temperature at 500 mb (degK)
September 5, 2014 Page A-2

Name Description

T700 Morning (1200 UTC) temperature at 700 mb (degK)

T850 Morning (1200 UTC) temperature at 850 mb (degK)

T925 Morning (1200 UTC) temperature at 925 mb (degK)

TAAVG Average apparent temperature (degK)

TAMAX Maximum apparent temperature (degK)

TAMIN Minimum apparent temperature (degK)

TAVGAM Morning (7-10 a.m. LST) average temperature (degK)

TAVGPM Afternoon (1-4 p.m. LST) average temperature (degK)

TD850 Morning (1200 UTC) dewpoint temperature at 850 mb (degK)

TDELTA Diurnal temperature change (degK)

TMAX Maximum surface temperature (degK)

TT1000 1000–500 mb thickness at 1200 UTC observed (m)

TT850 850–500 mb thickness at 1200 UTC observed (m)

UAVG Average daily wind U vector (m/s)

UAVGAM Average morning (7-10 a.m. LST) wind U vector (m/s)

UAVGPM Average afternoon (1-4 p.m. LST) wind U vector (m/s) 

VAVG Average daily wind V vector (m/s) 

VAVGAM Average morning (7-10 a.m. LST) wind V vector (m/s)

VAVGPM Average afternoon (1-4 p.m. LST) wind V vector (m/s) 

WD200 Upper-air wind direction at 200 mb 1200 UTC observed (deg)

WD250 Upper-air wind direction at 250 mb 1200 UTC observed (deg)

WD300 Upper-air wind direction at 300 mb 1200 UTC observed (deg)

WDAVG Average daily wind direction (deg)

WDAVGAM Average morning (7-10 a.m. LST) wind direction (deg)

WDAVGPM Average afternoon (1-4 p.m. LST) wind direction (deg)

WS200 Upper-air wind speed at 200 mb 1200 UTC observed (m/s)

WS250 Upper-air wind speed at 250 mb 1200 UTC observed (m/s)

WS300 Upper-air wind speed at 300 mb 1200 UTC observed (m/s)

WSAVG Average daily wind speed (m/s)

WSAVGAM Average morning (7-10 a.m. LST) wind speed (m/s)

WSAVGPM Average afternoon (1-4 p.m. LST) wind speed (m/s ) 
September 5, 2014 Page B-1

Appendix B

R Script to Perform Ozone Regression Analysis

To simplify the equation development and testing process for MetDat users, we created an R1

script that performs a regression analysis similar to the case study described in this technical 

memorandum. Our R script provides step-by-step instructions for building a linear regression 

model to estimate ozone concentrations using a data set from the MetDat database. The same 

ozone data for the Los Angeles Metropolitan case study described in this memorandum are 

discussed as an example throughout the script. However, users of the R script can easily import 

their own MetDat data file instead. Further, due to differences in the statistical techniques 

employed between Systat (the statistical package used for the case study) and R, the results for 

the case study analyses with each package are different. The differences are further discussed 

below. Regardless of the statistical technique or software used, there is not necessarily “one” 

ideal model, and variations in software may result in different equations that perform similarly 

well. 

Data Acquisition and Preparation

For the case study described in the script, ozone data for the Los Angeles Metropolitan Area 

were manually obtained from EPA AQS and pre-processed to compute the daily maximum 8-hr 

ozone concentration across all sites. Meteorological data were obtained from MetDat for the 

same time period and site location. The daily maximum ozone data (“MaxO3”) were then

manually combined with meteorological data from MetDat. The variable "TestCase" was created 

by first sorting the data in descending order by daily peak ozone concentration and then flagging 

every fourth point; this variable is used in the R script to create two data subsets, one for 

development and the second to test the regression equation. In addition, smoke from fires was 

also analyzed. A variable was manually added to the data set (“Smoke”), and each record in the 

data set was flagged as having been either impacted by smoke (Y) or not (N). The data 

preparation resulted in a csv file ready for use by the R script. The following fields are required 

in the file in order for the R script to run without modification:

“Date” – The date for each data record in month/day/year (m/d/yyyy) format.

“MaxO3” – The daily maximum ozone concentration.

“TestCase” – The Boolean flag for data records to be used in developing the regression 

equation (assigned a value of 0), and data records to be reserved for testing the regression 

equation (assigned a value of 1).

“Smoke” – The flag for data records with smoke impact (Y) and without smoke impact (N).

 1 R is a free and open-source statistical language and software package with many available libraries.
September 5, 2014 Page B-2

R Script Overview

The R script is designed to be run step-wise (rather than all at once) because the regression 

analysis requires review and decision-making throughout the process. There are ten basic steps 

included in the R script:

1. Set the working directory for your R project.

2. Import your data file from csv, and format the date variable.

3. Review the imported file and the data summary.

4. Use scatterplot matrices to explore relationships between variables.

5. Create two subsets of data, one for development of the regression equation and the 

other for independent testing of the regression equation.

6. Run forward step-wise regression.

7. Evaluate regression model performance using performance statistics and diagnostic 

plots.

8. Consider using the natural log of ozone concentration as the dependent variable.

9. Use the final regression equation to estimate ozone concentrations in the test data

set, and evaluate the prediction performance.

10. Export the results to a file (if desired).

Additional detail about each step and the results for the case study data are described in the 

next section.

Step-by-Step Results

1. Set the working directory for your R project.

The script uses a default location (C:/ProjectFiles/OzoneRegression/Output); this must be 

modified in the script to correspond to a user-specified output directory.

2. Import your data file from csv, and format the date variable.

Two R functions are invoked during this step. The first function,

(input.file <- read.csv(file.choose(), header = TRUE))

will launch a Windows Explorer dialogue to allow the user to select their data file. For this 

example, we are using the Los Angeles case study data obtained from MetDat and EPA AQS. 
September 5, 2014 Page B-3

The data include the column “MaxO3” which is the daily peak 8-hr ozone value. The second 

function specifies the format of the date/time field 

(input.file$Date1 <- as.POSIXct(strptime(input.file$Date1, format = 

"%m/%d/%Y", "GMT"))))

3. Review the imported file and the data summary.

The R summary() function can be used to review summary statistics about all of the variables 

included in the input data file. In particular, the summary statistics are useful to identify any 

variables that have significant missing data, as these may need to be excluded from the 

analysis. The regression will only be modeled for complete data records (e.g., records with a 

value for all variables); therefore, including variables with missing data might limit the sample 

set. In this example, three variables (e.g., “TT1000” shown in Figure 1) had significant missing 

data and will be excluded from the data set (see step 6).

Figure 1. Summary output for the example data set. The “NA’s” in the display provides

the count of missing data for each variable. In this example, the variable “TT1000” was 

missing 1315 records out of the total data set (1510 records).

4. Use scatterplot matrices to explore relationships between variables.

Due to the large number of variables in MetDat, it is useful to explore only a subset of variables 

at a time in a matrix. The R script creates a sample scatterplot matrix (Figure 2) using the initial 

recommended variables denoted with an asterisk at the top of the variable list in MetDat.

5. Create two subsets of data, one for development of the regression equation and the 

other for independent testing of the regression equation.

The example data set includes 1,510 observations for 63 variables. The variable "TestCase" 

was developed by first sorting the data in descending order by daily peak ozone concentration 

and then flagging every fourth point. The flagged data ("1") are reserved as the independent 

data set for later testing of the regression equation; the independent data set can include 

concentrations with smoke impacts. The nonflagged data ("0"), without smoke influence, are 

used as the development data set.
September 5, 2014 Page B-4

Figure 2. Scatterplot matrix showing the relationship between ozone (“MaxO3”) and the 

recommended initial meteorological parameters.

6. Run forward step-wise regression.

This script relies upon several built-in functions for regression modeling: linear model (lm), 

step(), and update(). These functions are available in the "stats" package, which is included in 

the basic R distribution. There are two basic steps to run step-wise regression using basic R. 

First, you must setup the range of the linear model. The null model includes only the intercept, 

while the full model includes all variables in the data set. Second, the step() function is invoked 

to perform stepwise forward regression between the null and full model formula. The linear 

model function in R uses the Akaike Information Criterion (AIC) to identify variables to include at 

each step in the forward regression. This is a different approach than Systat, which uses the pvalue and partial correlation between variables.
September 5, 2014 Page B-5

7. Evaluate regression model performance using performance statistics and diagnostic 

plots.

A summary of the regression is available using the R summary() function in Step 7(a). This 

provides the model formula, the coefficients for the variables included in the formula and their

standard errors, as well the residual standard error and multiple R-squared values. The 

summary information for the regression results using the case study data is shown in Figure 3.

The model performed well with an adjusted R-squared value of 0.72.

Additional diagnostic information can be obtained using the plot() function. This will produce a 

set of four plots: (1) the residuals versus fitted values, (2) a Q-Q plot of standardized residuals, 

(3) a scale-location plot (square roots of standardized residuals versus fitted values), and (4) a 

plot of residuals versus leverage. The diagnostic plots for the case study are provided in Figure 

4.

The model summary and diagnostic plots produced in Step 7(a) should be evaluated to consider 

the type and quantity of variables included in the model, the coefficients and their significance, 

and the R-squared and adjusted R-squared. In this example, there are some variables with 

insignificant slopes (p-value > 0.05), as well as some points with high leverage that might be 

excluded. Step 7(b) in the R script provides a useful function to review the regression results if 

outliers and data points with high leverage are excluded. In this example, data on 9/1/2009 

(case 906) and 6/3/2009 (case 745) were found to have large residuals and high leverage. Step 

7(c) in the R script provides the R function to update the linear model to exclude any 

insignificant predictor variables.

8. Consider using the natural log of ozone concentration as the dependent variable.

Because ozone concentrations are typically log-normally distributed, using the natural log of 

ozone concentrations as the dependent variable can improve model performance; predicting the 

natural log of ozone concentration can be done using the same linear model functions along 

with the R log() function. For the example data set, the linear model developed using the natural 

log of daily peak ozone concentrations performed better with an adjusted R-squared value of 

0.74. The regression summary and diagnostic plots were then further evaluated to examine 

outliers and insignificant variables. In total, 13 outlier cases were excluded, and then the model 

was re-generated; this improved the adjusted R-squared value to 0.7518. The model was then 

updated to exclude three variables that were not significant: average morning temperature, 

upper air wind speed at 200 mb, and the temperature at 925 mb. The exclusion of several 

additional variables was also evaluated but did not improve the model performance. The final 

model is summarized in Figure 5; the adjusted R-squared value was 0.7524.

9. Use the final regression equation to estimate ozone concentrations in the test data

set, and evaluate the prediction performance.

The predict() function in R is used to apply the model formula developed in Step 8 to the test 

data set. The function requires the same variable names in the test data set as in the model 

formula. The prediction performance for the example here is shown in Figure 6.
September 5, 2014 Page B-6
September 5, 2014 Page B-7

Figure 3. Summary of the linear regression results. The variables are defined in 

Appendix A.

Figure 4. Diagnostic plots showing outliers and influential cases.
September 5, 2014 Page B-8

Figure 5. Summary of the linear model using the natural log of daily peak ozone 

concentrations as the response variable.
September 5, 2014 Page B-9

Figure 6. Scatter plot of the natural logarithm of the raw and predicted daily peak 8-hr 

ozone concentrations using the independent test data set.

10. Export the results to a file (if desired).

In the final step, the results from the development of the model are merged with the original 

input file. The prediction results are then also merged into the results data frame. Finally, a text 

file is exported to the working directory.

Other Considerations

The approach described in this appendix resulted in a final regression equation with many 

predictor variables. However, you may desire to simplify the equation in order to remove 

variables that are not likely to be a major factor in ozone formation, and to avoid over-fitting in 

the model. R has many additional features and tools that could be used to fine tune the equation 

and evaluate regression diagnostics. In particular, given the large number of variables in the 

MetDat database, the analyst could explore an alternative to the basic step-wise linear 

regression approach described in this appendix. The R “leaps” package offers the ability to 

perform all-subsets regression in which the number of predictor variables to be included in the 

model can be defined as well the “nbest” models evaluated (e.g., the analyst might look at the 

three best regression models when 3, 5, or 10 predictor variables are allowed). For additional 

information, please visit the R documentation page.2

 2

 http://cran.r-project.org/web/packages/leaps/leaps.pdf.

y = 0.7392x - 0.6717

R2 = 0.7353

-4

-3

-2

-1

0

-4 -3 -2 -1 0

ln(Predicted Daily Peak Ozone)

ln(Daily Peak 8-hr Ozone)

Prediction Performance Using Test Dataset