Likelihood Methods for Basic Stratified Sampling, with Application to Von Bertalanffy Growth Model Estimation

Nan Zheng, Noel Cadigan

Centre for Fisheries Ecosystems Research, Fisheries and Marine Institute of Memorial University of Newfoundland, St. John’s, NL, Canada.

DOI: 10.4236/ojs.2019.96040                

Abstract

This paper mainly addresses maximum likelihood estimation for a response-selective stratified sampling scheme, the basic stratified sampling (BSS), in which the maximum subsample size in each stratum is fixed. We derived the complete-data likelihood for BSS, and extended it as a full-data likelihood by incorporating incomplete data. We also similarly extended the empirical proportion likelihood approach for consistent and efficient estimation. We conducted a simulation study to compare these two new approaches with the existing estimation methods in BSS. Our result indicates that they perform as well as the standard full information likelihood approach. Methods were illustrated using a growth model for fish size at age, including between-individual variability. One of our major conclusions is that the fully observed BSS data, the partially observed data used for stratification, and the sampling strategy are all important in constructing a consistent and efficient estimator.

Keywords

Length-Stratified Age Sampling, Response-Selective Sampling, Basic Stratified Sampling, Complete-Data Likelihood, Empirical Proportion Likelihood

1. Introduction

In stratified random sampling (SRS), the population or a random sample of the population is partitioned into relatively homogeneous subgroups, or strata, and then random samples are taken independently in each stratum for full observation. Such sampling design may also be regarded as a kind of two-phase sampling, with the population or the large sample before partitioning being the first phase sample, and the smaller and more extensive subsamples after partitioning being the second phase samples.

Practical implementations of SRS frequently fall into two categories as classified by [1]: 1) basic stratified sampling (BSS) where the maximum second phase subsample size (BSS1) or subsampling fraction (BSS2) in each stratum is prefixed, and 2) variable probability sampling (VPS) in which sequential units are independently generated from a model and then classified into strata where they are selected for full observation with pre-specified probabilities. [2] classified BSS2 as VPS, and hence all the inference methods for VPS are also suitable for BSS2.

Assume that there are a total of N sampling units on which the stratified sampling is conducted. Let $y_i$ and $x_i$, $i=1,\cdots ,N$, denote respectively the vectors of responses and covariates of the ith individual generated from the joint distribution $f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)=g\left(y\mathrm{<mo>\; |\; </mo>}x\mathrm{<mo>\; ;\; </mo>}\theta \right)h\left(x\right)$, with $\theta$ being a vector of all the parameters describing the conditional distribution of $y$ given $x$. In SRS $\left(y\mathrm{,}x\right)$ are fully observed only for a subset of size n of the N units, which are called complete data in this paper, and only a subset $z$ of $\left(y\mathrm{,}x\right)$ is observed for the other $N-n$ units, which are called incomplete data.

In SRS the unobserved elements of $y$ and/or $x$ are missing data, and missingness can be fully accounted for by variable $z$ which is observed for all the N units; that is, the unsampled variables are missing at random (MAR) in the terminology of [3]. In addition, for BSS and VPS, given the observed data, the missing probability for all the missing data is a constant involving no parameters $\theta$. As a result, the likelihood, which is called full information likelihood, is given by (see e.g. [1])

$L_F\left(\theta \right)=\left[{\displaystyle \underset{i=1}{\overset{n}{\prod }}}f\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}\theta \right)\right]\left[{\displaystyle \underset{i=n+1}{\overset{N}{\prod }}}u\left(z_i\mathrm{<mo>\; |\; </mo>}\theta \right)\right]\mathrm{,}$ (1)

where $u\left(z\mathrm{<mo>\; |\; </mo>}\theta \right)$ is the density function of $z$, $i=1,\cdots ,n$ enumerates the second phase complete data, and $i=n+1,\cdots ,N$ enumerates the first phase incomplete data.

If the response $y$ is not involved in the stratification, namely, vector $z$ contains no elements of vector $y$, $u\left(z\mathrm{<mo>\; |\; </mo>}\theta \right)=u\left(z\right)$ is independent of parameters $\theta$, and the full likelihood $L_F\left(\theta \right)$ reduces to ${\displaystyle {\prod }_{i=1}^n}g\left(y_i\mathrm{<mo>\; |\; </mo>}{x}_i\mathrm{<mo>\; ;\; </mo>}\theta \right)$, which is trivial since neither the sampling scheme nor the covariate distribution $h\left(x\right)$ needs to be taken into account. In this paper we consider only the SRS where the response $y$ is involved in stratification, which is often referred to as response-selective stratified sampling (RSSS).

In fisheries surveys, length-stratified age sampling (LSAS) is one of the most popular strategies for sampling the age distribution of a fish population. In the first phase of LSAS a large amount of caught fish of a certain species is measured for length, and classified into length strata (e.g. two centimeters, five centimeters). In the second phase a pre-specified small number of fish are randomly selected from each stratum for age measurement. LSAS is BSS, and since growth models generally describe how length increases as a function of age (i.e. length is the response and age is the covariate), it is response-selective. LSAS has been conducted world-wide for several decades. For example, the Canadian Department of Fisheries Oceans (DFO) conducts annual surveys since the 1970’s and uses LSAS for age sampling for many species such as cod and American plaice. Millions of length-at-age data have been accumulated for each species, which are invaluable for fisheries stock assessment and ocean ecosystem studies. In this paper we focus on BSS, with some of the methods and conclusions also applicable to VPS.

[4] suggested to model the age distribution of fish in a survey using the Gamma distribution. [5] also assumed a Gamma age distribution in their hierarchical model of growth for many fish populations, and they showed that parameter estimates did not change much when a more flexible parametric age distribution was used. [6] showed that a flexible Normal mixture distribution for age distribution is more robust to misspecification of the age distribution. Following these studies, in this paper we focus on the case where a valid parametric covariate distribution model is available. For the examples in the simulation studies and real data analysis, we simply use a Gamma distribution for age so that our comparison among various inference approaches is less influenced by numerical issues related to integrating over a complicated covariate distribution.

This is the motivation of this paper. [1] and [7] gave the complete-data likelihood for VPS, which is based solely on the second phase complete data and can be used when the information is not retained for units not selected for full observation. In this study we would like to derive a likelihood function for BSS requiring only the second phase complete data and the total sample size N, which can be used when the first phase BSS data is not available. Some authors ([e.g. [8]) applied a pseudoconditional likelihood approach [1] to LSAS data. We improved this approach to accommodate the first phase incomplete data and the complexities in fisheries LSAS. We conducted simulation studies to compare the new and existing likelihood and pseudolikelihood approaches that have been used or are conveniently applicable to fisheries LSAS. Our purpose is to identify the approaches with the best performance.

The outline of this paper is as follows. In Section 2 we define notations and review the likelihood and pseudolikelihood approaches relevant to this study. In Section 3 we derive the complete-data density function, complete-data likelihood and full-data likelihood for BSS. Application of an empirical proportion approach, which is an improved version of the pseudoconditional likelihood approach, to BSS is explored in Section 4. Results from simulation studies based on a linear model with between-individual (BI) variation and a Von Bertalanffy growth model with BI variation are presented in Section 5 to compare the performance of all these new and existing estimators discussed in this paper. The most promising estimators are then further illustrated in Section 6 by fitting the VonB model with BI variation to growth data for American plaice (Hippoglossoides platessoides) collected by DFO. Some further discussions are provided in Section 7.

2. Notation, Likelihoods and Pseudolikelihoods

Suppose that N units $\left(y_i\mathrm{,}{x}_i\right)\mathrm{,}i=\mathrm{1,2,}\cdots \mathrm{,}N$, are generated from the joint distribution $f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)$. As mentioned previously we always assume that an appropriate parametric covariate distribution is available, then $\theta$ here includes not only the parameters describing conditional distribution of response $y$ given covariate $x$, but also the parameters defining the covariate distribution. The range of $\left(y\mathrm{,}x\right)$ is divided into H exhaustive and mutually exclusive strata $S_1\mathrm{,}{S}_2\mathrm{,}\cdots \mathrm{,}{S}_H$. Denote the probability for $\left(y\mathrm{,}x\right)$ to fall into the hth stratum as $Q_h\left(\theta \right)$, namely,

$Q_h\left(\theta \right)=\text{Pr}\left\{\left(y\mathrm{,}x\right)\in S_h\right\}\mathrm{.}$ (2)

Define the indicator variable

$R_i=(\begin{array}{ll}\mathrm{1,}\hfill & \text{if}\left(y_i\mathrm{,}{x}_i\right)\text{is}\text{fully}\text{observed}\mathrm{,}\hfill \\ \mathrm{0,}\hfill & \text{if}\text{some}\text{information}\text{on}\left(y_i\mathrm{,}{x}_i\right)\text{is}\text{missing}\mathrm{.}\hfill \end{array}$ (3)

Because BSS2 can be classified as VPS [2], in the following we use BSS specially for BSS1. For BSS we assume that in each stratum $S_h$ there are $N_h$ units from which $n_h\le m_h$ units are randomly selected for full observation of $\left(y\mathrm{,}x\right)$, with $N={\displaystyle {\sum }_{h=1}^H}{N}_h$ and $n={\displaystyle {\sum }_{h=1}^H}{n}_h$. For the remaining $N_h-n_h$ units the values of $\left(y\mathrm{,}x\right)$ are only partially observed for a subset $z$. Here $m_h$ is the maximum sample size for full observation and

$n_h=(\begin{array}{ll}{N}_h,\hfill & \text{if}{N}_h<m_h,\hfill \\ m_h,\hfill & \text{if}{N}_h\ge m_h.\hfill \end{array}$ (4)

Although the likelihood for BSS (4) is given by (1), several published studies use other likelihoods, and some of these are described as follows.

[9] studied maximizing the likelihood function ${\displaystyle {\prod }_{i\mathrm{<mo>\; :\; </mo>}{R}_i=1}}{f}_c\left(y_i\mathrm{<mo>\; |\; </mo>}{x}_i\mathrm{,}{R}_i=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)$ for fitting regression models, and called this approach the conditional maximum likelihood. Under the assumption that a valid parametric covariate distribution is available, and the randomness in $n_h$ can be neglected for all the strata so that the $n_h$ in (4) are always equal to $m_h$ in all strata, in the Appendix we show that

$f_c\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}{R}_i=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\propto \frac{f\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}\theta \right)}{Q_{h_i}\left(\theta \right)}\text{if}\left(y_i\mathrm{,}{x}_i\right)\in S_{h_i}\mathrm{,}$ (5)

and the constant of proportionality does not depend on $\theta$. The conditional likelihood then becomes

$L_c\left(\theta \right)={\displaystyle \underset{i\mathrm{<mo>\; :\; </mo>}{R}_i=1}{\prod }}\frac{f\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}\theta \right)}{Q_{h_i}\left(\theta \right)}\mathrm{,}$ (6)

which is adopted in [10] and [11].

Weighted pseudo-likelihood estimators have been studied extensively since the 1980’s for problems involving response-selective sampling. For this topic we refer to [12] - [18]. In the most basic and popular version of this approach, the log-likelihood function if all the N units were fully observed, ${\displaystyle {\sum }_{i=1}^N}\ln f\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}\theta \right)$ with ln denoting natural logarithm, is estimated by the Horvitz-Thompson (HT) method based on the n units that are actually observed in full,

$l_w\left(\theta \right)={\displaystyle \underset{h=1}{\overset{H}{\sum }}}\frac{N_h}{n_h}{\displaystyle \underset{i=1}{\overset{n_h}{\sum }}}\ln f\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}\theta \right)\mathrm{.}$ (7)

Although this weighted log-pseudo-likelihood (7) may provide an unbiased parameter estimating equation, the HT approach is known to be inefficient, and can be seriously so in some situations such as when the sample unit values are not approximately proportional to the inclusion probabilities ([19]; [20], page 103-104).

An approach for addressing this inefficiency issue is to adjust the standard HT weights by using the whole set of incomplete data, namely, those with only a subset of $\left(y\mathrm{,}x\right)$ measured but available for all the N sample units (see e.g. [17] [18] [21]). We call this method the calibrated weighted likelihood approach. As an implementation of the calibrated weighted likelihood approach, in this study we modified the traditional Horvitz-Thompson weights by minimizing the chi-squared distance (see Equation (1.1) in [21]) between the original and modified weights subject to the constraint

${\displaystyle \underset{i=1}{\overset{N}{\sum }}}{R}_i{w}_i{y}_i={\displaystyle \underset{i=1}{\overset{N}{\sum }}}{y}_i,$ (8)

where $w_i$ are the modified weights. Similarly one can also calibrate up to higher order moments or calibrate the empirical distributions by imposing the constraints

${\displaystyle \underset{j=1}{\overset{N}{\sum }}}{R}_j{w}_j{1}_{y_j\le y_i}={\displaystyle \underset{j=1}{\overset{N}{\sum }}}{1}_{y_j\le y_i},$ (9)

where i enumerates all the subjects selected for full observation, and $1_{y_j\le y_i}=1$ if $y_j\le y_i$ and 0 otherwise. Nevertheless, these calibration strategies may not produce better estimates than (8) does, according to our simulation studies. Hence, in this paper we only report results with constraint (8). The calibrated weighted likelihood approach under all these constraints can be conveniently implemented with Equation (9) in [17].

In some applications (e.g. [8]) researchers use an approximate density based on variable probability sampling (VPS). In VPS, units are randomly selected for full observation from the $N_h$ partially observed units, with subsampling probabilities ${\gamma }_h$ that vary for each stratum h. The density approximation is based on the empirical subsampling probability ${\hat{\gamma }}_h=n_h/N_h$ (see Equation (2) in [2]),

$f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\text{BSS}\mathrm{<mo>\; ;\; </mo>}\theta \right)\approx f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\text{VPS}\mathrm{<mo>\; ;\; </mo>}\theta \right)=\frac{\frac{n_h}{N_h}f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)}{{\displaystyle \underset{h^{\prime }=1}{\overset{H}{\sum }}}\frac{n_{h^{\prime }}}{N_{h^{\prime }}}{Q}_{h^{\prime }}}\mathrm{.}$ (10)

Parameters are estimated based on the likelihood function defined from (10). Note that with the availability of a valid covariate distribution, a density function similar to (10) can also be constructed for the $N-n$ incomplete observations (i.e. those partially observed units). In LSAS there are always some empty strata with $N_h=0$ but non-negligible occupation probability $Q_h$, which are missing in the denominator of (10). We will address these issues in Section 4 and call the improved likelihood the “empirical proportion (EP) likelihood”.

3. Complete-Data Likelihood for BSS

As mentioned previously, the methods for VPS are applicable to BSS2, and the complete-data likelihood for VPS is given in [1] and [7]. Therefore in this section we only consider BSS1 and refer to BSS1 as BSS for convenience.

We denote $\text{dbin}\left(x\mathrm{,}N\mathrm{,}p\right)$ and $\text{pbin}\left(x\mathrm{,}N\mathrm{,}p\right)$ respectively as the binomial probability mass function and cumulative probability, with number of successes x, total number of events N and success probability p. The density function for a unit selected for full observation in BSS is denoted as $f_{BC}(\cdot )$ with “BC” indicating “BSS complete data”.

Theorem 1. In BSS the density function of a unit $\left(y\mathrm{,}x\right)$ selected for full observation is given by

$\begin{array}{l}{f}_{BC}\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)=\frac{f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)}{Q_h}\\ \times \frac{\left[{\displaystyle \underset{N_h=1}{\overset{m_h-1}{\sum }}}{N}_h\text{dbin}\left(N_h\mathrm{,}N\mathrm{,}{Q}_h\right)\right]+m_h\left[1-\text{pbin}\left(m_h-\mathrm{1,}N\mathrm{,}{Q}_h\right)\right]}{{\displaystyle \underset{h^{\prime }=1}{\overset{H}{\sum }}}\left\{\left[{\displaystyle \underset{N_{h^{\prime }}=1}{\overset{m_{h^{\prime }}-1}{\sum }}}{N}_{h^{\prime }}\text{dbin}\left(N_{h^{\prime }}\mathrm{,}N\mathrm{,}{Q}_{h^{\prime }}\right)\right]+m_{h^{\prime }}\left[1-\text{pbin}\left(m_{h^{\prime }}-\mathrm{1,}N\mathrm{,}{Q}_{h^{\prime }}\right)\right]\right\}}\end{array}$ (11)

if $\left(y\mathrm{,}x\right)\in S_h$.

The proof of Theorem 1 is given in the Appendix. As suggested in [9], [10], [22] and [23], the BSS complete-data (BC) likelihood can be constructed as

$L_{BC}={\displaystyle \underset{i\mathrm{<mo>\; :\; </mo>}{R}_i=1}{\prod }}{f}_{BC}\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}{R}_i=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\mathrm{.}$ (12)

With the same arguments for deriving (11), the density function for the partially observed units is

$f_{BI}\left(z\mathrm{<mo>\; |\; </mo>}R=\mathrm{0\; <mo>\; ;\; </mo>}\theta \right)=\frac{f\left(z\mathrm{<mo>\; |\; </mo>}\theta \right)}{Q_h}\times \frac{{\displaystyle \underset{N_h=m_h+1}{\overset{N}{\sum }}}\left(N_h-m_h\right)\text{dbin}\left(N_h\mathrm{,}N\mathrm{,}{Q}_h\right)}{{\displaystyle \underset{h^{\prime }=1}{\overset{H}{\sum }}}\left\{{\displaystyle \underset{N_{h^{\prime }}=m_{h^{\prime }}+1}{\overset{N}{\sum }}}\left(N_{h^{\prime }}-m_{h^{\prime }}\right)\text{dbin}\left(N_{h^{\prime }}\mathrm{,}N\mathrm{,}{Q}_{h^{\prime }}\right)\right\}}\mathrm{,}$ (13)

where the subscript “BI” denotes “BSS incomplete data”. The summations in (13) may be calculated more efficiently using

$\begin{array}{l}{\displaystyle \underset{N_h=m_h+1}{\overset{N}{\sum }}}\left(N_h-m_h\right)\text{dbin}\left(N_h\mathrm{,}N\mathrm{,}{Q}_h\right)\\ =N{Q}_h-{\displaystyle \underset{N_h\mathrm{<mo>\; =\; </mo>1}}{\overset{m_h}{\sum }}}{N}_h\text{dbin}\left(N_h\mathrm{,}N\mathrm{,}{Q}_h\right)-m_h\left[1-\text{pbin}\left(m_h\mathrm{,}N\mathrm{,}{Q}_h\right)\right]\mathrm{.}\end{array}$ (14)

Densities (11) and (13) incorporate respectively the information of complete data and incomplete data. We anticipate that they together can lead to better inference than using only complete data. The BSS full-data (BF) likelihood is

$L_{BF}\left(\theta \right)={\displaystyle \underset{h=1}{\overset{H}{\prod }}}\left[{\displaystyle \underset{i=1}{\overset{n_h}{\prod }}}{f}_{BC}\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}{R}_i=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\right]\left[{\displaystyle \underset{i=n_h+1}{\overset{N_h}{\prod }}}{f}_{BI}\left(z_i\mathrm{<mo>\; |\; </mo>}{R}_i=\mathrm{0\; <mo>\; ;\; </mo>}\theta \right)\right]\mathrm{.}$ (15)

Here and in the remainder of this paper, we enumerate the fully observed units in the hth stratum as $\mathrm{1,}\cdots \mathrm{,}{n}_h$, and the partially observed units in the same stratum as $n_h+\mathrm{1,}\cdots \mathrm{,}{N}_h$.

In some cases only the number of incomplete measurements, $\left(N_h-n_h\right)$, in each stratum are known, instead of the measured values of all $z_i$ ’s. In this situation we need to integrate out $z$ in (13) and rewrite the likelihood function (15) as

$\begin{array}{l}{L}_{BF}\left(\theta \right)={\displaystyle \underset{h=1}{\overset{H}{\prod }}}\left[{\displaystyle \underset{i=1}{\overset{n_h}{\prod }}}{f}_{BC}\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}{R}_i=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\right]\\ \times {\left[\frac{{\displaystyle \underset{N_h=m_h+1}{\overset{N}{\sum }}}\left(N_h-m_h\right)\text{dbin}\left(N_h,N,Q_h\right)}{{\displaystyle \underset{h^{\prime }=1}{\overset{H}{\sum }}}\left\{{\displaystyle \underset{N_{h^{\prime }}=m_{h^{\prime }}+1}{\overset{N}{\sum }}}\left(N_{h^{\prime }}-m_{h^{\prime }}\right)\text{dbin}\left(N_{h^{\prime }},N,Q_{h^{\prime }}\right)\right\}}\right]}^{N_h-n_h}.\end{array}$ (16)

In real data analysis it is important to examine residuals for the fitted model to assess the validity of assumptions. Equation (11) gives the density function for BSS complete data, and can be used to calculate residuals. For simplicity we assume response $y$ to be univariate y. Define the density function of $x$ conditional on $R=1$ as

$h_{BC}\left(x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)={\displaystyle \int f_{BC}\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\text{d}y}$.

$\text{E}\left(y\mathrm{<mo>\; |\; </mo>}x\mathrm{,}R=1\right)=\frac{{\displaystyle \int y{f}_{BC}\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\text{d}y}}{h_{BC}\left(x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)}\mathrm{,}$

$\text{E}\left(y^2\mathrm{<mo>\; |\; </mo>}x\mathrm{,}R=1\right)=\frac{{\displaystyle \int y^2{f}_{BC}\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\text{d}y}}{h_{BC}\left(x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)}\mathrm{,}$

$\text{Var}\left(y\mathrm{<mo>\; |\; </mo>}x\mathrm{,}R=1\right)=\text{E}\left(y^2\mathrm{<mo>\; |\; </mo>}x\mathrm{,}R=1\right)-{\left[\text{E}\left(y\mathrm{<mo>\; |\; </mo>}x\mathrm{,}R=1\right)\right]}^2\mathrm{.}$

The standardized residual for the ith observation $\left(y_i\mathrm{,}{x}_i\right)$ is

$\frac{y_i-\text{E}\left(y\mathrm{<mo>\; |\; </mo>}{x}_i\mathrm{,}R=1\right)}{\sqrt{\text{Var}\left(y\mathrm{<mo>\; |\; </mo>}{x}_i\mathrm{,}R=1\right)}}\mathrm{.}$ (17)

The measured data such as length and age are usually discrete, and the above integrations become summations, which are easier to evaluate.

4. Application of Empirical Proportion Approach to BSS

In this section we expand density (10) for application in BSS and especially in LSAS.

Empty strata ( $N_h=0$ ) always happen with LSAS. For the empty strata in (10), the empirical selection proportions $n_h/N_h\left(=0/0\right)$ are not defined. We need to assign selection probabilities for full and incomplete observations to those unobserved strata. In VPS these selection probabilities may be determined by the maximum likelihood method [10]. For sampling model (4), when $N_h\le m_h$, all the individuals in the hth stratum are selected for full measurement; hence, logically the empirical selection probability is 1 when $N_h=0<m_h$. We assume that in unobserved strata the probability for full observation is 1, and the probability for incomplete observation is 0. Hence, the empirical proportion (EP) density of the complete data with $\left(y\mathrm{,}x\right)$ fully measured is given by

$f_{\text{EP}}\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\mathrm{<mo>\; =\; </mo>}\frac{\frac{n_h}{N_h}f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)}{{\displaystyle \underset{h^{\prime }=1}{\overset{H_{obs}}{\sum }}}\frac{n_{h^{\prime }}}{N_{h^{\prime }}}{Q}_{h^{\prime }}+{\displaystyle \underset{h^{\prime }=H_{obs}+1}{\overset{H_{total}}{\sum }}}{Q}_{h^{\prime }}}\mathrm{.}$ (18)

Here $h=\mathrm{1,}\cdots \mathrm{,}{H}_{obs}$ enumerate the strata with data observed, and $h=H_{obs}+\mathrm{1,}\cdots \mathrm{,}{H}_{total}$ enumerate the strata without data. $H_{total}$ is the total number of strata with nonnegligible occupation probabilities $Q_h$ (see Equation (2)).

Similarly, we can include information from the incomplete observations using their EP density,

$f_{\text{EP}}\left(z\mathrm{<mo>\; |\; </mo>}R=\mathrm{0\; <mo>\; ;\; </mo>}\theta \right)=\frac{\frac{N_h-n_h}{N_h}f\left(z\mathrm{<mo>\; |\; </mo>}\theta \right)}{{\displaystyle \underset{h^{\prime }=1}{\overset{H_{obs}}{\sum }}}\frac{N_{h^{\prime }}-n_{h^{\prime }}}{N_{h^{\prime }}}{Q}_{h^{\prime }}}\mathrm{.}$ (19)

Here, without loss of generality, we assume that $z$ falls in the hth stratum. For an unobserved stratum h, since we have defined its proportion for full observation $n_h/N_h=1$, its proportion for partial observation $\left(N_h-n_h\right)/N_h=0$. The EP likelihood function then has the form

$L_{\text{EP}}\left(\theta \right)={\displaystyle \underset{h=1}{\overset{H}{\prod }}}\left[{\displaystyle \underset{i=1}{\overset{n_h}{\prod }}}{f}_{\text{EP}}\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}{R}_i=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\right]\left[{\displaystyle \underset{i=n_h+1}{\overset{N_h}{\prod }}}{f}_{\text{EP}}\left(z_i\mathrm{<mo>\; |\; </mo>}{R}_i=\mathrm{0\; <mo>\; ;\; </mo>}\theta \right)\right]\mathrm{.}$ (20)

If only the number of incomplete observations in each stratum is reported without knowing the $z$ values, $z$ in (19) needs to be integrated out and the likelihood (20) becomes

$L_{\text{EP}}\left(\theta \right)={\displaystyle \underset{h=1}{\overset{H}{\prod }}}\left[{\displaystyle \underset{i=1}{\overset{n_h}{\prod }}}{f}_{\text{EP}}\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}{R}_i=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\right]{\left[\frac{Q_h}{{\displaystyle \underset{h^{\prime }=1}{\overset{H_{obs}}{\sum }}}\frac{N_{h^{\prime }}-n_{h^{\prime }}}{N_{h^{\prime }}}{Q}_{h^{\prime }}}\right]}^{N_h-n_h}\mathrm{.}$ (21)

5. Simulation Study

In this section we examine the performance of the inference approaches for BSS described in the previous sections. We use two simple examples: a linear model with between individual (BI) variation, and a nonlinear Von Bertalanffy (VonB) growth model with BI variation. The simulation setup is as follows.

5.1. Linear Model with BI Variation

The linear model with BI variation is

$Y=a+BX+\epsilon ,$ (22)

where $B\sim N\left({\mu }_b\mathrm{,}{\sigma }_b^2\right)$, $X\sim N\left({\mu }_x\mathrm{,}{\sigma }_x^2\right)$ and $\epsilon \sim N\left(\mathrm{0,}{\sigma }_{\epsilon }^2\right)$. Capital letter B denotes the random effect of BI variation. We randomly generated $N=5000$ $\left(x_i\mathrm{,}{y}_i\right)$ pairs, $i=1,\cdots ,N$, from model (22). The parameters of the model were chosen as $a=0.5$, ${\mu }_b=0.2,0.5$ and 1.0, ${\sigma }_b=1.0$, ${\mu }_x=1.0$, ${\sigma }_x=5.0$, and ${\sigma }_{\epsilon }=0.7$. Here we selected a small intercept $a$ so that the issues with the relative performance in its estimation as defined by (25) can be clearly seen. Slope is an important parameter in linear model. Hence we selected small, moderate and large values for its mean ${\mu }_b$ and a relatively large standard deviation (SD) ${\sigma }_b$ to test different approaches in identifying the slope under various situations. The mean ${\mu }_x$ and SD ${\sigma }_x$ for covariate X are chosen so that the spread of the covariate allows reasonable estimates of the model parameters. We adopted a moderate error SD ( ${\sigma }_{\epsilon }$ ) relative to the other parameters. We stratified the data by length (Y) bins of size 2 and randomly selected a maximum of 15 units per length stratum to keep their X values, and dropped the X values of the other units not selected. This sampling design is close to the LSAS of fishery surveys that we would like to address in this study.

5.2. VonB Growth Model with BI Variation

The VonB model is a commonly used growth model in fisheries science (e.g. [24]). The basic VonB model is given by

$y\left(a\right)=l_{\infty }\left(1-{\text{e}}^{-k\left(a-a_0\right)}\right)\mathrm{,}$ (23)

where $y\left(a\right)$ denotes length at age $a$, $l_{\infty }$ is the maximum possible size (as $a\to \infty$ ), k is the growth rate parameter, and $a_0\left(<0\right)$ is the theoretical age at which the fish would have had zero length. Variation in growth is also important for population and community dynamics (e.g. [25]). Not accounting for individual variation in growth may lead to bias in estimating the population mean growth parameters and length at age, as noted by [26] and [27]. The VonB model with BI variation follows [11],

$Y=\mu \left(A\right)+\epsilon ,$ (24)

where Y is the measured length, $\mu \left(A\right)=l_{\infty }\left(1-{\text{e}}^{-k\left(A-a_0\right)}\right)$, $A\sim Gamma\left(\alpha \mathrm{,}\beta \right)$ and $\epsilon \sim N\left(\mathrm{0,}{\left[\text{CV}\times \mu \left(A\right)\right]}^2\right)$. The error $\epsilon$ here in fact includes both BI variation and Y observation error.

We randomly generated $N=5000$ ages from a gamma distribution with Case 1: $\left(\alpha ,\beta \right)=\left(3.643,1.225\right)$, and Case 2: $\left(\alpha ,\beta \right)=\left(11.227,0.641\right)$. $\alpha$ and $\beta$ are determined by matching the $\text{mean}=\alpha \beta$ and $\text{variance}=\alpha {\beta }^2$ with those of the age data for American plaice that we have been investigating. Case 1 represents a younger population with mean age = 4.46 and variance = 5.47, while case 2 represents an older population with mean = 7.20 and variance = 4.61. Case 1 has a broad age distribution close to the origin, and case 2 has a narrower distribution of ages. Lengths were then generated from model (24) with $l_{\infty }=70$, $k=0.2$, $a_0=-0.07$ and $\text{CV}=0.2$. We stratified the data by length classes of size 2 and randomly sampled a maximum of 15 units per length stratum to keep their ages and dropped all the other ages not selected.

5.3. Estimation Performance

Relative biases (RBias), relative standard errors (RSE), and relative square root mean squared errors (RRMSE) are defined as

$\begin{array}{l}\text{RBias}=100\times \frac{\text{Estimate}-\text{Truevalue}}{\left|\text{Truevalue}\right|}\mathrm{,}\text{RSE}=100\times \frac{\text{Standarderror}}{\left|\text{Truevalue}\right|}\mathrm{,}\\ \text{and}\text{RRMSE}=100\times \frac{\sqrt{\text{MSE}}}{\left|\text{Truevalue}\right|}\mathrm{.}\end{array}$ (25)

We derived these values using 500 simulations for the full information likelihood (1), conditional likelihood (6), weighted likelihood (7), calibrated weighted likelihood, complete-data likelihood (12), full-data likelihood (15), and EP likelihood (20) (see Tables 1-4). We also include the “random approach” based on maximizing the likelihood

Table 1. Relative bias (RBias), relative standard error (RSE) and relative square root mean squared error (RRMSE) of the estimates from various approaches for the parameters in the linear model with BI variation (22). ${\mu }_b=0.2$.

Table 2. Relative bias (RBias), relative standard error (RSE) and relative square root mean squared error (RRMSE) of the estimates from various approaches for the parameters in the linear model with BI variation (22). ${\mu }_b=0.5$.

Table 3. Relative bias (RBias), relative standard error (RSE) and relative square root mean squared error (RRMSE) of the estimates from various approaches for the parameters in the linear model with BI variation (22). ${\mu }_b=1.0$.

$L_R={\displaystyle \underset{i=1}{\overset{n}{\prod }}}f\left(y_i\mathrm{,}{x}_i\mathrm{<mo>\; |\; </mo>}\theta \right)$ (26)

as a reference point to see the difference between considering BSS and totally ignoring BSS.

For the linear model with BI variation (22), Tables 1-3 indicate that the full information, full-data and EP likelihood approaches have quite close performance, and in general they perform substantially better than all the other approaches in terms of RBias, RSE and RRMSE for all estimated parameters. The weighted likelihood (WL) and calibrated WL approaches have close performance, and there is no evidence that calibration improves the estimation; that is,

Table 4. Relative bias (RBias), relative standard error (RSE) and relative square root mean squared error (RRMSE) of the estimates from various approaches for the parameters in the VonB model with BI variation (24). Case 1: $\left(\alpha ,\beta \right)=\left(3.643,1.225\right)$. Case 2: $\left(\alpha ,\beta \right)=\left(11.227,0.641\right)$.

in some cases the calibrated WL has a little smaller RRMSEs than WL, and in the other cases the reverse happens, but the differences have no clear pattern, and are too small to draw reliable conclusions. Similarly, even though there is some difference in performance between the complete-data likelihood approach and the two WL approaches, it is not clear which method performs better. The two WL approaches have smaller RRMSEs for $a$ and ${\sigma }_{\epsilon }$ estimation, while the complete-data likelihood approach has smaller RRMSEs for other parameter estimation. The conditional likelihood approach based on (6) performs the worst among all the approaches in this study except the random approach. Especially for ${\mu }_x$, ${\sigma }_x$, ${\mu }_b$ and ${\sigma }_b$ estimation, its RRMSEs are more than twice of those from the complete-data likelihood approach. Nevertheless, the conditional likelihood approach performs substantially better than the random approach.

Simulation results presented in Table 4 provide a comparison of the various estimation approaches for a nonlinear VonB model with BI variation. The outcomes for this nonlinear case are similar to the linear case just described. The full information, full-data and EP likelihood approaches have only tiny differences in performance, and in general perform better than the other approaches. The WL and calibrated WL approaches have almost identical performance. The complete-data likelihood approach has close performance as the two WL approaches and it is not clear which method is better. In Case 1 the conditional likelihood approach performs better than the random approach, but worse than all the other approaches including the complete-data approach. In Case 2 its performance is much worse than all the approaches including the random approach. Actually, the conditional likelihood approach failed for this case because it did not converge in 107 of the 500 simulations. All the methods in this study cannot estimate $a_0$ well, with large RBias, large RSEs, and hence large RRMSEs. In practice we suggest to borrow information from other studies such as larvae studies to fix $a_0$, or equivalently to fix length at age 0, for the VonB model.

6. Real Data Analysis

The simulation study indicates that the full information likelihood (1), full-data likelihood (15) and EP likelihood (20) approaches perform better than the other estimation methods. In this section we apply these three approaches to fit the VonB model (24) using a dataset collected by DFO in NAFO Division 3N during the spring of 2011. Here we consider only female American plaice because males and females follow different growth models.

The LSAS within each Division involved measuring the length of all fish caught in research trawl tows, classifying them into 2 cm length strata, and subsampling a few or no otoliths from each length stratum. The sampling goal in each Division was to obtain about 25 age measurements per 2 cm length stratum by sex if $\text{length}\ge 10\text{cm}$, and about 15 age measurements per stratum without sex distinguishment if $\text{length}<10\text{cm}$.

Parameter estimates (ESTs) and the corresponding standard errors (SEs) are provided in Table 5. The three estimation approaches give similar values for all the parameters and SEs, which agrees with their close performance in the simulation study. For comparison, we also included estimates from the random approach (26), which result in a substantially larger value for $l_{\infty }$ and a smaller value for k. The standard errors of the estimates from the random approach are also larger, especially for $l_{\infty }$.

Applying (17), we obtained the standardized residuals of the second phase complete data for all approaches, whose box-and-whisker plots by age are shown in Figure 1. The standardized residuals from the full information likelihood, EP

Figure 1. Box-and-whisker plots of standardized residuals vs. age from fitting the VonB model with BI variation (24) to the American plaice data from DFO 2011 Spring survey in NAFO Division 3N by the four likelihood approaches: full information likelihood (Full information), empirical proportion likelihood (EP), full-data likelihood (Full-data) and random sample assumption based likelihood (Random). The black dots are the medians. The boxes indicate the lower and upper quartiles. The ends of the whiskers represent the lowest datum still within 1.5 IQR (interquartile range) of the lower quartile, and the highest datum still within 1.5 IQR of the upper quartile.

Table 5. Parameter estimates (EST) and standard errors (SE) for the VonB model with between-individual variation (24).

likelihood and full-data likelihood approaches do not indicate bias in fitted mean length at age from the data mean along the full range of age. The standardized residuals from the random approach (26) exhibit clear bias to negative values at ages larger than about 12, indicating over-estimation of $l_{\infty }$. In Figure 1 the interquartile range (IQR, the box) of the residuals is much larger at younger ages (≤4) compared to older ages (≥9). The standard deviation (SD) of the standardized residuals at each age is supposed to be 1. However, the calculated SDs (results not shown) transfers from being greater than 1 at younger ages (≤4) to being mainly smaller than about 0.6 at older ages (≥9). These suggest two problems with the model: 1) the BI variation model in (24) under-estimates the variation at shorter lengths and vice-versa at longer lengths for this data, and 2) due to reproduction, the juvenile female American plaice follows a different growth model from the adult female American plaice, which is neglected by the current model.

7. Discussion

We derived the density function (11) for BSS (basic stratified sampling) complete data, and constructed the complete-data likelihood (12), which allows statistical inference when the incomplete data are not well retained. The complete-data density can also be used for standardized residual calculation as discussed in Section 3. Residuals are important for validation of fitted models.

Both the complete-data likelihood approach and the random approach make use of only the complete data. The complete approach performs substantially better than the random approach in the simulation studies, indicating the importance of correctly incorporating the sampling scheme in the inference methods. The conditional likelihood (6) accounts for the sampling scheme approximately by ignoring the randomness in $n_h$ in all the strata. Therefore its performance lies between the random and the complete-data likelihood approaches in almost all the cases in the simulation study. However in some BSS sampling projects where the number of strata is small and the maximum subsample size $m_h$ for each stratum can usually be obtained, then the conditional likelihood (6) is appropriate.

Another method to incorporate the sampling scheme is to use the count information of the incomplete data in each stratum, as in the weighted likelihood (WL) and calibrated WL approaches. Even though in the simulation study the two methods of accounting for the sampling scheme, namely the complete-data likelihood and the (calibrated) WL approaches, have comparable performance, the complete-data likelihood requires an appropriate distribution model for covariates, which can limit its application. The WL and calibrated WL approaches are not subject to this restriction, and hence can be more practical.

A full utilization of the information in incomplete data is to incorporate the density function of the incomplete data in the likelihood. In this regard, we proposed two new likelihoods for BSS, namely, the full-data likelihood and the empirical proportion (EP) likelihood. If the covariate distribution can be properly modeled, the two new approaches perform as well as the standard full information likelihood approach, and they all perform substantially better than the other methods covered in this study. This result suggests the significance of the information in the incomplete data.

On the whole this study indicates that the complete data, the incomplete data, and the sampling scheme are all important for a consistent and efficient statistical inference from BSS data.

In this work we found that the EP likelihood approach, which was originally proposed for the variable probability sampling (VPS), works well (or the best together with the full-data and full information likelihood approaches) for BSS data. Its merits will further show up when covariates cannot be modeled effectively. This work is under the condition that a valid covariate distribution model is available, which may be a strong assumption in practice. We will explore the case when no appropriate covariate distribution model is available in another paper.

Acknowledgements

Research funding to Nan Zheng was provided by the Ocean Frontier Institute, through an award from the Canada First Research Excellence Fund. Funding was also provided by a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery grant to NC and NC’s Ocean Choice International Industry Research Chair program at the Marine Institute of Memorial University of Newfoundland.

Appendix: Proof of Theorem 1

Without loss of generality, we assume that $\left(y\mathrm{,}x\right)\in S_h$, then

$f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)=\text{Pr}\left(\left(y\mathrm{,}x\right)\in S_h\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\text{Pr}\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\left(y\mathrm{,}x\right)\in S_h\mathrm{,}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\mathrm{.}$

Since the selection for full observation is random given $\left(y\mathrm{,}x\right)\in S_h$,

$\text{Pr}\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\left(y\mathrm{,}x\right)\in S_h\mathrm{,}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)=\text{Pr}\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\left(y\mathrm{,}x\right)\in S_h\mathrm{<mo>\; ;\; </mo>}\theta \right)=\frac{f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)}{Q_h}\mathrm{,}$

and we have

$\begin{array}{l}f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)=\text{Pr}\left(\left(y\mathrm{,}x\right)\in S_h\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\frac{f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)}{Q_h}\\ ={\displaystyle \underset{n_h=0}{\overset{m_h}{\sum }}}\text{Pr}\left(n_h\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\text{Pr}\left(\left(y\mathrm{,}x\right)\in S_h\mathrm{<mo>\; |\; </mo>}{n}_h\mathrm{,}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\frac{f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)}{Q_h}\mathrm{,}\end{array}$ (27)

where $n_h$ is the sample size in the hth stratum as defined by (4).

$\text{Pr}\left(\left(y\mathrm{,}x\right)\in S_h\mathrm{<mo>\; |\; </mo>}{n}_h\mathrm{,}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\propto n_h$, that is, the probability for a selected unit to be in a stratum h is proportional to the number of vacancies in the stratum h. Also, $\text{Pr}\left(n_h\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)=\text{Pr}\left(n_h\mathrm{<mo>\; |\; </mo>}\theta \right)$, namely, the event {a unit is selected without any further information about its $\left(y\mathrm{,}x\right)$ } is independent of the event {there are $n_h$ units that are selected in the stratum h}.

$\text{Pr}\left(n_h\mathrm{<mo>\; |\; </mo>}\theta \right)=(\begin{array}{ll}\text{dbin}\left(N_h\mathrm{,}N\mathrm{,}{Q}_h\right)\mathrm{,}\hfill & \text{if}{N}_h<m_h\text{and}\text{hence}{n}_h=N_h\mathrm{,}\hfill \\ 1-\text{pbin}\left(m_h-\mathrm{1,}N\mathrm{,}{Q}_h\right)\mathrm{,}\hfill & \text{if}{N}_h\ge m_h\text{and}\text{hence}{n}_h=m_h\mathrm{.}\hfill \end{array}$

Hence, when $\left(y\mathrm{,}x\right)\in S_h$,

$\begin{array}{l}f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\propto \frac{f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)}{Q_h}\\ \times \left\{\left[{\displaystyle \underset{N_h=1}{\overset{m_h-1}{\sum }}}{N}_h\text{dbin}\left(N_h\mathrm{,}N\mathrm{,}{Q}_h\right)\right]+m_h\left[1-\text{pbin}\left(m_h-\mathrm{1,}N\mathrm{,}{Q}_h\right)\right]\right\}\mathrm{,}\end{array}$

which can be normalized into (11).

Note that in the case $\text{Pr}\left(n_h=m_h\mathrm{<mo>\; |\; </mo>}\theta \right)=1$ for all the strata $h=1,\cdots ,H$, $\text{Pr}\left(\left(y\mathrm{,}x\right)\in S_h\mathrm{<mo>\; |\; </mo>}{m}_h\mathrm{,}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)=m_h/{\displaystyle {\sum }_{h=1}^H}{m}_h$, which is a constant independent of $\theta$. Then (27) leads to $f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}R=\mathrm{1\; <mo>\; ;\; </mo>}\theta \right)\propto f\left(y\mathrm{,}x\mathrm{<mo>\; |\; </mo>}\theta \right)/Q_h$, which proved (5).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Lawless, J.F., Kalbfleisch, J.D. and Wild, C.J. (1999) Semiparametric Methods for Response-Selective and Missing Data Problems in Regression. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61, 413-438. https://doi.org/10.1111/1467-9868.00185
[2]	Jewell, N.P. (1985) Least Squares Regression with Data Arising from Stratified Samples of the Dependent Variable. Biometrika, 72, 11-21. https://doi.org/10.1093/biomet/72.1.11
[3]	Rubin, D.B. (1976) Inference and Missing Data. Biometrika, 63, 581-592. https://doi.org/10.1093/biomet/63.3.581
[4]	Cope, J.M. and Punt, A.E. (2007) Admitting Ageing Error When Fitting Growth Curves: An Example Using the Von Bertalanffy Growth Function with Random Effects. Canadian Journal of Fisheries and Aquatic Sciences, 64, 205-218. https://doi.org/10.1139/f06-179
[5]	Cadigan, N.G. and Campana, S.E. (2017) Hierarchical Model-Based Estimation of Population Growth Curves for Redfish (Sebastes mentella and Sebastes fasciatus) off the Eastern Coast of Canada. ICES Journal of Marine Science, 74, 687-697. https://doi.org/10.1093/icesjms/fsw195
[6]	Dey, R., Cadigan, N. and Zheng, N. (2019) Estimation of the von Bertalanffy Growth Model When Ages Are Measured with Error. Journal of the Royal Statistical Society: Series C (Applied Statistics), 68, 1131-1147. https://doi.org/10.1111/rssc.12340
[7]	Carroll, R.J., Ruppert, D., Stefanski, L.A. and Crainiceanu, C.M. (2006) Measurement Error in Nonlinear Models: A Modern Perspective. CRC Press, Boca Raton. https://doi.org/10.1201/9781420010138
[8]	Candy, S.G., Constable, A.J., Lamb, T. and Williams, R. (2007) A von Bertalanffy Growth Model for Toothfish at Heard Island Fitted to Length-at-Age Data and Compared to Observed Growth from Mark-Recapture Studies. CCAMLR Science, 14, 43-66.
[9]	Scott, A.J. and Wild, C.J. (2011) Fitting Regression Models with Response-Biased Samples. Canadian Journal of Statistics, 39, 519-536. https://doi.org/10.1002/cjs.10114
[10]	Hausman, J.A. and Wise, D.A. (1982) Stratification on Endogenous Variables and Estimation: The Gray Income Maintenance Experiment. In: Manski, C. and McFadden, D., Eds., Structural Analysis of Discrete Data: With Econometric Applications, Chapter 10, MIT Press, Cambridge, 365-391.
[11]	Piner, K.R., Lee, H.-H. and Maunder, M.N. (2016) Evaluation of Using Random-at-Length Observations and an Equilibrium Approximation of the Population Age Structure in Fitting the von Bertalanffy Growth Function. Fisheries Research, 180, 128-137. https://doi.org/10.1016/j.fishres.2015.05.024
[12]	Hsieh, D.A., Manski, C.F. and McFadden, D. (1985) Estimation of Response Probabilities from Augmented Retrospective Observations. Journal of the American Statistical Association, 80, 651-662. https://doi.org/10.1080/01621459.1985.10478165
[13]	Scott, A.J. and Wild, C.J. (1986) Fitting Logistic Models under Case-Control or Choice Based Sampling. Journal of the Royal Statistical Society. Series B (Methodological), 48, 170-182. https://doi.org/10.1111/j.2517-6161.1986.tb01400.x
[14]	Kalbeisch, J.D. and Lawless, J.F. (1988) Estimation of Reliability in Field-Performance Studies. Technometrics, 30, 365-378. https://doi.org/10.1080/00401706.1988.10488429
[15]	Kalbeisch, J.D. and Lawless, J.F. (1988) Likelihood Analysis of Multi-State Models for Disease Incidence and Mortality. Statistics in Medicine, 7, 149-160. https://doi.org/10.1002/sim.4780070116
[16]	Whittemore, A.S. (1997) Multistage Sampling Designs and Estimating Equations. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59, 589-602. https://doi.org/10.1111/1467-9868.00084
[17]	Breslow, N.E., Lumley, T., Ballantyne, C.M., Chambless, L.E. and Kulich, M. (2009) Improved Horvitz-Thompson Estimation of Model Parameters from Two-Phase Stratified Samples: Applications in Epidemiology. Statistics in Biosciences, 1, 32-49. https://doi.org/10.1007/s12561-009-9001-6
[18]	Saegusa, T. and Wellner, J.A. (2013) Weighted Likelihood Estimation under Two-Phase Sampling. Annals of Statistics, 41, 269. https://doi.org/10.1214/12-AOS1073
[19]	Robins, J.M., Rotnitzky, A. and Zhao, L.P. (1994) Estimation of Regression Coeficients When Some Regressors Are Not Always Observed. Journal of the American Statistical Association, 89, 846-866. https://doi.org/10.1080/01621459.1994.10476818
[20]	Thompson, S. (2012) Sampling. John Wiley & Sons, Hoboken. https://doi.org/10.1002/9781118162934
[21]	Wu, C.B. (2003) Optimal Calibration Estimators in Survey Sampling. Biometrika, 90, 937-951. https://doi.org/10.1093/biomet/90.4.937
[22]	Breslow, N.E. and Cain, K.C. (1988) Logistic Regression for Two-Stage Case-Control Data. Biometrika, 75, 11-20. https://doi.org/10.1093/biomet/75.1.11
[23]	Pfeffermann, D. and Sverchkov, M. (1999) Parametric and Semi-Parametric Estimation of Regression Models Fitted to Survey Data. Sankhyā: The Indian Journal of Statistics, Series B, 61, 166-186.
[24]	Quist, M.C., Pegg, M.A. and DeVries, D.R. (2012) Age and Growth. Fisheries Techniques. 3rd Edition, American Fisheries Society, Bethesda, 677-731.
[25]	Shelton, A.O., Satterthwaite, W.H., Beakes, M.P., Munch, S.B., Sogard, S.M. and Mangel, M. (2013) Separating Intrinsic and Environmental Contributions to Growth and Their Population Consequences. The American Naturalist, 181, 799-814. https://doi.org/10.1086/670198
[26]	Sainsbury, K.J. (1980) Effect of Individual Variability on the von Bertalanffy Growth Equation. Canadian Journal of Fisheries and Aquatic Sciences, 37, 241-247. https://doi.org/10.1139/f80-031
[27]	Wang, Y.-G., Thomas, M.R. and Somers, I.F. (1995) A Maximum Likelihood Approach for Estimating Growth from Tag-Recapture Data. Canadian Journal of Fisheries and Aquatic Sciences, 52, 252-259. https://doi.org/10.1139/f95-025