install.packages("devtools")
devtools::install_github("PRIMER-e/senlm")The models (mean functions and error distributions) described here are for non-negative data. If it is possible to simulate negative data from these models (e.g. Gaussian error) then it is understood that those data points will be set to zero. Also, if a mean function takes a negative value, then the mean function should be set to zero at that location. Any data simulated from a mean function equal to zero, must be zero.
If the data is Bernoulli or percentage data then the data is bounded inclusively between 0 and 1.
If the data is binomial count data, then then the data is bounded inclusively between 0 and \(n\).
The parameters for the mean function are stored in \(\theta_M\), the parameters for the error distribution are stored in \(\theta_E\). The binomial \(n\) value, and the \(\delta\) parameter of the tail-adjusted beta are stored in \(\theta_C\).
The available mean functions can be printed via the mean_function() function. The mean functions are further classified as being of type test, main, or sub. The test class functions are mainly used for testing code, and the sub class functions are special cases of the main functions.
| mean_fun | mean_class |
|---|---|
| constant | test |
| uniform | test |
| beta | main |
| sech | main |
| modskurt | main |
| gaussian | main |
| mixgaussian | main |
| hofV | main |
| sech.p1 | sub |
| sech.r0p1 | sub |
| mixgaussian.equal | sub |
| hofII | sub |
| hofIV | sub |
| hofIVb | sub |
| hofVb | sub |
All mean functions have an \(H\) parameter determining the maximum height of the mean function. The range of this parameter depends on the the error distribution (see next section). If the error distribution class is binary or percentage, then \(H\) is between 0 and 1; if the error distribution class binomial, then \(H\) is between 0 and \(n\); otherwise \(H\) must be positive.
The mean functions, \(\mu_{\theta_M}(x)\), are listed below in detail, with parameter descriptions, and their classification in parentheses in the title. The get_mean_fun_parnames() output displays the name of the parameters used in the package - sometimes the code uses a simpler roman character rather than a greek spelling e.g. w for \(\omega\) (omega).
The constant mean function.
\[ \begin{aligned} \theta_M &= \{H\} \\ H &= \mbox{Mean response (}H > 0, ~0 < H < 1, ~0 < H < n) \\ \mu_{\theta_M}(x) &= H \\ \end{aligned} \]
This mean function is uniform between \(c\) and \(d\), with constant value between \(c\) and \(d\) of \(H\), and zero elsewhere.
\[ \begin{aligned} \theta_M &= \{H, c, d\} \\ H &= \mbox{Mean response} ~(H > 0 ~~\mbox{or}~~ 0 < H < 1 ~~\mbox{or}~~ 0 < H < n) \\ c &= \mbox{Lower bound} ~(-\infty < c < \infty) \\ d &= \mbox{Upper bound} ~(-\infty < d < \infty, ~c < d) \\ \mu_{\theta_M}(x) &= H \times 1_{[c,d]}(x) \\ &= \left\{ \begin{array}{cl} 0, & x < c \\ H, & c \leq x \leq d \\ 0, & x > d \\ \end{array} \right . \end{aligned} \]
This mean curve is a modified version of the p.d.f. of a beta distribution. A standard beta distribution p.d.f. has the form:
\[ f(x) = \frac{x^{a-1} (1-x)^{b-1}}{B(a,b)}, \]
where \(a,b > 0\), \(0 \leq x \leq 1\), and \(B(a,b)\) is the beta function (normalising constant). The density curve takes different shapes as follows:
| \(a\) | \(b\) | Shape |
|---|---|---|
| \(< 1\) | \(~~~< 1\) | \(U\)-shaped |
| \(1\) | \(~~~1\) | Uniform |
| \(2\) | \(~~~1\) | Triangular (positive) |
| \(1\) | \(~~~2\) | Triangular (negative) |
| \(> 1\) | \(~~~< 1\) | \(J\)-shaped (positive) |
| \(< 1\) | \(~~~> 1\) | \(J\)-shaped (negative) |
| \(> 1\) | \(~~~> 1\) | Unimodal |
We modify this by translating the \(x\)-domain from \((0,1)\) to \((c,d)\), and multiplying by a height parameter \(H\). We also transform the \(a\) and \(b\) parameters to the \(u\) and \(v\) as follows:
\[ u = \frac{a}{a+b}, ~~~v=\frac{1}{a+b}. \]
Under this parameterisation, \(u\) (\(0<u<1\)) is the mean of the beta distribution, and \(v\) (\(v>0\)) is a shape parameter. We restrict ourselves to functions that are not \(U\) or \(J\) shaped (i.e \(a \geq 1\), and \(b \geq 1\)).
\[ \begin{aligned} \theta_M &= \{H, c, d, u, v\} \\ H &= \mbox{Mean response} ~(H > 0 ~~\mbox{or}~~ 0 < H < 1 ~~\mbox{or}~~ 0 < H < n) \\ c &= \mbox{Lower bound} ~(-\infty < c < \infty) \\ d &= \mbox{Upper bound} ~(-\infty < d < \infty, ~c < d) \\ u &= \mbox{Mean of the beta distribution} ~(0 < u < 1) \\ v &= \mbox{Shape parameter} ~(v > 0) \\ ~ & ~ \\ \mu_{\theta_M}(x) &= \frac{H}{H_m}\frac{1}{B(a,b)} \left(\frac{x-c}{d-c}\right)^{a-1} \left(\frac{d-x}{d-c}\right)^{b-1} \end{aligned} \] where \(a=\frac{u}{v}\) and \(b=\frac{1-u}{v}\), and \(H_m\) is the maximum value of the beta distribution.
This mean function is a modified sech function. The basic sech function has the following form:
\[ \mbox{sech} (x) = \frac{2}{e^x + e^{-x}}. \]
This is a symmetric, unimodal function centered around zero. This function is modified as follows.
\[ \begin{aligned} \theta_M &= \{H, m, s, r, p\} \\ H &= \mbox{Mean response} ~(H > 0 ~~\mbox{or}~~ 0 < H < 1 ~~\mbox{or}~~ 0 < H < n) \\ m &= \mbox{Mode} ~(-\infty < m < \infty) \\ s &= \mbox{Spread parameter} ~(s > \infty) \\ r &= \mbox{Skewness parameter} ~(-1 < r < 1) \\ p &= \mbox{Peakedness parameter} ~(p > 0) \\ ~ & ~ \\ \mu_{\theta_M}(x) &= \frac{H}{H_m} \exp\left(\frac{rp}{s}(x - (m-x_m))\right) \left(\mbox{sech}\left(\frac{x-(m-x_m)}{s}\right)\right)^p \end{aligned} \] where \(x_m = \frac{s}{2}\log\left(\frac{1+r}{1-r}\right)\) and \(H_m = \exp\left(\frac{rp}{s}x_m\right)\left(\sqrt{1 - r^2}\right)^p\) The values \(x_m\) and \(H_m\) are used to fix the curve so \(m\) is the location of the mode, and \(H\) is the maximum height of the curve (at \(m\)). If \(r\) is positive the curve is positively skewed.
The modskurt function is a modified sech function that has new parameter (\(b\)) that allows flattening of the peak of the previous \(sech\) functions. The some of the remaining parameters are renamed/reparameterised.
\[ \begin{aligned} \theta_M &= \{H, m, s, q, p, b\} \\ H &= \mbox{Mean response} ~(H > 0 ~~\mbox{or}~~ 0 < H < 1 ~~\mbox{or}~~ 0 < H < n) \\ m &= \mbox{Mode} ~(-\infty < m < \infty) \\ s &= \mbox{Spread parameter} ~(s > \infty) \\ q &= \mbox{Skewness parameter (modified $r$ from sech)} ~(0 < q < 1) \\ p &= \mbox{Peakedness parameter} ~(p > 0) \\ b &= \mbox{Flatness parameter} ~(b > 0) \\ ~ & ~ \\ \mu_{\theta_M}(x) &= H\left(\frac{1}{q\exp{\left(\frac{1}{p}\left(\frac{x-m}{qs}-b\right)\right)} + (1-q)\exp{ \left(-\frac{1}{p}\left(\frac{x-m}{(1-q)s}+b\right)\right) } - \exp{\left(-\frac{b}{p}\right)} + 1}\right)^p \end{aligned} \]
This mean function is a Gaussian curve with a maximum height \(H\), located at mean/mode \(m\), and a standard deviation of \(s\).
\[ \begin{aligned} \theta_M &= \{H, m, s\} \\ H &= \mbox{Maximum value of mean function} ~(H > 0, ~~\mbox{or}~~ 0 < H < 1, ~~\mbox{or}~~ 0 < H < n) \\ m &= \mbox{Location of maximum of mean function} ~(-\infty < m < \infty) \\ s &= \mbox{Standard deviation of mean function} ~(s > 0) \\ ~ & ~ \\ \mu_{\theta_M}(x) &= H \exp\left(-\frac{1}{2}\left(\frac{x-m}{s}\right)^2\right) \end{aligned} \]
This mean function is a mixture of two Gaussian curves. The first component has a maximum height located at mean/mode \(m_1\), with a standard deviation of \(s_1\). Similarly the second component has a maximum height located at mean/mode \(m_2\), with a standard deviation of \(s_2\). The mixture probability is given by \(\alpha\), and the maximum height of the curve is given by \(H\). The condition, \(m_1 \leq m_2\), ensures that component one is always on the left.
\[ \begin{aligned} \theta_M &= \{H, \alpha, m_1, m_2, s_1, s_2\} \\ H &= \mbox{Maximum value of mean function} ~(H > 0, ~~\mbox{or}~~ 0 < H < 1, ~~\mbox{or}~~ 0 < H < n) \\ \alpha &= \mbox{Mixture probability for first component} ~(0 < \alpha < 1) \\ m_i &= \mbox{Mean/mode of $i^{th}$ component} ~(-\infty < m_i< \infty) \\ s_i &= \mbox{Standard deviation of $i^{th}$ component} ~(s_i > 0) \\ ~ & ~ \\ \mu_{\theta_M}(x) &= \frac{H}{H_m}\left\{ \alpha\exp\left(-\frac{1}{2}\left(\frac{x-m_1}{s_1}\right)^2\right) + (1-\alpha)\exp\left(-\frac{1}{2}\left(\frac{x-m_2}{s_2}\right)^2\right) \right\} \end{aligned} \]
where \(H_m\) is the maximum value of the unscaled (\(H=1\)) gaussian mixture curve. This is determined numerically.
This mean function is a the same as the mixgaussian mean function, but with the constraint that \(s_1 = s_2 ~(=s)\).
This mean function is a scaled product of two sigmoid (logistic) growth curves (one increasing and one decreasing). The growth rates of each curve are \(\omega_1\) and \(\omega_2\) , and the points of inflection (half height locations) of each curve occur at \(m − k\) and \(m + k\) respectively.
\[ \begin{aligned} \theta_M &= \{H, m, \omega_1, \omega_2, k\} \\ H &= \mbox{Maximum value of mean function} ~(H > 0, ~~\mbox{or}~~ 0 < H < 1, ~~\mbox{or}~~ 0 < H < n) \\ m &= \mbox{Mid-point between points of inflection} ~(-\infty < m < \infty) \\ \omega_1 &= \mbox{Growth rate parameter of increasing growth curve} ~(\omega_1 > 0) \\ \omega_2 &= \mbox{Growth rate parameter of decreasing growth curve} ~(\omega_2 > 0) \\ k &= \mbox{Half the distance between points of inflections} ~(k > 0) \\ ~ & ~ \\ \mu_{\theta_M}(x) &= \frac{H}{H_m} \left(\frac{1}{1 + \exp{(-\omega_1(x-(m-k))}}\right) \left(\frac{1}{1 + \exp{(\omega_2(x-(m+k))}}\right) \end{aligned} \]
The variable \(H_m\) is the maximum value of the unscaled product of the sigmoidal functions. The value of \(H_m\) is determined numerically.
This mean function is a sigmoid (logistic) growth curve with a maximum height of \(H\), a point of inflection \(m\), and a growth rate parameter of \(\omega_0\).
\[ \begin{aligned} \theta_M &= \{H, m, \omega_0\} \\ H &= \mbox{Maximum value of mean function} ~(H > 0, ~~\mbox{or}~~ 0 < H < 1, ~~\mbox{or}~~ 0 < H < n) \\ m &= \mbox{Point of inflection} ~(-\infty < m < \infty) \\ \omega_0 &= \mbox{Growth rate parameter} ~(-\infty < \omega_0 < \infty) ~ & ~ \\ \mu_{\theta_M}(x) &= H\left(\frac{1}{1 + \exp{(-\omega_0(x-m))}}\right) \end{aligned} \]
This mean function is a scaled product of two sigmoid (logistic) growth curves (one increasing and one decreasing). The growth rates of each curve is \(\omega\), and points of inflection (half height locations) of each curve occur at \(m − k\) and \(m + k\) respectively.
\[ \begin{aligned} \theta_M &= \{H, m, \omega, k\} \\ H &= \mbox{Maximum value of mean function} ~(H > 0, ~~\mbox{or}~~ 0 < H < 1, ~~\mbox{or}~~ 0 < H < n) \\ m &= \mbox{Mid-point between points of inflection} ~(-\infty < m < \infty) \\ \omega &= \mbox{Growth rate parameter} ~(\omega > 0) \\ k &= \mbox{Half the distance between points of inflections} ~(k > 0) \\ ~ & ~ \\ \mu_{\theta_M}(x) &= \frac{H}{H_m} \left(\frac{1}{1 + \exp{(-\omega(x-(m-k))}}\right) \left(\frac{1}{1 + \exp{(\omega(x-(m+k))}}\right) \end{aligned} \]
The variable \(H_m\) is the maximum value of the unscaled product of the sigmoidal functions. \(H_m = (1 + \exp(−\omega k))^{−2}\).
This is the same mean function as the hofIV, except that \(k = 0\).
\[ \begin{aligned} \theta_M &= \{H, m, \omega\} \\ H &= \mbox{Maximum value of mean function} ~(H > 0, ~~\mbox{or}~~ 0 < H < 1, ~~\mbox{or}~~ 0 < H < n) \\ m &= \mbox{Location of points of inflection} ~(-\infty < m < \infty) \\ \omega &= \mbox{Growth rate parameter} ~(\omega > 0) \\ ~ & ~ \\ \mu_{\theta_M}(x) &= \frac{H}{0.25} \left(\frac{1}{1 + \exp{(-\omega(x-m)}}\right) \left(\frac{1}{1 + \exp{(\omega(x-m)}}\right) \end{aligned} \]
This is the same mean function as the hofV, except that \(k = 0\).
\[ \begin{aligned} \theta_M &= \{H, m, \omega_1, \omega_2\} \\ H &= \mbox{Maximum value of mean function} ~(H > 0, ~~\mbox{or}~~ 0 < H < 1, ~~\mbox{or}~~ 0 < H < n) \\ m &= \mbox{Mid-point between points of inflection} ~(-\infty < m < \infty) \\ \omega_1 &= \mbox{Growth rate parameter of increasing growth curve} ~(\omega_1 > 0) \\ \omega_2 &= \mbox{Growth rate parameter of decreasing growth curve} ~(\omega_2 > 0) \\ ~ & ~ \\ \mu_{\theta_M}(x) &= \frac{H}{H_m} \left(\frac{1}{1 + \exp{(-\omega_1(x-m)}}\right) \left(\frac{1}{1 + \exp{(\omega_2(x-m)}}\right) \end{aligned} \]
The error distributions are described in this section. Each error distribution is classified as to the type of data it is suited for. This information is shown in the table below.
| err_dist | err_class |
|---|---|
| bernoulli | binary |
| binomial.count | binomial |
| binomial.prop | binomial |
| poisson | count |
| negbin | count |
| zip | count |
| zinb | count |
| zipl | count |
| zinbl | count |
| zipl.mu | count |
| zinbl.mu | count |
| gaussian | abundance |
| tweedie | abundance |
| zig | abundance |
| zigl | abundance |
| zigl.mu | abundance |
| ziig | abundance |
| ziigl | abundance |
| ziigl.mu | abundance |
| tab | percentage |
| zitab | percentage |
In the sections that follow, the value of mean function at the \(i^{th}\) observation \(x_i\), is denoted by \(\mu_i\). In other words:
\[ \mu_i = \mu_{\theta_M} (x_i). \]
The Bernouli distribution is used for modelling binary data (\(0/1\)). The mean function is bounded between 0 and 1. There are no extra parameters for this distribution.
\[ \begin{aligned} y_i &\in \{0, 1\} \\ \theta_E &= \{\} \\ P (Y_i = y_i \mid \theta_M, ~\theta_E , ~x_i) &\sim \mbox{Bernoulli}(\mu_i) \\ \end{aligned} \]
The binomial distribution is used for modelling count data that is bounded between 0 and \(n\). The mean function is bounded between 0 and \(n\). There are no extra parameters for this distribution.
\[ \begin{aligned} y_i &\in \{0, 1, \ldots, n\} \\ \theta_E &= \{\} \\ \theta_C &= \{n\} \\ P (Y_i = y_i \mid ~\theta_M , ~\theta_C , ~x_i) &\sim \mbox{Binomial}(n, ~\mu_i) \\ \end{aligned} \]
This is that same model as the binomial.count distribution, except the data recorded as proportions (counts divided by \(n\)).
\[ \begin{aligned} y_i &\in \{\frac{0}{n}, \frac{1}{n}, \ldots, \frac{n}{n}\} \\ \theta_E &= \{\} \\ \theta_C &= \{n\} \\ P (Y_i = ny_i \mid ~\theta_M , ~\theta_C , ~x_i) &\sim \mbox{Binomial}(n, ~\mu_i) \\ \end{aligned} \]
The Poisson distribution has mean parameter \(\lambda=\mu_{\theta_M} (x))\). There are no extra parameters for this distribution.
\[ \begin{aligned} y_i &\in \{0, 1, 2, \ldots \} \\ \theta_E &= \{\} \\ P (Y_i = y_i \mid \theta_M , ~\theta_E , ~x_i) &\sim \mbox{Poisson}(\lambda_i=\mu_i) \\ \end{aligned} \]
The negative-binomial distribution we use is actually a Gamma-Poisson distribution with the following form:
\[ \begin{aligned} y_i &\in \{0, 1, 2, \ldots \} \\ \lambda_i &\sim \mbox{Gamma}(\mbox{shape}=\alpha, ~\mbox{rate}=\beta_i),\\ P(Y_i=y_i \mid \lambda_i) &\sim \mbox{Poisson}(\lambda_i) \\ ~ & ~ \\ \mbox{E}(\lambda_i) &= \frac{\alpha}{\beta_i} \\ \mbox{var}(\lambda_i) &= \frac{\alpha}{\beta_i^2}, \\ \mbox{E}(Y_i) &= \frac{\alpha}{\beta_i} + \frac{\alpha}{\beta_i^2} = \mu_i \\ \mbox{var}(Y_i) &= \frac{\alpha}{\beta_i^2} = \mu_i + \frac{\mu_i^2}{\alpha} = \mu_i + \phi \mu_i^2 \\ \end{aligned} \] where \(\phi = \dfrac{1}{\alpha}\) is the dispersion parameter (\(\phi=0\) corresponds to a Poisson distribution), and \(\mu_i\) is the mean parameter.
\[ \begin{aligned} y_i &\in \{0, 1, 2, \ldots \} \\ \theta_E &= \{\phi\} \\ \phi &= \mbox{Dispersion parameter}, ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{Gamma-Poisson}(\mu_i, ~\phi) \end{aligned} \]
The zero-inflated Poisson (ZIP) distribution has zero-spike parameter \(\pi\), and mean parameter \(\mu_i\).
\[ \begin{aligned} y_i &\in \{0, 1, 2, \ldots \} \\ \theta_E &= \{\pi\} \\ \pi &= \mbox{Zero-spike parameter}, (0 \leq \pi \leq 1) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZIP}(\mu_i, ~\pi) \\ &\sim \left\{\begin{array}{ll} \mbox{Poisson}(\mu_i), & \mbox{w.p.} ~~1-\pi, \\ 0, & \mbox{w.p.} ~~\pi. \end{array}\right. \end{aligned} \]
The zero-inflated negative-binomial (ZINB) distribution has zero-spike parameter \(\pi\), dispersion parameter \(\phi\), and mean parameter \(\mu_i\).
\[ \begin{aligned} y_i &\in \{0, 1, 2, \ldots \} \\ \theta_E &= \{\pi, \phi\} \\ \pi &= \mbox{Zero-spike parameter} ~(0 \leq \pi \leq 1) \\ \phi &= \mbox{Dispersion parameter} ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZINB}(\mu_i, ~\pi, ~\phi) \\ &\sim \left\{\begin{array}{ll} \mbox{Gamma-Poisson}(\mu_i, ~\phi), & \mbox{w.p.} ~~1-\pi, \\ 0, & \mbox{w.p.} ~~\pi. \end{array}\right. \end{aligned} \]
The zero-inflated Poisson linked (ZIPL) distribution is a Poisson distribution with mean parameter \(\mu_i\), with a zero spike that has parameters \(\gamma_0\) and \(\gamma_1\). The zero-spike parameters are linked to the mean function on the logit/log scale as follows:
\[ \mbox{logit}(\pi_i) = \gamma_0 + \gamma_1\log(\mu_i), \]
where \(\pi_i\) is the value of the zero-spike parameter at \(x_i\).
\[ \begin{aligned} y_i &\in \{0, 1, 2, \ldots \} \\ \theta_E &= \{\gamma_0, \gamma_1\} \\ \gamma_0 &= \mbox{Zero-spike parameter intercept} ~(-\infty \leq \gamma_0 \leq \infty) \\ \gamma_1 &= \mbox{Zero-spike parameter slope} ~(-\infty \leq \gamma_1 \leq \infty) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZIPL}(\mu_i, ~\gamma_0, ~\gamma_1) \\ &\sim \left\{\begin{array}{ll} \mbox{Poisson}(\mu_i), & \mbox{w.p.} ~~1-\pi_i, \\ 0, & \mbox{w.p.} ~~\pi_i. \end{array}\right. \end{aligned} \]
The zero-inflated negative-binomial linked (ZINBL) distribution is a negative-binomial distribution with mean parameter \(\mu_i\), dispersion parameter \(\phi\), and has zero-spike parameters \(\gamma_0\) and \(\gamma_1\). The zero-spike parameters are linked to the mean function on the logit/log scale as follows:
\[ \mbox{logit}(\pi_i) = \gamma_0 + \gamma_1\log(\mu_i), \]
where \(\pi_i\) is the value of the zero-spike parameter at \(x_i\).
\[ \begin{aligned} y_i &\in \{0, 1, 2, \ldots \} \\ \theta_E &= \{\gamma_0, \gamma_1, \phi\} \\ \gamma_0 &= \mbox{Zero-spike parameter intercept} ~(-\infty \leq \gamma_0 \leq \infty) \\ \gamma_1 &= \mbox{Zero-spike parameter slope} ~(-\infty \leq \gamma_1 \leq \infty) \\ \phi &= \mbox{Dispersion parameter} ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZINBL}(\mu_i, ~\gamma_0, ~\gamma_1, ~\phi) \\ &\sim \left\{\begin{array}{ll} \mbox{Gamma-Poisson}(\mu_i, ~\phi), & \mbox{w.p.} ~~1-\pi_i, \\ 0, & \mbox{w.p.} ~~\pi_i. \end{array}\right. \end{aligned} \]
The zero-inflated Poisson linked-\(\mu\) (ZIPL\(-\mu\)) distribution is a Poisson distribution with mean parameter \(\mu_i\), with a zero spike that has parameters \(\gamma_0\) and \(\gamma_1\). The zero-spike parameters are linked to the mean function on the logit/linear scale as follows:
\[ \mbox{logit}(\pi_i) = \gamma_0 + \gamma_1\mu_i, \]
where \(\pi_i\) is the value of the zero-spike parameter at \(x_i\).
\[ \begin{aligned} y_i &\in \{0, 1, 2, \ldots \} \\ \theta_E &= \{\gamma_0, \gamma_1\} \\ \gamma_0 &= \mbox{Zero-spike parameter intercept} ~(-\infty \leq \gamma_0 \leq \infty) \\ \gamma_1 &= \mbox{Zero-spike parameter slope} ~(-\infty \leq \gamma_1 \leq \infty) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZIPL-}\mu(\mu_i, ~\gamma_0, ~\gamma_1) \\ &\sim \left\{\begin{array}{ll} \mbox{Poisson}(\mu_i), & \mbox{w.p.} ~~1-\pi_i, \\ 0, & \mbox{w.p.} ~~\pi_i. \end{array}\right. \end{aligned} \]
The zero-inflated negative-binomial linked-\(\mu\) (ZINBL-\(\mu\)) distribution is a negative-binomial distribution with mean parameter \(\mu_i\), dispersion parameter \(\phi\), and has zero-spike parameters \(\gamma_0\) and \(\gamma_1\). The zero-spike parameters are linked to the mean function on the logit/linear scale as follows:
\[ \mbox{logit}(\pi_i) = \gamma_0 + \gamma_1\mu_i, \]
where \(\pi_i\) is the value of the zero-spike parameter at \(x_i\).
\[ \begin{aligned} y_i &\in \{0, 1, 2, \ldots \} \\ \theta_E &= \{\gamma_0, \gamma_1, \phi\} \\ \gamma_0 &= \mbox{Zero-spike parameter intercept} ~(-\infty \leq \gamma_0 \leq \infty) \\ \gamma_1 &= \mbox{Zero-spike parameter slope} ~(-\infty \leq \gamma_1 \leq \infty) \\ \phi &= \mbox{Dispersion parameter} ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZINB-}\mu(\mu_i, ~\gamma_0, ~\gamma_1, ~\phi) \\ &\sim \left\{\begin{array}{ll} \mbox{Gamma-Poisson}(\mu_i, ~\phi), & \mbox{w.p.} ~~1-\pi_i, \\ 0, & \mbox{w.p.} ~~\pi_i. \end{array}\right. \end{aligned} \]
The Gaussian distribution has mean parameter \(\mu_i\) and standard deviation \(\sigma\).
\[ \begin{aligned} y_i &\geq 0 \\ \theta_E &= \{\sigma\} \\ \sigma &= \mbox{Standard deviation} ~(\sigma > 0) \\ P(Y_i=y_i \mid ~\theta_M, ~\theta_E, ~x_i) &\sim \mbox{Normal}(\mu_i, ~\sigma^2) \end{aligned} \] Since our data is non-negative, any data simulated via this model will have negative values set to zero.
Tweedie distributions are a class of exponential dispersion models where the variance is a function of the mean of the form:
\[ \mbox{var}(Y_i) = \phi\mu_i^\rho,\]
where \(\rho \in (-\infty, 0] \cup [1, \infty)\). Special cases include the normal (\(\rho = 0\)), Poisson (\(\rho = 1\)), Gamma (\(\rho = 2\)) and inverse Gaussian (\(\rho = 3\)) distributions. For \(1 < \rho < 2\), the distributions are continuous for \(Y > 0\) and have a discrete mass at \(Y = 0\). We will restrict our attention to \(1 < \rho < 2\).
\[ \begin{aligned} y_i &\geq 0 \\ \theta_E &= \{\phi, ~\rho\} \\ \phi &= \mbox{Dispersion parameter} ~(\phi > 0) \\ \rho &= \mbox{Power parameter} ~(1 < \rho < 2) \\ P(Y_i=y_i \mid ~\theta_M, ~\theta_E, ~x_i) &\sim \mbox{Tweedie}(\mu_i, ~\phi, ~\rho) \\ \mbox{E}(Y_i) &= \mu_i \\ \mbox{var}(Y_i) &= \phi(\mu_i)^\rho \end{aligned} \]
The zero-inflated Gamma (ZIG) distribution has zero-spike parameter \(\pi\), dispersion parameter \(\phi\), and mean parameter \(\mu_i\).
\[ \begin{aligned} y_i &\geq 0 \\ \theta_E &= \{\pi, ~\phi\} \\ \pi &= \mbox{Zero-spike parameter}, ~(0 \leq \pi \leq 1) \\ \phi &= \mbox{Dispersion parameter}, ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZIG}(\mu_i, ~\pi, ~\phi) \\ &\sim \left\{\begin{array}{ll} \mbox{Gamma}(\mbox{shape}=1/\phi, ~\mbox{scale} = \phi\mu_i) & \mbox{w.p.} ~~1-\pi, \\ 0, & \mbox{w.p.} ~~\pi. \end{array}\right.\\ \mbox{E}(Y_i) &= (1 - \pi)\mu_i \\ \mbox{var}(Y_i) &= \pi(1-\pi)\mu^2 + (1-\pi)\phi\mu_i^2 \end{aligned} \]
The zero-inflated gamma linked (ZIGL) distribution is a Gamma distribution with mean parameter \(\mu_i\), a dispersion parameter \(\phi\), and a zero spike that has parameters \(\gamma_0\) and \(\gamma_1\). The zero-spike parameters are linked to the mean function on the logit/log scale as follows:
\[ \mbox{logit}(\pi_i) = \gamma_0 + \gamma_1\log(\mu_i), \]
where \(\pi_i\) is the value of the zero-spike parameter at \(x_i\).
\[ \begin{aligned} y_i &\geq 0 \\ \theta_E &= \{\gamma_0, ~\gamma_1, ~\phi\} \\ \gamma_0 &= \mbox{Zero-spike parameter intercept} ~(-\infty \leq \gamma_0 \leq \infty) \\ \gamma_1 &= \mbox{Zero-spike parameter slope} ~(-\infty \leq \gamma_1 \leq \infty) \\ \phi &= \mbox{Dispersion parameter}, ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZIGL}(\mu_i, ~\gamma_0, ~\gamma_1, ~\phi) \\ &\sim \left\{\begin{array}{ll} \mbox{Gamma}(\mbox{shape}=1/\phi, ~\mbox{scale} = \phi\mu_i) & \mbox{w.p.} ~~1-\pi_i, \\ 0, & \mbox{w.p.} ~~\pi_i. \end{array}\right.\\ \mbox{E}(Y_i) &= (1 - \pi)\mu_i \\ \mbox{var}(Y_i) &= \pi(1-\pi)\mu^2 + (1-\pi)\phi\mu_i^2 \end{aligned} \]
The zero-inflated gamma linked-\(\mu\) (ZIGL-\(\mu\)) distribution is a Gamma distribution with mean parameter \(\mu_i\), a dispersion parameter \(\phi\), and a zero spike that has parameters \(\gamma_0\) and \(\gamma_1\). The zero-spike parameters are linked to the mean function on the logit/linear scale as follows:
\[ \mbox{logit}(\pi_i) = \gamma_0 + \gamma_1\mu_i, \]
where \(\pi_i\) is the value of the zero-spike parameter at \(x_i\).
\[ \begin{aligned} y_i &\geq 0 \\ \theta_E &= \{\gamma_0, ~\gamma_1, ~\phi\} \\ \gamma_0 &= \mbox{Zero-spike parameter intercept} ~(-\infty \leq \gamma_0 \leq \infty) \\ \gamma_1 &= \mbox{Zero-spike parameter slope} ~(-\infty \leq \gamma_1 \leq \infty) \\ \phi &= \mbox{Dispersion parameter}, ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZIGL-}\mu(\mu_i, ~\gamma_0, ~\gamma_1, ~\phi) \\ &\sim \left\{\begin{array}{ll} \mbox{Gamma}(\mbox{shape}=1/\phi, ~\mbox{scale} = \phi\mu_i) & \mbox{w.p.} ~~1-\pi_i, \\ 0, & \mbox{w.p.} ~~\pi_i. \end{array}\right.\\ \mbox{E}(Y_i) &= (1 - \pi)\mu_i \\ \mbox{var}(Y_i) &= \pi(1-\pi)\mu^2 + (1-\pi)\phi\mu_i^2 \end{aligned} \]
The zero-inflated inverse-Gaussian (ZIIG) distribution has zero-spike parameter \(\pi\), dispersion parameter \(\phi\), and mean parameter \(\mu_i\).
\[ \begin{aligned} y_i &\geq 0 \\ \theta_E &= \{\pi, ~\phi\} \\ \pi &= \mbox{Zero-spike parameter}, ~(0 \leq \pi \leq 1) \\ \phi &= \mbox{Dispersion parameter}, ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZIIG}(\mu_i, ~\pi, ~\phi) \\ &\sim \left\{\begin{array}{ll} \mbox{Inverse-Gaussian}(\mu_i, \phi) & \mbox{w.p.} ~~1-\pi, \\ 0, & \mbox{w.p.} ~~\pi. \end{array}\right.\\ \mbox{E}(Y_i) &= (1 - \pi)\mu_i \\ \mbox{var}(Y_i) &= \pi(1 - \pi)\mu_i^2 + (1-\pi)\phi\mu_i^3 \end{aligned} \]
The zero-inflated inverse-Gaussian linked (ZIIGL) distribution is an inverse-Gaussian distribution with mean parameter \(\mu_i\), a dispersion parameter \(\phi\), and a zero spike that has parameters \(\gamma_0\) and \(\gamma_1\). The zero-spike parameters are linked to the mean function on the logit/log scale as follows:
\[ \mbox{logit}(\pi_i) = \gamma_0 + \gamma_1\log(\mu_i), \]
where \(\pi_i\) is the value of the zero-spike parameter at \(x_i\).
\[ \begin{aligned} y_i &\geq 0 \\ \theta_E &= \{\gamma_0, ~\gamma_1, ~\phi\} \\ \gamma_0 &= \mbox{Zero-spike parameter intercept} ~(-\infty \leq \gamma_0 \leq \infty) \\ \gamma_1 &= \mbox{Zero-spike parameter slope} ~(-\infty \leq \gamma_1 \leq \infty) \\ \phi &= \mbox{Dispersion parameter}, ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZIIGL}(\mu_i, ~\gamma_0, ~\gamma_1, ~\phi) \\ &\sim \left\{\begin{array}{ll} \mbox{Gamma}(\mbox{shape}=1/\phi, ~\mbox{scale} = \phi\mu_i) & \mbox{w.p.} ~~1-\pi_i, \\ 0, & \mbox{w.p.} ~~\pi_i. \end{array}\right.\\ \mbox{E}(Y_i) &= (1 - \pi)\mu_i \\ \mbox{var}(Y_i) &= \pi(1-\pi)\mu^2 + (1-\pi)\phi\mu_i^3 \end{aligned} \]
The zero-inflated inverse-Gaussian linked-\(\mu\) (ZIIGL-\(\mu\)) distribution is an inverse-Gaussian distribution with mean parameter \(\mu_i\), a dispersion parameter \(\phi\), and a zero spike that has parameters \(\gamma_0\) and \(\gamma_1\). The zero-spike parameters are linked to the mean function on the logit/linear scale as follows:
\[ \mbox{logit}(\pi_i) = \gamma_0 + \gamma_1\mu_i, \]
where \(\pi_i\) is the value of the zero-spike parameter at \(x_i\).
\[ \begin{aligned} y_i &\geq 0 \\ \theta_E &= \{\gamma_0, ~\gamma_1, ~\phi\} \\ \gamma_0 &= \mbox{Zero-spike parameter intercept} ~(-\infty \leq \gamma_0 \leq \infty) \\ \gamma_1 &= \mbox{Zero-spike parameter slope} ~(-\infty \leq \gamma_1 \leq \infty) \\ \phi &= \mbox{Dispersion parameter}, ~(\phi > 0) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~x_i) &\sim \mbox{ZIIGL-}\mu(\mu_i, ~\gamma_0, ~\gamma_1, ~\phi) \\ &\sim \left\{\begin{array}{ll} \mbox{Gamma}(\mbox{shape}=1/\phi, ~\mbox{scale} = \phi\mu_i) & \mbox{w.p.} ~~1-\pi_i, \\ 0, & \mbox{w.p.} ~~\pi_i. \end{array}\right.\\ \mbox{E}(Y_i) &= (1 - \pi)\mu_i \\ \mbox{var}(Y_i) &= \pi(1-\pi)\mu^2 + (1-\pi)\phi\mu_i^3 \end{aligned} \]
The tail-adjusted beta (TAB) distribution is a beta distribution where values less than \(\delta\) are rounded down to zero, and values greater than \(1-\delta\) are rounded up to one. This is used to model data that is between 0 and 1 (inclusive).
\[ \begin{aligned} y_i &\in [0,1] \\ \theta_E &= \{\sigma\} \\ \theta_C &= \{\delta\} \\ \sigma &= \mbox{Shape parameter}, (\sigma > 0) \\ \delta &= \mbox{Rounding parameter}, (0 < \delta < \frac{1}{2}) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~\theta_C, ~x_i) &\sim \mbox{TAB}(\mu_i, ~\sigma, ~\delta) \\ &\sim \left\{\begin{array}{ll} \frac{1}{\delta}p^{(0)}_i & y_i < \delta \\ \dfrac{y_i^{\alpha_i-1}(1-y_i)^{\beta-1}}{B(\alpha,\beta)} & \delta < y_i < 1 - \delta \\ \frac{1}{\delta}p^{(1)}_i & y_i > 1 - \delta \\ \end{array}\right.\\ \end{aligned} \]
where \(p^{(0)}_i\) is the probability the Beta(\(\alpha_i, ~\beta)\) distribution is less than \(\delta\), and \(p^{(1)}_i\) is the probability the Beta(\(\alpha_i, ~\beta)\) distribution is greater than \(1 - \delta\). The parameters of the TAB (\(\mu_i\), \(\sigma\)) are converted to beta distribution parameters (\(\alpha_i\), \(\beta\)) and back again as follows: \[ \begin{aligned} \mu_i &= \dfrac{\alpha_i}{\alpha_i + \beta}, &~ \sigma &= \dfrac{1}{\alpha_i + \beta} \\ \alpha & = \dfrac{\mu_i}{\sigma}, &~ \beta &= \dfrac{1 - \mu_i}{\sigma} \end{aligned} \]
We restrict ourselves to \(\alpha_i \geq 1\) and \(\beta \geq 1\) (uniform, triangular, and unimodal shapes).
The \(\delta\) parameter is not estimated via maximum likelihood like other parameters, it is estimated as the minimum of: the smallest non-zero \(y_i\), and one minus the largest non-one \(y_i\). In other words:
\[ \hat{\delta} = \min( \min_i(y_i : y_i>0), ~\min_i(1-y_i : y_i < 1) ). \]
The zero-inflated tail-adjusted Beta (ZITAB) distribution is the tail-adjusted beta distribution with a zero-inflated spike added.
\[ \begin{aligned} y_i &\in [0,1] \\ \theta_E &= \{\pi, ~\sigma\} \\ \theta_C &= \{\delta\} \\ \pi &= \mbox{Zero-spike parameter}, ~(0 \leq \pi \leq 1) \\ \sigma &= \mbox{Shape parameter}, (\sigma > 0) \\ \delta &= \mbox{Rounding parameter}, (0 < \delta < \frac{1}{2}) \\ P(Y_i=y_i \mid \theta_M, ~\theta_E, ~\theta_C, ~x_i) &\sim \mbox{ZITAB}(\mu_i, ~\pi, ~\sigma, ~\delta) \\ &\sim \left\{\begin{array}{ll} \frac{1}{\delta}\left[\pi + (1-\pi)p^{(0)}_i\right] & y_i < \delta \\ (1-\pi)\dfrac{y_i^{\alpha_i-1}(1-y_i)^{\beta-1}}{B(\alpha,\beta)} & \delta < y_i < 1 - \delta \\ \frac{1}{\delta}(1-\pi)p^{(1)}_i & y_i > 1 - \delta \\ \end{array}\right.\\ \end{aligned} \]
The set_models() function is used to create a data frame that contains the mean functions, error distributions, and if appropriate the binomial \(n\) parameter, of all models that one would wish to fit or simulate data from.
## All possible combinations of mean functions and error distrbiutions
set_models (mean_fun=c("gaussian", "beta"), err_dist=c("zinb", "zinbl"),
method="crossed")
mean_fun err_dist binomial_n
1 beta zinb NA
2 beta zinbl NA
3 gaussian zinb NA
4 gaussian zinbl NA
## Fit specific models
set_models (mean_fun=c("gaussian", "beta"), err_dist=c("zinb", "zinbl"),
method="paired")
mean_fun err_dist binomial_n
1 gaussian zinb NA
2 beta zinbl NA
## Fit binomial model - with n=40
set_models (mean_fun=c("gaussian"),
err_dist=c("poisson", "binomial.prop", "binomial.count"),
binomial_n=40, method="paired")
mean_fun err_dist binomial_n
1 gaussian poisson NA
2 gaussian binomial.prop 40
3 gaussian binomial.count 40
## Can also specify via mean and error class
set_models (mean_class="main", err_class="percentage")
mean_fun err_dist binomial_n
1 beta tab NA
2 beta zitab NA
3 gaussian tab NA
4 gaussian zitab NA
5 hofV tab NA
6 hofV zitab NA
7 mixgaussian tab NA
8 mixgaussian zitab NA
9 modskurt tab NA
10 modskurt zitab NA
11 sech tab NA
12 sech zitab NAThis section describes how to simulate data sets for any possible model. The create_default_par_list() function set default parameter values for each selected model.
## --- Set models
Models <- set_models (mean_fun=c("gaussian"), err_dist=c("zip", "zinb"))
## --- Set default parameters list
Pars <- create_default_par_list (Models)
## --- Simulate data
Data <- create_simulated_datasets (Pars, N=200, xmin=0, xmax=100, seed=12345)This code creates simulates data for the gaussian-zip, and gaussian-zinb models. The objects Pars is a list containing parameters for each model. The object Data is a list containing the data for each model; Data[[1]] is the data for the first model, etc.
The senlm packages contains the following data sets:
The first few rows of data are shown below.
data(haul, package="senlm")
head (haul)
year depth lattitude Albatrossia.pectoralis Sebastolobus.altivelis
1717 2003 527.5357 46.30425 0 37
1718 2003 673.7497 46.73217 0 397
1719 2003 276.1650 46.71926 0 0
1720 2003 823.6656 46.92419 3 711
1721 2003 172.6837 47.01439 0 0
1722 2003 985.7053 47.56385 18 534To fit one model to a data set we use the senlm() function. We illustrate this with the haul dataset.
Fit <- senlm (data=haul, xvar="depth", yvar="Sebastolobus.altivelis",
mean_fun="gaussian", err_dist="zip", conf.level=0.95)
summary (Fit)
$model
[1] "gaussian_zip"
$theta
H m s pi
817.5307540 917.8020537 225.5627466 0.5954255
$IC
npar nll AIC AICc BIC
4.00 36950.05 73908.10 73908.14 73928.01 Confidence intervals for the estimated parameters can be display as follows. The default confidence level is 95% which is set when fitting the model. Occasionly due to numerical instability of the estimated Hessian matrix the standard errors are not calculated.
rbind (lower=Fit$lb, upper=Fit$ub)
H m s pi
lower 812.7322 916.2962 224.2579 0.5615009
upper 822.3576 919.3079 226.8752 0.6284621To fit multiple models to multiple data sets, one can use the msenlm() function. The models to fit are specified with the set_models() function. The data should be stored in a list. It could be simulated data create using the create_simulated_datasets() function.
Models <- set_models (mean_class="test", err_class="binary")
Pars <- create_default_par_list (models=Models)
Data <- create_simulated_datasets (Pars, seed=12345)
Fits <- msenlm (models=Models, data=Data, xvar="x", yvar=names(Data)[-1])
mean_fun err_dist binomial_n
1 constant bernoulli NA
mean_fun err_dist binomial_n
2 uniform bernoulli NA
mean_fun err_dist binomial_n
1 constant bernoulli NA
mean_fun err_dist binomial_n
2 uniform bernoulli NA
print (Models)
mean_fun err_dist binomial_n
1 constant bernoulli NA
2 uniform bernoulli NA
summary (Fits)
convergence y x model mean_fun err_dist H c
1 0 constant_bernoulli x constant_bernoulli constant bernoulli 0.9140000 NA
2 0 constant_bernoulli x uniform_bernoulli uniform bernoulli 0.9140000 -4.092677
3 0 uniform_bernoulli x constant_bernoulli constant bernoulli 0.3580000 NA
4 0 uniform_bernoulli x uniform_bernoulli uniform bernoulli 0.9132653 30.345761
d npar nll AIC AICc BIC
1 NA 1 146.59213 295.1843 295.1923 299.3989
2 108.35812 3 146.59213 299.1843 299.2327 311.8281
3 NA 1 326.12939 654.2588 654.2668 658.4734
4 69.92516 3 57.80379 121.6076 121.6560 134.2514Alternatively, one can fit multiple models to real data. The msenlm.best function picks the best model for each response variable using the given criteria. The summary function summarises the fitted models.
Models <- set_models (mean_fun="gaussian", err_dist=c("zip", "zinb"))
Fits <- msenlm (models=Models, data=haul, xvar="depth",
yvar=c("Albatrossia.pectoralis", "Sebastolobus.altivelis"))
mean_fun err_dist binomial_n
1 gaussian zip NA
mean_fun err_dist binomial_n
2 gaussian zinb NA
mean_fun err_dist binomial_n
1 gaussian zip NA
mean_fun err_dist binomial_n
2 gaussian zinb NA
summary(Fits)
convergence y x model mean_fun err_dist H m
1 0 Albatrossia.pectoralis depth gaussian_zip gaussian zip 58.11781 1537.8453
2 0 Albatrossia.pectoralis depth gaussian_zinb gaussian zinb 29.74407 1145.7279
3 0 Sebastolobus.altivelis depth gaussian_zip gaussian zip 817.53123 917.8020
4 0 Sebastolobus.altivelis depth gaussian_zinb gaussian zinb 1448.37043 900.6961
s pi phi npar nll AIC AICc BIC
1 369.7516 4.303492e-01 NA 4 3022.4148 6052.830 6052.867 6072.739
2 213.8882 1.378740e-13 3.044749 5 927.2072 1864.414 1864.471 1889.301
3 225.5626 5.954216e-01 NA 4 36950.0517 73908.103 73908.141 73928.013
4 164.9688 1.548443e-01 1.575696 5 2910.8488 5831.698 5831.754 5856.584
summary(msenlm.best(Fits, best="AICc"))
convergence y x model mean_fun err_dist H m
1 0 Albatrossia.pectoralis depth gaussian_zinb gaussian zinb 29.74407 1145.7279
2 0 Sebastolobus.altivelis depth gaussian_zinb gaussian zinb 1448.37043 900.6961
s pi phi npar nll AIC AICc BIC
1 213.8882 1.378740e-13 3.044749 5 927.2072 1864.414 1864.471 1889.301
2 164.9688 1.548443e-01 1.575696 5 2910.8488 5831.698 5831.754 5856.584The maximum likelihood estimation needs good initial parameter estimates to work well. In this section we describe the methods used to obtain these estimates.
Given data \({x_i, y_i}\), we first approximate the mean function, using a smoothing spline (stats::smooth.spline) resulting in a set of points \({x_i,\mu_i}\). Since our data are non-negative, any negative values of the spline (\(\mu_i\) values) are set to zero.
From this we extract the following estimates.
The table below list the parameters of each error model.
| Model | Parameters | |||||
|---|---|---|---|---|---|---|
| bernoulli | ||||||
| binomial.count | ||||||
| binomial.prop | ||||||
| poisson | ||||||
| negbin | \(\phi\) | |||||
| zip | \(\pi\) | |||||
| zinb | \(\pi\) | \(\phi\) | ||||
| zipl | \(\gamma_0\) | \(\gamma_1\) | ||||
| zinbl | \(\gamma_0\) | \(\gamma_1\) | \(\phi\) | |||
| zipl.mu | \(\gamma_0\) | \(\gamma_1\) | ||||
| zinbl.mu | \(\gamma_0\) | \(\gamma_1\) | \(\phi\) | |||
| gaussian | \(\sigma\) | |||||
| tweedie | \(\phi\) | \(\rho\) | ||||
| zig | \(\pi\) | \(\phi\) | ||||
| zigl | \(\gamma_0\) | \(\gamma_1\) | \(\phi\) | |||
| zigl.mu | \(\gamma_0\) | \(\gamma_1\) | \(\phi\) | |||
| ziig | \(\pi\) | \(\phi\) | ||||
| ziigl | \(\gamma_0\) | \(\gamma_1\) | \(\phi\) | |||
| tab | \(\sigma\) | |||||
| zitab | \(\pi\) | \(\sigma\) |
The dispersion parameter, \(\phi\), is estimated as follows.
For a guassian error distribution, \(\sigma\) is a standard deviation and is estimated as the standard deviation of \(y_i - \mu_i\).
For a tail-adjusted beta, and a zero-inflated tail adjusted beta, \(\sigma\) is a shape parameter which is estimated using the method of moments.
The parameter \(\pi\) is usually estimated to be the proportion of zeros in the data (\(y\)) minus the expected number of zeros due to the spike. For the zitab model, \(\pi\) was estimated using the method of moments.
If the final estimate of \(\pi\) is less than 0.01 then it is set to 0.01.
The table below list the parameters of each mean function.
| Model | Parameters | ||||||
|---|---|---|---|---|---|---|---|
| constant | \(H\) | ||||||
| uniform | \(H\) | \(c\) | \(d\) | ||||
| beta | \(H\) | \(c\) | \(d\) | \(u\) | \(v\) | ||
| gaussian | \(H\) | \(m\) | \(s\) | ||||
| mixgaussian | \(H\) | \(m_1\) \(m_2\) | \(s_1\) \(s_2\) | \(\alpha\) | |||
| mixgaussian.equal | \(H\) | \(m_1\) \(m_2\) | \(s\) | \(\alpha\) | |||
| sech | \(H\) | \(m\) | \(s\) | \(r\) | \(p\) | ||
| sech.p1 | \(H\) | \(m\) | \(s\) | \(r\) | |||
| sech.r0p1 | \(H\) | \(m\) | \(s\) | ||||
| modskurt | \(H\) | \(m\) | \(s\) | \(q\) | \(p\) | \(b\) | |
| hofV | \(H\) | \(m\) | \(w_1\) \(w_2\) | \(k\) | |||
| hofII | \(H\) | \(m\) | \(\omega_0\) | ||||
| hofIV | \(H\) | \(m\) | \(\omega\) | \(k\) | |||
| hofIVb | \(H\) | \(m\) | \(\omega\) | ||||
| hofVb | \(H\) | \(m\) | \(\omega_1\) \(\omega_2\) |
The estimate of \(H\) is the average of \(y_i\). When the error distribution is bernoulli, tab, or zitab (i.e. \(y_i\) is 0-1), \(H\) is thresholded to a maximum value of 0.9.
The boundary parameters \(c\) and \(d\) are estimated to be the minimum and maximum values of \(x_i\) respectively, conditional on the \(y_i\) being positive.
The beta parameters \(u\) and \(v\) are estimated using the standard method of moments.
The estimates of the means \(m_1\) and \(m_2\) are taken to be the 35th and 65th percentiles of the \(x\) values (when \(y\) is positive). This is largely independent of the \(y\) data but is a simple method for keeping the two modes from converging too early in the fitting algorithm.
The standard deviations \(s\), or \(s_1\) and \(s_2\), are taken to be the estimate of \(s\) from the smoothing spline.
The mixing parameter, \(\alpha\), is set to 0.5.
The mode location parameter, \(m\), is estimated from the smoothing spline.
The skewness parameter, \(r\), is set to 0.
The peakedness parameter, \(p\), is set to 1.
For the spread parameter, \(s\): let \(s_x\) be the standard deviation of the \(x\) values weighted by the \(y\) values, and let \(s_0\) be the smallest non-zero distance between the \(x\) values. The initial estimate of \(s\) is taken to me the maximum of \(s_x\) and \(s_0\).
The parameters \(m\), and \(s\) are estimated as for the Sech function.
The skewness parameter, \(q\) is set to 0.5 (symmetric).
The peakedness parameter, \(p\), is set to 1.
We initially set the flatness parameter to \(2\). The parameter estimation works better starting from a flatter curve than needed, rather than a more peaked curve.