Family functions in singleRcapture package
Package singleRcapture utilizes various family type functions that specify variable parts of population size estimation, regression, diagnostics and other necessary information that depends on the model. These functions are used as model argument in estimatePopsize function.
chao(lambdaLink = "loghalf", ...) Hurdleztgeom( lambdaLink = c("log", "neglog"), piLink = c("logit", "cloglog", "probit"), ... ) Hurdleztnegbin( nSim = 1000, epsSim = 1e-08, eimStep = 6, lambdaLink = c("log", "neglog"), alphaLink = c("log", "neglog"), piLink = c("logit", "cloglog", "probit"), ... ) Hurdleztpoisson( lambdaLink = c("log", "neglog"), piLink = c("logit", "cloglog", "probit"), ... ) oiztgeom( lambdaLink = c("log", "neglog"), omegaLink = c("logit", "cloglog", "probit"), ... ) oiztnegbin( nSim = 1000, epsSim = 1e-08, eimStep = 6, lambdaLink = c("log", "neglog"), alphaLink = c("log", "neglog"), omegaLink = c("logit", "cloglog", "probit"), ... ) oiztpoisson( lambdaLink = c("log", "neglog"), omegaLink = c("logit", "cloglog", "probit"), ... ) zelterman(lambdaLink = "loghalf", ...) zotgeom(lambdaLink = c("log", "neglog"), ...) zotnegbin( nSim = 1000, epsSim = 1e-08, eimStep = 6, lambdaLink = c("log", "neglog"), alphaLink = c("log", "neglog"), ... ) zotpoisson(lambdaLink = c("log", "neglog"), ...) ztHurdlegeom( lambdaLink = c("log", "neglog"), piLink = c("logit", "cloglog", "probit"), ... ) ztHurdlenegbin( nSim = 1000, epsSim = 1e-08, eimStep = 6, lambdaLink = c("log", "neglog"), alphaLink = c("log", "neglog"), piLink = c("logit", "cloglog", "probit"), ... ) ztHurdlepoisson( lambdaLink = c("log", "neglog"), piLink = c("logit", "cloglog", "probit"), ... ) ztgeom(lambdaLink = c("log", "neglog"), ...) ztnegbin( nSim = 1000, epsSim = 1e-08, eimStep = 6, lambdaLink = c("log", "neglog"), alphaLink = c("log", "neglog"), ... ) ztoigeom( lambdaLink = c("log", "neglog"), omegaLink = c("logit", "cloglog", "probit"), ... ) ztoinegbin( nSim = 1000, epsSim = 1e-08, eimStep = 6, lambdaLink = c("log", "neglog"), alphaLink = c("log", "neglog"), omegaLink = c("logit", "cloglog", "probit"), ... ) ztoipoisson( lambdaLink = c("log", "neglog"), omegaLink = c("logit", "cloglog", "probit"), ... ) ztpoisson(lambdaLink = c("log", "neglog"), ...)
lambdaLink: a link for Poisson parameter, "log"
by default except for zelterman's and chao's models where only \mjseqn \ln \left (\frac x2\right ) is possible.
...: Additional arguments, not used for now.
piLink: a link for probability parameter, "logit" by default
nSim, epsSim: if working weights cannot be computed analytically these arguments specify maximum number of simulations allowed and precision level for finding them numerically respectively.
eimStep: a non negative integer describing how many values should be used at each step of approximation of information matrixes when no analytic solution is available (e.g. "ztnegbin"), default varies depending on a function. Higher value usually means faster convergence but may potentially cause issues with convergence.
alphaLink: a link for dispersion parameter, "log" by default
omegaLink: a link for inflation parameter, "logit" by default
A object of class family containing objects:
makeMinusLogLike -- A factory function for creating the following functions: \mjseqn \ell (\boldsymbol \beta), \frac \partial \ell\partial \boldsymbol \beta, \frac \partial ^2\ell\partial \boldsymbol \beta^T\partial \boldsymbol \beta functions from the \mjseqn \boldsymbol y vector and the \mjseqn \boldsymbol X_vlm matrix (or just \mjseqn \boldsymbol X if applied to model with single linear predictor)which has the deriv argument with possible values in c(0, 1, 2) that determine which derivative to return; the default value is 0, which represents the minus log-likelihood.
links -- A list with link functions.
mu.eta, variance -- Functions of linear predictors that return expected value and variance. The type argument with 2 possible values ("trunc" and "nontrunc") that specifies whether to return \mjseqn \mathbb E(Y|Y>0), \text var(Y|Y>0) or \mjseqn \mathbb E(Y), \text var(Y) respectively; the deriv argument with values in c(0, 1, 2) is used for indicating the derivative with respect to the linear predictors, which is used for providing standard errors in the predict method.
family -- A string that specifies name of the model.
valideta, validmu -- For now it only returns TRUE. In the near future, it will be used to check whether applied linear predictors are valid (i.e. are transformed into some elements of the parameter space subjected to the inverse link function).
funcZ, Wfun -- Functions that create pseudo residuals and working weights used in IRLS algorithm.
devResids -- A function that given the linear predictors prior weights vector and response vector returns deviance residuals. Not all family functions have these functions implemented yet.
pointEst, popVar -- Functions that given prior weights linear predictors and in the latter case also estimate of \mjseqn \text cov(\hat \boldsymbol \beta) and \mjseqn \boldsymbol X_vlm
matrix return point estimate for population size and analytic estimation of its variance.There is a additional boolean parameter contr in the former function that if set to true returns contribution of each unit.
etaNames -- Names of linear predictors.
densityFunction -- A function that given linear predictors returns value of PMF at values x. Additional argument type
specifies whether to return \mjseqn \mathbb P(Y|Y>0) or \mjseqn \mathbb P(Y).
simulate -- A function that generates values of a dependent vector given linear predictors.
getStart -- An expression for generating starting points.
\loadmathjax
Most of these functions are based on some "base" distribution with support \mjseqn \mathbb N_0=\mathbb N\cup \lbrace 0\rbrace that describe distribution of \mjseqn Y before truncation. Currently they include: \mjsdeqn \mathbb P(Y=y|\lambda ,\alpha )=\left \lbrace
\begin arraycc
\frac \lambda ^ye^-\lambday! & \text Poisson distribution
\frac \Gamma (y+\alpha ^-1)\Gamma (\alpha ^-1)y!
\left (\frac \alpha ^-1\alpha ^-1+\lambda\right )^\alpha ^-1
\left (\frac \lambda\alpha ^-1+\lambda\right )^y & \text negative binomial distribution
\frac \lambda ^y(1+\lambda )^y+1 & \text geometric distribution
\end array
\right .
where \mjseqn \lambda is the Poisson parameter and \mjseqn \alpha is the dispersion parameter. Geometric distribution is a special case of negative binomial distribution when \mjseqn \alpha =1 it is included because negative binomial distribution is quite troublesome numerical regression in fitting. It is important to know that PMF of negative binomial distribution approaches the PMF of Poisson distribution when \mjseqn \alpha \rightarrow 0^+.
Note in literature on single source capture recapture models the dispersion parameter which introduces greater variability in negative binomial distribution compared to Poisson distribution is generally interpreted as explaining the unobserved heterogeneity i.e. presence of important unobserved independent variables. All these methods for estimating population size are tied to Poisson processes hence we use \mjseqn \lambda as parameter symbol instead of \mjseqn \mu to emphasize this connection. Also will not be hard to see that all estimators derived from modifying the "base" distribution are unbiased if assumptions made by respective models are not violated.
The zero-truncated models corresponding to "base" distributions are characterized by relation: \mjsdeqn \mathbb P(Y=y|Y>0)=\left \lbrace
\begin arraycc
\frac \mathbb P(Y=y)1-\mathbb P(Y=0) & \text wheny\neq 0
0 & \text wheny=0 \end array\right .
which allows us to estimate parameter values using only observed part of population. These models lead to the following estimates, respectively: \mjsdeqn \begin aligned
\hat N &= \sum _k=1^N_obs\frac 11-\exp (-\lambda _k) & \text For Poisson distribution
\hat N &= \sum _k=1^N_obs\frac 11-(1+\alpha_k\lambda _k)^-\alpha _k^-1 & \text For negative binomial distribution
\hat N &= \sum _k=1^N_obs\frac 1+\lambda _k\lambda _k & \text For geometric distribution
\end aligned
One common way in which assumptions of zero-truncated models are violated is presence of one-inflation the presence of which is somewhat similar in single source capture-recapture models to zero-inflation in usual count data analysis. There are two ways in which one-inflation may be understood, they relate to whether \mjseqn \mathbb P(Y=0) is modified by inflation. The first approach is inflate (\mjseqn \omega parameter) zero-truncated distribution as: \mjsdeqn \mathbb P_new(Y=y|Y>0) = \left \lbrace \begin arraycc
\omega + (1 - \omega )\mathbb P_old(Y=1|Y>0)& \text when: y = 1
(1 - \omega ) \mathbb P_old(Y=y|Y>0) & \text when: y \neq 1 \end array\right .
which corresponds to: \mjsdeqn \mathbb P_new(Y=y) = \left \lbrace \begin arraycc
\mathbb P_old(Y=0) & \text when: y = 0
\omega (1 - \mathbb P(Y=0)) + (1 - \omega )\mathbb P_old(Y=1) & \text when: y = 1
(1 - \omega ) \mathbb P_old(Y=y) & \text when: y > 1 \end array\right .
before zero-truncation. Models that utilize this approach are commonly referred to as zero-truncated one-inflated models. Another way of accommodating one-inflation in SSCR is by putting inflation parameter on base distribution as: \mjsdeqn \mathbb P_new(Y=y) = \left \lbrace \begin arraycc
\omega + (1 - \omega )\mathbb P_old(Y=1)& \text when: y = 1
(1 - \omega ) \mathbb P_old(Y=y) & \text when: y \neq 1 \end array\right .
which then becomes: \mjsdeqn \mathbb P_new(Y=y|Y>0) = \left \lbrace \begin arraycc
\frac \omega1 - (1-\omega )\mathbb P_old(Y=0) + \frac (1 - \omega )1 - (1-\omega )\mathbb P_old(Y=0)\mathbb P_old(Y=1)& \text when: y = 1
\frac (1 - \omega )1 - (1-\omega )\mathbb P_old(Y=0)\mathbb P_old(Y=y) & \text when: y > 1 \end array\right .
after truncation. It was shown by Böhning in 2022 paper that these approaches are equivalent in terms of maximizing likelihoods if we do not put formula on \mjseqn \omega. They can however lead to different population size estimates.
For zero-truncated one-inflated models the formula for population size estimate \mjseqn \hat N does not change since \mjseqn \mathbb P(y=0) remains the same but estimation of parameters changes all calculations.
For one-inflated zero-truncated models population size estimates are expressed, respectively by: \mjsdeqn \begin aligned
\hat N &= \sum _k=1^N_obs\frac 11-(1-\omega _k)\exp (-\lambda _k) &\text For base Poisson distribution
\hat N &= \sum _k=1^N_obs\frac 11-(1-\omega _k)(1+\alpha_k\lambda _k)^-\alpha _k^-1 &\text For base negative binomial distribution
\hat N &= \sum _k=1^N_obs\frac 1+\lambda _k\lambda _k + \omega _k &\text For base geometric distribution
\end aligned
Zero-one-truncated models ignore one counts instead of accommodating one-inflation by utilizing the identity \mjsdeqn \ell _\text ztoi=\boldsymbol f_1\ln \frac \boldsymbol f_1N_obs
+(N_obs-\boldsymbol f_1)\ln \left (1-\frac \boldsymbol f_1N_obs
\right ) + \ell _\text zot
where \mjseqn \ell _\text zot is the log likelihood of zero-one-truncated distribution characterized by probability mass function: \mjsdeqn \mathbb P(Y=y|Y>1)=\left \lbrace
\begin arraycc
\frac \mathbb P(Y=y)1-\mathbb P(Y=0)-\mathbb P(Y=1) & \text wheny > 1
0 & \text wheny\in \lbrace 0, 1\rbrace
\end array\right .
where \mjseqn \mathbb P(Y) is the probability mass function of the "base" distribution. The identity above justifies use of zero-one-truncated, unfortunately it was only proven for intercept only models, however numerical simulations seem to indicate that even if the theorem cannot be extended for (non trivial) regression population size estimation is still possible.
For zero-one-truncated models population size estimates are expressed by: \mjsdeqn \begin aligned
\hat N &= \boldsymbol f_1 + \sum _k=1^N_obs
\frac 1-\lambda _k\exp (-\lambda _k)1-\exp (-\lambda _k)-\lambda _k\exp (-\lambda _k)
&\text For base Poisson distribution
\hat N &= \boldsymbol f_1 + \sum _k=1^N_obs
\frac 1-\lambda _k(1+\alpha _k\lambda _k)^-1-\alpha _k^-11-(1+\alpha _k\lambda _k)^-\alpha _k^-1-\lambda _k(1+\alpha _k\lambda _k)^-1-\alpha _k^-1
&\text For base negative binomial distribution
\hat N &= \boldsymbol f_1 + \sum _k=1^N_obs
\frac \lambda _k^2+\lambda _k+1\lambda _k^2
&\text For base geometric distribution
\end aligned
Pseudo hurdle models are experimental and not yet described in literature.
Lastly there are chao and zelterman models which are based on logistic regression on the dummy variable \mjsdeqn Z = \left \lbrace \begin arraycc
0 & \text ifY = 1
1 & \text ifY = 2 \end array\right .
based on the equation: \mjsdeqn \text logit(p_k)= \ln \left (\frac \lambda _k2\right )= \boldsymbol \beta\mathbf x_k=\eta _k
where \mjseqn \lambda _k is the Poisson parameter.
The zelterman estimator of population size is expressed as: \mjsdeqn \hat N=\sum _k=1^N_obs1-\exp \left (-\lambda _k\right )
and chao estimator has the form: \mjsdeqn \hat N=N_obs+\sum _k=1^\boldsymbol f_1+\boldsymbol f_2
\frac 1\lambda _k+ \frac \lambda _k^22
estimatePopsize()
Piotr Chlebicki, Maciej Beręsewicz
Useful links