Moonshot

Air quality in the United States at the county level

Most of us go about our daily lives without giving much thought to the air we breathe. But there are times when we must be mindful of the air where we live, because air pollution may be high enough that our health could be affected.

The US Environmental Protection Agency, or EPA, measures the “cleanliness” of the air across the country. Thanks to these measurements, we can make informed plans regarding our outdoor activities.

The EPA uses the Air Quality Index ( $AQI$ ) to determine how benign the local air is at a given time. The $AQI$ is a number between 0 and 500. The lower the number, the cleaner the air. According to the EPA booklet “Air Quality Index – A Guide to Air Quality and Your Health”, the $AQI$ is calculated for four major air pollutants regulated by the Clean Air Act: ground-level ozone, airborne particles, carbon monoxide, and sulfur dioxide. The first two of these pose the greatest hazard to human health. Ozone can cause coughing, throat irritation or a burning sensation in the airways. Particle pollutants, which are microscopic solids or liquid droplets that can reach the lungs, may cause irritation of the nose, eyes and throat. Those at risk from ground-level ozone and particulate matter include people with asthma, children, and the elderly.

The $AQI$ associated with a specific pollutant is obtained from measurements of the pollutant concentration $C$ . The EPA uses a table of concentration ranges to which different $AQI$ values are pre-assigned (see the booklet link above). The desired $AQI$ is computed using a piece-wise linear function of $C$ :

$AQI = \frac{AQI_{\mathrm{high}} - AQI_\mathrm{low}}{C_{\mathrm{high}} - C_{\mathrm{low}}}\left(C - C_{\mathrm{low}}\right) + AQI_\mathrm{low}$ ,

where $C_{\mathrm{low}}$ and $C_{\mathrm{high}}$ are the lower and upper ends, respectively, of the tabulated pollutant concentration range such that $C_\mathrm{low} \leq C \leq C_\mathrm{high}$ , and $AQI_\mathrm{low}$ and $AQI_\mathrm{high}$ are the endpoints of the $AQI$ range assigned to the interval $\left[C_\mathrm{low}, C_\mathrm{high}\right]$ . Again, the lower the $AQI$ , the better the air quality.

The $AQI$ is formally computed by the above equation for each pollutant. Each measuring site might measure the concentration of all pollutants, of some of them, or of only one. When the concentration for more than one pollutant is measured, the final reported $AQI$ corresponds to the largest $AQI$ .

Which areas in the United States have the cleanest air, on average? More specifically, what are the time averages of the $AQI$ in different counties? To answer this question, one can obtain county-level $AQI$ data from this EPA webpage (scroll to the middle of the page to the table called “Daily AQI”; the csv files containing those data, as well as the R code to analyze the data, can be found in this GitHub repository). The reported daily $AQI$ values correspond to only 1,316 US counties (less than half of all US counties), because reports for areas with a population of less than 350,000 are not required by law. Note also that $AQI$ s for each sampled county are not reported strictly on a daily basis (for example, many measurements are reported every three days).

I used data from the year 2000 up to April 2020. The following map in Figure 1 shows the $AQI$ averaged over this period:

US_county_meanAQI — **Figure 1**: $AQI$ by county averaged between January 2000 and April 2020. Different colors correspond to different ranges of $AQI$ values, as indicated in the figure legend. Data for counties in gray are unavailable.

**Figure 1**: $AQI$ by county averaged between January 2000 and April 2020. Different colors correspond to different ranges of $AQI$ values, as indicated in the figure legend. Data for counties in gray are unavailable.

As described by the figure legend, the counties in green have the best air quality (remember that the lower the $AQI$ , the cleaner the air). In EPA terminology, the green counties have “Good” air quality. The counties in yellow have “Moderate” air quality. The air in the two regions in orange is “Unhealthy for Sensitive Groups”. A distribution of these mean $AQI$ s is shown in Figure 2:

hist_AQI — **Figure 2**: Distribution of the time-averaged $AQI$ values shown in Fig. 1.

**Figure 2**: Distribution of the time-averaged $AQI$ values shown in Fig. 1.

Figure 2 shows that most (84.1%) of the US counties for which $AQI$ data exist have had clean air (i.e., “Good” air quality) since 2000. The air quality was moderate, on average, for 15.7% of counties. Only two counties (0.15% of the total measured) had an average air quality that was unhealthy for certain people. These two counties are Kings, in California, and Ada, in Idaho.

During the 2000 – 2020 period examined here, what were the best and worst recorded air qualities in each of the studied counties? In the case of the best air quality, all 1,316 examined counties reported minimum $AQI$ s that were below 51. In other words, at different points in that period, each of those counties enjoyed “Good” air quality. For that reason, the minimum- $AQI$ map shows all of those counties in green, and it is therefore not included here.

For the case of the worst air quality reached in each county in the same period, Fig. 3 shows a US map with the maximum measured $AQI$ :

US_county_maxAQI — **Figure 3**: Maximum $AQI$ recorded in each available county in the period between January 2000 and April 2020.

**Figure 3**: Maximum $AQI$ recorded in each available county in the period between January 2000 and April 2020.

Remarkably, there are a few counties in which the air quality remained “Good” throughout. A list of those counties is shown in Fig. 4:

**Figure 4**: US counties with a maximum $AQI$ of at most 50, in the period between January 2000 and April 2020. The air quality in these counties was always “Good”.

**Figure 4**: US counties with a maximum $AQI$ of at most 50, in the period between January 2000 and April 2020. The air quality in these counties was always “Good”.

The distribution of maximum $AQI$ s is shown in Fig. 5. Note that there were several $AQI$ values (25 of them, corresponding to 2% of the total number of reported counties) well above 500, extending the tail of the distribution to the right into the thousands. Those extreme values were not included in the distribution of Fig. 5 for visual convenience. Nevertheless, it would be worthwhile to inquire further with EPA experts as to the meaning of such outliers.

**Figure 5**: Distribution of the maximum $AQI$ values shown in Fig. 3.

**Figure 5**: Distribution of the maximum $AQI$ values shown in Fig. 3.

From Fig. 5, 42% of the studied counties reported “Unhealthy” air conditions at some point during the 2000 – 2020 period. “Very unhealthy” conditions occurred in 29% of the counties, and “Hazardous” conditions in 4% of the counties.

Finally, calculating the standard deviation of each county’s $AQI$ ( $AQI_{\mathrm{RMS}}$ , where RMS stands for “root-mean-square”) gives an idea of the variation in air quality over the same period. This variation is shown in Fig. 6. In this figure, the numbers shown correspond to those counties for which the $AQI_{\mathrm{RMS}}$ is less than 100. Counties in blue have a low standard deviation, while counties in magenta have a high standard deviation. There were 2 counties in which $AQI_{\mathrm{RMS}}$ exceeded 100, both in California (neither of which is colored): Kings ( $AQI_{\mathrm{RMS}} = 533$ ) and Napa ( $AQI_{\mathrm{RMS}}=380$ ).

**Figure 6**: Standard deviation of $AQI$ recorded in each available county in the period between January 2000 and April 2020. Counties in which the $AQI$ varied little are shown in blue, while those in which the $AQI$ varied much are shown in magenta. Two California counties in which the root-mean-square (RMS) $AQI$ is greater than 100 are not shown.

**Figure 6**: Standard deviation of $AQI$ recorded in each available county in the period between January 2000 and April 2020. Counties in which the $AQI$ varied little are shown in blue, while those in which the $AQI$ varied much are shown in magenta. Two California counties in which the root-mean-square (RMS) $AQI$ is greater than 100 are not shown.

In conclusion, the EPA data on daily county $AQI$ values show that

on average, most of the counties for which $AQI$ values are available enjoyed clean air between January 2000 and April 2020;
most of those counties experienced at least one instance of unhealthy air conditions in the same period; and
all but two counties had a root-mean-square $AQI$ ( $AQI_{\mathrm{RMS}}$ ) of less than 100, with a median $AQI_{\mathrm{RMS}}$ of 20.

The top 50 astronomy keywords of all time

In a previous post I took a look at the collaboration numbers among astronomers, and verified that groups of more than five coauthors have become prevalent in the last five years or so. A natural question that arises from this is “What are the most popular topics in astronomy publications?”.

To answer this question, one can look for the “Keywords” section that appears in most “recent” (i.e., from 1970 onwards) astronomy papers, and count the number of actual keywords listed. Once again, the NASA ADS Astronomy Abstract Service proves invaluable for this task. So without further ado, here are the results:

These keywords correspond to refereed articles only, and they have all been converted to lower case. The search gave 128, 394 different keywords. The 3 least-frequent keywords were “zz ceti”, “zz psc” and “zz uma”, with only one instance each. The keyword “dark energy” was number 311 in the list, “mars” was number 881, and “extrasolar planets” was number 1,410.

It is also interesting to look at the frequency of keywords over time. I picked four keywords (“stars: atmospheres”, “brown dwarfs”, “gamma rays: observations”, and “astrochemistry”), and plotted their available yearly counts since 1990. I did this for each of three journals: The Astrophysical Journal (ApJ, shown in red below), Monthly Notices of the Royal Astronomical Society (MNRAS, blue), and Astronomy & Astrophysics (A & A, green). I also included the yearly counts for all the journals combined (shown in black):

Can you explain the two peaks in the black curve for “gamma rays: observations”?

Check out this cool interactive website built by Stefano Meschiari, in which you can choose a keyword and see its behavior in time.

A quick look at salaries in Austin

The city of Austin, Texas, is well known for being a center for technology and business. Many large companies, such as Dell, IBM, Google and National Instruments, to name a few, are headquartered or have regional offices here. According to Wikipedia, Austin was ranked number 1 by Forbes among all big cities for jobs in 2012, and it was also ranked number 1 by the Wall Street Journal Marketwatch for growing businesses, naming it the “San Francisco of Texas”.

I wanted to learn a bit more about the type of industries that have presence in Austin, as well as how much they pay (useful to know if you are in the job market), and I turned to the website Glassdoor.com, which contains first-hand information on companies, salaries associated with different job positions, and even employee reviews.

I did a search of all the companies listed on the website, and there were a little over 2,400. Not all of the listed companies included information about their positions or salaries (especially true for small companies). The first thing I scraped off the website was the distribution of companies over different industries. The results are summarized in this graph:

Austin’s character as a tech hub is obvious from this bar chart, where IT companies stand out. Keep in mind that this information is changing constantly.

How are salaries distributed across industries? Gathering the available salary figures for all positions in each company, and for every industry, gave this result:

Note that the salaries taken to produce this plot are averages for each job title. The black dot associated with each industry represents the median of the corresponding average salaries. The left end of each dashed bar corresponds to the minimum salary for that industry, while the right end is the maximum value. The left and right edges of the rectangular boxes denote the 1st and 3rd quartiles, respectively. There are several outliers indicated by the empty circles. Two industries in this plot (Non-Profit and Insurance) do not appear in the first plot, because not enough companies were listed in those industries on the website, which of course doesn’t mean that there aren’t many such companies in Austin.

Next, I counted the total number of different job titles. Unsurprisingly, a tech-related title stood out. I include the 15 most numerous job titles below:

The final count of different job titles turned out to be 13, 525. This number includes those cases in which there exist several variations of a job title. For example, similarly to Software Engineer, there are titles such as Software Development Engineer I, Software Development Senior Engineer, Software Applications Engineer, etc.

Salaries for some job titles are shown in this figure:

It’s also interesting to look at how salaries for one job title vary across different companies. Take a peek at this graph for Electrical Engineer:

The highest salary is approximately twice the amount of the smallest one. Now take a look at the case of Account Manager:

The abrupt changes shown in the bars are due to my selection of companies to include in the plot; since there are many more companies that can fit in one graph (over 100 of them), I picked three groups of companies from the ones available. The highest-paying company for Account Manager offers more than five times the salary of the one that pays the least.

In the two preceding plots, it’s worth bearing in mind that the same job title may involve different duties, depending on the company, so one has to be cautious about interpreting this raw information.

How do physicists and astronomers team up to write research papers?

The way in which physicists and astronomers team up to write technical papers has changed over the years, and not only is it interesting to look at this behavior for its own sake, but by analyzing the data it may be possible to better understand what role, if any, does the number of authors have on the scientific impact of a paper. Likewise, such an analysis can allow physics and astronomy journals to make decisions about their publishing policies.

I was curious about the trends in the number of authors per refereed astronomy paper, so I set out to write an R script that would read in data from the NASA Astrophysics Data System, an online database of both refereed and non-refereed academic papers in astronomy and physics. The script counts the monthly number of refereed astronomy and physics papers between January 1967 and September 2013, as well as the number of authors in each of those papers, and plots the results. It was a good exercise in manipulation of regular expressions. It is possible that the script can be made more efficient by using different R structures than the ones I used, but I wanted to look at the results first.

Here’s a plot of the monthly mean and maximum number of authors:

The upward trend in both sets of data is readily apparent. The maximum is much more variable, though. The record maximum number of authors for a single paper is 3062, for May 2011, which belongs to a study of proton-proton collisions performed at the Large Hadron Collider. The variability of the mean number has increased in the last three years or so.

Let’s decompose both series into their constituent parts, namely trend, seasonal and random:

In both cases, there is very little seasonality, and the random component is more stable for the series of mean number of authors.

As it turns out, there was a nice study by Schulman et al. (1997) in which it is reported that the annual mean number of authors per astronomy paper, published in different journals, increased between the years 1975 and 1995 by a factor of between 1.5 and 2. We see something similar here. In their article, those authors showed that the fraction of papers written by a single author decreased steadily in the same period by a factor of about 3. To check that out, I looked at the monthly data of fraction of papers written by N authors, with N=1, 2, 3, 4, 5, and also N>5. This is what I found:

Single-author papers have decreased by a factor of almost 4 since 1967 (black line). The period covered by Schulman et al. shows a similar trend here, even though I’m also including physics articles. Papers written by two authors (red line) grew in quantity until the mid-1980s, and then slowly decayed to their present fraction of 20%.

Perhaps the most noteworthy behavior is that of papers written by more than five authors (orange line). From making up only 1% of the total papers in 1967, they are now the most numerous, at 28%. There was an abrupt rise in the trend of these N>5 publications in 1990, which coincides with the launch of the Hubble Space Telescope, and it hasn’t changed significantly ever since. Performing a Holt-Winters fit, which is a method that uses exponentially weighted moving averages to update estimates of the seasonally adjusted mean, the slope and seasonals, we can make a forecast of the fraction of papers that will be written by more than five authors in the next 24 months. R has functions called HoltWinters and forecast.HoltWinters that do this. The result can be seen in this figure:

The blue line represents the forecast of relative number of papers until December 2014. The dark and light shaded areas are the 85% and 95% confidence regions, respectively.

Schulman et al. (1997) have pointed out possible explanations for the decrease in the number of single-author papers and increasing numbers of multi-author articles. One is the “growth of multiwavelength astrophysics… which requires astronomers to be proficient in multiple wavebands or to collaborate with experts in other wavebands”. Another possible cause is that the increasing competition for jobs and grants forces applicants to look for collaborations that will result in a longer publication list for their CV. Also, the development and availability of rapid communication through the Internet has boosted the interactions among colleagues.

It would be interesting to look at economic indicators side by side with the data shown above, to try to elucidate if there’s a connection between number of authors and economic trends.

Reference: Schulman, E. , French, J. C., Powell, A. L., Eichhorn, G., Kurtz, M. J. and Murray, S. S. 1997, Publications of the Astronomical Society of the Pacific, 109, 1278

Analyzing the ups and downs of planet factories

Planets, asteroids and comets are assembled inside mysterious gaseous structures called protoplanetary disks, which form around a newborn star (see Figure 1 below). The gas in these disks is turbulent due to a combination of magnetic, chemical and kinematic effects, and this turbulence is the driving agent behind a disk’s evolution. The stronger the turbulence, the faster the disk gas spirals into the central star.

Artist's impression of a protoplanetary disk. Credit: NASA — *Figure 1: Artist’s impression of a protoplanetary disk. Credit: NASA*

One way in which astrophysicists can understand these systems better is through numerical simulations of small disk patches . One of the simulation outputs is a very important quantity, known as “α” (alpha) in the planet-formation field, which is basically a dimensionless measure of the intensity of the turbulence. If scientists are able to decipher how disk turbulence varies during the lifetime of a protoplanetary disk (which can be on the order of 1 to 10 million years), we could potentially explain many fascinating features of our own solar system and of other planetary systems.

Up until now, astronomical observations of disks have provided only indirect, order-of-magnitude estimates of α. In spite of this, simulations that include the relevant physics allow us to glean the random variability of a disk’s α, and this provides theorists with a vital piece of information that can be used as input for their disk models.

However, thorough numerical descriptions of disks require large amounts of computational resources, particularly computing time, as well as knowledge of the complex codes that are used. What if it was possible to obtain α-values without having to carry out full 3-D simulations of disks every time someone needs them?

One way to do this is by performing a time series analysis of α-data from previously-run disk simulations, and obtaining a stochastic model that can be used in subsequent disk calculations. The statistical language R contains many packages and functions ideally suited for this task, and the analysis below uses those. Excellent introductions to the use of R in time series analysis can be found in “Introductory Time Series with R”, by P. S. P. Cowpertwait and A. V. Metcalfe (Springer), and in “Time Series Analysis WIth Applications in R”, by J. D. Cryer and K. Chan (also published by Springer). What follows is not intended to be a rigorous or comprehensive study, and I would welcome comments and suggestions to improve these results.

Figure 2 shows a realization of α from a typical protoplanetary disk calculation, using an astrophysical-fluid-dynamics code. The unit of time is one orbit, which is one revolution of the disk patch being simulated around the central star (this can easily be converted to years, or any other desired unit). Data points are spaced by 0.5 orbits. Here we are discarding the first 70 orbits, because it takes the turbulence some time to set in after the beginning of the simulation, and we want to make sure that the data are truly capturing the turbulent state.

alfa1 — *Figure 2: Time series of the turbulent intensity of a protoplanetary disk, as calculated from a numerical simulation.*

The correlation between α-values at time t and values at times t-k, with k=1, 10 and 200 representing the number of time steps or “lag” for each data point, are shown in Figure 3. Values that are separated by one lag (half an orbit) are highly correlated, as evidenced by the black symbols. On the other hand, α-values that are 200 lags (100 orbits) apart show almost no correlation, as can be seen from the blue symbols. The corresponding correlation coefficients are indicated in the figure.

alfacor — *Figure 3: Scatter plot showing the correlation of turbulence intensity values between time t and times t-1 (black symbols), t-10 (red symbols) and t-200 (blue symbols).*

Another way of examining correlations between pairs of data points is to plot the sample autocorrelation function (or sample ACF). A correlogram is a helpful visual tool that can be used for this purpose. Figure 4 presents the sample ACF of the data in Figure 2; for each k, the correlogram shows the correlation between α₁ and α_1+k, α₂ and α_2+k , α₃ and α_3+k, etc.

*Figure 4: Sample autocorrelation function (ACF) of data from Figure 2. See text for details.*

The blue horizontal lines in Figure 4 are drawn at ±2/√N, where N is the number of data points. That is, the lines represent two approximate standard errors above and below zero. The sample ACF has a sampling distribution from which we can infer whether the data come from a completely random series. If the sample ACF falls outside the lines (as it does for most lags in Figure 4), we have evidence against the null hypothesis that the autocorrelation function equals zero, at the 5% level. In other words, the correlations are statistically significant. If the sample ACF remained within the 95% confidence band, it would be reasonable to infer that no correlations exist.

The text annotation in Figure 4 reports the result of an augmented Dickey-Fuller (ADF) test on the α-data of Figure 2. I won’t go into the details of the test (an exposition can be found in the books mentioned above), but generally speaking, it can be used to quantify the certainty about whether a time series is non-stationary (strictly speaking, whether it contains a unit root, which is the null hypothesis). A non-stationary series has a non-constant mean and/or variance, and shows significant correlations in its correlogram. Additionally, in that case the Dickey-Fuller statistic will be a negative number with a relatively small absolute value. That seems to be the case for our data, and Figure 4 supports the notion that the disk turbulence intensity from the simulation is non-stationary.

Now, to proceed further we need to transform the original time series into a stationary one, that is, into a process which is, in a sense, in statistical equilibrium. We can achieve that by taking the natural logarithm of, and then differencing, the series. This is shown in Figure 5, where the left panel corresponds to the transformed series, and the right panel to its sample ACF.

alfalogdiff1_acf — *Figure 5: Differences of the natural logarithm of turbulent intensity (left panel), and their associated ACF (right panel). The results of an ADF test are shown in the right panel.*

The ACF shows a slow decay, but after some 350 lags no correlations remain. At this point, I would be tempted to say that even after this transformation (and some others that I’m not showing here) the series still remains non-stationary, given what we see in the correlogram. However, the augmented Dickey-Fuller test gives us reason to reject the null hypothesis of a unit root, or in other words, non-stationarity. Moreover, the transformed series shows a relatively stable mean compared to the original series in Figure 2, and the variance exhibits less change, except at a few times around 250 orbits. These larger changes would need to be incorporated into a more elaborate model, but for the time being I’m going with the result of the ADF test and take the transformed series to be stationary.

So, knowing that it was necessary to take the differences of the logarithms of α, we can use the “arima” function in R to obtain a model of that particular realization of α. “ARIMA” stands for Auto Regressive Integrated Moving Average, which describes the structure of a time series model consisting of recursive terms (the AR part) and terms of independent and identically distributed random variables (the MA part). The I in ARIMA refers to the fact that we had to difference the series to make it stationary. Exactly how many AR and MA terms are needed, and how many times we need to take differences, lies at the heart of the time series analysis problem. In this case, by trial and error, 5 AR terms and 5 MA terms are found to most accurately describe our stochastic process.

Using arima on the natural logarithm of the α-series, and once we specify that we differenced it once, R computes the coefficients of the ARMA components and their associated standard errors:

	AR1	AR2	AR3	AR4	AR5	MA1	MA2	MA3	MA4	MA5
Coefficient	0.295	-0.058	-0.058	0.942	-0.352	-0.555	0.017	-0.019	-0.982	0.570
S.E.	0.237	0.067	0.062	0.067	0.179	0.232	0.023	0.023	0.025	0.210

We can see that the coefficients for terms AR2, AR3, MA2 and MA3 are not significant. The resulting model is then

where the w_t‘s are white noise variables. This expression could, potentially, be used in relatively simple numerical models of disks that share physical characteristics with the one that produced α in Figure 2.

It is necessary to perform some diagnostics before choosing the model as a statistical descriptor of our data. To that end, we examine the standardized residuals obtained from our model fit. If the model is adequately specified and the estimates of the parameters are reasonably close to the true values, then the residuals should behave roughly like independent, identically distributed normal variables. If we plot their behavior over time, they should scatter approximately equally above and below zero, and not show any trends. This appears to be the case, as shown in Figure 6:

*Figure 6: Standardized residuals of the ARIMA fit to our data.*

We can assess the normality of these residuals by means of a quantile-quantile, or Q-Q, plot. This graphical tool can be used to display the quantiles of the data versus the quantiles of a normal distribution. If the data are normally distributed, they will form an approximately straight line. They are close to doing that in Figure 7, except at the edges:

*Figure 7: Normal Q-Q plot of the residuals from our ARIMA model. The results of a Shapiro-Wilk test are shown. See text for details.*

A Shapiro-Wilk test for normality can be performed to calculate the correlation between the residuals and the normal quantiles. A p-value of 0.69 is obtained, and we therefore can not reject the null hypothesis that the residuals are normally distributed.

As a final diagnostic for our model in this exercise, we look at the ACF of the residuals to establish whether any autocorrelations remain.:

*Figure 8: ACF of residuals obtained from the ARIMA fit to the turbulence intensity data.*

Figure 8 shows that is not the case, and our model could be a good candidate to be used as a proxy of the turbulent activity in a protoplanetary disk, provided its physical features are similar to those of the 3-D numerical simulation.