Spatio-temporal smoothing of large ecological datasets describing species distributions can be made challenging by high computational costs and deficiencies in the available data. We present an application of a GAM-based smoothing method to a large ordinal categorical dataset on the distribution of wildebeest in the Serengeti ecosystem

Ferguson, Elaine A.

Husmeier, Dirk

Matthiopoulos, Jason

English

Enlighten

Constructing wildebeest density distributionsby spatio-temporal smoothing of ordinalcategorical data using GAMsElaine A. Ferguson1, Jason Matthiopoulos1, Dirk Husmeier21 Institute of Biodiversity, Animal Health and Comparative Medicine, College ofMedical, Veterinary and Life Sciences, University of Glasgow, Glasgow, UK2 School of Mathematics and Statistics, College of Science and Engineering, Uni-versity of Glasgow, Glasgow, UKE-mail for correspondence: e.ferguson.2@research.gla.ac.ukAbstract: Spatio-temporal smoothing of large ecological datasets describingspecies distributions can be made challenging by high computational costs anddeficiencies in the available data. We present an application of a GAM-basedsmoothing method to a large ordinal categorical dataset on the distribution ofwildebeest in the Serengeti ecosystem.Keywords: GAMs; Ordinal categorical; Smoothing; Spatio-temporal; Wilde-beest.1 IntroductionSpatio-temporal smoothing of species distribution data has many potentialuses in ecology; for example, to provide a smooth density function thatcan be used with gradient matching approaches (Xun et al. 2013) to fitpartial differential equation (PDE) models of animal movement. A rangeof smoothing methods (kernel density estimation, splines, Gaussian pro-cesses, etc.) have been developed in the statistical literature. However, thepracticalities and expense involved in collecting species distribution dataover large areas in the field can mean that the data are not in a formthat these methods can readily be applied to. Ordinal categorical data, forexample, may be collected when it is infeasible to accurately count all indi-viduals in a population, so that the abundance at each point in space andtime is instead estimated as belonging to a broader abundance category. AThis paper was published as a part of the proceedings of the 32nd Interna-tional Workshop on Statistical Modelling (IWSM), Johann Bernoulli Institute,Rijksuniversiteit Groningen, Netherlands, 3–7 July 2017. The copyright remainswith the author(s). Permission to reproduce or extract any parts of this abstractshould be requested from the author(s).2 Smoothing ordinal categorical data using GAMsrelatively small number of approaches have been developed for smoothingdata of this type, where we need to recover the underlying true density ofindividuals from the categories (Chu and Ghahramani 2005, Wood et al.2016). Smoothing large datasets in multiple dimensions can also be madechallenging by high computational costs. Methods that allow smoothingof these datasets even when computational resources are limited wouldtherefore be very useful. Here we present an application of a method forapplying spatio-temporal smoothing to a large ordinal categorical dataseton the distribution of wildebeest in the Serengeti ecosystem of Tanzaniaand Kenya.2 MethodsThe wildebeest distribution data, which have been described and utilisedin a number of previous studies (Norton-Griffiths 1973, Maddock 1979,Boone et al. 2006, Holdo et al. 2009), were obtained from monthly aerialsurveys of the Serengeti ecosystem during the period from August 1969to August 1972. Each cell in a grid of 25km2 cells was assigned to one offive wildebeest abundance categories: 0, 1-25, 26-250, 251-2,500 and >2,500individuals per 25km2. There were 2,576 cells making up the spatial grid,all of which were sampled on 33 occasions during the time period, resultingin a large dataset with a total 85,008 data points.To smooth the data in time, t, and the two spatial dimensions (x, y) we fit-ted GAMs (generalised additive models) with a tensor product (composedof cubic regression spline smooths, where overfitting was prevented by pe-nalisation of the integral of the squared second derivatives) between thesethree variables using the mgcv package (Wood 2011) in R (R Core Team2015). We used the ordinal categorical GAM method described in Wood etal. (2016), where the linear predictor gives the value of a latent variable,here representing the wildebeest density underlying the ordinal categories.The cut-off points that demarcate the five ordinal categories were speci-fied, and the probability that a point in space and time belongs to a givencategory equals the probability that the latent variable lies between thecorresponding category cut-offs at that point.In Wood et al. (2016), the latent function can range from −∞ to ∞, butwe know that wildebeest density has a minimum 0 and a finite maximumWmax. We can introduce these constraints by applying a sigmoidal trans-formation to the latent function L after the GAM has been fitted, giving apreliminary wildebeest density Wˆ as follows:Wˆ (x, y, t) =Wmax1 + exp (−L (x, y, t)) (1)Note that this also required that an inverse sigmoid transform be appliedto the category cut-offs c prior to the GAM fitting:Ferguson et al. 3c¯ = − log(Wmaxc− 1)(2)Wmax was estimated by first assuming that the wildebeest densities in thegrid cells assigned to the lower four ordinal categories, which had knownupper and lower bounds, were equal to the mid-points of those categories.The sum of the densities in these lower category cells for each month wasthen subtracted from the total number of wildebeest WT known to be inthe region from a population count in 1971 (Norton-Griffiths 1973). Theremaining wildebeest for each month were assumed to be divided evenlybetween the cells in the highest ordinal category (which was unboundedabove) for that month. We took Wmax to be the largest wildebeest densityestimated for cells in the highest abundance category over all months.Even after applying sensible upper and lower bounds to the latent function,large fluctuations in the area under Wˆ (which represents the total numberof wildebeest in the region) can occur over time. This is undesirable, sincewe expect wildebeest numbers to remain relatively stable at WT over thetime period of interest. We therefore consider the normalised wildebeestdensity W¯ , where the total number of animals is maintained at WT bynormalising Wˆ as follows:W¯ (x, y, t) =Wˆ (x, y, t)WT∫Wˆ (x, y, t)dxdy(3)Due to computational time and memory constraints, a sufficiently flexibleGAM could not be fitted to the entire large dataset simultaneously. Wetherefore divided the time series into three contiguous intervals and fitteda GAM in (x, y, t) to each interval separately. Each GAM had 20 knotsin the marginal smooth in each spatial dimension, and a number of knotsin the marginal smooth in time that was equal to the number of timepoints present in the data subset to which the GAM was fitted (11 or 12).This resulted in the effective degrees of freedom, which are determined bythe degree of penalization (selected during fitting) applied to the integral ofthe squared second derivatives, being considerably lower than the maximumnumber available, suggesting that the number of knots was sufficient (Wood2006). The three GAMs were joined together by averaging at the link timesli (i ∈ 1, 2), with smoothness being maintained by allowing the influence ofeach GAM on the others to decline smoothly, according to the parameterσ, as distance from the point of joining increased. For a given point (x¯, y¯, t¯),therefore, we obtain a final estimate of wildebeest density W byW (x¯, y¯, t¯) = W¯GAMj (x¯, y¯, t¯) +2∑i=1ai exp(−(t¯− li)22σ2)mi (t¯) (4)Here W¯GAMj is the normalised wildebeest density obtained from the GAMfitted to time interval j, where4 Smoothing ordinal categorical data using GAMsj = 1 if t¯ ≤ l12 if l1 < t¯ ≤ l23 if t¯ > l2(5)The ai are given byai (x¯, y¯, li) =W¯GAMi (x¯, y¯, li)− W¯GAMi+1 (x¯, y¯, li)2(6)FIGURE 1. Model fit in space at three different time points. A-C: The wildebeestspatial distribution data for months 1, 18 and 35. D-F: The smooth wildebeestdensity distribution estimated in space by the model for months 1, 18 and 35.The two contours indicate the boundaries between abundance categories 0, 1 and2. G-I: Estimated wildebeest abundance categories based on D-F.Ferguson et al. 5and the mi, which ensure that the adjustments are made in the correctdirection on either side of each link point, aremi (t) ={ −1 if t¯ ≤ li1 if t¯ > li(7)If the influence of the adjoining GAMs declines too slowly with distancefrom the link points, relative to the rate at which changes occur in W¯GAMi(i.e. σ is too large), unrealistic negative values of W can occur. We thereforetuned σ by starting with a relatively large value and gradually decreasingit until no negative values of W occurred.3 Results and ConclusionThe method described was found to successfully produce a smooth functionin space that resembles the original data (Figure 1). The resulting functionis also observed to be smooth in time, with no evidence that the wilde-beest density changes either more slowly or more rapidly around the GAMlink times than it does elsewhere in the time period (Figure 2). This sug-gests that our approach of linking models that have been fitted to subsetsof a larger dataset is a promising means of reducing the high computa-tional costs of smoothing large datasets in multiple dimensions. Using thismethod, we have recovered realistically bounded wildebeest abundance es-timates from coarse ordinal categories; an ability that could be useful inthe field of ecology where such imperfect data are common. By producinga smooth surface from which spatial and temporal gradients in density canFIGURE 2. Changes in the estimated wildebeest density in six grid cells (indi-cated by different colours/line types) over the time period of interest. The linktimes between the three GAMs are indicated by dashed vertical lines.6 Smoothing ordinal categorical data using GAMsbe calculated, our method also promises to enable statistical inference forPDE models of animal movement using the gradient matching approach ofXun et al. (2013), which we will investigate in future work.Acknowledgments: Special thanks to Ricardo M. Holdo for providingaccess to the wildebeest distribution data. E.A.F. is funded by a Universityof Glasgow Lord Kelvin/Adam Smith PhD scholarship.ReferencesBoone, R.B., Thirgood, S.J. and Hopcraft J.G.C. (2006). Serengeti wilde-beest migratory patterns modeled from rainfall and new vegetationgrowth. Ecology, 87, 1987 – 1994.Chu, W. and, Ghahramani, Z. (2005). Gaussian processes for ordinal re-gression. Journal of Machine Learning Research, 6, 1019 – 1041.Holdo, R. M., Holt, R. D. and Fryxell, J. M. (2009). Opposing rainfall andplant nutritional gradients best explain the wildebeest migration inthe Serengeti. The American Naturalist, 173, 431 – 445.Maddock, L. (1979). The migration and grazing succession. In: Serengeti:dynamics of an ecosystem, Sinclair, A. R. E., and Norton-Griffiths,M. (editors), Chicago: University of Chicago Press, 104 – 129.Norton-Griffiths, M. (1973). Counting the Serengeti migratory wildebeestusing two-stage sampling. East African Wildlife Journal, 11, 135 –149.R Core Team (2015). R: A language and environment for statistical com-puting. R Foundation for Statistical Computing, Vienna, Austria.Wood, S.N. (2006). Generalized Additive Models: An Introduction with R.CRC Press.Wood, S.N. (2011). Fast stable restricted maximum likelihood and marginallikelihood estimation of semiparametric generalized linear models.Journal of the Royal Statistical Society, Series B, 73, 3 – 36.Wood, S.N., Pya, N., and Sfken, B. (2016). Smoothing parameter and modelselection for general smooth models. arXiv:1511.03864v2.Xun, X., Cao, J., Mallick, B., Maity, A., and Carroll, R.J. (2013). Param-eter Estimation of Partial Differential Equation Models. Journal ofthe American Statistical Association, 108, 1009 – 1020.

Constructing Wildebeest Density Distributions by Spatio-temporal Smoothing of Ordinal Categorical Data Using GAMs

Enlighten: Publications

Constructing wildebeest density distributions
by spatio-temporal smoothing of ordinal
categorical data using GAMs
Elaine A. Ferguson1, Jason Matthiopoulos1, Dirk Husmeier2
1 Institute of Biodiversity, Animal Health and Comparative Medicine, College of
Medical, Veterinary and Life Sciences, University of Glasgow, Glasgow, UK
2 School of Mathematics and Statistics, College of Science and Engineering, Uni-
versity of Glasgow, Glasgow, UK
E-mail for correspondence: e.ferguson.2@research.gla.ac.uk
Abstract: Spatio-temporal smoothing of large ecological datasets describing
species distributions can be made challenging by high computational costs and
deficiencies in the available data. We present an application of a GAM-based
smoothing method to a large ordinal categorical dataset on the distribution of
wildebeest in the Serengeti ecosystem.
Keywords: GAMs; Ordinal categorical; Smoothing; Spatio-temporal; Wilde-
beest.
1 Introduction
Spatio-temporal smoothing of species distribution data has many potential
uses in ecology; for example, to provide a smooth density function that
can be used with gradient matching approaches (Xun et al. 2013) to fit
partial differential equation (PDE) models of animal movement. A range
of smoothing methods (kernel density estimation, splines, Gaussian pro-
cesses, etc.) have been developed in the statistical literature. However, the
practicalities and expense involved in collecting species distribution data
over large areas in the field can mean that the data are not in a form
that these methods can readily be applied to. Ordinal categorical data, for
example, may be collected when it is infeasible to accurately count all indi-
viduals in a population, so that the abundance at each point in space and
time is instead estimated as belonging to a broader abundance category. A
This paper was published as a part of the proceedings of the 32nd Interna-
tional Workshop on Statistical Modelling (IWSM), Johann Bernoulli Institute,
Rijksuniversiteit Groningen, Netherlands, 3–7 July 2017. The copyright remains
with the author(s). Permission to reproduce or extract any parts of this abstract
should be requested from the author(s).
2 Smoothing ordinal categorical data using GAMs
relatively small number of approaches have been developed for smoothing
data of this type, where we need to recover the underlying true density of
individuals from the categories (Chu and Ghahramani 2005, Wood et al.
2016). Smoothing large datasets in multiple dimensions can also be made
challenging by high computational costs. Methods that allow smoothing
of these datasets even when computational resources are limited would
therefore be very useful. Here we present an application of a method for
applying spatio-temporal smoothing to a large ordinal categorical dataset
on the distribution of wildebeest in the Serengeti ecosystem of Tanzania
and Kenya.
2 Methods
The wildebeest distribution data, which have been described and utilised
in a number of previous studies (Norton-Griffiths 1973, Maddock 1979,
Boone et al. 2006, Holdo et al. 2009), were obtained from monthly aerial
surveys of the Serengeti ecosystem during the period from August 1969
to August 1972. Each cell in a grid of 25km2 cells was assigned to one of
five wildebeest abundance categories: 0, 1-25, 26-250, 251-2,500 and >2,500
individuals per 25km2. There were 2,576 cells making up the spatial grid,
all of which were sampled on 33 occasions during the time period, resulting
in a large dataset with a total 85,008 data points.
To smooth the data in time, t, and the two spatial dimensions (x, y) we fit-
ted GAMs (generalised additive models) with a tensor product (composed
of cubic regression spline smooths, where overfitting was prevented by pe-
nalisation of the integral of the squared second derivatives) between these
three variables using the mgcv package (Wood 2011) in R (R Core Team
2015). We used the ordinal categorical GAM method described in Wood et
al. (2016), where the linear predictor gives the value of a latent variable,
here representing the wildebeest density underlying the ordinal categories.
The cut-off points that demarcate the five ordinal categories were speci-
fied, and the probability that a point in space and time belongs to a given
category equals the probability that the latent variable lies between the
corresponding category cut-offs at that point.
In Wood et al. (2016), the latent function can range from −∞ to ∞, but
we know that wildebeest density has a minimum 0 and a finite maximum
Wmax. We can introduce these constraints by applying a sigmoidal trans-
formation to the latent function L after the GAM has been fitted, giving a
preliminary wildebeest density Wˆ as follows:
Wˆ (x, y, t) =
Wmax
1 + exp (−L (x, y, t)) (1)
Note that this also required that an inverse sigmoid transform be applied
to the category cut-offs c prior to the GAM fitting:
Ferguson et al. 3
c¯ = − log
(
Wmax
c
− 1
)
(2)
Wmax was estimated by first assuming that the wildebeest densities in the
grid cells assigned to the lower four ordinal categories, which had known
upper and lower bounds, were equal to the mid-points of those categories.
The sum of the densities in these lower category cells for each month was
then subtracted from the total number of wildebeest WT known to be in
the region from a population count in 1971 (Norton-Griffiths 1973). The
remaining wildebeest for each month were assumed to be divided evenly
between the cells in the highest ordinal category (which was unbounded
above) for that month. We took Wmax to be the largest wildebeest density
estimated for cells in the highest abundance category over all months.
Even after applying sensible upper and lower bounds to the latent function,
large fluctuations in the area under Wˆ (which represents the total number
of wildebeest in the region) can occur over time. This is undesirable, since
we expect wildebeest numbers to remain relatively stable at WT over the
time period of interest. We therefore consider the normalised wildebeest
density W¯ , where the total number of animals is maintained at WT by
normalising Wˆ as follows:
W¯ (x, y, t) =
Wˆ (x, y, t)WT∫
Wˆ (x, y, t)dxdy
(3)
Due to computational time and memory constraints, a sufficiently flexible
GAM could not be fitted to the entire large dataset simultaneously. We
therefore divided the time series into three contiguous intervals and fitted
a GAM in (x, y, t) to each interval separately. Each GAM had 20 knots
in the marginal smooth in each spatial dimension, and a number of knots
in the marginal smooth in time that was equal to the number of time
points present in the data subset to which the GAM was fitted (11 or 12).
This resulted in the effective degrees of freedom, which are determined by
the degree of penalization (selected during fitting) applied to the integral of
the squared second derivatives, being considerably lower than the maximum
number available, suggesting that the number of knots was sufficient (Wood
2006). The three GAMs were joined together by averaging at the link times
li (i ∈ 1, 2), with smoothness being maintained by allowing the influence of
each GAM on the others to decline smoothly, according to the parameter
σ, as distance from the point of joining increased. For a given point (x¯, y¯, t¯),
therefore, we obtain a final estimate of wildebeest density W by
W (x¯, y¯, t¯) = W¯GAMj (x¯, y¯, t¯) +
2∑
i=1
ai exp
(
−(t¯− li)2
2σ2
)
mi (t¯) (4)
Here W¯GAMj is the normalised wildebeest density obtained from the GAM
fitted to time interval j, where
4 Smoothing ordinal categorical data using GAMs
j =
 1 if t¯ ≤ l12 if l1 < t¯ ≤ l2
3 if t¯ > l2
(5)
The ai are given by
ai (x¯, y¯, li) =
W¯GAMi (x¯, y¯, li)− W¯GAMi+1 (x¯, y¯, li)
2
(6)
FIGURE 1. Model fit in space at three different time points. A-C: The wildebeest
spatial distribution data for months 1, 18 and 35. D-F: The smooth wildebeest
density distribution estimated in space by the model for months 1, 18 and 35.
The two contours indicate the boundaries between abundance categories 0, 1 and
2. G-I: Estimated wildebeest abundance categories based on D-F.
Ferguson et al. 5
and the mi, which ensure that the adjustments are made in the correct
direction on either side of each link point, are
mi (t) =
{ −1 if t¯ ≤ li
1 if t¯ > li
(7)
If the influence of the adjoining GAMs declines too slowly with distance
from the link points, relative to the rate at which changes occur in W¯GAMi
(i.e. σ is too large), unrealistic negative values of W can occur. We therefore
tuned σ by starting with a relatively large value and gradually decreasing
it until no negative values of W occurred.
3 Results and Conclusion
The method described was found to successfully produce a smooth function
in space that resembles the original data (Figure 1). The resulting function
is also observed to be smooth in time, with no evidence that the wilde-
beest density changes either more slowly or more rapidly around the GAM
link times than it does elsewhere in the time period (Figure 2). This sug-
gests that our approach of linking models that have been fitted to subsets
of a larger dataset is a promising means of reducing the high computa-
tional costs of smoothing large datasets in multiple dimensions. Using this
method, we have recovered realistically bounded wildebeest abundance es-
timates from coarse ordinal categories; an ability that could be useful in
the field of ecology where such imperfect data are common. By producing
a smooth surface from which spatial and temporal gradients in density can
FIGURE 2. Changes in the estimated wildebeest density in six grid cells (indi-
cated by different colours/line types) over the time period of interest. The link
times between the three GAMs are indicated by dashed vertical lines.
6 Smoothing ordinal categorical data using GAMs
be calculated, our method also promises to enable statistical inference for
PDE models of animal movement using the gradient matching approach of
Xun et al. (2013), which we will investigate in future work.
Acknowledgments: Special thanks to Ricardo M. Holdo for providing
access to the wildebeest distribution data. E.A.F. is funded by a University
of Glasgow Lord Kelvin/Adam Smith PhD scholarship.
References
Boone, R.B., Thirgood, S.J. and Hopcraft J.G.C. (2006). Serengeti wilde-
beest migratory patterns modeled from rainfall and new vegetation
growth. Ecology, 87, 1987 – 1994.
Chu, W. and, Ghahramani, Z. (2005). Gaussian processes for ordinal re-
gression. Journal of Machine Learning Research, 6, 1019 – 1041.
Holdo, R. M., Holt, R. D. and Fryxell, J. M. (2009). Opposing rainfall and
plant nutritional gradients best explain the wildebeest migration in
the Serengeti. The American Naturalist, 173, 431 – 445.
Maddock, L. (1979). The migration and grazing succession. In: Serengeti:
dynamics of an ecosystem, Sinclair, A. R. E., and Norton-Griffiths,
M. (editors), Chicago: University of Chicago Press, 104 – 129.
Norton-Griffiths, M. (1973). Counting the Serengeti migratory wildebeest
using two-stage sampling. East African Wildlife Journal, 11, 135 –
149.
R Core Team (2015). R: A language and environment for statistical com-
puting. R Foundation for Statistical Computing, Vienna, Austria.
Wood, S.N. (2006). Generalized Additive Models: An Introduction with R.
CRC Press.
Wood, S.N. (2011). Fast stable restricted maximum likelihood and marginal
likelihood estimation of semiparametric generalized linear models.
Journal of the Royal Statistical Society, Series B, 73, 3 – 36.
Wood, S.N., Pya, N., and Sfken, B. (2016). Smoothing parameter and model
selection for general smooth models. arXiv:1511.03864v2.
Xun, X., Cao, J., Mallick, B., Maity, A., and Carroll, R.J. (2013). Param-
eter Estimation of Partial Differential Equation Models. Journal of
the American Statistical Association, 108, 1009 – 1020.


http://eprints.gla.ac.uk/144327/7/144327.pdf

Constructing Wildebeest Density Distributions by Spatio-temporal Smoothing of Ordinal Categorical Data Using GAMs

Abstract

Similar works

Full text

Available Versions

Enlighten

Enlighten: Publications