The effects of spatial variability of nitrous oxide emissions from grazed pastures on the sampling distribution of chamber measurements

Summary

Chamber sampling is a common method for measuring nitrous oxide (N₂O) emissions from agricultural soils. However, for grazed pastures, the patchy nature of urine deposition results in very high levels of spatial variability in N₂O emissions. In the present study, the behaviour of the sample mean was examined by simulating a large number (9999) of random N₂O chamber samples under different assumptions regarding the underlying N₂O distribution. Using sample sizes of up to 100 chambers, the Central Limit Theorem did not apply. The distribution of the sample mean was always right-skewed with a standard deviation varying between 12·5 and 135% of the true mean. However, the arithmetic mean was an unbiased estimator and the mean of the sample mean distribution was close to the true mean of the simulated N₂O distribution. The properties of the sample mean distribution (variance, skewness) were affected significantly by the assumed distribution of the emission factor, but not by distribution of the urine patch concentration. The geometric mean was also investigated as a potential alternative estimator. However, although its distribution had lower variance, it was also biased. Two methods for bias correcting the mean were investigated. These methods reduced the bias, but at the cost of increasing the variance. Neither of the bias-corrected estimators were consistently better than the arithmetic mean in terms of skewness and variance. To improve the estimation of N₂O emissions from a grazed pasture using chambers, techniques need to be developed to identify urine patch and non-urine patch areas before sampling.

INTRODUCTION

In New Zealand, nitrous oxide (N₂O) emissions from agricultural soils account for 13·6% of anthropogenic greenhouse gas emissions on a CO₂-equivalent basis (Ministry for the Environment 2014). Of these, 56·3% are direct emissions from grazing animals’ excreta, 18·4% direct emissions from other sources of nitrogen (N) to soil (e.g. fertilizer, animal waste and crop residues) and 25·4% from indirect emissions. The present study focuses on direct emissions from urine patches, as the emissions from dung are comparatively small (van der Weerden et al. 2011).

Nitrous oxide emissions have a high degree of spatial and temporal variability due to the complexity of the various interacting processes and variability of the driving factors (Mosier et al. 1996). A number of different factors, such as topography and soil properties, have been linked to spatial variability of N₂O emissions (Parkin 1987; Pennock et al. 1992; Velthof et al. 1996, 2000; Ball et al. 1997; Grant & Pattey 2003; Yanai et al. 2003; Turner et al. 2008). In grazed pasture systems there is additional spatial variability due to the highly non-uniform nature of urine deposits (i.e. small areas of high N-concentration) (Bolan et al. 2004). Previous studies have looked at the distribution of urine patches both experimentally (Moir et al. 2011) and theoretically (Pleasants et al. 2007); however, these studies did not investigate N₂O emissions. Giltrap et al. (2014) investigated the spatial variability of N₂O emissions within a field and found differences in the emission response to applied urine at different locations.

Closed chambers are a relatively inexpensive method for measuring N₂O emissions. In this method a number of chambers are inserted into the soil. The tops are then sealed to form closed chambers and the air within the headspaces sampled at regular intervals. The N₂O flux is calculated from the change in N₂O concentration over time. Sealing and sampling of the chambers can be automated or performed manually. Manual chambers in particular are relatively cheap to deploy, although the labour requirements increase as more chambers are used. International researchers have published guidelines on the use of chambers for N₂O measurements (de Klein & Harvey 2012). These guidelines do not specify a minimum number of chambers, but do state that chambers should cover the largest total area practicable while providing information at the smallest scale required.

Integrative, micrometeorological methods are an alternative to chamber measurements. These methods measure emissions over a wider spatial scale and are therefore less affected by spatial variability; however, there can be limitations due to weather conditions and topography. Also atmospheric N₂O concentrations tend to be quite low, and may be close to the detection limit of many sensors.

Nitrous oxide emissions are frequently expressed in terms of an emission factor (EF). That is, the amount of N lost as N₂O (in excess of any background emission rate) as a fraction of the amount of N input. However, EF values vary with soil properties, weather conditions and other factors. Almost 90% of the uncertainty in New Zealand's inventory value for N₂O emissions from agricultural soil in 2002 (the relative error of the 95% confidence interval is –42% to +74%) was attributed to uncertainty in the value of EF (Ministry for the Environment 2014).

Studies for determining EF for animal urine in New Zealand used closed chambers to which known amounts of urine-N were applied (de Klein et al. 2003). This eliminated the variability due to uneven urine deposition, although there was still high variability in N₂O emissions due to variability in soil properties (CV ranged from 22 to 55% for different sites). However, these types of experiments may not be completely representative of what occurs during an actual grazing event (e.g. effects such as compaction caused by animal treading and differences in urine concentration between animals are neglected). Some studies have been conducted to measure N₂O emissions from actual grazed pastures (Saggar et al. 2004, 2007; Hoeft et al. 2012; Luo et al. 2013). These studies generally (although not always) used a larger number of chambers than the applied urine studies to account for the additional spatial variability of urine patch distribution. The question is how well can such a sample estimate the mean emissions from a field?

Many common statistical techniques assume that the data come from a normal distribution. This is clearly not the case for sampling a grazed pasture as there is both the probability of any given chamber landing on a urine patch combined with the variability of emissions from a urine patch. Neither is there an obvious transformation to normalize the data. However, according to the Central Limit Theorem (CLT), the arithmetic mean of a number (n) of independent random variables, each with a finite mean (μ) and variance (σ ²), will be approximately normally distributed with mean = μ and variance = σ ²/n as n → ∞ (Bain & Engelhardt 1992).

In the case of N₂O emissions from grazed pastures, the distribution is highly skewed and it cannot be taken for granted that the CLT can be applied for practical sample sizes. To estimate the probability density function of the sampling distribution, computer simulations were used to generate a large number of samples from the underlying N₂O emission distribution using a number of different sample sizes and assumptions about the properties of the underlying distribution. As the number of samples taken was large, the distribution of these samples should closely approximate the probability density function of the sampling distribution, and this distribution shall be referred to as the empirical sampling distribution (ESD). The properties of the ESDs were then examined to observe how they changed with different sample sizes and with changes to the underlying emission distribution. The arithmetic mean of each sample was initially used as the estimator of the mean of the underlying N₂O emissions; some alternative estimators were then examined for comparison.

While the current research has been framed in terms of N₂O measurements, the results may also be applicable to other situations where point measurements are made on a grazed field and the results are likely to be influenced by the presence or absence of a urine patch (e.g. nitrate leaching).

MATERIALS AND METHODS

Nitrous oxide emission model

For the simulations it was assumed that a number of chambers, n, were used to estimate the N₂O emissions from a field following a grazing event. Each chamber had a probability, p, of landing on a urine patch, where p represents the proportion of the field covered by urine patches. This assumed that the total area of chambers was sufficiently small that p could be considered the same for all chambers. The possibility of overlapping urine patches or chambers having only part of their area on a urine patch was not considered.

As the present study was concerned with spatial rather than temporal variation, it was assumed that each chamber measured the total N₂O emitted over time accurately and that the duration of the measurement campaign was long enough to capture all the additional N₂O emissions due to the grazing event.

The following equations describe the N₂O emissions measured in a single chamber:

(1)$$\eqalign{{\rm N}_2 &{\rm O}-{\rm N} = B\comma\,\cr & \hbox{if the chamber was not on a urine patch}}$$

(2)$$\eqalign{{\rm N}_2 &{\rm O} - {\rm N} = B + {\rm EF} \times U\comma\ \cr & \hbox{if the chamber was on a urine patch}}$$

where N₂O–N is the total N emitted as N₂O (kg N/ha), B is the background N₂O emissions over the measurement period (kg N/ha), EF is the emission factor (which is specific to a given location and time period, in kg/kg) and U is the N concentration of the urine patch (kg N/ha).

This formulation assumed that EF and U were independent. Some studies have suggested that EF may increase with increasing N concentration (Kim et al. 2013); however, de Klein et al. (2014) suggested that this dependency only occurs at low values, while Kelliher et al. (2014) found no dependency.

Simulations

Empirical sample distributions were generated by simulating 9999 random samples in R (R Foundation for Statistical Computing, http://www.r-project.org/) and were generated for different sample sizes, urine coverages and assumptions about the distribution of U and EF. Table 1 shows the parameters investigated and the range of values considered. The maximum value of n considered was 100 chambers, which is higher than generally used in practice but still potentially feasible. Increasing the maximum number of chambers will eventually result in the sample mean distribution approaching a normal distribution (by the CLT). However, in this case it would no longer be necessary to use the methods described in the present study.

Table 1. Parameters varied to investigate the effect on the sample mean distribution

The background emissions (B) were considered to have a constant value of 1 kg N/ha. While there is likely to be some variability in background emissions, this was not considered likely to be a major contributor to overall spatial variability.

Different underlying distributions were investigated for the urine concentration and EF. For urine, only the constant and normal distributions were considered, whereas for the EF, gradient and log-normal distributions were also examined.

In all cases the distributions were parameterized to have the same expected value (0·01 for EF and 1000 kg N/ha for urine concentration). The expected value of 0·01 for EF corresponds to the New Zealand specific value used in the national inventory for calculating N₂O emissions from urine applied to soil (Ministry for the Environment 2014). However, measured EF values have a high degree of variability with a positive skew (de Klein et al. 2003). Therefore, for the normal distribution a standard deviation of 0·005 was selected to give the largest variability possible while keeping the EF positive within two standard deviations of the mean.

Few studies have measured the N loading in urine patches in terms of g N/ha (Selbie et al. 2015). Haynes & Williams (1993) give a range of 8–15 g/l for urine N concentration, with urination volumes in the range 1·6–2·2 l and patch size 0·16–0·49 m² for dairy cattle. The expected value of 1000 kg N/ha for urine patch N concentration in the present study is consistent with this and is a value often used in artificial urine patch experiments (e.g. Luo et al. 2008). The distributions considered are described briefly below:

• Constant distribution (Const) – assumes there is no variability in the underlying quantity. While this is unrealistic, it may be a reasonable approximation when the variability of the quantity is negligible compared with other sources of variability. For the urine concentration the constant value was 1000 kg N/ha, while for EF the constant value was 0·01.

• Normal distribution (Norm) – commonly occurs in nature and many standard statistical methods assume a normal distribution. For the urine concentration the normal distribution was assumed to have mean = 1000 kg N/ha and standard deviation = 200 kg N/ha. For the EF distribution the mean was 0·01 and the standard deviation 0·005.

• Gradient distribution (Grad) – this distribution simulates the case where the variation in EF is due to correlation with some property (e.g. soil moisture) that varies linearly across the field. Sampling a random location is therefore equivalent to drawing a sample from a uniform distribution. In these simulations, the minimum and maximum EF values were 0 and 0·02, respectively.

• Log-normal distribution (LNorm) – the skew in the chamber samples may not be due solely to the patchy nature of the urine deposits, but may also reflect a skewness in the EF distribution. Studies looking at N₂O emissions from known applications of urine N often find positively skewed results, suggesting this may be the case (e.g. van der Weerden et al. 2011). Log-normal distributions frequently occur when several factors act in a multiplicative fashion. For example, N₂O emission rates could be affected by soil moisture, temperature and availability of substrates. The expected value of the log-normal distribution is given by:

  (3)$$E(X) = \exp \left( {\mu + \displaystyle{{\sigma ^2} \over 2}} \right)$$
  where μ is the log-mean and σ is the log-standard deviation. For the log-normal EF simulations the log-mean was set to –5·105 and the log-standard deviation to 1. There are an infinite number of combinations of μ and σ that could be selected that satisfy E(X) = 0·01. The particular combination selected here resulted in a 95% confidence interval for EF of 0·00085–0·043. Note that due to rounding, the expected value of this distribution was 0·017% higher than the other distributions.

Properties of the empirical distribution

Nitrous oxide emissions from grazed pastures are highly skewed and it cannot be taken for granted that the CLT can be applied for practical sample sizes. Empirical sample distributions were generated for a number of estimators (arithmetic mean, geometric mean and two methods of biased-corrected geometric mean of a sample) using different sample sizes and assumptions about the distributions of EF and U (e.g. constant, normal, gradient and log-normal). The influence of these changes on the properties of the sampling distribution was then investigated. The properties of the ESDs investigated were the mean, variance, relative bias (rBias), probability that given a sample was an underestimate of the true mean (pUnder) and skewness.

The true mean of the N₂O emissions distribution is the expected value given by:

(4)$$E({\rm N}_2 {\rm O}) = E(B + p{\rm EF} \times U) = B + pE({\rm EF})E(U)$$

assuming that the EF is independent of the urine concentration. All the distributions examined were parameterized, so that E(EF) = 0·01 and E(U) = 1000 kg N/ha. Therefore, for an unbiased estimator the mean of the ESD should = 1 + 10p.

When the CLT applies, the arithmetic mean of a sample should have variance = σ ²/n. However, even when the CLT does not apply it is desirable for an estimator to have a low variance.

rBias was defined as:

(5)$$r{\rm Bias} = ({\rm Mean(ESD)} - E({\rm N}_2 {\rm O}))/E({\rm N}_2 {\rm O})$$

where Mean(ESD) is the mean of the empirical sampling distribution and E(N₂O) = 1 + 10p as shown above.

The pUnder is the proportion of samples that underestimate the true mean. In an unbiased, symmetrical distribution this should be 0·50. A pUnder that is much higher or lower than 0·50 indicates either asymmetry (skew) and/or bias. Note that in the absence of any variability in the urine concentration or EF, pUnder could be calculated using the binomial distribution.

Finally, skewness is defined as the third standardized moment of a distribution. Skewness should be 0 for a symmetric distribution (such as the Normal), positive or negative values indicate positive (right) or negative (left) skew, respectively.

It should be noted that a skewed distribution is not the same as a biased distribution. A biased distribution is one in which the expected value (average over a large number of samples) differs from the true mean. However, a skewed distribution could still produce the correct average over a large number of samples, even though there is a high probability of any individual sample being an under- (or over-) estimate.

Modelling the distribution properties

Having generated ESDs, the next step was to model how the properties of the sample distributions varied with sample size and with the properties of the underlying population distribution. The obvious explanatory variables were p and n, as well as the distributions assumed for the urine patch concentration (U _dist) factors and the emission factor (EF_dist). However, some non-linear terms also suggested themselves. The sampling distribution was related to the binomial distribution as the number of urine patches sampled would follow a binomial distribution. The variance of a binomial distribution is npq = np(1 − p) = np–np ². Therefore np and p ² were included. When the CLT applies the variance of the sample mean depends on 1/n, so 1/n was added to the list of explanatory variables. The terms n ², 1/p, n/p and p/n were also considered.

For each variable the interactions with EF_dist, U _dist, and the EF_dist : U _dist interaction were considered. Thus the full model was:

(6)$$\eqalign{{\rm Property} \sim & {\rm EF}_{\rm dist} \,{\ast}\, {\it U}_{\rm dist} \,{\ast} ({\it n} + {\it p} + {\it np} + {\it n}^2 + {\it p}^2 \cr & + 1/{\it n} + 1/{\it p} + {\it p}/{\it n} + {\it n}/{\it p})}$$

where ‘*’ indicates that both main and interaction effects were considered.

The best model was initially selected using stepwise selection (with the Akaike Information Criterion [AIC] as the selection criterion), in both the forward and backward directions. However, in many cases the selected model still contained a large number of explanatory variables. These models were replaced with simplified models when it could be demonstrated that explanatory variables could be dropped while retaining satisfactory model performance (i.e. good explanatory power, lack of structure and constant variance in the residuals).

A standard linear model was applied for all properties except pUnder, for which a binomial generalized linear model was used.

Alternative estimators

Nitrous oxide emissions from grazed pastures are highly skewed and therefore sample sizes will probably need to be large before the CLT can be applied. However, it is sometimes possible to use data transformations to improve the normality of the data. In particular, log-transformations are commonly used to reduce positive skew.

The geometric mean (g) is the mathematical equivalent of back-transforming the mean of the log-transformed data and is given by:

(7)$$g = \left( {\mathop \prod \limits_{i = 1}^n {x_i}} \right)^{1/n} = \exp \left( {\displaystyle{{\mathop \sum \nolimits_{i = 1}^n {\log(x_i)} \over n}}} \right)$$

where x _i is the ith observation and n is the number of observations.

Therefore, g could potentially be used as an alternative estimator to the arithmetic mean (a). Unfortunately g is a biased estimator, that is, its expected value is less than the true mean of the N₂O distribution.

Assuming that the underlying N₂O distribution was actually log-normal, then from Eqn (3) a reasonable correction factor might be exp(σ ²/2). Therefore the second alternative estimator investigated was:

(8)$$g_{c1} = g \times {\exp(s^2/2)}$$

where s ² is the sample variance.

Neyman & Scott (1960) gave bias-correction formulae for a number of variable transformations. However, for the log-transformation this requires solving a Bessel function of an imaginary argument. Given that the N₂O distribution is not precisely log-normal a slightly simpler formula was used to calculate the correction factor: exp((1–1/n)s ²/2). When the Taylor expansion is applied to this correction factor it only differs from that of Neyman & Scott (1960) for the third and higher terms. However, it does have the property that it approaches exp(s ²/2) as n → ∞. Therefore the final estimator considered was:

(9)$$g_{c2} = g \times {\exp}\left( {\left( {1 - \displaystyle{1 \over n}} \right)\displaystyle{{s^2} \over 2}} \right)$$

All simulations and analyses were performed using the R statistical software package version 3.0.2

RESULTS

Sample mean

Using the arithmetic mean (a) as the estimator for each simulated N₂O sample and then fitting the sample mean (the mean of the 9999 randomly generated samples) to the proportion of urine coverage gave the following relationship:

(10)$$\eqalign{&{\rm Mean}(a) = 1 {\cdot} 01( \pm 0 {\cdot} 001) + 9 {\cdot} 9( \pm 0 {\cdot} 01)p, \cr & \quad R_{{\rm adj}}^2 = 0 {\cdot} 997}$$

where Mean(a) is the sample mean (using the arithmetic mean as the estimator for each sample).

This agrees fairly well with the true mean as expressed in Eqn (4), indicating that the arithmetic mean is largely unbiased irrespective of the number of chambers used or the underlying distributions of the urine concentration or EF.

Figure 1 shows the modelled means (values calculated from Eqn (10)) plotted against the observed values (simulated sample means) as well as the residuals plotted against fitted values, n and p. There was an obvious decrease in the variability of the residuals as n increased. It was also apparent that the largest residuals correspond to the log-normal EF distribution. However, the largest relative residual was still quite small at only 3% of the fitted value.

Fig. 1. Plots comparing values fitted using model (10) with the observed mean of the ESD for a (arithmetic mean of a sample): (a) Fitted value v. observed mean, (b) residuals v. fitted values, (c) residuals v. number of chambers (n), (d) residuals v. urine coverage (p).

Variance

Figure 2 shows a plot of Var(a) × n v. n. When the CLT applies Var(a) → σ ²/n, in which case Var(a) × n should be constant with respect to n. From Fig. 2(a) this appears to be the case for all but the log-normal EF distribution. The log-normal EF distribution shows variability in Var(a) × n for different n (but not in any easily to discernible pattern, Fig. 2(b)).

Fig. 2. The variance of the ESD for a × n (arithmetic mean times the number of chambers) plotted against the number of chambers (should be constant if CLT applies). (a) All EF distributions and (b) log-normal EF only.

The variance could be quite well explained by a model involving p/n, EF _DIST and their interaction (Eqn (11)).

(11)$$\eqalign{&{\rm Var}(a) \sim p/n + {\rm EF}_{{\rm DIST}} + {\rm EF}_{{\rm DIST}} :p/n, \cr & \quad R_{{\rm adj}}^2 = 0 {\cdot} 9753}$$

According to Eqn (11) the variance depends on the EF distribution and increases with p/n. The model is a reasonable approximation for the constant, normal and gradient EF distributions, but there are some large residuals for the log-normal EF distribution (Fig. 3). The variance was always highest for the log-normal EF distribution. For low p/n values the constant, normal, gradient EF distributions have similar variances. As p/n increases the variances of the different EF distributions diverge.

Fig. 3. Var(a) v. p/n. Straight lines indicate model (11) estimates.

Equation (11) is useful for determining the level of precision with which the mean can be determined. For example, suppose P = 0·05, n = 20. The variance of the sample mean would be 0·24 for a constant EF, corresponding to a standard deviation of 0·49 kg N/ha or 33% of the true mean. If the EF was normally distributed instead, then the standard deviation would be 38% of the true mean.

Relationship between mean and variance

The mean and variance appear to behave in a manner consistent with the CLT (with the exception of the variance when using the log-normal EF distribution). However, another requirement for a normal distribution is that the variance and mean should be independent. This is also an important assumption of linear regression.

From inspection of models (10) and (11) it can be seen that both the mean and variance depend on p. Therefore, it seems likely that the mean and variance will be correlated. This is confirmed in Fig. 4. However, the slope of the relationship between variance and mean decreases as n increases, and for n = 100 the variance is almost constant with respect to the mean.

Fig. 4. Var(a) v. mean(a) for (a) constant, (b) normal, (c) gradient and (d) log-normal EF distributions.

The significance of the slopes of the relationships between variance and mean were determined by fitting a linear regression model of the variance of the ESD. The predictors were the mean of the ESD, a 24-level factor (EF_n) with each level corresponding to a combination of EF distribution and number of chambers, and their interactions (Eqn (12)). The following relationship was fitted:

(12)$${\rm Var}(a) \sim {\rm EF}\_n + {\rm EF}\_n:{\rm Mean}(a) - 1$$

With the exception of the levels corresponding to n = 100, all the coefficients for the EF_n interaction with the ESD mean were significantly different from 0 (at a 0·01 confidence level). None of the EF_n:Mean interaction coefficients were significant for n = 100 (at a 0·1 confidence level). Therefore the variance of the sample mean depends on the sample mean for the chamber numbers used in practice. However, at 100 chambers the variance was independent of the mean. This relationship held across all EF distributions.

As a consequence of this dependence of the variance on the mean (for n < 100) ordinary linear regression methods (including ANOVA) should not be performed on samples of this type. However, this can be corrected using a weighted regression provided the relationship between the mean and the variance is known. If the variance can be written as:

(13)$${\rm Var}(Y_i) = \sigma ^2 f\,(\mu _i )$$

then the weighted regression can be performed using 1/f(μ _i) as weights. There are several options for f(μ _i). Using model (12) gives f(μ _i) = $\bar x$. This assumes the same n and EF for all samples and will be slightly inaccurate due to neglecting the constant term. If the EF distribution is known, then the more precise form of f(μ _i) = $a + b\bar x$ can be used. Alternatively, model (11) gives f(μi) = p/n. The constant terms in model (11) are not statistically significant, so can justifiably be neglected. If the urine coverage can be assumed to be constant, then it can be incorporated into σ ² leaving f(μi) = 1/n.

Relative bias

Earlier it was found that the mean of the ESDs were very close to the theoretical population means (Eqn (10)). Therefore bias should be small. Figure 5 shows the bias (relative to the mean) plotted against p and n.

Fig. 5. Relative bias v. (a) number of chambers and (b) urine coverage.

It is difficult to discern any patterns from these plots, although the variability in the rBias reduces with increasing n. However, the most extreme rBias values are only –2% and +3%, both of which occur for log-normal distributions.

Bias becomes important when it is high relative to the variability. One measure of the relative importance of the bias would be to take the ratio of the bias to the standard deviation. This is equivalent to the ratio of the rBias to the CV. Figure 6 shows that the magnitude of the rBias is at most 2% of the CV. That is, the absolute bias is at most 2% of the standard deviation of the sample mean and therefore not a major concern.

Fig. 6. Ratio of relative bias to CV plotted against CV.

Probability that a given sample was an underestimate of the true mean (pUnder)

With the exception of some constant EF distribution points, pUnder was > 0·5 (Fig. 7). Both pUnder and the spread of pUnder decreased with increasing number of chambers and urine coverage. Even with 100 chambers the probability of underestimating the mean can be much larger than 0·5. Given that no evidence of bias was found in the sample mean distribution, these high pUnder values indicate a right-skewed distribution.

Fig. 7. Probability of a sample mean being less than the true mean by (a) number of chambers and (b) urine coverage.

pUnder is a probability and therefore bounded by 0 and 1. Therefore a binomial generalized linear model with a logit† link function (Eqn (14)) was used for modelling:

(14)$$\qquad\quad\eqalign{&{\rm logit} (\,p {\rm Under}(a))\sim 1/n \pm p/n \cr &+ {\rm EF}_{\rm dist} (\hbox{83\% Deviance explained})}$$

Model (14) produced a reasonable fit for all but the constant EF distribution (Fig. 8). As an example, a sample from a normal EF distribution with n = 20 and P = 0·05 would have log-odds of 0·42 of being an underestimate. This corresponds to a pUnder of 60%, indicating a slight right-skew in the sample mean distribution.

Fig. 8. (a) Modelled v. observed pUnder and (b) deviance residuals v. fitted values for model (14).

Skewness

Skewness is the third standardized moment and relates to the lack of symmetry in a distribution. Skewness should be 0 for a symmetric distribution (such as the normal). Figure 9 shows the observed skewness of the sample mean distribution plotted as a function of number of chambers and urine coverage. Unsurprisingly, skewness was largest for the right-skewed log-normal EF distribution. Skewness reduced with increasing number of chambers, but remained positive over the range of conditions tested. There did not seem to be an obvious pattern to skewness with increasing urine coverage. Model (15) below produced a reasonable fit to the data, except for the log-normal EF distribution (Fig. 10.):

(15)$$\eqalign{&{\rm Skewness}\,(a) \sim {\rm EF}_{{\rm dist}} + 1/n + p/n + {\rm EF}_{{\rm dist}} \cr & \quad \quad \quad :1/n + {\rm EF}_{{\rm dist}} :p/n, \cr & \quad R_{{\rm adj}}^2 = 0 {\cdot} 8075}$$

Fig. 9. ESD Skewness as a function of (a) number of chambers and (b) urine coverage.

Fig. 10. (a) Modelled v. observed skewness and (b) residuals v. fitted values for model (15).

In model (15) only the log-normal EF distribution terms were significantly different from the constant EF distribution. The variability of the residuals increased with increasing skewness (Fig. 10.(b)). As the highest skewness occurred in the log-normal EF distribution, this distribution resulted in some large residuals and indicated that a better model of skewness may be required for log-normal EF distributions.

For the other EF distributions, model (15) indicates that for low values of p there will always be positive skewness regardless of the size of n. This threshold for p ranges from 0·168 to 0·178 depending on the EF distribution. This is higher than the maximum urine coverage considered in the present study (0·15), indicating that skewness cannot be eliminated for the range of urine coverages likely to be found in the field. Increasing p decreased skewness, but as n increases the skewness approaches a non-zero value (which depends on the EF distribution). By this model it is not possible to eliminate the skewness by increasing the sample size. While the CLT means the skewness should approach 0 for sufficiently high n, this must happen at n values well beyond those simulated in the present study, for which model (15) is valid.

Model (15) can be applied to constant, normal or uniform EF distributions, but is not recommended for log-normal EF-distributions. It also requires knowledge of the number of chambers and urine coverage. For 20 chambers and P = 0·05, a 20% increase in p resulted in a ~4% reduction in skewness, while a 20% increase in n resulted in a 7–8% decrease in skewness for the constant, normal and gradient EF distributions.

ALTERNATIVE ESTIMATORS

Mean

Fitting a linear regression of the alternative estimators on p gave the following relationships:

(16a)$$\eqalign{{\rm Mean}\,(g) &= 0 {\cdot} 984\,( \pm \,0 {\cdot} 0051) + 3 {\cdot} 02\,( \pm \,0 {\cdot} 054)p,\cr R_{{\rm adj}}^2 &= 0 {\cdot} 94} $$

(16b)$$\eqalign{&{\rm Mean}\,(g_{c1} ) = 0 {\cdot} 94\,( \pm \,0 {\cdot} 031) + 9 {\cdot} 0\,( \pm \,0 {\cdot} 33)p,\cr &\quad R_{{\rm adj}}^2 = 0 {\cdot} 80}$$

(16c)$${\rm Mean}\,(g_{c2} ) = 0 {\cdot} 95\,( \pm \,0 {\cdot} 016) + 8 {\cdot} 0\,( \pm \,0 {\cdot} 16)p,\quad R_{{\rm adj}}^2 = 0 {\cdot} 92$$

The geometric mean seriously underestimates the N₂O emissions, with the coefficient of p being less than one-third of what it should be. The first corrected geometric mean g _c1 gets much closer to the true mean, but not quite as close as the arithmetic mean (Eqn (10)). The second corrected mean (g _c2) also performs much better than the uncorrected g but still underestimates by more than g _c1.

For all estimators the variability increased as the mean increased (Fig. 11.).

Fig. 11. Mean of the ESD plotted against the true mean for (a) geometric mean (g), (b) first corrected geometric mean g _c1 and (c) second corrected geometric mean g _c2. The 1–1 line is also indicated.

rBias

Figure 12 shows rBias plotted against n and p for each of the three alternative estimators. The three estimators generally had negative bias, although g _c1 actually had a strong positive bias for n = 4. The bias was strongest in the uncorrected geometric mean. The rBias for g _c2 was always more negative than g _c1 but the difference approached 0 as the number of chambers increased. However, in terms of the magnitude of the bias, g _c2 was actually less biased than g _c1 for n = 4, but more biased at higher n. rBias decreased slightly (became more negative) then stabilized with increasing n for all three estimators. g showed a strong negative trend (i.e. became more biased) with increasing p, whereas the two bias-corrected geometric means had increasing variability in rBias as urine coverage increased. The range of rBias was from –13 to –46% for g, –26 to +29% for g _c1, and –26 to +1·3% for g _c2. These are large compared with the relative biases of a.

Fig. 12. rBias plotted against number of chambers (a, c, e) and urine coverage (b, d, f) for geometric mean (a, b), first corrected geometric mean (c, d), second corrected geometric mean (e, f).

Variance

The uncorrected geometric mean consistently had a lower variance than the arithmetic mean (variance ratio < 1, Fig. 13.). Given that low variance is a desirable property, the relative advantage (from a variance perspective) of g relative to a increases as n increases but decreases with increasing p. In Fig. 2, it was seen that for most cases Var(a) was inversely proportional n. Therefore the variance of g must decrease even more rapidly with n.

Fig. 13. Ratio of the variance of g, g _c1 and g _c2 to the variance of a plotted against the number of chambers (a, c, e) and urine coverage (b, d, f). Note log scale on the y-axis for the corrected geometric means.

The bias corrections for the geometric mean use the sample standard deviation of the log-transformed data to estimate the population log-standard deviation. Unfortunately, this increases the variance relative to the uncorrected geometric mean. For low sample sizes the corrected geometric means can sometimes have larger variances than the arithmetic mean. Both the corrected geometric means show similar trends in variance relative to Var(a) with increasing n and p. The second bias correction method always produced a lower variance than the first, although the difference decreased as the number of chambers increased.

For all three alternative estimators the relative improvement in variance over the arithmetic mean was greatest when the EF distribution was log-normal.

Skewness

None of the estimators consistently reduced (or increased) the skew relative to the arithmetic mean (Fig. 14.). In general, for log-normal EF the geometric mean (and the corrected versions) tended to have lower skew than the arithmetic mean. For the other EF distributions the skew tended to be higher than the arithmetic mean.

Fig. 14. Ratio of the skewness of g, g _c1 and g _c2 to the skewness of a plotted against number of chambers (a, c, e) and urine coverage (b, d, f). Note log scale on the y-axis for the corrected geometric means.

DISCUSSION

The arithmetic mean of chamber measurements was an unbiased estimator of the N₂O emissions from a grazed pasture. However, the CLT did not apply over the range of conditions examined (even with 100 chambers) and the sample mean distribution was right-skewed. The variance was relatively high (with standard deviation ranging from 12 to 135% of the actual mean) but did appear to be proportional to 1/n except when the EF distribution was log-normal. For n < 100 the variance was proportional to the mean, therefore attempts to use regression modelling on chamber data should include weights to account for the non-constant variance. The lack of normality means that statistical tests that assume normality of the sample mean (e.g. t-test and F-test) should not be used for this type of data. With knowledge of the urine coverage and EF distribution, equivalent tests could be created by generating empirical sampling distributions as we have done here. However these tests would probably have lower power due to the high probability of not sampling any urine patches.

Formulae were derived that gave reasonable estimates for variance, pUnder and skewness of the sample distribution of the arithmetic mean. These formulae required knowledge of the number of chambers, the urine coverage and the EF distribution. Whether the urine concentration was constant or normally distributed had little impact on the sample mean distribution. An experimenter will know the number of chambers being used, but will not know the urine coverage or the EF distribution. However, there are methods for estimating the urine coverage if the number of animals grazed and the duration of the grazing event are known (e.g. Pleasants et al. 2007) and there have recently been some studies looking at measuring urine patch concentration experimentally (Moir et al. 2011). The EF distribution is trickier. In particular, the formulae for variance and skewness are poor predictors when the EF distribution is log-normal. Therefore determining the nature of the EF distribution should be a priority.

Some alternative estimators were also examined. Log-transformations are commonly used to make right-skewed data more symmetrical. The geometric mean is the mathematical equivalent of back-transforming the mean of the log-transformed data. Unfortunately, the geometric mean was a biased estimator that consistently underestimated the true mean, and although the geometric mean did have a lower variance than the arithmetic mean, it did not necessarily have less skew. Still, the lower variance of the geometric mean could make it a useful estimator for situations where the aim is to compare N₂O emissions from different treatments, rather than to estimate the actual emissions.

Two methods for bias-correcting the geometric mean were also investigated. These methods both used the sample standard deviation of the log-transformed data to calculate the bias correction. While these methods did reduce the bias they also had the unfortunate side effect of increasing the variance. None of the alternative estimators universally outperformed the arithmetic mean. However, in some circumstances (e.g. when the EF distribution is log-normal) a bias-corrected geometric mean might make a better estimator, but this would need to be determined on a case-by-case basis.

A number of simplifying assumptions were made for the current simulations, which included that the urine patches did not overlap and that the EF was independent of the urine concentration. The first assumption is more likely to be satisfied when stocking rates are low and grazing durations short. According to the method of Pleasants et al. (2007) with the mean patch area (0·37 m²) and urination rate (0·53 events/cow/hour) measured by Moir et al. (2011), one day grazing at a stocking density of 180 cows/ha would result in 0·081 of the area being covered by urine patches, but 0·96 of the urine patch area would be single urine patches. So overlapping urine patches would account for only 0·003 of the total area. However, this assumes that the distribution of urine patches is completely random. In reality, animal behaviour is likely to result in non-random urine patch distribution. For example, urine patches are likely to be more frequent in areas where animals congregate (e.g. near water troughs or shelter), but there could also be a tendency to avoid grazing near existing urine patches. The measurements of Moir et al. (2011) found urine patch coverage in a season ranged from 0·035 to 0·076 with cow-days/ha/season ranging from 120 to 450. This lower coverage indicates a higher rate of overlap than estimated using the Pleasants method.

The second assumption requires that the relationship between N₂O emissions and urine N concentration be linear. Kim et al. (2013) hypothesized that the relationship between N₂O emissions and applied N could have linear and non-linear phases depending on the limiting process. Kelliher et al. (2014) found a strong linear relationship between N₂O emissions and applied N, but no significant quadratic term for urea applied at rates of 25 to 1500 kg N/ha. However, in a field study, de Klein et al. (2014) found an increasing trend of EF with urine N concentration in a free-draining soil, but not on a poorly draining soil. In addition to a direct relationship between urine concentration and EF, there could also be an indirect effect if the spatial variability of EF was driven by some factor (such as soil moisture or slope) that also affects the amount of time the animals spend in a given part of the field.

The present study also did not consider the effect of using different means and standard deviations for the EF and U concentrations. This is an area that could be further investigated if more data regarding these distributions were available. However, completely random sampling on a grazed pasture is always going to be problematic due to the large difference between the emissions from urine and non-urine patches and the lack of knowledge of the number of urine patches included in a given sample. If non-destructive methods for determining the urine patch locations shortly after a grazing event can be developed then a stratified sampling scheme could be used to give more robust estimates of emissions. Alternatively, if the presence or absence of a urine patch in each chamber could be conclusively determined after the fact (e.g. from soil properties or grass growth), then the relative contribution from the urine and non-urine containing chambers could be scaled according to the overall urine coverage. The method of applying a known amount of urine-N to each chamber avoids the errors introduced by not knowing the location of the urine patches, but requires the total amount of urine-N excreted during grazing to be estimated, which may miss effects caused by animal behaviour (e.g. treading).

CONCLUSIONS

In the present study formulae were derived that can be used to calculate the mean, bias, variance and skewness of the sample distribution when measuring N₂O emissions using chambers on a grazed pasture. These formulae require knowledge of the fractional area covered by urine patches, number of chambers, and the nature of the EF distribution.

The arithmetic mean was an unbiased estimator, and over a large number of trials the mean of all the samples will approach the true value. Unfortunately, the sample mean distribution had a high variance and skew. This meant for any particular trial the sample mean could be a long way from the true value and there was a >50% chance that any given sample would be an underestimate. The CLT did not apply for realistic numbers of chambers, therefore many common statistical techniques requiring normal data (e.g. t-tests and F-tests) should not be applied to data of this type. The behaviour of the sample mean distribution did not depend on whether the urine concentration was constant or normally distributed, but did depend on the distribution of the EF.

The geometric mean had lower variance, but high bias compared with the arithmetic mean. Two bias corrections for the geometric mean were tried. These methods did reduce the bias of the geometric mean, but as the corrections used estimates of the standard deviation, the variances were high. The corrected geometric means were not a universal improvement on the arithmetic mean, although in some circumstances they had lower variances and skewness (e.g. when the EF was log-normally distributed and n > 4). However, all these alternative estimators still suffered from non-normality.

The best way to improve the estimation of N₂O emissions from a grazed pasture would be to sample separately from known urine and non-urine patch areas. This, however, requires the development of a non-disruptive method for identifying urine patches shortly after a grazing event. In the absence of such a method there will always be some inherent problems (e.g. skewness) in the sample mean distribution. The distribution of the EF was also an important factor in the properties of the sample mean distribution. Knowledge of how the EF varies spatially is therefore useful for determining the properties of the sample mean distribution.