A Statistical Investigation to Monitor and Understand Atmospheric CFC Decline with the Spatial-Longitudinal Bent-Cable Model

Concerns about the concentrations of chlorofluorocarbons (CFCs) in the atmosphere are based on their effects on the ozone layer by catalytically destroying ozone. The recent steady decline in atmospheric concentration of CFC could be a direct result of the Montréal Protocol’s ban on CFC products, in effect since 1989. However, CFCs have long atmospheric residence times because of their low chemical reactivity, and as a consequence have already been distributed globally. To study the spatial effects and extent of the decline, we apply the proposed spatial-longitudinal bent-cable model to CFC data observed over a global detection network. The bent cable is a parametric regression model to study data that exhibits a trend change. It comprises two linear segments to describe the incoming and the outgoing phases, joined by a quadratic bend to model the transition period. For spatial longitudinal data, measurements taken over time are nested within spatially dependent locations. Here, it is useful to extend the existing longitudinal bent-cable regression to handle spatial effects. We do so in a hierarchical Bayesian framework by allowing the error terms to be correlated across space. The methodology is illustrated with applications to CFC-11 and 12 data. Our analysis reveals that: (a) there is a strong spatial relationship among all the monitoring locations across the globe; (b) both CFC-11 and CFC-12 increased significantly before entering into a transition zone; (c) after completing the transition, CFC-11 has been decreasing significantly from the atmosphere, but a slow (insignificant) decrease for CFC-12 is observed; and (d) it may take almost 5 times longer to diminish CFC-12 from the atmosphere compared to CFC-11.


Introduction
Atmospheric ozone layer protects all living things from the detrimental effects of ultraviolet (UV) radiation from the sun.In particular, about 90% of the total atmospheric ozone lies in the stratosphere, which is the major source of blocking the harmful UV radiation to reach in the earth's surface (Struijs et al., 2010).The impact of thinning ozone layer caused by chlorouorocarbons (CFCs) on human health and environment is now globally recognized (Han, Kennedy, Mackie, & Dlugogorski, 2010;Struijs et al., 2010); negative effects of increased UV exposure include skin cancer and cataracts to humans, and irreversible damages to plants and drifting organisms (e.g.animals, plants, archaea, bacteria) in the ocean's photic-zone (Khan, Chiu, & Dubin, 2009).
CFCs comprise about 90% of the Ozone Depleting Substances (ODSs) in the atmosphere, and more than 80% of these are accounted by CFC-11 and CFC-12, the most abundant CFCs in the atmosphere (Moulijn, Makkee, Wiersma, & van de Sandt, 2000).These CFCs are sources of chlorine to the atmosphere, which catalytically destroy stratospheric ozone (Molina & Rowland, 1974).In particular, each chlorine atom can break down an average of 100, 000 ozone molecules during it's 1-2 years atmospheric lifetime.Furthermore, they are strong infrared absorbers in the atmospheric window region and are very potent greenhouse gases (Zhang et al., 2010).The 1987 Montréal Protocol on Substances That Deplete Ozone Layer (World Meteorological Organisation, 2011) formally recognized the significant threats of the ODSs, and provided a mechanism to phase-out global production and consumption of CFCs along with other ODSs.
Since the Montréal protocol came into effect, emissions of CFCs have changed dramatically, and their atmospheric concentrations have either leveled off or decreased.Relative to the corresponding peak values in the late 1980s, the annual aggregated production and consumption of ODSs have declined by 95%, and annual emissions by 82% of which around 80% are accounted by CFCs (Ravishankara et al., 2008).Although the production of CFCs has decreased substantially and is expected to continue to decrease in the future, the amount of CFCs in existing products and equipments has the potential to make an important contribution to future CFC emissions (Zhang et al., 2010).In addition, due to their extended lifetimes (45 years for CFC-11 and 100 years for CFC-12 (World Meteorological Organisation, 2011)), CFCs can stay in the upper atmosphere for a substantially long period of time.During that time, winds spread them across the world.Consequently, the atmospheric concentration of CFCs is now a global concern, and it is hypothesized that they have already been distributed globally (i.e., all areas around the globe are neighbours with respect to the CFC concentration in the atmosphere).So, it is important to understand the spatial distribution of the CFC concentration around the globe by fitting an appropriate statistical model.
CFCs are monitored from different stations across the globe.Each station constitutes an individual curve for CFC trend, which is different from the others due to the differences in instruments, sampling techniques, exposure to wind, actual levels of the CFC during measurement, and so on.Therefore, the appropriate statistical method to analyze these data should involve a longitudinal model which can characterize a trend change over time.Chiu, Lockhart and Routledge (2006) and Chiu and Lockhart (2010) developed a special changepoint model to analyze a single profile that exhibits a trend change.The model is parsimonious with respect to the number of parameters, and appealing due to its simple structure, flexibility and interpretability.It comprises two linear segments (incoming and outgoing) joined by a quadratic bend.There are five parameters: an intercept and a slope for the incoming phase, two parameters to characterize the center and the half-width of the bend, and a slope for the outgoing phase.Khan et al. (2009) proposed an extension of the bent-cable model for longitudinal data, but with virtually no methodological details.They addressed the global concern of CFC-11 by fitting their longitudinal bent-cable regression model.However, they did not take into account spatial effects, an important aspect to understand the atmospheric distribution of CFCs around the globe, and did not investigate CFC-12, not to mention had considered a shorter period of follow-up for CFC-11.
Accounting for the spatial effects, we develop a statistical framework for the longitudinal bent-cable model to analyze both the CFC-11 and 12 data.We will address the following questions statistically by fitting our proposed spatial-longitudinal bent-cable model: (1) What is the spatial distribution of the CFC concentration around the globe?(2) What was the time point at which the CFC concentration took a downturn from an increasing trend?(3) Has the concentration of CFC been decreasing significantly after the transition?(4) Are there any differences and/or similarities in the atmospheric concentration of CFCs in different parts of the world?Addressing these questions could be useful not only to policy makers, but also for human awareness.For example, (a) through the estimation of the strength and extent of the distribution of CFCs around the globe, our analysis will reveal whether the concern about CFCs is global or regional, and (b) through the estimation of the transition and the rate by which CFCs have been declining from the atmosphere, our analysis will provide insight about the extent of the ozone destruction process, which, in turn, will convey information regarding precautions against the UV radiation.In a broader sense, we will also comment on the effectiveness of the Montreal Protocol in reducing the use of CFCs.
In Section 2, we describe the CFC data.In Section 3, we present our modeling framework to account for spatial effects in the longitudinal bent-cable model.This is all constructed under a Bayesian framework (Section 4).We then apply our method to the CFC-11 and CFC-12 data to address the aforementioned questions (Sections 5).In Section 6, we present a simulation study to demonstrate the performance of our proposed methodology.We summarize our findings in Section 7.

CFC Data
CFCs are monitored from different stations across the globe by the Global Monitoring Division of the National Oceanic and Atmospheric Administration (NOAA/ESRL halocarbons group, http://www.esrl.noaa.gov/gmd/),and Atmospheric Lifetime Experiment/Global Atmospheric Gases Experiment/Advanced Global Atmospheric Gases Vol. 1, No. 2; 2012 Experiment (ALE/GAGE/AGAGE, http://agage.eas.gatech.edu/dataarchive/) program sponsored by NASA.Both these programs currently use gas chromatographic (GC) technique (known as in situ program) to measure CFCs.Below we summarize the process (also see Figure 1) as described in the NOAA/ESRL website http://www.esrl.noaa.gov/gmd/hats/insitu/insitu.html.
The process of measuring CFCs begins by drawing outside air through a sampling line using a pump.The air then enters the gas sampling valve, where it flushes the valve for about five minutes.This process leads the sample to bleed down to ambient pressure.The GC analysis begins by injecting the sample into a separation column.Carrier gas (usually nitrogen or a mixture of methane and argon) is then used to sweep the sample through the separation process.Each type of molecule has a different rate of progression.Therefore, by examining the time to reach the end of the column, different components of the sample can be separated.An electron capture detector is used to produce chromatogram, the dip of which is considered as the measure of the amount of chemical present.The ALE/GAGE/AGAGE program (Prinn et al., 2000) consists of three phases corresponding to advances and upgrades in GC techniques: ALE (began in 1978), GAGE (began between 1981and 1985), and AGAGE (began between 1993and 1996).Under this program, CFCs are measured from stations located in (6) Mace Head, Ireland, (7) Cape Grim, Tasmania, (8) Ragged Point, Barbados, (9) Cape Matatula, American Samoa, and (10) Trinidad Head, California.
We take into account monthly mean CFC data based upon two considerations: (a) the availability of complete data, and (b) the necessity of a sufficiently long study period to capture the trend change over time.For (a), we choose eight monitoring stations (stations 1-8) for which data are available until September of 2010.Although the CFC monitoring process began in the late 1970's for some stations, complete data for all stations are not available until January of 1988.Therefore, we choose a study period beginning from January of 1988 to September of 2010 (273 months), which is also long enough to observe the trend change over time.Table 2 summarizes the geographical locations of the stations and the instrumentations used to record data.

Spatial-longitudinal Bent-cable Model
We index the stations by i = 1, 2, . . ., m and the months by j = 1, 2, . . ., n.We model the CFC measurement for the i th station at time t i j , denoted by y i j , by the bent-cable regression model where the bent-cable function f (t i j , θ i ) is defined by with θ i = (β i , α i ) , β i = (β 0i , β 1i , β 2i ) and α i = (γ i , τ i ) .As illustrated in Figure 2, β 0i and β 1i are, respectively, the intercept and slope of the incoming phase; β 1i + β 2i is the slope of the outgoing phase; and γ i and τ i , the halfwidth and center of the bend, respectively, characterize the transition period.Note that the beginning and end of a transition can be represented by τ i − γ i and τ i + γ i , respectively.Chiu & Lockhart (2010) also defined the critical time point (CTP) at which the bent cable changes direction (i.e., takes a downturn from an increasing trend, and vice versa) by To take into account the spatial effects, we model the error component as follows: where δ is the spatial autocorrelation coefficient that measures the strength of spatial association, W = (w ik ) m×m is the spatial weight matrix that describes the spatial arrangement of the locations (stations), and v i j ∼ N(0, σ 2 i ) is the random error component.Here, w ik is assumed known representing the a priori assumption about the spatial relationship between stations i and k, with all the diagonal elements w ii set to zeros.The matrix W is typically normalized in an attempt to constrain δ in the interval (−1, +1) (Elhorst, 2010).One such normalization technique is to divide the elements of W by its largest characteristic root.Note that W = I, where I is the identity matrix, implies no spatial effects, and constitutes the longitudinal bent-cable model of Khan et al. (2009).
Equations ( 1)-( 3) constitute our spatial-longitudinal bent-cable model.The Bayesian hierarchical formulation of the model involves three levels.The first level (Level 1 or within-individual level) characterizes an individual trajectory by taking into account correlation and variation among the repeated measurements over time.The second level (Level 2 or between-individual level) models the station-specific regression coefficients θ i 's to take into account different patterns manifested by different trajectories.The third level (Level 3) quantifies prior information for the random quantities in Levels 1 and 2. Let y i = (y i1 , y i2 , . . ., y in ), y = (y 1 , y 2 , . . ., y m ) , and f i j ≡ f (t i j , θ i ).Also, let Θ be the vector of all the model parameters collectively.Then, based on our choices of the distributions for the relevant quantities, the hierarchical formulation of the model can be expressed as [ where I m − δW n in (4) represents the Jacobian term taking into account the endogeneity of m k=1 w ik y k j (Anselin, 1988, p. 63), and g(•) is used to denote a density function.In the hierarchical formulation of the model, we use N p , LN p , W, G and U to represent p-variate normal, p-variate lognormal, Wishart, gamma and uniform distributions, respectively, with all the hyperparameters assumed known.Note that μ β = (μ 0 , μ 1 , μ 2 ) and Σ β are, respectively, the mean and covariance of β i , which characterize the global (population) trends in the incoming and outgoing phases: μ 0 and μ 1 are, respectively, the population intercept and slope for the incoming phase, μ 1 + μ 2 is the population slope for the outgoing phase, and Σ β provides information about the variability of the linear parameters β i .Similarly, μ α = (μ γ , μ τ ) characterizes the global transition.Since our assumption for α i involves a lognormal distribution, the medians M γ = exp {μ γ } and M τ = exp {μ τ } can be regarded as the half-width and the center of the bend, respectively, for the population.We can also use the standard deviations of γ i and τ i to describe the between-station variability for the transition parameters.

Bayesian Inference and Implementation
Bayesian inference is based on the posterior distribution of a parameter given the data: the density function provides the behavior of the parameter, and the posterior mean (median if the posterior density is noticeably asymmetric) and variance can be considered a point estimate of the parameter and the uncertainty of the estimate, respectively.The posterior distribution of a parameter, say μ β , denoted by π(μ β |y), can also be used to determine points p 1 and p 2 such that p 2 p 1 π(μ β |y)dμ β = 0.95.This gives a 95% Bayesian credible interval, (p1, p2) for μ β .
Bayesian inference is carried out here by Markov chain Monte Carlo (MCMC) technique.We employ the Metropolis within Gibbs algorithm (Smith & Roberts, 1993) to approximate the posterior densities of the parameters.The algorithm is based on drawing random samples from the the full conditionals of each of the components of Θ; we work out the full conditionals for each component to implement this algorithm.We have written our own code in C to generate MCMC samples, which we subsequently analyze using the "coda" (Plummer, Best, Cowles, & Vines, 2006) and "MCMCpack" (Martin, Quinn, & Park, 2011) packages in R (R Development Core Team, 2012).

The Spatial Matrix and Model Building
The elements of the spatial matrix W is defined by where i ∼ k denotes the contiguity of station i with station k.The specification of the spatial weights is an important problem and is often chosen on an ad hoc basis (Anselin, LeGallo, & Jayet, 2008).
For the CFC data, since winds can spread CFCs across the world, there is no clearly defined boundary by which the locations of the monitoring stations can be differentiated.Therefore, we consider several models by defining the spatial configurations based on the distances between the stations, and choose the one for which the estimate of the deviance information criterion (DIC) is minimum (Spiegelhalter, Best, Carlin, & Van der Linde, 2002).Specifically, letting d ik be the distance between stations i and k, we consider the following spatial configurations: Model 1: W = I (i.e., no spatial effects), Note that for Model 4, w ik = 1 ∀ i k, hypothesizing that CFCs are distributed around the globe.Note also that we normalize W by dividing all of its elements by its largest characteristic root.
For MCMC, we construct two Markov chains each of 5, 000, 000 iterations to approximate the posterior density.The initial 100, 000 iterations are discarded as burn-in, and the inferences are based on every 200 th iteration of the chains (thinning), resulting in a total of 24, 500 iterations per chain.In Table , the DICs are given for each model.The smallest DIC results from Model 4 for both CFC-11 and 12. Hence, according to the DIC values, Model 4 is estimated to be the model that would best fit the data among the models under consideration.This also supports the notion of the global distribution of the CFCs (i.e., each station is a neighbour of the remaining ones).Now, we report the results for Model 4 in Section 5.2.

Results
The trace and density plots for δ for both types of CFCs are presented in Figure 3.The lack of any trend in the trace plot indicates good mixing.The density plots display no signs of multimodality.These plots show stationarity of the Markov chains.We also test the stationarity of the chains for δ by the Gelman-Rubin statistic R (Gelman & Rubin, 1992) and the Geweke Z score (Geweke, 1992).Values of R are 1 for both CFC-11 and 12, and the |Z| scores are 1.28 and 0.53 for CFC-11 and 12, respectively, none of which are greater than the critical value 1.96 at 5% level of significance.So, these two diagnostic criteria also provide no evidence against stationarity of the chains for δ.For CFC-12, the posterior mean for δ is 0.78 with 95% credible interval (0.76, 0.80), and that for CFC-11 is 0.52 Vol. 1, No. 2; 2012 with 95% credible interval (0.48, 0.58).These indicate significant spatial autocorrelations for both types of CFCs, with stronger association for CFC-12.The stronger association for CFC-12 is also an evidence of its extended lifetime in the atmosphere: 45 years for CFC-11 and 100 years for CFC-12.We now turn to the behaviour of the station-specific trajectories described by the regression coefficients β i and α i , and the between-station variability characterized by Σ β and Σ α .The trace plots and convergence diagnostic criteria (not shown) provide no evidence against the stationarity of the chains for β i , α i , Σ β and Σ α .Table 4 summarizes the posterior characteristics of the station-specific fits, and Figure 7 5.We see small variation in the slope parameters β 1i and β 2i across the stations.This suggests that the rates at which CFCs increased and decreased in the incoming and outgoing phases are similar across stations.However, we see some variation for the transition parameters γ i and τ i (i.e., in the times to transition zones).
www.ccsenet.org/ijspInternational Journal of Statistics and Probability Vol. 1, No. 2; 2012 increase of CFC-11; (3) on average, CFC-11 completed the transition between 08/1988 and 09/1995, whereas CFC-12 between 05/1988 and 04/2001; (4) CFC-11 has been decreasing significantly after completing the transition, but the rate at which CFC-12 has been decreasing is not significant; and (5) although the rates of increase and decrease are similar across stations for each type of CFC, there are considerable variations in the times to transition zone.Though the two types of CFCs may pose potentially differential health effects, the above findings also indicate a much more severe global concern for CFC-12 with respect to its concentration in the atmosphere.

Simulation
In this section, we present a simulation study to demonstrate the efficacy of our spatial-longitudinal bent-cable methodology.In particular, we supplement the motivation of our methodology by illustrating the importance of properly taking into account the spatial component.To this end, we present a scenario, where, in reality, spatial dependency exists among the units.We then analyze the data using the true model, as well as models with misspecified spatial configuration.
For each simulation, 500 data sets are generated, and 10, 000 MCMC iterations are used to approximate posterior distributions per set.DICs and posterior summaries for each parameter are then averaged over the 500 sets.In addition to the DICs, relative mean square error (RMSE) of the posterior summaries is used to evaluate the performance of the analytic procedures.RMSE takes into account both bias (measured by the relative bias, RB) and variability (gauged by the relative standard deviation, RSD), according to the formula RMSE = RB 2 + RSD 2 , where RB = bias/true parameter value, RSD = SD/true parameter value, and SD is the standard deviation of the posterior means/medians across the 500 simulated data sets.Note that the smaller the RMSE, the more accurate an estimator is.The above simulation results demonstrate the importance of modelling the spatial component when, in reality, spatial dependency exists among the units.Moreover, smaller bias and highly accurate (smaller RMSE) estimates indicate that our proposed methodology performs well to analyze data that resemble those from the CFC study.

Conclusion
Under a spatially dependent process, each unit is typically affected by those of the neighbouring units.Therefore, the assumption of independence across the neighbouring units might be unrealistic, especially when they represent geographical locations.Moreover, addressing the spatial effects provide important insights about such a process.In this article, we propose an extension of the longitudinal bent-cable model by taking into account spatial effects.We have tailored our work especially for the scientific context of the CFC data.Due to their extended lifetimes, CFCs persist long enough in the atmosphere, and consequently are believed to have spread across the world.Therefore, CFCs monitored from one station may depend on those from another station, giving rise to a presumed spatially dependent longitudinal process.Our methodology provides understanding not only of the spatial distribution, but also of the global threat that CFCs may pose to all living organisms.It also reveals useful information regarding the atmospheric CFC decline throughout the globe.
Since the Montréal protocol came into effect, a global decrease in the CFCs is monitored and confirmed by our analysis.Note that the Montréal Protocol contains an extended CFC phase-out schedules-1996 for developed countries and 2010 for developing countries.Thus, many countries at various geographical locations continued to contribute CFCs to the atmosphere during the 273 months in our study period.In our analysis, this fact is reflected with a slow decrease in CFC concentrations from the atmosphere.In fact, our analysis does not reveal a

Figure 1 .
Figure 1.The analytic technique of separation and detection of CFCs using gas chromatograph

Figure 2 .
Figure 2. A graphical representation of the bent-cable function.The trend change is modeled by a quadratic bend characterizing the transition period (τ − γ, τ + γ).Two linear models with slopes β 1 and β 1 + β 2 characterize the incoming and outgoing phases.The time point at which the cable takes a downturn from the increasing trend is defined by CTP = τ − γ − 2βγ/β, and is marked by the dotted vertical line

Figure 3 .
Figure 3. Trace and density plots for the posteriors of the spatial autocorrelation parameter δ from two chains

Figure 4 .
Figure 4. CFC-12 analysis-trace and density plots for the posteriors of the population parameters from two chains displays one representative fitted curve for each of CFC-11 and CFC-12.The figure shows that our model fits the data well, with the observed data and individual fits agreeing quite closely.All the station sites exhibit significant increase of both types of CFCs before the transition took place.The rate of increase was maximum in South Pole ( β1 = 1.49 and 0.98 for CFC-12 and CFC-11, respectively).The transition ended in around 2001 for CFC-12, and between 07/1994 and 12/1996 for CFC-11.After completing the transition, CFC-11 decreased at a higher rate across all stations.The most rapid decrease took place in Niwot Ridge ( β 1 + β 2 = −0.11and −0.178 for CFC-12 and CFC-11, respectively).Posterior characteristics of the standard deviations describing the between-station variability are summarized in Table

Table 1 .
Geographical locations (latitude, longitude, and elevation in meter above sea level (masl)) of the stations and the instrumentations used to record data

Table 2 .
Model comparison for the CFC data using DICs

Table 6 .
Simulation results for data generated from Model 4 with n = 273 and m = 8: absolute relative bias (RB), absolute relative standard deviation (RSD) and relative mean square error (RMSE) of the population regression coefficients and the spatial autocorrelation coefficient, calculated by considering average of 500 posterior means Numerical results for the population regression coefficients and spatial autocorrelation parameter are summarized in Table6.For μ 0 , both |RB| and RMSE are very close to zero regardless of the assumption about the spatial configuration.Smallest |RB| and RMSE for each of the other parameters are observed when we analyze using the true model (i.e., Model 4), whereas the largest |RB| and RMSE are observed if we completely ignore spatial dependency (i.e., Model 1).Similar results are obtained for σ 2 i 's (Table7), and Σ β and Σ α (Table7), that is, the smallest |RB| and RMSE occur when we analyze using the true model, and the largest |RB| and RMSE if we ignore spatial dependency.The only exception occurs for (Σ α ) 12 (Table7), for which the smallest |RB| and RMSE are observed for Model 2. The average DICs are 7893, 7720, and 7379 for Models 1, 2 and 4, respectively.This indicates that the best fit is achieved by Model 4 (smallest DIC), and a comparatively poor fit if we completely ignore the spatial dependency (largest DIC).

Table 7 .
Simulation results for data generated from Model 4 with n = 273 and m = 8: absolute relative bias (RB), absolute relative standard deviation (RSD) and relative mean square error (RMSE) of σ 2 i 's, calculated by considering average of 500 posterior medians

Table 8 .
Simulation results for data generated from Model 4 with n = 273 and m = 8: absolute relative bias (RB), absolute relative standard deviation (RSD) and relative mean square error (RMSE) of the covariance matrices Σ β and Σ α , calculated by considering average of 500 posterior means for the covariances (Σ β ) i j and (Σ α ) i j for i j, and average of 500 posterior medians for the variances (Σ β ) ii and (Σ α ) ii