3418
Ecology,85(12), 2004, pp. 3418±3427
q2004 by the Ecological Society of America
NONADDITIVE EFFECTS OF THE ENVIRONMENT ON THE SURVIVAL OF
A LARGE MARINEFISH POPULATION
L
ORENZO
C
IANNELLI,1,4
K
UNG
-SIKCHAN
,
2
K
EVIN
M. B
AILEY
,
3ANDN
ILSCHR.STENSETH
4,5
1
University of Washington, Joint Institute for the Study of Atmosphere and Ocean, Seattle, Washington 98105 USA
2
Department of Statistics and Actuarial Sciences, University of Iowa, Iowa City, Iowa 52242 USA
3Alaska Fisheries Science Center, NOAA, Seattle, Washington 98115 USA
4
Centre for Ecological and Evolutionary Synthesis (CEES), Department of Biology, University of Oslo,
Oslo N-0316 Norway
Abstract.Climate can affect population dynamics in indirect ways via nonadditive
forcing by external variables on internal demographic rates. Current analytical techniques,
employed in population ecology, fail to explicitly include nonadditive interactions between
internal and external variables, and therefore cannot ef®ciently address indirect climate
effects. Here, we present the results of an analysis, employing speci®cally developed sta-
tistical methodology, on density-dependent survival of walleye pollock (
Theragra chal-
cogramma
) prerecruitment stages in relation to background environmental variables in the
Gulf of Alaska. We found that spring winds and water temperature mediate the intensity
of density-dependent survival from the eggs to the age-0 stage. Fall water temperature and
juvenile pollock predator abundance mediate density dependence from the age-0 to the age-
1 stage. The inclusion of such nonadditive and nonlinear effects in a population dynamics
model improved our ability to simulate pollock recruitment. Our results point to the im-
portance of understanding nonadditive and nonlinear interactions between external (climate)
and internal factors in the presence of underlying environmental variation. These topics
are discussed in the context of current research priorities in population ecology and con-
servation biology.
Key words: density dependence; environmental change; Gulf of Alaska; nonadditivity; phase
dependence;
Theragra chalcogramma;walleye pollock.
INTRODUCTION
The ecological effects of environmental change on
population dynamics have received much attention,
particularly with respect to climate variation (Stenseth
et al. 2002). It is well known that self-regulating and
internal forces, such as density dependence (e.g., Sten-
seth et al. 1999), and external environmental variables,
such as climate (e.g., Chavez et al. 2003), can directly
affect individual survival. However, in the presence of
both internal and external control, an understanding of
their joint effect is less clear. A survey of the literature
indicates that internal and external variables may in-
directly and nonadditively in¯uence individual surviv-
al. For example, through in situ experimental studies
it was found that the coral reef ®shes
Chromis cyanea
andThalassoma hardwickeshow different levels of
density-dependent mortality in different patch reefs,
each characterized by a unique set of external envi-
ronmental variables (Hixon and Carr 1997, Shima and
Osenberg 2003). Also, blue petrels
Halobalena caeru-
lea
in the southern Indian Ocean exhibit sharp declines
in abundance during periods of high population density
5
Corresponding author.E-mail: [email protected]
and unfavorable oceanographic conditions (Barbraud
and Weimerskirch 2003).
Most analytical approaches used in population ecol-
ogy assume additive effects of external variables on
population demographic parameters (e.g., Grenfell et
al. 1992, Cury et al. 1995, Stenseth 1999, Stenseth et
al. 1999, Dennis and Otten 2000), and are therefore
inadequate to address nonadditive interactions. In the
present study, we focus on the structure of density-
dependent survival of prerecruitment Gulf of Alaska
(GOA) walleye pollock (
Theragra chalcogramma)in
relation to environmental variables known to affect pol-
lock early life survival. We develop and test a new
model approach where the shape of the density-depen-
dent function may change according to whether some
external environmental covariates are below or above
a threshold value; hence our model explicitly includes
nonadditive interactions. We compare the results from
the threshold formulation with those from a nonaddi-
tive but continuous (rather than threshold) formulation
and with those from a fully additive model. Through
these comparisons, we address not only the effect of
internal and external variables on pollock survival but
also the type of dynamics involved (i.e., additive, non-
additive threshold, nonadditive continuous).
Walleye pollock currently constitutes the second
largest single species ®shery in the world (FAO 2002;

December 2004 3419 NONADDITIVE ENVIRONMENTAL CONTROL
F
IG
. 1. The Gulf of Alaska (GOA) and Shelikof Strait
region, where the bulk of the GOA pollock spawns during
early April. ACC stands for Alaska Coastal Current.
FIG. 2. Time series of the environmental and population variables used in our analysis.TS
5springseasurfacetemperature;TF5fall sea surface temperature;
P
G5
ground®sh predation index.
second only to the Peruvian anchoveta). Most of the
landings come from the Bering Sea, but a large pollock
stock of about 800 000 Mg (averaged over the last 20±
30 yr) is also present in the GOA (Dorn et al. 2002).
The GOA stock has served as a model population for
research, as the early life stages are aggregated in
smaller areas than in those of the Bering Sea popula-
tion. Typically, the bulk of the GOA pollock population
spawns at the beginning of April in Shelikof Strait, a
limited area between Kodiak Island and the Alaska
Peninsula (Fig. 1). After spawning, eggs and larvae are
advected in the Alaska Coastal Current along the Alas-
ka Peninsula. The egg stage lasts approximately two
weeks and the larval stage lasts from 40 to 60 days,
after which larvae metamorphose into the juvenile pe-
lagic form (age-0). Pollock recruit to the mature pop-
ulation and to the ®shery at about age-4. In this study
we focus on the prerecruitment survival (i.e.,
,
age-2)
of the GOA pollock population spawning in Shelikof
Strait.
M
ETHODS
The data
We analyzed the density-dependent structure of pol-
lock survival in the Gulf of Alaska (GOA) during the
period 1975±1996 between the following prerecruit-
ment stages: from eggs to age-0, from age-0 to age-1,
and from age-1 to age-2. The analysis was done in
relation to environmental variables known to have an
effect on pollock early life survival (Megrey et al.
1995, Bailey 2000). These are wind speed cubed (
W;
proportional to water column turbulence), spring (
T
S)
and fall (
T
F) mean sea surface temperature, and ground-
®sh predation intensity (
P
G
; Fig. 2). The time series
from 1975 to 1996 was chosen to insure the quality of
data for the study. Information on
P
Glarval and juvenile
abundance is not available for earlier time periods, and
estimates of recruitment strength are less reliable. In
addition, more recent time series on age-speci®c larval
abundance are not available, nor are recent estimates
of age-0 abundance.
Pollock prerecruitment abundance estimates were
taken from Bailey et al. (1996) and Bailey (2000). In
brief, egg abundance (Fig. 2) was calculated from the
spawning biomass, using fertility and sex ratio esti-
mates. Estimates of age-0 abundance were available
from late summer/autumn shrimp trawl and juvenile

3420
LORENZO CIANNELLI ET AL. Ecology, Vol. 85, No. 12
pollock surveys (described in Bailey and Spring 1992),
with additional data from Wilson et al. (1996). Values
from years with missing data were estimated by inter-
polating across stages and years from a life table anal-
ysis (Bailey et al. 1996). In the total life table, 12 of
112 possible cells were ®lled by this interpolation
method. In years when data were available, estimated
values were compared to observed values, showing that
the method is fairly robust (
r
5
0.97;
P,0.001). For
the age-0 data alone, comparison of observed vs. es-
timated values was also robust (
r50.82;P,0.05).
Estimates of age-1 abundance are from acoustic sur-
veys done in early spring, as described in Bailey
(2000). Since the survey covers only part of the full
distribution range of the age-1 juveniles, it was con-
sidered a relative index of abundance. The estimated
number of age-1 juveniles was proportional to the
abundance index, scaled to fall within the range of age-
0 and age-2 abundance from life tables. Age-2 recruits
were available from an age-structured stock assessment
model tuned statistically to the commercial harvests
and acoustic and bottom trawl surveys (Dorn et al.
2002).
The environmental data were obtained from a variety
of databases (Fig. 2).
T
S
was the mean value over the
April±June earlylifeperiod,and
T
F
was the mean value
between October and November. Both were centered
on 56
8N, 156
8W, and were obtained from reanalysis
of the National Center for Environmental Prediction
data set. The
Wvariable (m
3
/s
3
) was derived from gra-
dient winds and averaged over the month of May (Bai-
ley and Macklin 1994). The ground®sh predation index
(
P
G) is proportional to the amount eaten by the major
predators of age-0 and age-1 pollock, including arrow-
tooth ¯ounder (
Atherestes stomias
), adult pollock, cod
(
Gadus macrocephalus), halibut (Hippoglossus steno-
lepis
), and ¯athead sole (Hippoglossoides elassodon).
P
G
wascalculatedfromamodelaccountingforpredator
biomass,dailyration,overlappingofpredatorandprey,
and percent pollock in the diet (Bailey 2000).
The models
We hypothesize that internal (i.e., demographic) and
external (i.e., environmental) variables have nonaddi-
tive effects on pollock survival. Further, we hypothe-
size that nonadditive interactions, if present, act to
change the level of density-dependent survival. To
study the density-dependent survival of young pollock
across different environmental conditions we used var-
iations of generalized additive model (GAM) formu-
lations, as implemented in the mgcv library of R (Wood
2000, Wood 2001). Speci®cally, let be the natural
a
X
t
logarithm of the abundance of the prerecruitment pop-
ulation stage
a
at time
t.Let be the survival at time
a
St
tbetween two consecutive prerecruitment stages (cal-
culated as the difference between the logarithmic abun-
dances). LetE
t
be a vector of environmental variables
at time
t
whose single components are identi®ed by the
superscript
j; letfand
g
jbe nonparametric, smoothing
functions, respectively, specifying the effect of popu-
lation abundance (i.e., density dependence) and envi-
ronmental forcing of the covariateEon the demo-
j
graphic variableS. By assuming no changes of density
dependence across the entire modeled phase (i.e., ad-
ditivity), the formulation becomes:
j
aa
S5
b
1f(X)
1g
(E)
1«
(1)Ottjtt
j
wherebis an intercept and«is a noise term; for the
model identi®cation, all additive function components
are henceforth assumed to be centered, i.e., having zero
mean over the data.
In contrast, if there is a change of density dependence
across contrasting environmental phases (i.e., nonad-
ditivity) the model assumes either a continuous (if
changes of density dependence are gradual) or a thresh-
old (if changes are abrupt) formulation. In the contin-
uous case, Model 1 is enriched by the inclusion of
second order interaction terms, whose effect is regu-
lated by the nonparametric functions
s(Hastie and Tib-
shirani 1990). Speci®cally,
aajak ÃÃS5b1f(X)1g(E)1s(XE)1«. (2)OOttjtkttt
jk
We refer to this model as CGAM, where C stands for
``continuous interaction.'' The variables and are
akà ÃX E
tt
the standardized and shifted (a constant is added in
order to make the variable
$0) population and envi-
ronmental covariates, respectively. The superscript
k,
uniquely identi®es the interacting environmental co-
variates among all those present in the vectorE.
The threshold nonadditive formulation (hereon re-
ferred to as TGAM, where T stands for ``threshold
interaction'') is composed by two additive formula-
tions. In particular, assuming a change of density de-
pendence (i.e., nonadditivity) as a function of a linear
combination of the environmental vector (
aE; where
ais a row vector of coef®cients) the effect of internal
and external variables on pollock survival can be mod-
eled as follows:
ajìb1f(X)1g(E)1«ifaE#rO11 tjtt
ï j
a
(3)S5ít
aj
ïb1f(X)1g(E)1«otherwiseO22 tjtt
î j
whereris an environmental threshold across which the
density-dependent function switches from
f
1tof
2, with
possible changes in the intercept as well (from
b
1to
b
2), with the additive environmental effects otherwise
unaltered.
In the case of Model 3, if the environmental covariate
is two dimensional,
aE#randaE.rde®ne two
half spaces between a set of two external variables,
separated by a threshold straight line,
aE5r,tobe
estimated from the data. The estimation procedure of
the threshold line for Model 3 consists of the following

December 2004 3421 NONADDITIVE ENVIRONMENTAL CONTROL
steps. Each pair of distinct observed covariate vectors
implicitly de®nes a threshold line
aE
5r. For a ®xed
threshold line, Model 3 may be written as a GAM with
known multiplicative factors, as follows:
aa
S5
bI(
aE
#
r)1f(X)I(aE
#
r)1bI(aE
.
r)
t 1 t 1tt 2 t
aj
1f(X)I(aE
.r)
1
g
(E)
1«
(4)O2tt jtt
j
where the operator
Idesignates indicator variables (1,
0). Model 4 is an extended GAM and can be estimated
by the method of penalized least squares (Wahba 1990,
Green and Silverman 1994; see also Appendix A for
explanation of penalized least squares), with the cor-
responding generalized cross validation (GCV) denot-
ed by GCV(
a
,
r). We estimate
a
andrby minimizing
GCV(
a,r) over all lines de®ned by pairs of observed
covariate vectors, except those lines whose associated
half spaces contains
,20% of the data cases; this is
done to ensure adequate data for estimation. The min-
imum GCV(
a,r) is then de®ned as the GCV of the
estimated TGAM.
However, the GCV so de®ned has not properly ac-
counted for the parameters de®ning the threshold line;
hence it cannot be directly compared with that of the
best ®tted additive environmental-effect Model 1 or
nonadditive continuous effect Model 2. To handle this
problem, we compute the ``genuine'' cross-validatory
squared prediction error (genuine CV), as follows.
Each data case is deleted one at a time with the response
of the deleted case predicted by the examined model
(GAM, CGAM, or TGAM) ®tted to all the remaining
cases. A square prediction error is then calculated, and
the same routine is repeated across all data cases. The
mean of all squared prediction error is then the (gen-
uine) cross-validatory squared prediction error for the
examined model. While this procedure is rather com-
putationally intensive, it properly accounts for the es-
timation of the threshold line and the estimation of the
degrees of freedom for the functions appearing in all
additive and nonadditive formulations (see Efron and
Gong 1982, Li 1986, 1987).
Our model building strategy consisted of four steps.
First, we used a TGAM formulation to assess, based
on the model GCV, the effect (i.e., additive or non-
additive) of the environmental covariates under scru-
tiny. We only examined covariates that from previous
studies were known to in¯uence pollock recruitment.
In the survival from egg to age-0 these were
W
and
T
S,
and in the survival from age-0 to age-1 (and from age-
1 to age-2) the environmental covariates were
P
Gand
T
F
.
TF
and
P
Gwere not included in the analysis of egg
to age-0 survival because the ®rst refers to a time of
the year (October±November) that does not include the
passage from the egg (April) to the age-0 (August)
stage, while the second refers to the predation on stages
older than the age-0. Likewise,
T
Sand
W
were not
considered in the survival analysis of the juvenile stag-
es because we assume that most of the juvenile pollock
mortality, which is not accounted for by predation, oc-
curs during fall and winter, a view largely shared in
recruitment studies of marine ®sh populations residing
in a seasonal environment (Schultz et al. 1998). The
second step of the model building strategy was to in-
clude the same covariates used in the TGAM (both
additive and interacting) in an additive formulation
(GAM). The third step was to add a continuous inter-
action term to the GAM, to obtain a CGAM formu-
lation. We only examined second order interaction
terms (third order interactions are not properly ad-
dressed using the adopted scheme [Hastie and Tib-
shirani 1990]). One of the two covariates of the CGAM
interaction term was always initial population size (
X
t),
and the other was one of the environmental variables
found to signi®cantly interact in the TGAM formula-
tion. Because the TGAM formulation always had two
interacting covariates (besides
X), there were two sec-
ond order interaction terms that could potentially be
added to the CGAM formulation. Initially, we included
all interacting terms; however, in the ®nal formulation,
we only retained those that minimized the model GCV.
The fourth and ®nal step of the model building strategy
was to compare all inspected formulations using the
genuine CV criterion.
To further assess the signi®cance of the additive
(Model 1) vs. the nonadditive (Models 2 and 3) model
formulations we simulated the dynamics of pollock
prerecruitment stages using the estimated survival rates
from all ®nal models. This was done for the calibration
period (1975±1996), and for a set of three independent
data points (1972±1974) for which we only had esti-
mates of age-2 abundance and environmental variables.
In the simulation, each model formulation was used to
®rst predict age-0 from the observed egg abundance
and environmental variables, then age-1 from the pre-
dicted age-0 (plus error), and ®nally age-2 from the
predicted age-1 (plus error). All estimated abundances
(age-0, age-1, and age-2) were the means of 1000 sim-
ulations, each resulting from the addition of an error
term (randomly sampled from the model residuals) to
the model estimate.
R
ESULTS
The survival from egg to age-0 was marginally best
modeled by the nonadditive continuous (CGAM) for-
mulation, followed by the threshold (TGAM) formu-
lation, and lastly by the fully additive (GAM) formu-
lations (Table 1). The survival from age-0 to age-1 was
best modeled by the TGAM formulation, followed by
the GAM, and ®nally by the CGAM formulations. In
the egg to age-0 and the age-0 to age-1 survival phases,
the interaction of the CGAM was best modeled (based
on the GCV criteria) by a linear term (Table 1). In
contrast to earlier stages, initial density in the survival
from age-1 to age-2 did not enter as a signi®cant factor
in an additive model (
P50.11). Based on this result

3422
LORENZO CIANNELLI ET AL. Ecology, Vol. 85, No. 12
T
ABLE
1. Final formulations and genuine cross validation (CV) of all inspection models in the analysis of pollock pre-
recruitment survival.
Model type Stages Formulation CV
GAM (Model 1) egg to age-0
age-0 to age-1
age-1 to age-2
S
5 b
1
f (
X
)
1 g1
( T
S,t
)
1g
2
(
W
t
)
1«
t
aa
tt
S5
b
1
f(
X
)
1
g
1(T
F,t
)1
g
2
(
P
G,t)1«
t
aa
tt
S
5b
1
g
1
(
T
F,t)
1
g
2(P
G,t)
1«
t
a
t
0.728
0.627
0.279
CGAM (Model 2) egg to age-0
age-0 to age-1
age-1 to age-2
S
5b
1
f(
X
)1g
1(T
S,t)1g
2(
W
t)
1a
TÃ
S,tXÃ
1«
t
aa a
tt t
S5
b1f
(X)
1
g
1(T
F,t
)
1g
2(P
G,t)1a
TÃ
F,tXÃ1«
t
aa a
tt t
NA
0.602
0.754
NA
TGAM (Model 3) egg to age-0
age-0 to age-1
age-1 to age-2S5
b
1
1f
1(X)
1
g
1
(W
t
)
1«
t
ifa
{
T
S,t
,
W,
t}#r
aa
tt
S5b
21
f
2(
X)1
g
1
(W
t
)1«
tif
a
{T
S,t
,
W,
t}.r
aa
tt
S5
b
11f
1
(X
)
1g
1(T
F,t)1«
t
if
a{T
F,t,P
G,t}#r
aa
tt
S5b
21
f2
(X)1g
1(T
F,t)1«
tifa{T
F,t,P
G,t
}
.r
aa
tt
NA
0.617
0.329
NA
Notes:
Initial population size was not found to signi®cantly affect the survival from age-1 to age-2; therefore, such passage
was only inspected with an additive model formulation.
T
S5spring sea surface temperature;
T
F5fall sea surface temperature;
W
5wind speed cubed;P
G5ground®sh predation index on juvenile pollock;a5linear coef®cient;r5environmental
threshold. See Eqs. 1±3 (
Methods: The models
) for explanation of additional model parameters.
FIG. 3. Scatter plot of the interacting en-
vironmental covariate affecting the density-de-
pendent survivalofpollock, as estimated from
the TGAM formulation. Also shown are the
threshold lines (detected by the GCV criteria)
that divide the plot of environmental variables
in two regions, one above and one below the
line. Regions above the threshold line corre-
spond to environmental regimes with negative
density dependence, while no signi®cant den-
sity-dependent structure was detected in regions
below the threshold line. Plots are shown for
covariates affecting egg to age-0 survival (
W5
wind speed cubed;
T
S
5
spring sea surface tem-
perature) and age-0 to age-1 (
T
F5fall sea sur-
face temperature;
P
G5ground®sh predation in-
dex) survival.
we inferred that there was no evidence of density-de-
pendent survival from age-1 to age-2. Therefore, no
conclusions could be drawn on age-1 survival, regard-
ing nonadditivity of environmental and demographic
variables.
The TGAM formulation indicated that during the egg
to age-0 stage, high
W, and to a lesser extent lowT
S,
induced negative density dependence (i.e., compensa-
tion). Under the inverse conditions there was no no-
ticeable effect of density on survival, suggesting an
overall relaxation of density-dependent effects (Fig. 3).
High
P
GandT
Finduced compensation in pollock sur-
vival from the age-0 to the age-1 stage, while no par-
ticular density-dependent structure was detected during
the opposite environmental conditions (Fig. 3).
The resulting qualitative picture obtained for the
CGAM formulation was similar to that of the TGAM
(Fig. 4). One exception was that water temperature,
rather than wind speed, was the most important inter-
acting covariate in the survival from egg to age-0 (Ta-
ble 1, Fig. 4). However, in our study region,
WandT
S
are to some degree negatively correlated (R
2
50.28;
high wind induces deeper water mixing and cooling of
surface temperature; Fig. 3), so their ®nal effect on egg
survival is also correlated. Survival from age-1 to age-
2 was not signi®cantly affected by the initial population
density, while it was affected by
P
GandT
F,inanad-
ditive fashion (Fig. 4).
The nonadditive threshold formulation (Model 3)
was always found to best reproduce the observed abun-
dance of pollock prerecruitment stages from initial con-
ditions (i.e., egg density and environmental covariates;
Fig. 5). At the recruitment level of age-2, for example,
the TGAM captured 69% of the observed variability,
while the CGAM captured 64% and the GAM captured
59%. In addition, during the independent period 1972±
1974, the TGAM closely predicted two out of three
age-2 abundance estimates, while both CGAM and GAM
grossly over predicted all three data points (Fig. 5).
From visual inspection of normal probability plots
of the model residuals, we did not observe any major
departure from normality, given the constraint of the
small sample size (see Appendix B for residual anal-
ysis). We also noticed the presence of a potential outlier
in the CGAM formulation of the egg to age-0 survival
(Year 1991) and in the TGAM formulation of the age-
0 to age-1 survival (Year 1982). However, the removal
of the offending case, or the addition of an outlier effect

December 2004 3423
NONADDITIVE ENVIRONMENTAL CONTROL
FIG. 4. Surface plots of survival for each of the inspected pollock prerecruitment stages and model formulations. For
egg to age-0 survival and age-0 to age-1 survival (®rst and second row, respectively), the three-dimensional plots only show
the effect of the variables that were signi®cantly interacting in a CGAM formulation (see Table 1), while all covariates
affecting age-1 survival (third row) were used in the surface plot (see Table 1). The three-dimensional plots were derived
by predicting survival from a regular grid in which the plotted variables ranged from their minimum to maximum values,
and the nonplotted variables were held constant as their average values. The only exception was for
W
in the egg survival
plots (®rst row), which was varied as a function of
T
S(W
andT
Sare negatively correlated). The environmental regimes for
the contrasting phases of the TGAM formulations are identi®ed in Fig. 3.
(indicator variable) to the model formulation, further
normalized the residual plot, and did not change any
of the model conclusions (see Appendix B). Residuals
from a GAM formulation that only contained the de-
mographic variable (i.e., population size), against the
distance from the threshold line of the corresponding
TGAM formulation (shown in Fig. 3), revealed the
presence of a linear trend. This pattern was observed
for both the egg to age-0 and the age-0 to age-1 sur-
vivals (see residual analysis in Appendix B). The pres-
ence of linear correlation between the GAM residuals
and the interacting covariate (i.e., distance from the
threshold line of Fig. 3) further supports the claim of
nonadditivity.
D
ISCUSSION
We have demonstrated that external environmental
variables can affect pollock survival in a nonadditive
and interactive fashion, and that the inclusion of these
nonadditive effects in a population dynamics model
improves our ability to simulate recruitment. We fur-
ther demonstrate that nonadditive interactions are the
result of changes in density-dependent survival in re-
lation to background environmental variables.
Previous studies have simulated pollock recruitment
in the Gulf of Alaska and eastern Bering Sea (Megrey
et al. 1995, Quinn and Niebauer 1995, respectively)
with linear regression models, solely based on abiotic
variables. These earlier analyses have been very rele-
vant in isolating key variables that affect pollock re-
cruitment (e.g., air temperature and ice cover [Quinn
and Neibauer 1995], and precipitation and wind mixing
[Megrey et al. 1995]). However, due to the direct link
between the examined covariates and the recruitment
level, these studies could not provide extensive expla-
nations for the mechanisms involved in the regulation

3424
LORENZO CIANNELLI ET AL. Ecology, Vol. 85, No. 12
F
IG. 5. Observed and simulated abundances of pollock
prerecruitment stages during the calibration (clear areas) and
validation period (shaded areas). Correlations for the ln(age-
0) simulations (
R
2
values) are: GAM, 0.43; CGAM, 0.58; and
TGAM, 0.65. For the ln(age-1) simulations,
R
2
values are:
GAM, 0.36; CGAM, 0.44; and TGAM, 0.51. For the ln(age-
2) simulations,
R2values are: GAM, 0.59; CGAM, 0.64; and
TGAM, 0.69. Only the values included in the calibration pe-
riod were used to calculate
R
2
.
of pollock survival. In our analysis, we considered den-
sity-dependent forces and biotic factors (e.g., predation
index) in addition to abiotic variables, and compared
a variety of contrasting survival dynamics (i.e., addi-
tive vs. nonadditive, continuous vs. threshold). More-
over, we focused on pollock stage-speci®c survival
(i.e., from eggs to age-0, from age-0 to age-1, and from
age-1 to age-2) over a longer time series than those
available to earlier studies. Consequently, we also pro-
vide a more comprehensive treatment of the mecha-
nismsinvolved in regulating pollock prerecruitment
survival.
The threshold nonadditive statistical formulation
(TGAM) that we developed can be regarded as a non-
parametric generalization of the Threshold Autore-
gressive Models(Tong 1990). This linkage may be
helpful for investigating the probabilistic properties of
the TGAM and large sample properties of the proposed
estimation scheme. A continuous (CGAM) model for-
mulation was marginally superior (in terms of genuine
cross validation) to the TGAM only from egg to age-
0, while survival from age-0 to age-1 was best modeled
with a TGAM. However, the CGAM was not ef®cient
in simulating three independent recruitment points
(1972±1974), while the TGAM closely predicted two
out of three points (Fig. 5). The out-of-sample forecast
re¯ects the general problem of modeling interactions
by product terms. A product term is of second degree
and, as such, provides adequate within-sample forecast
but may give a worse prediction upon extrapolation
due to the stronger curvature of second degree terms.
In general, piecewise linear models such as the TGAM
can be expected to perform better in extrapolation.
Hence, on the basis of model selection criteria (genuine
CV) and out-of-sample forecast performance, we prefer
the TGAM to the continuous interaction model for both
egg to age-0 and the age-0 to age-1 stages. In addition,
as we show in our analysis, the TGAM methodology
can be easily adapted to explore interactions of order
.
2(e.g.,P
G3T
F3X
a
), not easily addressed by CGAM
formulations (Hastie and Tibshirani 1990).
The observed changes of pollock density-dependent
survival are in agreement with what is expected from
previous studies. Increased negative density depen-
dence of egg and larval survival during conditions of
high wind and low temperature may be partly induced
by water column turbulence. For pollock it has been
shown that high turbulence may reduce the foraging
success of young larvae (Bailey and Macklin 1994).
Strong winds can also have a direct effect on survival
by increasing transport and advection of larvae and
juveniles into offshore and unfavorable nursery areas
(Bailey and Macklin 1994). Similar negative effects of
high wind on recruitment were also observed in small
pelagic ®shes (e.g., sardine and anchovies) residing in
Ekman-type upwelling systems (Cury and Roy 1989).
The effect of
TS
on egg and larval survival is not as
straightforward to interpret since temperature may af-
fect survival across multiple substages both directly
(e.g., through growth) and indirectly (e.g., through food
availability and competition with other species). How-
ever, in our study region,
WandT
Sare to some degree
negatively correlated, which explains their opposite
outcome on the regulation of pollock density depen-
dence.
From age-0 to age-1, several mechanisms may ex-
plain the observed increase in compensation with in-
creased predation and water temperature. One possi-
bility is an increase of interspeci®c competition for
food. In fact, ground®sh predators, besides feeding on
young pollock, also consume a number of zooplankton
prey common in the diet of juvenile pollock (e.g., eu-
phausiids and large copepods; see Brodeur and Wilson
1996 and references therein). In addition, interspeci®c

December 2004 3425 NONADDITIVE ENVIRONMENTAL CONTROL
competition for food may increase during high water
temperature regimes, due to increased metabolic de-
mands of young pollock (Ciannelli et al. 1998). An-
other possible explanation is habitat reduction, as warm
sea surface temperatures force juvenile pollock to re-
side in a restricted region of the water column (Brodeur
and Wilson 1996). The absence of a signi®cant effect
of initial density on the survival of pollock between
age-1 and age-2 may also be expected, due to a change
of ecological niche between the two stages (Duffy-
Anderson et al. 2003). However, care should be placed
in interpreting this result, given the narrow margin by
which the initial density term failed to signi®cantly
differ from zero (
P
50.11, rejection at
P,
0.05), and
the uncertainty of the demographic variables used in
our study.
In evaluating the results of our analysis it bears not-
ing that a great deal of uncertainty is inherent in each
of the data sources, such as sampling errors or esti-
mation of missing values (Bailey et al. 1996, Bailey
2000). Also, the application of the nonadditive models
to the data in hand, while extremely interesting, is also
limited by the rather brief time series. For example, we
found that in the CGAM formulations, linear interac-
tion terms were favored against nonlinear ones; how-
ever, nonlinear relationships are not easily identi®able
under the constraint of a small sample size.
In ®sheriesscience,ithaslongbeenrecognizedthat
prerecruitment dynamics are overwhelmingly impor-
tant in shaping the abundance of the entire population
(Hjort 1914). Recruitment of marine ®shes is a com-
plexprocess,beingtheemergentpropertyofnonlinear
interactions among external and demographic vari-
ables,aswellaslowfrequencyphysicalandbiological
events(Baileyetal.2003).Insuchcontext,ourresults
offer a template whereby the joint effect of demo-
graphic and external variables may be better under-
stood. This level of comprehension is particularlyrel-
evant in studying the outcome of climate variability
on the population ecology of marine species. Our
study may, in fact, provide insights in understanding
theeffectofthemid1970sclimatechangeontheGOA
pollock dynamics. Following the climate regime shift
in the North Paci®c, there was an increase in both
watertemperature(Mantuaetal.1997)andground®sh
biomass (Anderson and Piatt 1999) in the Gulf of
Alaska. Bailey (2000) showed that in concomitance
with such changes, recruitment control of pollock
shifted from the larval to the juvenile stage. We have
demonstrated that both
T
F
andP
Gtrigger changes in
pollock density dependence and notably affect sur-
vival. Therefore, it is possible that the shift in re-
cruitment control was the result of increased density
dependence during the juvenile stage.
A mechanistic understanding of the interaction be-
tween external and internal population variables is rel-
evant to the management of renewable marine resourc-
es in the face of underlying environmental variability.
Density dependence is, in fact, widely assumed in many
®sheries management models (Hilborn and Walters
1992). However, it also is one of the most elusive pro-
cess to detect in marine populations (Rose et al. 2001).
Failure to detect density dependence can arise in part
because, as we demonstrate, its intensity can change
over temporally distinct environmental phases. Hence,
analyses that integrate over time, without allowing for
explicit changes of density-dependent survival, will in-
evitably integrate across periods of contrasting popu-
lation phases.
The statistical methodology developed here can be
easily applied to other marine and terrestrial popula-
tions. In terrestrial ecology, for example, many studies
point to the importance of understanding how internal
and external variables jointly affect population abun-
dance (Leirs et al. 1997, Lima et al. 2001, Stenseth et
al. 2003). Previous efforts to model joint effects of
these two forcing mechanisms include linear (Barbraud
and Weimerskirch 2003), and a priori speci®ed non-
linear formulations (e.g., parametric threshold models;
Tong 1990, Stenseth et al. 1998). However, these ap-
proaches fall short in the presence of nonlinear dynam-
ics or poorly understood functional relationships. In
such a context, the phase-dependent and nonparametric
nature of the TGAM and CGAM methodologies will
advance our ability to address internal and external
mechanisms of population control.
A
CKNOWLEDGMENTS
We thank Geir Ottersen and Kyrre Lekve for technical
assistance during the early stage of this study. Comments
from Meta Landys, Mauricio Lima, two anonymous re-
viewers, and the subject-matter editor improved earlier
versions of this paper. Karna McKinney helped with the
®gures. LC was funded by the NOAA Steller Sea Lion
Research and the Pollock Conservation Cooperative Re-
search Center. KSC was partially supported by the Uni-
versity of Iowa BSI program. NCS acknowledged support
from the Norwegian Science Council to the EcoClim-proj-
ect as well as from the University of Oslo. This research is
contribution FOCI-0492 to NOAA's Fisheries±Oceanography
Coordinated Investigations.
L
ITERATURECITED
Anderson, P. J., and J. F. Piatt. 1999. Community reorgani-
zation in the Gulf of Alaska following ocean climate regime
shift. Marine Ecology Progress Series189:117±123.
Bailey, K. M. 2000. Shifting control of recruitment of walleye
pollock
Theragra chalcogrammaafter a major climatic and
ecosystem change. Marine Ecology Progress Series198:
215±224.
Bailey, K. M., R. D. Brodeur, and A. B. Hollow ed. 1996.
Cohort survival patterns of walleye pollock,
Theragra chal-
cogramma
, in the Shelikof Strait, Alaska: a critical factor
analysis. Fisheries Oceanography5(S1):179±188.
Bailey, K. M., L. Ciannelli, and V. N. Agostini. 2003. Com-
plexity and constraints combined in hybrid models of re-
cruitment. Pages 293±301
inH. I. Browman and A. Skif-
tesvik, editors. The big ®sh bang, proceedings of the 26
th
annual larval ®sh conference. Institute of Marine Research,
Bergen, Norway.
Bailey, K. M., and A. S. Macklin. 1994. Analysis of patterns
in larval walleye pollock
Theragra chalcogrammasurvival

3426
LORENZO CIANNELLI ET AL. Ecology, Vol. 85, No. 12
and wind mixing events in Shelikof Strait, Gulf of Alaska.
Marine Ecology Progress Series113:1±12.
Bailey, K. M., and S. M. Spring. 1992. Comparison of lar-
val, age-0 juvenile and age-2 recruit abundance indexes
of walleye pollock,
Theragra chalcogramma, in the west-
ern Gulf of Alaska. ICES Journal of Marine Science49:
297±304.
Barbraud, C., and H. Weimerskirch. 2003. Climate and den-
sity shape population dynamics of a marine top predator.
Proceedings of the Royal Society of London B270:2111±
2116.
Brodeur, R. D., and M. T. Wilson. 1996. A review of the
distribution, ecology and population dynamics of age-0
walleye pollock in the Gulf of Alaska. Fisheries Ocean-
ography5(S1):148±166.
Chavez, F. P., J. Ryan, S. E. Lluch-Cota, and M. C. NÄiquen.
2003. From anchovies to sardines and back: multidecadal
change in the Paci®c Ocean. Science299:217±222.
Ciannelli, L., R. D. Brodeur, and T. W. Buckley. 1998. De-
velopment and application of a bioenergetics model for
juvenile walleye pollock. Journal of Fish Biology52:879±
898.
Cury, P., and C. Roy. 1989. Optimal environmental window
and pelagic ®sh recruitment success in upwelling areas.
Canadian JournalofFisheriesandAquaticScience46:670±
680.
Cury, P., C. Roy, R. Mendelssohn, A. Bakun, D. M. Husby,
and R. H. Parrish. 1995. Moderate is better: exploringnon-
linear climatic effects on the Californian northern anchovy
(
Engraulis mordax). Pages 417±424inR. J. Beamish, ed-
itor. Climate-change and northern ®sh populations. Cana-
dian Special Publication of Fisheries and Aquatic Sciences
121.
Dennis, B., and M. R. M. Otten. 2000. Joint effects of the
density dependence and rainfall on abundance of San Joa-
quin kit fox. Journal of Wildlife and Management64:388±
400.
Dorn, M., S. Barbeaux, M. Guttormsen, B. Mergrey, A. Hol-
lowed, E. Brown, and K. Splalinger. 2002. Assessment of
walleye pollock in the Gulf of Alaska.
InStock assessment
and ®shery evaluation report for the ground®sh resources
of the Gulf of Alaska. North Paci®c Fishery Management
Council, Anchorage, Alaska, USA.
Duffy-Anderson, J. T., L. Ciannelli, T. Honkalehto, K. M.
Bailey, S. Sogard, A. Springer, and T. Buckley. 2003.
Distribution of age-1 and age-2 walleye pollock in the
Gulf of Alaska and eastern Bering Sea: sources of var-
iation and implications for higher trophic levels. Pages
381±394
inH. I. Browman and A. Skiftesvik, editors.
The big ®sh bang, proceedings of the 26th annual larval
®sh conference. Institute of Marine Research, Bergen,
Norway.
Efron, B., and G. Gong. 1982. A leisurely look at the boot-
strap, the jackknife and the cross-validation. American
Statistician37:36±48.
FAO (Food and Agricultural Organization). 2002. The state
of world ®sheries and aquaculture. FAO, Rome,Italy.
^http:
//www.fao.org
&
Green, P. J., and B. W. Silverman. 1994. Nonparametric
regression and generalized linear models: a roughness
penalty approach. Chapman and Hall, New York, New
York, USA.
Grenfell, B. T., O. F. Price, S. D. Albon, and T. H. Clutton-
Brock. 1992. Overcompensation and population cycles in
an ungulate. Nature355:823±826.
Hastie, T. J., and R. J. Tibshirani. 1990. Generalized additive
models. Chapman and Hall, New York, New York, USA.
Hilborn, R., and C. J. Walters. 1992. Quantitative ®sheries
stock assessment: choice, dynamics and uncertainty. Chap-
man and Hall, New York, New York, USA.
Hixon, M. A., and M. H. Carr. 1997. Synergistic predation,
density dependence, and population regulation in marine
®shes. Science277:946±949.
Hjort, J. 1914. Fluctuations in the great ®sheries of Northern
Europe. Viewed in the light of biological research. Rapports
et ProceÁs-Verbaux des ReÂunions, Conseil International pour
l'Exploration de la Mer20:1±228.
Leirs, H., N. C. Stenseth, J. D. Nichols, J. E. Hines, R. Ver-
hagen, and W. Verheyen. 1997. Stochastic seasonality and
nonlinear density-dependent factors regulate population
size in an African rodent. Nature389:176±180.
Li, K. C. 1986. Asymptotic optimality of
C
Land generalized
cross validation in the ridge regression with application to
spline smoothing. Annals of Statistics14:1101±1112.
Li, K. C. 1987. Asymptotic optimality for
C
p,C
L, cross-
validation and generalized cross-validation: discrete index
set. Annals of Statistics15:958±975.
Lima, M., N. C. Stenseth, N. G. Yoccoz, and F. M. Jaksic.
2001. Demography and population dynamics of the mouse
opossum (
Thylamys elegans) in semi-arid Chile: season-
ality, feedback structure and climate. Proceedings of the
Royal Society of London B268:2053±2064.
Mantua, N. J., S. R. Hare, Y. Zhang, J. M. Wallace, and R.
C. Francis. 1997. A Paci®c interdecadal climate oscillation
with impacts on salmon production. Bulletin of the Amer-
ican Meteorological Society78:1069±1079.
Megrey, B. A., S. J. Bograd, W. J. Rugen, A. B. Hollowed,
P. Stabeno, S. A. Macklin, J. D. Schumacher, and W. J.
Ingraham. 1995. An exploratory analysis of associations
between biotic and abiotic factors and year-class strength
of Gulf of Alaska walleye pollock. Pages 227±243
inR. J.
Beamish, editor. Climate-change and northern ®sh popu-
lations. Canadian Special Publication of Fisheries and
Aquatic Sciences121.
Quinn, T. J., II, and H. J. Niebauer. 1995. Relation of
eastern Bering Sea walleye pollock (
Theragra chalco-
gramma
) recruitment to environmental and oceano-
graphic variables. Pages 497±507
inR. J. Beamish, ed-
itor. Climate-change and northern ®sh populations. Ca-
nadian Special Publication of Fisheries and Aquatic Sci-
ences121.
Rose, K. A., J. H. Cowan, K. O. Winemiller, R. M. Myers,
and R. Hilborn. 2001. Compensatory density dependence
in ®sh populations: importance, controversy, understanding
and prognosis. Fish and Fisheries2:293±327.
Schultz, E. T., D. O. Conover, and A. Ehtisham. 1998. The
dead of winter: size-dependent variation and genetic dif-
ferences in seasonal mortality among Atlantic silverside
(Atherinidae:
Menidia menidia) from different latitudes.
Canadian Journal of Fisheries and Aquatic Sciences55:
1149±1157.
Shima, J. S., and C. W. Osenberg. 2003. Cryptic density
dependence: effects of covariation between density and site
quality in reef ®shes. Ecology84:46±52.
Stenseth, N. C. 1999. Population cycles in voles and lem-
mings: density dependence and phase dependence in a sto-
chastic world. Oikos87:427±461.
Stenseth, N. C., O. N. Bjùrnstad, W. Falck, J. M. Fromentin,
J. Gjùsaeter, and J. S. Gray. 1999. Dynamics of coastal
cod populations: intra- and intercohort density dependence
and stochastic processes. Proceedings of the Royal Society
of London B266:1645±1654.
Stenseth, N. C., K.-S. Chan, E. Framstad, and H. Tong. 1998.
Phase-dependent population dynamics in Norwegian lem-
mings: interaction between deterministic and stochastic
processes. Proceedings of the Royal Society of London B
265:1957±1968.
Stenseth, N. C., A. Mysterrud, G. Ottersen, J. W. Hurrel,
K.-S. Chan, and M. Lima. 2002. Ecological effects of cli-
mate ¯uctuations. Science297:1292±1296.

December 2004 3427
NONADDITIVE ENVIRONMENTAL CONTROL
Stenseth, N. C., G. Ottersen, J. W. Hurrell, A. Mysterud, M.
Lima, K.-S. Chan, N. G. Yoccoz, and B. AÊdlandsvik. 2003.
Studying climate effects on ecology through the use of
climate indices: the North Atlantic Oscillation, El NinÄo
Southern Oscillation and beyond. Proceedings of the Royal
Society of London B270:2087±2096.
Tong, H. 1990. Non-linear time series: a dynamical system
approach. Oxford University Press, Oxford, UK.
Wahba, G. 1990. Spline models for observational data,
SIAM, Philadelphia. CBMS-NSF Regional Conference Se-
ries in Applied Mathematics59.
Wilson, M. T., R. D. Brodeur, and S. Hinckley. 1996. Dis-
tribution and abundance of age-0 walleye pollock,
Thera-
gra chalcogramma
, in the western Gulf of Alaska during
September 1990. NOAA Technical Report NMFS126:
11±24.
Wood, S. J. R. 2000. Modelling and smoothing parameter
estimation with multiple quadratic penalties. Journal of the
Royal Statistical Society B62(2):413±428.
Wood, S. J. R. 2001. mgcv: GAMs and Generalized Ridge
Regression for R. R News1(2):20±25.
APPENDIX A
A description of the method of penalized least squares is available in ESA's Electronic Data Archive:
Ecological Archives
E085-119-A1.
APPENDIX B
A description of residual analysis and ®gures showing normal probability plots and correlation between residuals
and distance from the threshold line of Fig. 3 are available in ESA's Electronic Data Archive:
Ecological Archives
E085-119-A2.