Geostatistical tools for assessing sampling designs applied to a Portuguese bottom trawl survey field experience

We present a bottom trawl survey (bts) field experience carried out off the Portuguese Continental shelf to test two sampling designs proposals previously analysed by simulation which implement a hybrid random-systematic and a systematic sampling strategy. We used a common base regular grid covering the survey area and overlapped it with the existent random design to build the hybrid design while the systematic design added a set of regular locations at smaller distances creating four denser sampling areas. We use hake (Merluccius merluccius) abundance and model-based geostatistics to compute measures such as mean abundance, μ, and the 95% percentile, p95, which summarise the areal behaviour; coverage of the prediction confidence interval, ξ, to assess the adequacy of the model; and a modified generalised cross validation index, p, to evaluate prediction precision. the hybrid design showed a lower coefficient of variation for μ (11.89% against 13.25%); a slightly higher coefficient of variation for p95 (11.31% against 11.09%); similar ξ (0.94); and lower p (16.32 against 18.82). We conclude that the hybrid design performs better, our procedure for building it can be used to adjust bts designs to modern geostatistical techniques, and the statistics used constitute valuable tools for assessing bts performance.

introduCtion designs for bottom trawl surveys (bts) rely on previous knowledge of the spatial distribution and population structure of the target species combined with statistical analysis of preliminary data (e.g. ault et al., 1999;Hata and berkson, 2004) or simulation procedures (e.g. schnute and Haigh, 2003;anon., 2005b). in defining the bts sampling design, one is confronted with operational constraints such as the existence of trawlable grounds and vessel time available.the survey design is typically reviewed from time to time to adjust the stratification (e.g.smith and Gavaris, 1993;Folmer and Pennington, 2000), tow duration (e.g.Cerviño and saborido-rey, 2006;Wieland and storr-Paulsen, 2006), technical issues such as gear changes (e.g.Zimmermann et al., 2003;Cooper et al., 2004) and other factors which may change over the years.
several authors have discussed the advantages of systematic designs over random designs for sampling spatial correlated variables such as fish abundance (Cochran, 1960;ripley, 1981;thompson, 1992;Cressie, 1993;Chiles and delfiner, 1999;Kimura and somerton, 2006;diggle and ribeiro Jr, 2007).nevertheless, in the case of spatial correlated variables there are two conflicting objectives that cannot be combined in a single criterion: estimation of the covariance function parameters and prediction (müller, 2001). in the first situation it is important to have locations at short distances to inspect the behaviour of the correlation function close to the origin, and locations at distances close to the limit of spatial correlation to estimate the correlation range (müller, 2001). in the second situation the best predictions will result from the design with the highest covariance with the locations to be predicted (thompson, 1992). in the case of predicting fish abundance it is common to require a complete map of the study area and the best choice will be a design that covers the area evenly.However, when the covariance function is unknown (a common characteristic of fish abundance analysis), it must be estimated from the data before predicting and the two objectives must be combined.several authors suggest designs that mix a set of locations covering the area with additional locations at short distances (müller, 2001;diggle and lophaven, 2006;Zhu and stein, 2006) to balance the two objectives.such designs applied to bottom trawl surveys have received limited attention (selzenmuller et al., 2005), although fish abundance characteristics fit well in the assumptions of these proposals.
our analysis adopts a model-based geostatistical method (diggle et al., 1998;diggle and ribeiro Jr, 2007) to explicitly take into account spatial patterns of abundance and provide a flexible modelling framework.it is worth mention previous work by simard et al. (1992) and Petitgas (2001) wich also addressed the issue of sampling design, although not using model-based geostatistcs.
the designs are accessed by a set of statistics to provide information about aspects of the data that are relevant for modelling fish abundance.in a global perspective, referring to the entire study region, we use mean abundance and the 95% percentile to summarise the areal abundance, commonly used for studying time trends and building abundance indices for stock assessment.in a local perspective, referring to particular locations within the study area, we use the observed values to assess the suitability of the model by computing the coverage of the prediction confidence interval, and the prediction precision by computing a modified generalised cross-validation index.note that the assessment of the model suitability and the prediction precision are extremely valuable statistics, once that kriging is a linear predictor and the maps produced with it will be used to estimate the spatial distribution of abundance and the abundance index mentioned above.With relation to the analysis reported here we rely on our experience with bottom trawl surveys (anon., 2002, 2003, 2004, 2005a, 2006; sousa et al., 2005; mendes et al., 2007; sousa et al., 2007) to provide contextual information which supports the adoption of a particular class of model, and avoid model mis-specification as far as possible.
the work described on this paper aims to: (i) report a bts field experience for testing sampling designs, and (ii) describe geostatistical tools for assessing the performance of sampling designs.matErial the Portuguese bts started in June 1979, covering the continental shelf and following a stratified random design.in 1989 the stratification was defined by 12 sectors along the coast subdivided into 4 depth ranges: 20-100 m, 101-200 m, 201-500 m and 501-750 m, with a total of 48 strata.due to constraints in the vessel time available the sample size was set to 97 locations evenly allocated to each stratum.the coordinates of the sampling locations were selected randomly, albeit constrained by the historical records of clear tow positions and other information about the sea floor, thus avoiding places where trawling was not possible.during this period haul duration was set to one hour but recent experiments proved that half-hour hauls provide the same information about length distributions (Cardador, pers. comm.). in light of this findings haul duration was reduced to half an hour and an additional set of hauls were available which led to a revision of the sampling design.the revision was split into a preliminary phase using simulations and geostatistical analysis (Jardim and ribeiro Jr, 2007) and a second phase during which a field test was performed to provide real information about the proposed sampling designs. in the third phases the decision will be based on the scientific data provided and the existing financial and administrative constraints.
the field experience was carried out during the summer of 2001, with r/V noruega off the southwest of the Portuguese Continental shelf (Fig. 1) using a norwegian Campbell trawl 1800/96 (nCt) with a codend of 20 mm, a mean vertical opening of 4.8 m and a mean horizontal opening between wings of 15.5 m. the survey performed two sampling designs selected from the simulation study reported by Jardim and ribeiro Jr (2007).the survey area was limited in the south by the cape of s.Vicente (37.00º north), in the north by setubal's Canyon (38.30º north), in the east by the 20 m depth isoline and in the west by the 500 m isoline.the survey area was approximately 4300 km 2 in size and the maximum distance within the area was approximately 150 km. the data collected in both designs and considered here consisted of date/time, geographical location and hake (Merluccius merluccius) catch in weight (kg).Geographical coordinates were transformed into utm units and hake abundance was computed in kg∕km and assigned to the haul starting coordinates.the area swept was computed using the haul start and ending positions to correct haul speed variations.mEtHods this section describes the sampling designs to be tested and how they were built.it also describes the geostatistical modelling framework and the ad-justments considered to cope with the small dataset available-a common characteristic of bts due to its high price.Finally, we describe the technical details of the performance statistics chosen.

Sampling designs
our sampling designs were built by mixing a set of operational constraints with the geostatistical principles described above and the need to keep the continuity of the survey history. in particular, the two designs tested were built to distinguish between a hybrid random-systematic sampling strategy and a systematic strategy.
the sampling effort available for the candidate design was 36 locations.We built two candidate de- signs using as a basis a regular grid with 19 locations covering the survey area, and added seventeen additional ones at shorter distances.two candidate designs were built, a hybrid design that allocated the additional locations randomly and a systematic design that allocated them at regular locations.the hybrid design overlaps the regular grid with the existent random design, keeping some continuity with the survey historical records (top-left plot in Fig. 2). the systematic design includes regular locations at smaller distances, creating 4 denser sampling areas (bottom-left plot in Fig. 2). the vessel time available did not allow other possibilities to be tested.

Geostatistical model
Geostatistical observations consist of pairs (x,y) with elements (x i ,y i ) : i = 1,…,n, where x i denotes the coordinates of each of the n spatial locations within a study region A ⊂ r 2 and y i the measurement of the corresponding observable study variable.We adopted the box-Cox transformed Gaussian model with transformation parameter λ as presented in Christensen et al. (2001).denoting by z i the transformed values, such that g λ (y i ) = z i , the model for the vector of variables Z observed at locations x can be written as a linear model Z(x) = S(x) + ε, where S is a stationary Gaussian stochastic process, with E[S(x)] = µ, Var[S(x)] = σ 2 and an isotropic correlation function ρ(h) = Corr[S(x),S(x′)], where h = ||x-x′|| is the Euclidean distance between locations x and x′. the terms ε are assumed to be mutually independent and identically distributed, ε ∼ Gau(0,τ 2 ).For the correlation function ρ(h) we adopt the exponential function with algebraic form ρ(h) = exp{-h∕φ} where φ is the range parameter such that ρ(h) ∼ 0.05 when h = 3φ.Following usual geostatistical terminology (isaaks and srivastava, 1989) we call σ τ 2 = τ 2 + σ 2 the total sill, σ 2 the partial sill, τ 2 the nugget effect and 3φ the practical range.Geometric anisotropy (isaaks and srivastava, 1989;Cressie, 1993) is considered an extension of this model with an extra parameter ψ = {ψ A ,ψ R } , where ψ A is the anisotropic angle and ψ R is the anisotropic ratio.
Hereafter we use [•] to denote the distribution of the quantity indicated within brackets.Following the adopted model, [g λ (Y)] ∼ mVGau(µ1,Σ), i.e. [Y ] is multivariate trans-Gaussian with expected value µ and covariance matrix Σ parametrised by {σ 2 ,φ,τ 2 }.Parameter estimates can be obtained by maximum likelihood (Cressie, 1993;diggle et al., 1998;diggle and ribeiro Jr, 2007) and used for spatial prediction. in its simplest format, spatial prediction given by the kriging predictor consists in obtaining expected values and associated variances at unsampled locations.more generally, the predictive distribution of quantities of interest can be obtained analytically, if possible, or by sampling from this distribution.Considering a prediction target T(x 0 ) = g λ -1 (S(x 0 )), the realised value of the process in the original measurement scale at spatial locations x 0 .simulations from the conditional distribution [T(x 0 )|Y (x)] are obtained by simulating from the multivariate Gaussian [S(x 0 )|Y(x)] and back transforming the simulated values to the original scale of measurement (Chiles and delfiner, 1999;diggle and ribeiro Jr, 2007).these simulations are called conditional simulations referring to the fact they are obtained from the distribution of the quantity of interest conditioned to the observed values Y (x).
We split inference into two steps.First the box-Cox transformation parameter λ and the anisotropy parameter ψ R are investigated by pooling all the observations in a single dataset and computing profile likelihoods (diggle and ribeiro Jr, 2007).We consider the north-south coastal orientation of the study region as the direction of greatest spatial continuity and fix ψ A in 0 degrees azimuthal angle.afterward, having estimated these two parameters we regard their point estimates as constants in the model and proceed by computing, for each design, the maximum likelihood estimates of the remaining model parameters.the reasoning for the two-step procedures is twofold.Pragmatically, it overcomes the difficulty of identifying all parameters with a small dataset.With regards to modelling assumptions it considers the transformation and anisotropy parameters as part of the model specification reflecting the nature of the data and contextual information, therefore not requiring to be identified by the designs.thereafter, we compute kriging predictions on a 2 × 2 km grid within the study area, x 0 , with a total of 1070 locations, and obtain 1,000 conditional simulations from [Y (x 0 )|Y ] for each design.

Performance statistics
Consider E[Z(x i )] and σ z 2 (x i ) the kriging predictor and its variance on the Gaussian scale at location x i ∈ x 0 .For a value of 0.5 of the transformation parameter λ, the back transformation to the original scale gives . the variance of µ ˆ, denoted by σ ˆµ2 , is computed by the mean of all terms in the covariance matrix , where Σ Z (x 0 ) is the covariance matrix of [S(x 0 )|Z(x)].more generally, inferences on other quantities of interest T(x 0 ) are obtained from the conditional simulations.denote by t s (x 0 ), s = 1,…,S = 1,000 conditional simulations from [T(x 0 )|Y (x)].For example, an α-th percentile is estimated by p ˆ = S -1 Σ s p ˆs where p ˆs = p α (t s (x 0 )), the average of the empirical distribution p ˆ obtained from the conditional simulations.the variance of p ˆ is given by σ ˆp2 = (S-1) -1 Σ s (p ˆs -p ˆ)2 .
the coverage of the prediction confidence interval, ξ, and the generalised cross validation index, p, were computed using cross-validation statistics (Hastie et al., 2001) combined with conditional simulations as follows.First, create a new data set by leaving one observation out at a location x i , simulate 1,000 values of the variable at that location, and repeat this procedure visiting all data locations.subsequently, consider y(x i ) an observation of the process Y at location x i , i = 1,…,n; y(x (i) ) the observed data set without the observation y(x i ) and t s (x i ) a conditional simulation s = 1,…,S of [T(x i )|Y = y(x (i) )] at location x i .the predictive confidence interval is given by CI(x i ) = [p 2.5 (t s (x i )),p 97.5 (t s (x i ))] and the proportion of observations lying inside the intervals ξ = n -1 Σ i (y(x i ) ∈ CI(x i )) provides the coverage of the prediction confidence interval.the cross validation index is given by p = n -1 Σ i (S -1 Σ s (t s (x i )y(x i )) 2 ), the average of the mean quadratic error at each location estimated using the full set of conditional simulations.rEsults the two sampling designs and the observations of hake abundance are presented in the leftmost panels of Figure 2, where the base regular design is represented by the black triangles.the abundance of hake observed showed that the distribution of abundance was spread over the area, with lower values in the north and a small number of zeros.
the 95% confidence interval obtained for the box-Cox transformation parameter was [0.12,0.55]and we set λ ˆ = 0.5, corresponding to a square root transformation.the profiled log-likelihood of the anisotropy ratio showed no evidence of anisotropy.nevertheless, we carried out analysis using different values of ψ R to check the sensitivity of the results, which proved negligible.
Covariance parameter estimates showed higher values for the hybrid design than the corresponding ones given by the systematic design (table 1). the total variance σ ˆT2 was 3.75, with σ ˆ2 = 0.75 and τ ˆ2 = 3.00; and φ ˆ = 16.64 .the systematic design estimates were σ ˆT2 = 3.20, with σ ˆ2 = 0.61 and τ ˆ2 = 2.59; and φ ˆ = 10.21.Considering τ ˆ2REL , which computes the relationship between the random variability and the spatial structure variability, and σ 2 φ -1 , which give information about the "size" of the spatial process.the two designs showed similar relative nuggets.However, the hybrid design showed a lower ratio between sill and range, reflecting a higher spatial structure of the stochastic process.notice that the practical range, 3φ, was ≈ 50 km for the hybrid and ≈ 30 km for the systematic design.
the rightmost panels of Figure 2 show the abundance maps predicted and their variance for each design.both predictions are similar and the spatial pattern of variance reflects the influence of the observations, showing lower variability near the observed locations and higher variability in areas where extrapolation was further extended.the hybrid design had higher variance in the centre-east of the study area and lower variance in the north due to a better coverage in this area.
the estimates of µ and p 95 were similar, although the hybrid design showed slightly lower values.the hybrid design showed a lower coefficient of variation for µ, CV µ = 11.89%,than the systematic design, CV µ = 13.25%.sampling statistics computed for these designs showed a similar pattern (table 1). the p 95 variance was slightly lower for the systematic design, CV p95 = 11.09%,than for the hybrid design, CV p95 = 11.31%.the coverage of the prediction confidence intervals was 0.94 for both designs.these results reinforce our modelling choices given that if the model was wrong we would expect ξ to be different from the nominal value of the confidence interval.the generalised cross validation index showed a lower estimate with the hybrid design, 16.32, than with the systematic design, 18.82, indicating a higher prediction precision of the hybrid design.the above results show that the higher spatial structure of the stochastic process estimated for the hybrid design exceeded its higher total variability with relation to the estimation of these performance statistics.disCussion assessing sampling designs for bts raises interesting questions about appropriate methodologies for analysing data and obtaining statistics of interest, which are particularly relevant considering the multipurpose/multispecies nature of bts and the small sample sizes.the adoption of a formal criteria and loss function to find an optimum design seems unrealistic in practice due to the multidimensionality of the data and the conflicting objectives of inference and prediction.Here we followed a pragmatic approach to sampling design, starting with a design that joins a regular grid with the old random design, following with a design that uses the same regular grid but reallocates the random locations in a regular shape.Considering the wide literature that supports the use of systematic designs for spatial correlated variables, these designs implement the two most promising strategies and test the possibility of keeping the continuity with the historical time series.to compare these proposals we rely on spatial modelling to compute statistics of primary interest and look for consistency among them, exploring several aspects of the same dataset.We advocate that the approach described above will provide valuable information to support the decision-making process.
although the results obtained are constrained by the characteristics of the area and the species analysed, we believe that the methodology defined by our approach can be applied to other areas and species, providing an important source of information for revising sampling design.it would not be surprising if similar results are found for other species, because the principles behind the construction of the sampling designs tested are quite generic and can be applied to most fish species.
the performance statistics were selected to reflect relevant characteristics and different aspects of spatial prediction.the global mean is the most widely used index of abundance, often estimated by the sample average.We favour the geostatistical estimator presented and its variance as a measure of uncertainty, because it takes into account the spatial dependency within the area and insights about the spatial process.the 95th percentile estimated by conditional simulations can be used to identify areas of high abundance, giving information about candidate areas to protect. the coverage of the prediction confidence intervals is a diagnostic tool.a small coverage reflects an underestimation of the variance or the inadequacy of the model to explain the available data. the cross-validation index combined with conditional simulations incorporates the prediction precision in the index, which is not taken into account by the traditional cross-validation.For example, if a location has the same value predicted by different designs but with different prediction variances, this index will distinguish the two situations.
our results showed that the hybrid design performed better in all cases except for σ p 2 .a clear parallel can be established with the lattice plus closed pairs designs of diggle and lophaven ( 2006), the EK-optimal designs of Zimmerman (2006) and the D EA designs of Zhu and stein (2006).although following different construction, these designs cover the entire study areas in study, include a set of positions at small distance and performed better than their random or systematic competitors.Common to all these studies and our work is the fact that the analyses were carried out in situations in which the model parameters were considered unknown and needed to be estimated from the data, which made it clear that both parameter estimation and prediction are important for the precision of the prediction target.
in conclusion, we consider that our results indicate that keeping the old random design and adding  -, the variance of the sampling mean.model parameters are: τ 2 , the short distance variance or nugget effect; σ 2 the variance of the spatial process; σ 2 T the total variance; φ the correlation range parameter; and the transformation parameters λ, the box-Cox parameter and the anisotropy parameters {ψ A ,ψ R }. the relative nugget, τ 2 REL , and the ratio between relative sill and range σ 2 φ -1 , were computed to give more insights about the spatial process.Performance statistics are: µ ˆ and σ ˆµ2 , the mean and variance of the global abundance; p ˆ95 and σ ˆp2 , the mean and variance of the 95th percentile of the global abundance; p, the generalised cross validation index and ξ, the coverage of the prediction confidence interval with nominal level of 0.95.

Fig. 1 .
Fig. 1. -survey area on the southwest of the Portuguese Continental shelf between 20 m and 500 m.

Fig. 2 .
Fig. 2. -study area on the Portuguese southwest coast.the top panels show information about the hybrid random-systematic design and the bottom panels about the systematic design.the leftmost plots show the sampling designs locations, the black triangles represent the regular grid common to both designs, and the open circles the additional locations.Follows the observations of hake abundance (kg/km 2 ) and the predictions obtained by kriging, both on the square root scale.the rightmost plots present the kriging variance

Table 1 .
-sampling statistics, estimates of model parameters and performance statistics by design.sampling statistics are: n, the sample size; Y -, the sampling mean; s 2 regular grid to build a new design can be a good and pragmatic solution for adjusting bts designs to modern geostatistical techniques.secondly, the performance statistics described above seem to capture the most important features of the data with relation to abundance estimation, constituting good measures for assessing bts performance.