ICES Journal of Marine Science, 54: 213–225. 1997
Sole egg distributions in space and time characterised by a geostatistical model and its estimation variance Pierre Petitgas Petitgas, P. 1997. Sole egg distributions in space and time characterised by a geostatistical model and its estimation variance. – ICES Journal of Marine Science, 54: 213–225. Geostatistics is used to model the spatio-temporal variability in egg volumetric density measured by monitoring surveys performed on a major spawning ground of sole in the Bay of Biscay (France). An approximation formula is proposed for computing the estimation variance of the mean egg density estimated over space and time. The major structural feature is the time invariance in the spatial distribution of the relative density (density(x,t)/mean(t)). This leads to interpretation of the different samplings as different realisations of the same random function in space. Because of the timerepeated samplings at each station, drift (or trend) and residuals can be separated. The drift is defined as the spatial structural component that is constant over time. The spatio-temporal variability was described by a multiplicative model with two effects: a time effect corresponding to a proportionality to the spatial average at each time and a space effect corresponding to the time-invariant spatial distribution of the relative density. Residuals showed structure in space but not in time and no space–time interaction. Geostatistical aggregation curves were used to characterise invariance over time of the aggregation pattern, the drift was mapped using non-stationary IRF-k (Intrinsic Random Functions of order k) kriging and variograms and crossvariograms were used to study the structure in the residuals. ? 1997 International Council for the Exploration of the Sea
Key words: spatial distribution, geostatistics, survey precision, egg production. Received 25 October 1995; accepted 12 June 1996. P. Petitgas: ORSTOM, HEA, BP 5045, 34032 cedex 1, Montpellier, France.
Introduction Series of monitoring surveys provide a certain amount of information on the spatio-temporal variability of fish distributions. Here geostatistics is used to characterise this variability. Variability can be interpreted in terms of pure randomness, correlation or drift, depending on the degree of determinism. Repeatable spatial patterns over the years have been observed in groundfish distributions (Hilborn and Walters, 1992; Houghton, 1987; Myers and Stokes, 1989) and this consistency over time will be considered here as a deterministic drift (or trend). In geostatistics the sampled spatial distribution is interpreted as one realisation of an underlying process. But for one data set (one realisation) Matheron (1971) showed that the separation between drift and residuals is not possible because the drift estimation procedure affects the residual estimation. The residuals then denote what is left after the method of drift estimation has been applied and they cannot quantify the natural variability around the spatial drift component. To overcome this 1054–3139/97/020213+13 $25.00/0/jm960184
difficulty, Sullivan (1991) used ancillary information. In an example where fish density is related to depth, the author used the regression of fish density on depth to estimate the drift in fish density from the known bathymetry in the area. In this paper each survey from a time series of surveys is interpreted as a realisation of the same underlying process and fitted to an overall spatio-temporal model. Casting each of the surveys as repeated observations enables estimation of the drift and separation of this drift from the residuals. The drift is defined here as that component of the spatial distribution which is invariant over time and does not necessarily need to be a smooth function over space. This is the main difference between geostatistical and other mapping procedures (Foote and Stefansson, 1993). Then, because residuals in geostatistics characterise a component of process variability, the inference of their spatio-temporal structure enables one to derive the estimation variance for the spatio-temporal mean estimate. Such an approach is particularly interesting for repeated ichtyoplankton surveys performed regularly ? 1997 International Council for the Exploration of the Sea
214
P. Petitgas
over the spawning area during the spawning period to assess egg and larval production. In the case of systematic design in space and time, a spatio-temporal model characterising the entire series of surveys is needed to estimate the survey precision of the egg or larval mean density. This paper analyses the planktonic sole egg distributions in space and time on a major spawning ground of sole (Solea solea, L.) in the Bay of Biscay (France). First, a consistent component in the spatial distribution is characterised and modelled as a drift in space. It is estimated by averaging over time using the different surveys. Then residuals are shown to exhibit a proportional effect. Finally a multiplicative model is proposed where time and space are the two factors. An estimation variance formula is derived for the spatiotemporal egg estimation. To my knowledge, this is the first attempt in the fisheries literature to evaluate spatiotemporal estimation variance for ichtyoplankton monitoring surveys.
Table 1. Monitoring surveys on the spawning ground of sole off the Loire river (France). The date is the mid day of the survey. Data are egg densities: number of eggs (stage-I) per 10 m3.
Material and methods
average stage-I egg density at time t on the surveyed area V:
Sampling design and data A major spawning ground of sole (Solea solea, L.) in the Bay of Biscay is located off the mouth of the river Loire (Arbault et al., 1984) where artisanal trawlers operate on sandy sediments (Dardignac, 1988). From 1986 to 1989, Institut Français de Recherche pour l’Exploitation de la Mer (IFREMER) undertook a monitoring programme on this spawning ground to measure the spatio-temporal variability of the ichtyoplankton and the hydrodynamics in the general context of studies in sole recruitment fluctuations (Koutsikopoulos, 1991). Over the 4 years of the programme, 13 surveys were performed from the end of February to the beginning of April, 7 of which were in the peak of the spawning period (Table 1). All surveys sampled the same grid of 35 stations (Fig. 1). The same 11 stations were replicated at each survey. The eggs were sampled by vertical hauls, in which the net was hauled on a V trajectory, from surface to bottom to surface again. The average filtered volume per haul was 544 m3. Details on the net, the hauls and the sole egg vertical structure are given in Koutsikopoulos (1991). In the laboratory, sole eggs were identified and counted for each of the development stages of Riley (1974). The counts were made on the whole samples (no subsampling). The variable studied is the number of sole eggs at stage I per 10 m3 of sea water.
Modelling Let us denote Z(x,t) to be the stage-I egg density at point x and time t. In addition, ZV(t) denotes the
Year
Date
Spatial mean ZV(t)
1986 1986 1986 1986 1986 1987 1987 1987 1988 1988 1989 1989 1989
March 2 March 6 March 10 March 18 April 7 March 8 April 8 April 14 March 4 March 30 February 23 March 14 April 8
2.22 2.81 2.83 1.40 0.41 3.12 0.22 0.19 2.64 0.12 0.47 2.53 0.14
Maximum density 10.6 21.6 23.5 6.4 2.7 28.3 1.4 0.9 17.6 0.7 2.8 16.8 0.7
The relative density at x and t is the contribution of station x to ZV(t) and is denoted Zp(x,t):
The experimental relative density showed constant statistical characteristics over time. Thus, the relative density was interpreted with a Universal kriging approach (Matheron, 1971) which assumes that: (i) the drift is the expectation of the random function over many realisations, (ii) the scale of the drift is large in comparison to the surveyed area and (iii) the drift is separated from the residuals and the residuals are then analysed. The drift was defined by:
where the expectation is taken over time at point x. M(x) is the (not necessarily smooth) drift that is constant over time and R(x,t) are spatio-temporal residuals. Residuals R(x,t) showed a proportional effect to M(x) which was modelled as: R(x,t)=aM(x)U(x,t) where residuals U(x,t) have zero mean and variance one. Thus, the complete model used for characterising the data was:
A geostatistical model in space and time
215 47°30'
30m
34 Belle Ile 32
26
LO
28 24
30
IR E
100m
18 20
47° 16
14
22
8
10
Yeu
6
12
2 4 46°30' N
2°30'
2°W
Figure 1. The monitoring survey grid on the spawning ground of sole off the mouth of the river Loire (France). The stations are numbered. The replicated stations are figured by double crosses.
between the two replicates was used. The residuals R were estimated by: which was written as: Z(x,t)=ZV(t)M(x)Y(x,t) with Y(x,t)=1+aU(x,t). The residuals U(x,t) exhibited a correlation structure in space, pure randomness in time and no space–time interaction. The time process ZV(t) was assumed to be purely random and stationary during the period studied. For this model, we propose an approximation for the estimation variance of the spatial temporal mean estimate in the case of regularly spaced surveys in time. Drift and residuals were estimated on the data directly because all of the surveys used the same survey grid. For each survey, the zonal mean, ZV(t), was estimated by the simple average of the survey data. At each station the drift M(x) was estimated by a simple average in time:
where nt denotes the number of surveys and nx the number of stations. At the replicated stations, the mean
and at the replicated station by:
Residuals have a zero experimental mean in time by definition. They also have zero experimental mean in space because the experimental mean of the relative density equals one. Residuals U were estimated by:
The model was inferred conditional on the values of ZV(t). Estimated values were considered as true process values for all the process components of the model, ZV(t), Zp(x,t), M(x), U(x,t). Petitgas (1991) studied the influence of estimation errors on the time average variogram of the residuals U and showed that it was small.
216
P. Petitgas
Geostatistical tools Spatial distributions were compared between years by computing aggregation curves (‘‘selectivity’’ curves of Matheron, 1981). Structure in the residuals U was investigated by computing variograms and crossvariograms (Matheron, 1971). The estimation variance formula of the spatio-temporal mean estimate was based on the model used. For estimating this variance we needed the variogram model of the residuals and the map of the drift values. We interpolated the drift by kriging with an intrinsic function of order 2 (Matheron and Delfiner, 1980). Aggregation curves The curve Q(T) of Matheron (1981) (geostatistical ‘‘selectivity’’ curve) relates the abundance Q(z) to the surface T(z) occupied by densities greater than the threshold z, when the z values are ranked from maximum to minimum. In fisheries science the use of the term ‘‘selectivity’’ is confusing as it has another meaning. I suggest the use of the term ‘‘aggregation’’ curve because Q(z) measures the maximum abundance that can be on any surface T(z) in the survey area. The curve Q(T) is similar to a Lorenz curve (Dagum, 1985). The difference is that in the aggregation curve, densities z are ranked from maximum to minimum. The curve Q(T) is clearly related to the cumulative density distribution but because abundance is plotted against surface, the aggregation curve is a convenient tool for performing comparisons between years. T ranges from 0 to V and Q ranges from 0 to the total VZV. Here, V is considered constant and the surfaces T(z) are expressed as percentages of V. As such, T will range from 0 to 1 and Q from 0 to ZV. We computed aggregation curves for the relative density using the ratio P(z)=Q(z)/ZV, for each survey. P(z) measures the maximum percentage of the total abundance on any surface T(z). T and P range from 0 to 1. We computed the curve P(T) for each year as follows:
survey design (for a recent review of applications see Pelletier and Parma, 1994). The variogram quantifies mean square differences in value between pairs of sample points as a function of the distance, h, separating them. If the variogram curve increases with distance, then there is correlation structure as closer points in space are more alike than distant ones. When the variogram curve stays flat, there is no correlation structure. In general, the variogram curve will increase with distance then stabilise around a sill. The distance at which the sill is reached is the range and measures the range of the correlations. Different estimators of the variogram have been proposed (Cressie, 1991) and the classical experimental variogram was used, computed on the residuals U for each survey at time t as follows:
where D(h)={(i,j)Qxi "xjP=h} contains n(h) pairs and i and j are two sample points. The replicates were used to estimate the variogram for very small distances and a special distance class was considered. Cross-variograms The cross-variogram is a multivariate stationary structural tool characterising how two variables covary in space and is used in the cokriging equations (Matheron, 1971). Cross-variograms were computed between residuals of two surveys. Cross-variograms were computed for all possible pairs of surveys then averaged for each distance. We analysed how the differences between points x and x+h covaried at times t1 and t2 and cross-variograms were computed as follows:
where D(h)={(i,j)Qxi "xjP=h} contains n(h) pairs and i and j are two sample points. where ni is the number of samples in class i, zi is the mean in class i, n is the total number of classes and ZV is the data average that estimates the zonal mean. Variograms Variograms and cross-variograms were first introduced by Matheron (1971) for characterising spatial variation of natural processes and interpolating the density field by kriging. Fundamentals can be found in Matheron (1971), Journel and Huijbregts (1978) or Armstrong et al. (1992). Stationary, linear and univariate geostatistics have been applied in fisheries science for mapping and
Kriging in IRF-k Intrinsic random functions of order-k (IRF-k) were introduced by Matheron (1973) for mapping nonstationary spatial distributions. Here, we are interested in mapping the non-stationary relative density M(x) the values of which are known on the survey grid. For the purpose of mapping, the values M(x) have in the kriging procedure the status of random variables and will be denoted m(x). In all other parts they are considered as scalars, M(x). Kriging estimates E[m(x)]=M(x) where the expectation is taken over the different realisations. In the IRF-k methodology, m(x) is considered to be the
A geostatistical model in space and time
217
sum of a drift and residuals in a neighbourhood centred on x:
where the f l are mononomials of the coordinates. The IRF-k kriging estimator at point x0 is a linear combination of the known values at points xá:
The total estimation variance in space and time, ó2tot, was approximated by the sum of three major error variance terms: ó2tot =Var(ZVT "Z*VT)=Var(Err1)+Var(Err2)+Var(Err3) where
This linear combination is such that (i) the drift Ókl=0alf l(x0) is filtered (no bias) and (ii) the estimation variance depends only on a function K(h) refered to as the generalised covariance which is a polynomial in h. The unbiased condition results from k+1 contraints on the weights:
The estimation variance is expressed as:
where Káâ is the generalised covariance for the distance between points xá and xâ. The IRF-k kriging system is obtained by minimising the estimation variance with the k+1 contraints and is given by:
where ìl are k+1 Lagrange parameters. The values of the coefficients al do not need to be estimated as they are not used for kriging. Only the degree k (i.e. the shape of the drift) need be estimated. In comparison to Universal kriging, drift and residuals are not dissociated in the IRF-k analysis and they depend on the scale at which they are estimated. Estimation of the order k and of the generalised covariance K(h) were performed simultaneously by using the automatic procedure described in Delfiner (1975) and Matheron and Delfiner (1980).
Spatio-temporal estimation variance The spatio-temporal average to estimate is:
and refers to the error in time,
which is the average error in space and Err3 is referred to as a geometric error, which will be discussed in the next section. Geometric error Matheron (1971) defines geometric error as the error associated with the mean estimate due to the error on the limits of the area on which the mean is estimated. For a systematic survey design, the geometric error variance on the mean is (Matheron, 1971):
where D2 is the process variance in the field V and ó2V/V2 is the relative estimation variance of the area estimate. Here, we neglected the geometric error variance term in space as we considered the limits of the spawning ground known rather precisely because of the link between adults and sediments. For the seven surveys during the spawning period, the spawning area was not observed to change. Surveyed area and spawning area were considered to be the same. On the contrary, limits in time of the spawning season were not known precisely and this uncertainty will affect the spatio-temporal estimation. Thus, the formula above was computed in time. Let D2(ZV) be the variance of the process ZV(t) during the period T. The relative variance on T was computed as suggested by Matheron (1971, p. 46). For regularly spaced surveys in time, limits of the spawning season are considered to be random within the between survey interval. The spawning period T is: T=(n"1)l+å1 +å2
where T is the spawning period and V the spawning ground area. V is assumed not to change in time during T and to be known. On the contrary T must be estimated. The spatio-temporal mean is estimated by:
where n is the number of surveys, l the interval between surveys and å1 and å2 are two independent and uniform random variables on the interval [0,1]. Thus the variance of T is:
218
P. Petitgas
The third term of ó2tot was estimated by:
Now, we give an approximation for computing ó2spa. Recall that the relative density is: Zp(x,t)=M(x)+aM(x)U(x,t).
Error in time The process ZV(t) is assumed to be purely random over time and therefore we have:
The estimation error on the zonal mean of Zp decomposes to the sum of two errors:
where D2(ZV) is the variance of the process ZV(t) over T and nt is the number of surveys. Average error in space The estimation errors on the zonal means ZV(t) are not correlated between surveys as the model has no space– time interaction. Thus the second variance term of ó2tot is:
where the variance Var(ZV(t)"Z*V(t))=ó2 is the average error variance made in space for any value ZV(t) (non-conditional variance). The variance ó2 is not a function of time. It can be computed by conditioning on ZV(t):
The error e1 is not considered to be a random variable in space as M(x) is a scalar. The error e2 is considered random and its mean equals zero. The estimation variance on Zp is: ó2spa =e21 +E[e22]. Difference between true and estimated drift values at the sample points was neglected. The first term e21 was approximated by using the kriged map of M(x). The spatial integral was approximated by the discrete sum of the kriged gridded values. The second term is developed as:
ó2 =E[Var[(ZV(t)"Z*V(t))/ ZV(t)]]+Var[E[(ZV(t)"Z*V(t))/ZV(t)]] As the estimator Z*V(t) is unbiased the second term disappears and we only need to compute the first term. Conditional on ZV(t) we have:
The error variance on the zonal mean of Zp is time invariant because (i) its spatial structure is consistent in time (consistent drift and consistent variogram of residuals) and (ii) the grid of sampled stations is the same for each survey. Thus we have:
where CU(Pxi "yjP) is the covariance of residuals U for the distance Pxi "yjP. The first two terms can be approximated by replacing the integrals by discrete sums on the gridded kriged drift. The third term is a double summation on the sampled grid points and is easily computed. In the present study we computed the third term only. This may result in slightly overestimating E[e22]. Finally, the spatio-temporal estimation variance was approximated by:
ó2c(t)=ó2spaZ2V(t). Now, we take the average of ó2c(t): ó2 =E[ó2spaZ2V(t)]=ó2spaE[Z2V] The second variance term of ó2tot was estimated by:
where ó2spa was approximated by:
A geostatistical model in space and time
219 March 6, 1986
March 2, 1986 47°30'
47°
46°30' N 2°30'
2°30'
2°W
2°W
47°30'
47°
46°30' March 10, 1986
March 18, 1986
Figure 2. Proportional representation of sole egg densities in 1986. Circle radius is proportional to density value. Maximum value for each survey (see Table 1) is attributed the maximum radius and other radii are deduced proportionally.
Data analysis and model fitting Time process for spatial means The date and the spatial mean ZV(t) for each survey are presented in Table 1. Seven surveys with a spatial mean greater than one egg per 10 m3 were considered to fall inside the spawning season and were retained in the present analysis for the model inference. The other surveys were considered to fall outside the spawning season as their spatial average was low and the extension of spawning in space was reduced. They were not used in this analysis. The spawning was modelled by a binary process in time: either we are in the spawning period (ZV(t) greater than 1) and ZV(t) fluctuates around a mean value or we are not in the spawning period and
ZV(t) is null. Arbault et al. (1984) and Le Bec (1985) estimated the spawning season to last approximately 2 mo. However, this fact was not obvious from the information given in Table 1. Fluctuations of ZV(t) were assumed to be stationary and purely random inside the spawning season.
Comparing spatial pattern and aggregation between years Proportional representations of the data are given in Figures 2 and 3. The spatial distributions are similar: the high density area is always in the same location, it has an elliptic shape centred in the middle of the area around the 60 m isobath. This big patch was permanent over
220
P. Petitgas March 4, 1988
March 8, 1987 47°30'
47°
46°30' N 2°30'
2°30'
2°W 2°30'
2°W 2°W 47°30'
47°
46°30' N March 14, 1989 Figure 3. Proportional representation of sole egg densities from 1987 to 1989. Circle radius is proportional to density value. Maximum value for each survey (see Table 1) is attributed the maximum radius and other radii are deduced proportionally.
time during the spawning season in 1986 and is also present in the following three years, showing the same shape and location. The aggregation curves P(T) for the different surveys are presented in Figure 4. No obvious differences were observed between the curves and thus curves were not considered to be significantly different. At any time t of the spawning season, for any year, the eggs were considered to distribute in space in the same manner proportionally to the global mean at time t generating the same patch with the same gradients. The relative densities zp(x,t1) and zp(x,t2) sampled at times t1 and t2 were interpreted as two realisations of the same non-stationary random function M(x).
Proportional effect A proportional relationship between the variance of the residuals computed in time at point x, Et(r2(x,t)), and the drift value M(x) is shown in Figure 5. Thus, temporal variations around the drift are of higher amplitude at the centre than on the sides. This justified the use of a multiplicative model in space and time to interpret the data. The relation between the residual standard deviation and the drift is modelled by a least square linear regression which goes through the origin (Fig. 5). We have: óR(x)=aM(x) and R(x,t)=aM(x)U(x,t) with a=0.58
A geostatistical model in space and time
221
1.0 0.9 0.8 0.7
P(T)
0.6 0.5 0.4 0.3 0.2 0.1 0.0
0.1
0.2
0.3
0.4
0.5 T
0.6
0.7
0.8
0.9
Figure 4. Experimental aggregation curves P(T) for the different surveys.
8 7 6
Et(r2) (x)
5 4 3 2 1
0
1
2
3
4
5 6 M(x)
7
8
9
10
Figure 5. Relation between the drift value M(x) and the temporal standard deviation at point x of the residual R(x,t).
where U(x,t) are residuals with zero mean and unit variance, both in space and time. Plots of the variance of U against values of M(x) and ZV(t) indicated no correlation with M(x) or ZV(t) and are not shown here.
Structure in space and time of the residuals Non-directional variograms were computed on the residuals U(x,t), for each survey with a distance bin of
five nautical miles (5 nmi). For the different times t (seven surveys), the experimental variograms were similar and their mean was calculated for each distance bin. We thus infer a time-averaged spatial structure for the residuals. The variogram of the residuals is shown in Figure 6 where there appears to be a clear structuring with a correlation range of 11 nmi. The temporal variogram of U was computed at each point x and the average of these variograms was taken
222
P. Petitgas
1.25
Variogram
1.00
0.75
0.50
0.25
0
10
20 Distance (nm)
30
40
Figure 6. Variogram in space of the residuals U. The model fitted is a nugget of value 0.2 added to a spherical model with sill 1 and range 11 nmi.
over all the stations. The time distances were binned each year. The variogram stayed flat meaning no structure in time and is not shown. The intra-annual variability between surveys is of the same order of magnitude as the inter-annual variability. Non-directional cross-variograms were computed on the residuals U(x,t) between pairs of surveys. Distances were binned by classes of 5 nmi. There are 7 surveys, thus (72)=21 cross-variograms were calculated. Their average was taken over time. This time averaging was done in bins of one year. The spatio-temporal crossvariograms are shown in Figure 7. For the first distance (replicated stations) the cross-variogram values are approximate zero. After the first bin (5 nmi) all cross variogram values stabilise around a sill of "0.1. This small negative correlation is interpreted as being due to the inference procedure of the model because the drift M(x) is estimated by averaging in time. Thus, at point x, positive residuals at time t are necessarily counterbalanced at another time by negative ones because the temporal mean of the residuals at x equals zero. The cross-variograms were considered to be flat, indicating no space–time interaction. The cross-variograms between surveys and with translations in space were not computed. We assumed they were also flat. A more practical statement of the residual structure would be as follows: the residual highs and lows around the time invariant drift have on average an isotropic
dimension of 11 nmi, they move in space over time according to a pure random process and without memory of their former location.
Use of the model Mapping the drift by kriging with IRF-2 In Figure 8 the drift M(x) is presented as estimated by averaging in time at each station. The western margin of the spawning area was only sampled by one station but it was assumed that this station was close enough to the western limit of the survey grid for modelling purposes. In the north-east part, the spawning ground extends slightly to the coast. The drift M(x) is not a smooth surface in space. The neighbourhoods chosen for inference of the IRF-k model were circles with a diameter of 50 Nm with 850 neighbourhoods located at random on the surveyed area. The average number of samples in each neighbourhood was 12. The drift modelled in these neighbourhoods was a quadratic function (order k equals two). The generalised covariance estimated was: K(h)=0.539 h. More details on the different fitting trials can be found in Petitgas (1991). The interpolation grid had a mesh size of 2#2 nmi2 (average distance between nearest stations is 5 nmi). The interpolation grid extended the survey grid by one inter-station distance on all sides. The kriged map is presented in Figure 9.
A geostatistical model in space and time
223
0.1
Crossvariogram
0.0
–0.1
–0.2
–0.3
0
10
20 Distance (nm)
30
40
Figure 7. Cross-variograms between space and time for the residuals U. -=same year, /=one year lag, + =two year lag.
47°30' N
47°
46°30' 2°30' W
2°
Figure 8. Proportional representation of the time invariant drift M(x). Maximum drift value is 10.1 eggs per 10 m "3.
Estimation variance for the spatio-temporal mean estimate A numerical application of Equation (1) was performed using four surveys regularly spaced in time during the spawning season. The spawning season was assumed to last 2 mo from Le Bec (1985). Values of ZV(t) used were
those of year 1986 as four surveys were performed during the spawning season in this year. This approach is possible because in the model, (i) the spatial mean is a random variable in time during the spawning season and (ii) each survey gives one realisation which has no residual correlation with another one.
224
P. Petitgas 47°30' N
2
5
7
0
6
1
47°
8 4
1
2
1
3
3 2
0
4
0
46°30' 2°30' W
2°
Figure 9. Kriged map of the time invariant drift M(x) obtained by non-stationary kriging in IRF-k. Maximum estimated drift value is 9.2 eggs per 10 m "3.
Recall that three error terms were retained (in the section on spatio-temporal estimation variance): Err1 is the error in time, Err2 is the error in space and Err3 is the geometric error made in time. Using same notations as previously, we have:
nt =4; D2(ZV)=0.4522; E[Z2V]=5.6984
It takes 4 d to survey the present grid of 35 stations. In the 2 mo period of the spawning, eight surveys at the maximum are possible. We computed the formula of ó2tot with eight surveys (all other parameters unchanged). The relative error ótot/Z*VT reduces to 11.2%. This is considered the lowest obtainable precision for 1986 with the survey grid of 35 stations. A precision between 17.1% and 11.2% could be obtained by putting more surveys in time and less stations in space.
ó2tot =Var(Err1)+Var(Err2)+Var(Err3)=0.1130+ 0.0393+0.0047=0.157.
Discussion
ó2spa =e21 +E[e22]=0.0004+0.0272
The relative spatio-temporal error is:
The contribution of each error term to the total variance ó2tot are 72% for Err1, 25% for Err2 and 3% for Err3. The major source of imprecision lies in the estimation of the average value over time of the spatial means ZV(t). In the variance ó2spa, the term e21 is small meaning that the major error in the estimation of the drift surface M(x) comes from the residuals U. The relative error in space is: óspa =16.6%. In the studied surveys, more effort had been put in space than in time but final precision is good.
The major spatio-temporal structural characteristics of the sole egg distributions were analysed and a multiplicative model was used. Multiplicative models were introduced in fisheries science by Robson (1966) for calibrating fishing effort and generalised linear models (McCullagh and Nelder, 1989) have been used for analysing spatio-temporal survey data (Stefansson, 1991; Anon., 1992). The contribution of this paper was (i) to establish the model in terms of spatial structure inferred on untransformed data and (ii) to use this model to derive the corresponding estimation variance for a spatio-temporal estimation of the mean egg density. In the present approach, we took advantage of the fact that the same grid of stations in space was surveyed
A geostatistical model in space and time at different times. At a given time, the average density over the grid of stations was used to estimate a relative density at each station. Spatial distribution patterns of the relative density were then shown to be constant over the years and within a year. In the model inference, the data average was used as if it was the true zonal mean. This has simplified the computations a great deal in deriving the variance equations. The applicability of this approach clearly depends on the number of repetitions in time. If too few, estimation error in the drift and residuals will be studied as process variability making the approach useless. This was addressed by Petitgas (1991) for the variogram of the residuals which had little noise in them. The process of the zonal mean was assumed to be a purely random and stationary variable over the spawning season. During this period a time invariant drift was characterised on the relative variable. Increase and decrease of zonal mean density at the beginning and end of the spawning season were not modelled. MacCall (1990) has observed for anchovy larvae that large fluctuations of the zonal mean (abundance) were associated with different patterns in the spatial distribution. In this study, such situation was not observed between years. The estimation variance formula could be used for optimising survey effort in terms of number of stations in space vs. number of surveys in time. This requires to estimate precisely the statistics of the process ZV(t) and to generate the drift values at the sampling stations. In the present case study, the major source of variation left was in time on the process ZV(t) and thus it could be suggested to increase the effort in time and reduce it in space.
Acknowledgements This research was a contract between Institut Français de Recherche pour l’Exploitation de la Mer (IFREMER) and École des Mines de Paris (Centre de Géostatistique). Data were collected by IFREMER with the R/V ‘‘Thalassa’’. Drs C. Koutsikopoulos, J. Rivoirard, and P. Fréon, as well as two anonymous referees and assistant editor S. J. Smith are thanked for their insights and comments that helped improve the manuscript.
References Anon. 1992. Report of the workshop on the analysis of trawl survey data. ICES CM1992/D:6. Arbault, S., Camus, P., and Le Bec, C. 1984. Estimation du stock de soles (Solea vulgaris, Quensel 1806) dans le Golfe de Gascogne à partir de la production d’oeufs. Journal of Applied Ichtyology, 4: 145–156. Armstrong, M., Renard, D., Rivoirard, J., and Petitgas, P. 1992. Geostatistics for fish survey data, course publicised by ICES, Centre de Géostatistique, ENSMP, Fontainebleau.
225
Cressie, N. 1991. Statistics for spatial data. John Wiley and Sons, New York. Dagum, C. 1985. Lorenz curve. In Encyclopedia of Statistical Sciences, Vol. 5. Ed. by S. Kotz and N. Johnson. John Wiley and Sons, New York. 156–161 pp. Delfiner, P. 1975. Linear estimation of non-stationary spatial phenomena. In Proceedings of NATO ASI, Advanced Geostatistics in the Mining Industry, Rome 1975. D. Reidel and Co., Dordrecht. 49–68 pp. Dardignac, J. 1988. Les pêcheries du Golfe de Gascogne – bilan des connaissances, Rapports Scientifiques et Techniques de l’IFREMER, no. 9. Foote, K., and Stefansson, G. 1993. Definition of the problem of estimating fish abundance over an area from acoustic line-transects measurements of density. ICES Journal of Marine Science, 50: 369–381. Hilborn, R., and Walters, C. 1992. Quantitative fisheries stock assessment: choice dynamics and uncertainty. Chapman and H all, New York. Houghton, R. 1987. The consistency of the spatial distributions of young gadoids with time. ICES CM1987/D:15. Journel, A., and Huijbregts, Ch. 1978. Mining geostatistics. Academic Press, London. Koutsikopoulos, C. 1991. Recrutement de la sole (Solea solea, L.) du Golfe de Gascogne: influence de l’hydrologie et de l’hydrodynamisme. Thèse de Doctorat, Université de Bretagne Occidentale, Brest. Le Bec, G. 1985. Cycle sexuel et fécondité de la sole (Solea solea, L.) du Golfe de Gascogne. Revue des Travaux de l’Institut des Pêches Maritimes, 47: 179–189. MacCall, A. 1990. Dynamic geography of marine fish populations. University of Washington Press, Seattle. McCullagh, P., and Nelder, J. 1989. Generalised linear models. 2nd edit. Chapman and Hall, London. Matheron, G. 1971. The theory of regionalised variables and its applications, Les Cahiers du Centre de Morphologie Mathématique, Fascicule 5, Centre de Géostatistique, ENSMP, Fontainebleau. (English version.) Matheron, G. 1973. Intrinsic random functions and their application. Advances in Applied Probability, 5: 439–468. Matheron, G. 1981. La sélectivité des distributions, Note N-686, Centre de Géostatistique, ENSMP, Fontainebleau. Matheron, G., and Delfiner, P. 1980. The intrinsic random functions of order k, course C-80, Centre de Géostatistique, ENSMP, Fontainebleau. (English version.) Myers, R., and Stokes, T. 1989. Density-dependent habitat utilisation of groundfish and the improvement of research surveys. ICES CM1989/D:15. Pelletier, D., and Parma, A. 1994. Spatial distribution of pacific halibut (Hippoglossus stenolepis): an application of geostatistics to longline survey data. Canadian Journal of Fisheries and Aquatic Science, 51: 1506–1518. Petitgas, P. 1991. Contributions géostatistiques à la biologie des pêches maritimes. Thèse de Doctorat, Centre de Géostatistique, ENSMP, Paris. Riley, J. 1974. The distribution and mortality of sole eggs. In The early life history of fish. Ed. by J. Blaxter, Springer Verlag, Berlin. Robson, D. 1966. Estimation of the relative fishing of individuals ships. ICNAF Research Bulletin, 3: 5–14. Stefansson, G. 1991. Analysis of groundfish survey data: combining the GLM and delta approaches. ICES CM1991/D:9. Sullivan, P. 1991. Stock abundance estimation using depthdependent drifts and spatially correlated variation. Canadian Journal of Fisheries and Aquatic Science, 48: 1691–1703.
