26 
Fishery Bulletin 110(1) 
plex equation (having the most parameters) had three 
parameters and represented the response of LCPUE 
as a symmetrical dome-shaped function of the habitat 
variables, so that 
LCPUE h = a h + P h X h + <j) h Xf r (2) 
Here, X h = habitat variable h\ and a h , P h , and p h are 
parameters fitted to the data. 
The second equation describes LCPUE as an exponen- 
tial function of the habitat variables, so that 
LCPUE h - a h X li e~ b ‘ ,x '' . (3) 
In this case, only two parameters, a h and b h , are fitted. 
With the simplest equation (least parameters), the pre- 
dicted rockfish abundance was computed as proportional 
to the habitat variables X h , so that 
LCPUE h -c h X h , (4) 
where c h = the only parameter fitted in the equation. 
All components of LCPUE were combined before fit- 
ting the parameters. For example, the initial (full) mod- 
el for the shortspine thornyhead analyses estimated 18 
parameters plus one dummy parameter for each year, 
by using the 3-parameter equation ( Eq. 2) for each of 
the six variables so that 
LCPUE — (Xj-f + P^Xq + PdXq + (Xj, + PpXp + ppX p + 
& TD PtD^TD pTD^TD r Ps^S Ps^S T ® CS ^ ^ 
P CS X CS + p cs X^.g\- oc p + PpXp + (ppXp + Y + £. 
All 18 parameters, plus the year effects, were fitted 
simultaneously. 
Model parameters were estimated by minimizing the 
negative log-likelihood (Hilborn and Mangel, 1997), 
by using either a normal or a gamma distribution de- 
pendent on the characteristics of the model residuals. 
For Pacific ocean perch, shortspine thornyhead, short- 
raker rockfish, sharpchin rockfish, and rougheye and 
blackspotted rockfish, species for which there were no 
major departures from normality, the normal distribu- 
tion was used. In the cases of dusky rockfish, harlequin 
rockfish, and northern rockfish, a gamma distribution 
with shape = 0.5 and scale=l was used to fit the model 
parameters. 
Models were reduced by sequentially removing one 
parameter for a variable (e.g., the depth relationship 
was changed from Eq. 2 to Eq. 3) and parameters were 
refitted. Next another parameter for that variable was 
removed (e.g., the depth relationship was changed from 
Eq. 3 to Eq. 4) and the model was refitted. This was 
repeated until the variable was no longer included in 
the equation (all parameters were removed). Then these 
steps were repeated for the next five variables. The 
models were compared by using the Akaike information 
criterion (AIC) for non = nested models to determine the 
best fitting model: 
AlC = 2L + 2p, ( 8 ) 
where p = the number of parameters in the model: and 
L = the sum of the negative log-likelihood of the 
model given the data (Akaike, 1992). 
The best of this series of models was then evaluated 
against the full 18-parameter model. If there was a re- 
duction in the AIC from the full model, the new model 
was kept as the best fitting model. Further parameter 
reductions were implemented by repeating the steps 
above and comparing the results to the best-fitting 
model from the previous series, until the reduction in 
the number of parameters or elimination of variables 
resulted in no reduction in the AIC score, and this final 
model was deemed best for the data set analyzed. The 
r 2 (squared correlation coefficient) between the observed 
and predicted values was used to determine the per- 
centage of variance in the LCPUE data set explained 
by the model. 
To examine potential spatial correlation that was un- 
related to habitat variables we analyzed the residuals 
for the best fitting model. The residuals from the best- 
fitting model were kriged across the geographic area of 
the survey by fitting a generalized least squares trend 
surface (Venables and Ripley, 2002) and assuming a 
spherical covariance function for a range of distances 
from 0.01° to 1° of latitude and longitude. To predict the 
values at each bottom trawl survey station, the trend 
surface value at each station was added back to the 
model prediction (Hengl et al., 2007). The correlation 
between the observed and predicted values (plus the 
kriged surface value) at the range of distances was used 
to determine the scale of spatial correlation. The values 
from the best-fitting surface (with the highest corre- 
lation coefficient) were then compared to the results 
without the addition of the trend surface to determine 
whether a significant amount of residual variance could 
be explained by spatial autocorrelation in the residuals. 
Abundance indices 
The annual abundance index was estimated by the 
dummy variable for the year effect (Y). Because this was 
the variable of interest in the modeling (i.e., producing 
a time series of annual abundance from the model), 
it was included in all of the models and not tested for 
significance with the AIC method above. Errors for all 
parameter estimates in the best-fitting model (including 
the year effect) were estimated by bootstrapping. The 
data were resampled 500 times with replacement, the 
form of the best-fitting model was refitted to the resa- 
mpled data, and the parameters were recalculated for 
each bootstrap. Confidence intervals for the mean were 
estimated for each of the years from the bootstrap data. 
The year-effect estimate and confidence intervals were 
then back-transformed from log(CPUE) and the constant 
