Assessing the impact of environmental variability on harvest in a heterogeneous fishery: a case study of the Canadian lobster fishery

ABSTRACT

Global fisheries face significant challenges in the coming years due to climate change. Understanding and anticipating the impacts of climate change is a necessity for implementing appropriate fisheries management. This study uses a panel dataset of individual fishing vessels to examine how variation in ocean temperature affects fish harvest. Using the American lobster (Homarus americanus) fishery in the Maritimes region of Canada as a case study, this paper employs a generalised linear mixed model (GLMM) taking into account heterogeneity amongst fishers, gear, vessels, and fishing areas. The GLMM is found to have better performance and estimations when compared against alternative specifications. As expected, a significant and positive relationship was found, further contributing to the existing evidence of warming impacts on the lobster fishery. The implications of this study are twofold: first, it provides further evidence that environmental change does have a significant positive impact on harvest. This information should be considered by fishing industry and fisheries authorities when implementing appropriate adaptive management strategies and measures in their decision making. Second, it illustrates that allowing for mixed-effects using GLMMs is a valuable empirical tool when dealing with hierarchical data structures.

KEYWORDS:

• Generalised linear mixed model
• mixed-effects
• bioeconomics
• climate change
• environmental variable
• heterogeneity

JEL CODE:

• C18
• C33
• Q21
• D22
• Q50
• C52
• Q54

1. Introduction

Global fisheries face significant challenges in the coming years, many of which are caused either directly or indirectly by climate change. According to the Intergovernmental Panel on Climate Change (IPCC)’s special report on the ocean and cryosphere, the world’s oceans are becoming warmer and more acidic, impacting the productivity, abundance, and distribution of marine species (Bindoff et al., 2019). As waters warm beyond species’ optimal range, species that once inhabited certain areas may begin to move their distribution northward and into deeper waters, while some may begin to die out entirely. A global study revealed that climate change extensively affects the distribution of global catch potential leading to changes in fisheries productivity, with increase in the polar regions and a loss in the tropics (Cheung et al. 2010). Harvesters and communities that are heavily reliant on fisheries revenue are the most vulnerable to these changes. This is especially true for the ones who cannot diversify the species they catch or obtain alternative employment. Climate-induced changes in productivity and distribution pose challenges to fisheries management (Bryndum-Buchholz et al. 2020). Thus, anticipating the consequences of climate change on fishing operations and advancing our understanding of these impacts is crucial for developing appropriate mitigation and adaptation measures.

This paper seeks to empirically investigate the impacts of changing environmental conditions on a fishery resource using a panel dataset of harvest and fishing effort of commercial vessels. The inclusion of environmental variables in the harvest production function is a convenient way to estimate the relative effects on harvest by treating the environmental variable as an additional input in production. However, this is not without its limitations – there are several pitfalls that can arise when attempting this type of modelling exercise. For example, the model can suffer from omitted variables such as missing stock size data or unobserved characteristics between subjects. These latent variables may be correlated with the explanatory variables, violating the key assumption of independence that is required for ordinary linear regression. While traditional fisheries production models assume homogeneity among fishers, vessels, and fishing areas, more recent research confronts the issue of how to incorporate heterogeneity among these factors of production. This is a key issue that is explored in depth in this paper.

A generalised linear mixed model (GLMM) with both fixed and random effects is employed to account for heterogeneity among fishing vessels and fishing locations. It will attempt to isolate the effect of ocean bottom temperature on annual harvest from other production inputs; this effect is disentangled from the operational and effort-based measures such as number of fishing days, vessel size, and vessel power. Exploiting the hierarchical structure of the data, we are able to mitigate some of the common pitfalls that are often encountered, and to show that efficiency in estimation is improved. Compared against a model without mixed effects, we find that the mixed-effects model performs better in estimation. As a case study, we use the American lobster (Homarus americanus) fishery in Canada to explore how ocean bottom temperature affects harvest, but it is important to note that this method can also be applied to other fisheries and to other environmental variables. The main contribution of this paper is to provide a robust econometric framework for empirical estimation of environmental elasticities as an input in production of a natural resource.

2. Case study: the American lobster fishery in Atlantic Canada

The American lobster (Homarus americanus) is geographically distributed along the coast of the northwestern Atlantic, ranging from North Carolina to Newfoundland and Labrador (DFO 2020). The most abundant populations are found in the Gulf of Maine, the Nova Scotian shelf, and the southern Gulf of St. Lawrence. The lobster fishery is the most commercially important fishery in Atlantic Canada, with annual landed value exceeding CAD $1 billion (DFO 2021a).¹ The fishery in the Scotia-Fundy region alone provides employment for approximately 7,500 people and generates many other direct and indirect economic benefits (DFO 2020). It has become the backbone for the inshore commercial fisheries in the region. It is ecologically important to the biodiversity of the area, and it is an intrinsic part of the culture and identity of the East Coast of Canada (Greenan et al. 2019). Although the Canadian lobster fishery has seen record landings in recent years, the future of the fishery is uncertain due to risks posed by climate change. It is difficult to predict the net effects of climate change on lobster populations as stocks will be influenced by changes at both the regional level and larger scale changes in ocean conditions (DFO 2020).

Lobster fisheries in Canada are primarily located in inshore waters of the Maritimes region. This paper focuses on the lobster fisheries in the Scotia-Fundy region of Canada, which includes the eastern coast of Nova Scotia and the Bay of Fundy. These correspond to lobster fishing areas (LFAs) 27–34 and 35–38, respectively (Figure 1). This region encompasses the inshore waters from the northern tip of Cape Breton to the New Brunswick-Maine border and is one of the most productive regions for lobster fishing in Canada. This paper focuses only on the inshore lobster fishery, as the offshore lobster fishery (LFA 41) is managed separately and accounts for a small amount of annual landings.

Figure 1. Inshore lobster fishing areas (LFAs) in the Scotia-Fundy region of Canada. The shaded regions depict average annual landings in thousand tonnes for the years 2014–2018.

The inshore lobster fisheries are managed through effort-based controls with limits on number of licences and traps, delimited seasons and zones, and protection of juvenile and ovigerous females. The fishing areas (LFAs) are the primary management tool to control fishing effort for each designated area, and the opening and length of fishing season vary across different LFAs (Reid-Musson and Neis 2022). There are also restrictions on number of traps per licence holder that also varies by licence type and LFA. Fishing effort is fully competitive, and this type of management may potentially create a ‘race to fish’ scenario and intensify harvesters’ safety risks (Reid-Musson and Neis 2022).

In the North Atlantic, air temperatures are rising and ocean circulation patterns are changing, leading to higher temperatures both at the surface and in deeper waters of the ocean (DFO 2021b). It has been well-established that rising ocean temperature has an impact on lobsters’ habitat preference, coincided with the species’ range shifting further north (Le Bris et al. 2018; Greenan et al. 2019; Goode et al. 2019; Tanaka et al. 2020). As ocean temperatures are cooler in the Scotian Shelf and Gulf of Maine than in areas that have experienced stock depletion such as New England, warming ocean temperatures have increased lobster habitat suitability in these regions (ASMFC 2020). Higher temperatures also have an impact on lobsters’ growth, size at maturity, and reproduction. It affects molting phenology making them more vulnerable to predators, and it increases susceptibility to epizootic shell disease (Groner et al. 2018). These climate-induced changes in habitat and biological functions present challenges for the management of lobster fisheries in the region.

3. Literature review

There are two broad categories of economic production models to analyse fisheries: optimal or simulated bioeconomic production models (e.g. the classical Gordon-Schaefer model) and empirical econometric modelling. The former combines biological/ecological and economic components in the optimal or simulation setting for homogenous fisheries, while the latter applies econometric techniques based on cross-sectional, time-series, and panel data on individual fishing vessels. A detailed review of fisheries production models can be found in Squires and Walden (2021).

It has long been acknowledged that ocean biophysical conditions such as water temperature, salinity, waves, wind, and storms have effects on fish stocks (Bindoff et al., 2019). Given the uncertainty and complexity of the linkages between environmental conditions and fish stocks, the incorporation of environmental variables in fisheries production models has been sporadic. Population dynamics are influenced by environmental characteristics in many complex, and often nonlinear ways. The limitations in methodologies and availability of data have impeded researchers’ abilities to empirically assess the impacts of climate change on fishery resources and fishing sectors. The inherent complexity of ecological systems lends to interactions with the environment in bioeconomic modelling being simplistic. Despite this, there are some notable examples of economic models being applied to estimate production in the midst of a changing environment.

Environmental variables might be included as an input in production (Barbier 2000), but alternatively they might enter the biomass growth function, or alter consumers’ utility functions. Lynne, Conroy, and Prochaska (1981) examined the Florida Gulf Coast blue crab fishery, which relies on the threatened mangrove forests as habitat. In this case, the extent of mangrove forest is the environmental input and the relationship between catch and mangrove area is modelled. Similarly, Barbier and Strand (1998) modelled catch against mangrove area in the shrimp fishery in the Bay of Campeche, Mexico. Kahn and Kemp (1985) estimate the economic losses from the destruction of submerged aquatic vegetation on the commercial striped bass fishery in Chesapeake Bay. They estimated industry supply and demand functions where environmental degradation enters the supply function, and the equations are solved to find the bioeconomic equilibrium catch under different levels of degradation. Foley, Van Rensburg, and Armstrong (2010a) and Foley, Kahui, Armstrong (2010b) examine the bioeconomic interplay between cold water coral as habitat for fish species and the impacts of habitat reduction on these species. Cheung et al. (2010) used a dynamic bioclimate envelope model to project the maximum exploitable catch of a species under climate change scenarios, and the findings suggest that the polar region is benefiting while the tropics are losing from climate change.

Several empirical studies have linked environmental variables to harvest in the lobster fisheries of Canada and the United States. Bell and Fullenbaum (1972) were one of the first to include a variable for environmental quality in the analysis of the inshore lobster fishery in the United States. In this case study, seawater temperature appeared directly in the production function. The results indicated that seawater temperature has a positive effect on the growth of the lobster stock, citing trends that suggest that declining seawater temperature is partially responsible for declining coastal lobster catches. Henderson and Tugwell (1979) estimated a production function for two lobster fishing areas in Nova Scotia that included both current and lagged bottom temperature. The assumption is that temperature affects the catchability of lobster, as lobsters tend to move around more and cover more territory when temperature rises. Several other studies found correlations between lobster harvest and ocean temperature, such as McCleese and Wilder (1958), Dow (1961), and Flowers and Saila (1972). Hudon (1994) and Drinkwater et al. (2006) similarly found correlations between harvest and temperature, and in addition found wind to be a significant determinant of catch amounts.

As with some of the existing studies, this current paper incorporates environmental values by having the environmental variable enter the model as an input directly in the harvest production function. With the increasing availability of panel data and more advanced econometric techniques, empirical vessel-level analysis in response to environmental issues has garnered more interest. Huang, Smith, and Craig (2010) used a differenced bioeconomic framework and individual fishing data combined with oxygen monitoring data to quantify the economic effects of hypoxia on the brown shrimp fishery in North Carolina. Later, Huang and Smith (2014) applied a restricted Cobb–Douglas production function to model harvest which is determined by both fishing inputs and other inputs such as season closure, wind, waves, and stock size. Autoregressive models combined with Seemingly Unrelated Regression (SUR) models were applied to estimate an output production function. Nguyen (2022) used a partial equilibrium analysis with combined production, demand, and aggregate supply functions to project welfare impacts of climate change on fisheries in Vietnam. In this study, the production function included sea surface temperature, precipitation, number of typhoons, maximum wind speeds of typhoons, and the Southern Oscillation Index. An autoregressive distributed lag model (ARDL) was used to predict fishery yields.

When the panel data are complex and nested, e.g. fishers operating within fishing areas, multilevel modelling (otherwise known as mixed-effects or hierarchical modelling) is a powerful tool that can accommodate various data structures and improve inference (Gelman 2006; Gelman and Hill 2006). Generalised linear models (GLMs) are a generalisation of linear regression that allow a linear model to be related to the response variable via a link function (Nelder and Wedderburn 1972). Generalised linear mixed models (GLMMs) are an extension of GLMs that allow for both fixed and random effects to be estimated when data has a hierarchical structure (Breslow and Clayton 1993). The fixed effect is the population-averaged effect, while the random effects are the subject-specific effects that manifest through variances that represent individual-specific heterogeneity. Mixed-effects modelling can be a useful tool in fisheries research because there are often many unobserved characteristics that confound estimation (Hyun, Cadrin, and Roman 2014; Thorson and Minto 2015), but this method has also been widely applied in disciplines such as ecology (Venables and Dichmont 2004; Bolker et al. 2009; Harrison et al. 2018), psychology (Moscatelli, Mezzetti, and Lacquaniti 2012; Meteyard and Davies 2020; Bono, Alarcón, and Blanca 2021), and medicine (Dean and Nielsen 2007; Casals, Girabent-Farrés, and Carrasco 2014).

4. Material and methods

4.1. Methodology

A Cobb–Douglas harvest production function is specified, the goal of which is estimate and assess the relative effects of an environmental variable (ocean temperature) as well as technical and effort-based variables (number of days fished, length of vessel, and vessel tonnage) on lobster harvest. The hope is that the model yields reliable coefficient estimates that can allow us to glean something about the relationship between environmental variability and harvest without undue complexity.

First, we consider the following basic model, (1) $y_{\mathrm{ijt}}={\beta }_0+{\beta }_1{d}_{\mathrm{ijt}}+{\beta }_2{l}_{\mathrm{ijt}}+{\beta }_3{g}_{\mathrm{ijt}}+{\beta }_4{t}_{\mathrm{jt}}+{\epsilon }_{\mathrm{ijt}}$(1) where $y$ is harvest, d is days fished, l is length of vessel, g is gross tonnage, t is temperature anomaly, and ϵ is the idiosyncratic error term. The subscript i corresponds to each vessel, j to LFA, and t to year. Since the data have a hierarchical structure (vessels operating within fishing areas), we exploit this by using a generalised linear mixed model (GLMM) that allow for random effects to be specified. In matrix form, the specification is given by, (2) $g\{E(y|X,u)\}=\mathrm{X\beta }+\mathrm{Zu}$(2) where the dependent variable y is an $n\times 1$ vector. $X$ is an $n\times k$ design matrix for the fixed effects $\beta$. This contains the explanatory variables and their associated coefficients. $Z$ is an $n\times m$ design matrix for the random effects $u$. $g$ is the invertible link function which can take on many different functional forms. In our case, the outcome variable is annual harvest of individual fishing vessels, which is a continuous and repeated variable and is heavily right-skewed, not normally distributed (Figure 2). Different distributions can be used to deal with this, but log-normal is chosen as we suspect this most closely approximates the data-generating process. Therefore, we choose a log link function and a Gaussian distributional family.

Figure 2. Histogram of harvest observations (kg).

Taking the hierarchical structure of the data into consideration and allowing for random intercepts, we consider the following extended model, (3) $y_{\mathrm{ijt}}={\beta }_0+{\beta }_1{d}_{\mathrm{ijt}}+{\beta }_2{l}_{\mathrm{ijt}}+{\beta }_3{g}_{\mathrm{ijt}}+{\beta }_4{t}_{\mathrm{jt}}+u_{\mathrm{ij}}+u_j+{\epsilon }_{\mathrm{ijt}}$(3) where the first six terms represent the fixed component and are equivalent to Equation (1). This is the dependent variable, the explanatory variables, and their estimated coefficients. The last three terms represent the random component which are two additional random intercepts and the error term. To illustrate the intuition behind the three levels of effects, consider the following equations, (4) $y_{\mathrm{ijt}}={\gamma }_{0j}+{\gamma }_{0\mathrm{ij}}+{\beta }_1{d}_{\mathrm{ijt}}+{\beta }_2{l}_{\mathrm{ijt}}+{\beta }_3{g}_{\mathrm{ijt}}+{\beta }_4{t}_{\mathrm{jt}}+u_{\mathrm{ij}}+u_j+{\epsilon }_{\mathrm{ijt}}$(4) (5) ${\gamma }_{0j}={\beta }_{00}+u_{\mathrm{oj}}$(5) (6) ${\gamma }_{0\mathrm{ij}}={\beta }_{000}+u_{\mathrm{oij}}$(6)

The equation for the intercept ${\gamma }_{0j}$ consists of the LFA-level mean intercept ${\beta }_{00}$ and an LFA-specific random intercept $u_{\mathrm{oj}}$. The equation for the intercept ${\gamma }_{0\mathrm{i}j}$ consists of the overall mean intercept ${\beta }_{000}$ and a vessel-specific random intercept $u_{\mathrm{oij}}$. Rearranging and collapsing to one intercept for the fixed portion of the model, we are left with Equation (2) where $u_j$ represents the LFA-level unobserved effects and $u_{\mathrm{ij}}$ represents the vessel-within-LFA unobserved effects.

The reasoning behind this choice of model is that in addition to stochastic effects, there are unobserved effects that are particular to LFAs and vessels. At the vessel level, there may be technological inputs or crew skill that are not accounted for. At the LFA level, there are oceanographic, biological, and ecological influences that are difficult to observe or quantify. Accounting for multilevel structures in the data can improve statistical inferences. As opposed to ordinary linear regression which treats all explanatory variables as independent and calculates standard errors using only the residual variance, mixed-effects models allow subjects at each level to deviate by its own mean and calculates standard errors using both the residual variance and the variance between the higher levels of the hierarchy. By not accounting for the variances at the different levels when hierarchical structures exist, coefficient estimates are likely to be biased upward and can lead to a type I error, i.e. finding statistical significance when none exists.

This method enables us to obtain reliable estimates despite missing pertinent data. Perhaps most importantly, biomass data are not directly observed, and this creates difficulties for this type of modelling exercise. Biomass is a crucial component of the traditional harvest function, however, data on population size is often unavailable or spotty across space and time. Since all vessels operating in the same LFA share the same stock in each time period, the biomass can be treated as a random intercept. Although not a perfect solution, this method should still allow us to obtain coefficient estimates for the environmental variables, which is what we are interested in.

4.1.1. Estimation procedure

To justify the application of a multilevel model, we analyse the variance components by calculating the intraclass correlation coefficient (ICC). Assuming that our model is correctly specified, conditional on the explanatory variables (the fixed part of the model), the ICC calculates the dispersion of harvest weight around a mean value for each level of the hierarchy (the random intercepts$u_i$ and$u_{\mathrm{ij}}$). The intuition behind this is that the higher the correlation is within the clusters (the larger the ICC) the lower the variability is within the clusters and the higher the variability is between the clusters. Having a high degree of correlation within clusters can lead to biased estimates if regular pooled regression is used, as it violates the assumption of independence.

The ICC is calculated and the level-3 intraclass correlation, which is the correlation of observations between the same LFA, is estimated to be 0.21. The level-2 intraclass correlation, which is the correlation between yearly observations for each vessel, is estimated to be 0.68. Therefore, conditional on the covariates, we find that annual harvest is only weakly correlated within the same LFA, but strongly correlated across year classes for each vessel. A rule of thumb can be used to determine if the hierarchical structure should be taken into consideration. The literature defines the design effect of a sample statistic as the ratio of the actual variance of a sample to the variance of a simple random sample of the same number of elements (Kish 1965). In multilevel modelling, the design effect (DE) is estimated as a function of the ICC and the average size of the clusters (c) (Muthen and Satorra 1995): (7) $\mathrm{DE}=1+(c-1)\ast \mathrm{ICC}$(7) The rule of thumb is that when the design effect is less than 2, the multilevel structure can be ignored. The average number of vessels in each LFA is calculated and it equals 598.9. This leads to a design effect of ∼127, and thus the multilevel structure of the data cannot be ignored. The coefficients for the mixed-effects model with random intercepts are estimated using maximum likelihood.

In addition to the mixed-effects model, two other versions of the model are employed for comparison (Table 2). The first column is the model with only the fixed component (Eq. 1). The second column is the mixed-effects model but with random effects only at the vessel level, and the last column is the mixed-effects model with random effects at both the vessel and the LFA level (Eq. 3).

4.2. Data collection

4.2.1. Logbook data

The landings data used for this analysis are from commercial fisheries logbook data collected by Fisheries and Oceans Canada (DFO). The variables retrieved from the logbook data are catch in kilograms, date landed, lobster fishing area (LFA), vessel identifier, vessel length, and vessel tonnage. The data span the years 2006–2018. The logbook data is geographically delineated by NAFO division, subdivision, and LFA, but LFA was chosen as the fishery is managed at the LFA level. LFAs 28 and 29 were excluded as sampled temperature data points were sparse in this area, and LFA 37 is a small area shared by LFAs 36 and 38. All variables are continuous except for tonnage class which is a dummy variable that is coded 1 if the vessel is greater than 25 tons and 0 if the vessel is less than 25 tons. The reason for this is that it is a categorical variable in the logbook data. After compiling the logbook data, we are left with 10 fishing areas and 4,064 vessels operating over the 13 years analysed. The data have a hierarchical structure in that each LFA contains vessel observations, and each vessel contains repeated observations for each year. The panel is balanced in the spatial dimension in that there are observations for each year for each LFA, but highly unbalanced at the fleet level. Some vessels drop in and out, some are retired, and some are newly added. A total of 36,672 observations will be used for the analysis.

4.2.2. Bottom temperature data

Complete ocean temperature profiles at all depths were retrieved from DFO’s Marine Environmental Data Section (MEDS) (DFO 2021c). To determine which profiles reached the seafloor, these data were cross-referenced with the Canadian Hydrographic Service Non-Navigational (NONNA) Bathymetric Data (DFO 2021d). Although there are a large number of observations that span the study area, the locations sampled are not consistent from year to year. Therefore, bottom temperature data are expressed in terms of anomalies from their long-term mean. Present climate data were retrieved from DFO’s Bedford Institute of Oceanography North Atlantic model (BNAM) (Wang et al. 2018). Bottom temperature in the model are monthly averages for the years 1990–2015 on a spatial grid of 1/12 degrees. Each real temperature observation was matched with its nearest geodetic neighbour from the BNAM long-term averages, and the difference between these points were calculated by subtracting the long-term average observation from each real temperature observation. Shapefiles that delineate the LFAs’ polygons were used to determine which observations fall into which LFA and the observations were grouped accordingly. The average temperature anomalies for each LFA for each year were then calculated, and these anomalies are proxies for temperature change. A negative observation represents a colder-than-average anomaly, and a positive anomaly represents a warmer-than-average anomaly. The minimum temperature anomaly is – 2.9 degrees Celsius and the maximum is 6. Since logarithms cannot be taken for negative values, the temperature anomalies are rescaled so that the minimum is zero and the maximum is 100. All of the variables used in the analysis and their associated summary statistics are given in Table 1.

Table 1. Variables and associated summary statistics.

	Mean	S.D.	Min	Max
Bottom temperature anomaly
Degrees Celsius	1.62	1.61	−2.88	6.03
Rescaled	50.56	18.02	0.28	99.997
Vessel and effort variables
Weight landed (kg)	14,796	12,564	6	170,994
Number of days fished per vessel	42	23	1	624*
Length of vessel (feet)	38.08	6.08	0	64
Vessel tonnage ≥25 tons (dummy)	0.07	0.25	0	1

*Since number of logbook entries is used as a proxy for days fished, there are some instances where fishers record multiple entries per day. There are two instances in which the number of a vessel’s logbook entries exceeds the number of calendar days: 624 in 2006, and 397 in 2013.

5. Results and discussion

For all model specifications the coefficients are statistically significant at the 1% level, with the exception of vessel tonnage ≥25 tons (Table 2). The model that achieves the best fit according to the Akaike information criterion (AIC) is the mixed-effects model with random effects at both the vessel and LFA level. All variables except vessel tonnage are continuous, therefore the log-transformed variables are elasticities. Vessel tonnage ≥25 tons is a dummy variable, so the exponentiated coefficient is the ratio of the mean harvest weight for vessels ≥25 tons to the mean harvest weight for vessels <25 tons. Therefore, this would give an estimate of the expected percent increase in mean harvest weight that would result when going from vessels <25 tons to vessels ≥25 tons, holding other variables constant. However, the estimated coefficients are statistically insignificant for all four models. The negative coefficients for two of the three models suggest that an increase in tonnage would decrease output which is contrary to what would be expected. The statistical insignificance could be because harvest is inelastic to vessel tonnage, but it should be noted that the vast majority of vessels (93%) are < 25 tons.

Most coefficients have the expected signs, and reasonably plausible magnitudes. The results from the estimation of our chosen models using the historical data available suggest that a 1% increase in ocean bottom temperature anomaly results in a 4.8%–15.6% increase in harvest by weight, depending on model specification. The directional signal tracks with previous studies that found positive correlation between bottom temperature and harvest. However, it is important to note that the precision of the estimated coefficients rely heavily upon the transformation of the temperature data to anomalies. Although this was chosen as the least problematic method given inconsistencies in the sampling areas, this makes interpretation difficult and we would be better served by using raw temperature data if this was possible.

As expected, number of days fished and vessel length have a large and positive impact on harvest. Thus, controlling for technical efficiency of the fleet is crucial as changes in harvest can be attributed to changes in biomass or changes in the amount of fishing pressure that is applied. The increase in landings seen in recent years may be caused by an increase in abundance, which could be linked to climate-related factors, but could also be the result of increased fishing capacity. Although the number of licence holders has decreased and the number of traps allocated to each licence holder has stayed relatively constant, it is likely that effort capacity and fishing efficiency has effectively increased. According to the 2007 report by the Fisheries Resource Conservation Council (FRCC), there is general agreement that harvesters’ ability to catch lobster has improved due to improved gear, vessels, and technology (Fisheries Resource Conservation Council 2007).

Table 2. Coefficients and standard errors from maximum likelihood estimation: the model with the fixed component only, the mixed-effects model with only vessel-level random effects, and the mixed-effects model with random effects at both the vessel and fishing area-level.

	Fixed component only (Eq. 1)	Mixed-effects model (random effects only at vessel level)	Mixed-effects model (random effects at both vessel and LFA, Eq. 3)
Variable	Coeff. (St. error)	Coeff. (St. error)	Coeff. (St. error)
Intercept	−3.4821 (0.0807)***	−2.0762 (0.1806)***	0.1789 (0.2228)
Days fished	0.5354 (0.0045)***	0.5729 (0.0046)***	0.5726 (0.0047)***
Vessel length	2.8232 (0.0221)***	2.498 (0.0495)***	1.9138 (0.0526)***
Vessel tonnage ≥25 tons	−0.0115 (0.0147)	0.0261 (0.0334)	−0.0273 (0.0298)
Bottom temperature anomaly	0.1564 (0.0072)***	0.0558 (0.0054)***	0.0476 (0.0054)***
Random effects
LFA			0.1301 (0.0598)**
Vessel		0.3707 (0.0092)***	0.2994 (0.0075)***
Overall variance	0.4833 (0.0036)	0.2105 (0.0017)	0.2015 (0.0016)
AIC	76388.45	57931.53	56197.42
BIC	76439.42	57991.00	56265.40
Log-likelihood	−38188.22	−28958.76	−28090.71

*** p < 1%; **p < 5%; *p < 10%

Results from the mixed-effects models are compared with results from the model with only the fixed component. Since the portion of the variance attributed to LFA is fairly small, a mixed-effects model with random effects only at the vessel level was also compared. It is useful to compare the fixed coefficients that are estimated under the different specifications as it provides insight into what exactly mixed-effects modelling does, how it is different, and why it is useful. The coefficient for number of days fished does not change significantly with the different specifications. Meanwhile, the coefficient for vessel length decreases from 2.8 with the fixed model to 1.9 with the mixed-effects model. This suggests that not accounting for the variance at the vessel level results in an over-estimation of the coefficient for vessel length.

The difference in the coefficients between the fixed model and the mixed-effects model can be attributed to the concept of partial pooling, otherwise known as shrinkage. When the data are pooled, each observation has an equal chance of success. With partial pooling, each unit (vessel in our case) has a different chance of success and this is informed by the vessel-specific characteristics. This allows vessels with less observations and more extreme values to borrow strength from vessels with more observations and less extreme values, and therefore ‘shrinks’ the estimated coefficient back to a more reasonable value. It acknowledges that each unit has characteristics in common, while also placing less importance on extreme values within each unit (Clark 2019).

A notable takeaway is that the estimated bottom temperature anomaly elasticities vary substantially across the models. It decreases from 15.6% in the fixed-only model to 5.6% in the model with random effects at the vessel level only to 4.8% in the model with random effects at both the vessel and area levels. It makes sense that the coefficient does not differ significantly between the latter two, as the variance at the LFA level is not very large. The difference between the coefficients from the fixed and mixed-effects models can again be attributed to the idea of partial pooling. The influence of bottom temperature on harvest may be higher when considering each unit separately, but when imposing a normal distribution on each vessel’s observations this makes extreme values less probable, thus shrinking the coefficients back toward a more reasonable value.

This empirical exercise is useful for investigating the historical impact of bottom temperatures on harvest as well as illustrating the merits of mixed-effects models. However, it is important to note that this alone cannot be used to make predictions of future harvest. Lobsters’ tolerance to temperature exhibits a bell-shaped curve. Although lobster abundance has increased steadily as temperatures have increased, it is still at the ascending part of the curve. Once temperatures reach a certain point, productivity will start to decline and distribution will be shifted offshore and into deeper water (Oppenheim et al. 2019). Additionally, there are many complex interactions at play that are not accounted for by this convenient production framework. Although there is a well-established link between landings and temperature, there are different pathways through which this is realised. For example, temperature can impact landings through changes in recruitment, likelihood to enter traps, availability of food sources, etc. Models that contain these interactions are complex and are beyond the scope of this paper. Making predictions based on future environmental scenarios would involve a much more complex model, and this would be interesting to explore as future research.

5.1. Model validation

In addition to the coefficient estimates, best linear unbiased predictions (BLUPs) are retrieved for the random effects (Henderson 1950). Usually, random effects are only reported in terms of variance components, but BLUPs can also be estimated in addition to coefficients as another form of model selection. By inserting BLUPs into the estimation equation and solving, fitted values can be obtained. To assess goodness of fit, we plot the log-transformed harvest weight observations against the fitted values. We compare the fitted values of the model with the fixed component only with the mixed-effects model with random effects at both the vessel and area level (Figure 3). We also plot the residuals from these models against the fitted values (Figure 4).

Figure 3. Observed log-transformed harvest versus fitted values from the model with the fixed component only (left) and the mixed-effects model (right).

Figure 4. Residuals versus fitted values from the model with the fixed component only (left) and the mixed-effects model (right).

Visualising the fitted values and the residual errors allows us to see how the model performs with and without the random effects. Figure 3 clearly shows that the GLMM results in a better goodness-of-fit than the model without mixed effects. Figure 4 shows that the residuals are more tightly centred around zero with the GLMM; they are dispersed fairly evenly above and below zero, although they tend more toward negative values. This implies that our model tends to over-estimate, but this is made much less severe with the GLMM.

5.2. Random effects estimates

It is also helpful to look at how the estimated random effects change over time. The median vessel-level random effects hover around zero but there is a significant amount of dispersion (Figure 5, left panel). The estimated unobserved effects at the vessel level also appear to be increasing over the time period analysed. This suggests that indeed there are unobservable effort-related factors at play such as more crew or more advanced technology, and that these are increasing over time. This is consistent with the FRCC report’s suggestions that harvesters’ ability to catch lobster is improving. This can muddy the waters for estimation because the traditional effort-based measures will be under-estimating the operational influences on harvest. Consequently, this makes incorporating environmental variables into the analysis difficult as the true effect will be harder to isolate from the noise. The random effects at the fishing area level remain relatively constant over the time period, suggesting that unobserved effects at these larger spatial scales are time-invariant and less problematic for estimation (Figure 5, right panel). Perhaps most importantly, how much of the increase in harvest is attributed to increasing stock size rather than increased fishing pressure or technological change remains a mystery for now. Catch per unit of effort (CPUE) in terms of catch per fishing trip in the logbook data displays an increasing trend, but with the available data it is not possible to say conclusively what is driving this. Further research that uses more detailed effort data such as number of trap hauls may be able to uncover this. Unfortunately, the number of trap hauls per fishing trip is not available in the data that were available to us.

Figure 5. Box and whisker plots of vessel-level random effects estimates (left) and LFA-level random effects estimates (right).

6. Conclusions

The aim of this paper was to examine the economic consequences of ocean temperature increases on the lobster fishery in Atlantic Canada. This study pulled data from multiple sources and set up an empirical econometric framework that models annual lobster harvest as a function of environmental and operational variables. Coefficients were estimated using a maximum-likelihood mixed-effects estimator and these were compared with two other model specifications to see which model fits the best. The mixed-effects model was selected as the best fitting model according to the AIC, and goodness-of-fit plots visually substantiated this. All coefficients were statistically significant at the 5% level apart from vessel tonnage. The estimated coefficients suggest that a 1% increase in bottom temperature anomaly results in approximately a 5% increase in harvest weight with the chosen model specification, and this was the most conservative estimate of the three models. The estimated coefficients for the effort-related variables vessel size and number of fishing days are large and positive, as expected. Best linear unbiased predictions (BLUPs) were estimated for the random effects at the vessel level. The unobserved effects at the vessel level are significant and have increased over the years 2006–2018. This complicates matters when trying to isolate the relative effects of environmental variables, and suggests that finding a way to incorporate technological progress would improve estimation. There is a vast body of research that aims to measure technical efficiency in fisheries (i.e. stochastic frontier analysis) and incorporating this would be a boon to future research.

The implications of this study are twofold: first, it provides further evidence that environmental change does have a significant impact on harvest. Although this empirical framework does not capture all of the intricate ecological systems at play, it is highly likely that the significant and positive effect of temperature on harvest is reflective of real-world phenomena. As fisheries are confronted with considerable environmental uncertainty in the coming years, there is a fear that warming waters will exacerbate other issues such as excess fishing pressure and competition. This must be considered by policymakers when implementing management measures, as distributional changes may tempt increases in allowances or access. In the case of the Canadian lobster fishery, the current management measures include designated fishing areas, fishing seasons, and limits on the number of traps per licence holder. This management regime incentivizes a race to fish the most productive areas, and fishers have the incentive to outcompete others by increasing vessel power, size, speed, number of crew, or make other capital investments aimed at maximising catch (Pfeiffer and Gratz 2016). As lobsters’ preferred habitat shifts, the existing management measures may not be adequate to mitigate excess capacity. In the colder parts of lobsters’ range, increasing temperatures may lead to greater abundance, and these areas might become more attractive for fishing. In warmer areas that are approaching the upper limit of lobsters’ thermal range, the negative impacts on the species’ physiology will begin to manifest. The strong relationship between harvest and temperature underscores that conservation measures should be taken while temperatures are still at a manageable level. This analysis is more exploratory than prescriptive, and caution must be taken when making assumptions about the future of the fishery based on historical trends. However, it calls attention to something that must be delved into deeper.

The second implication of this analysis is one of a more technical nature: that mixed-effects modelling can be a useful part of the natural resource economist’s toolbox when the data are hierarchically structured. Given that fisheries management is often area-based, mixed-effects models are surprisingly under-utilized in fisheries economics research. Models that combine both fixed and random effects provide a more flexible approach for analysing the data that are not normally distributed. In this case, we find that the model that most closely resembles reality is the model with random effects at the vessel level and the fishing area level. On the other hand, a pooled model with only fixed effects that ignores the hierarchical structure of the data resulted in a poorer model fit. This is a lesson in the importance of accounting for heterogeneity when these data structures exist.

Acknowledgements

We thank Barbara Best and Gisele Magnusson at Fisheries and Oceans Canada (DFO) for their gracious assistance in data collection and manuscript review. We also thank Luc Bujold at DFO’s Marine Environmental Data Section for providing bottom temperature data, Zeliang Wang at DFO’s Bedford Institute of Oceanography for providing BNAM model output, and Christine Irwin at DFO’s Economic Analysis and Statistics group for providing shapefiles.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1 Annual landed values in the last five years were consistently above $1 billion with the exception of 2020, which was valued at $761 million.