Let us consider again the data from Roux et al. (2001). This data set
is built into the package and can be accessed from the
Psurveys object. P and S probability terms arise in the
activity-level interpretation of trace evidence; Evett and Buckleton
(1990) give the foundational glass-evidence formulation, while Coulson
et al. (2001) describe the clothing-survey setting in which zeta
modelling became useful. That is, we can type:
Printing the survey object gives a compact summary of the observed frequencies:
Number of Groups
n rn
--- ----
0 754
1 9
2 8
3 4
4 1
Roux C, Kirk R, Benson S, Van Haren T, Petterd C (2001).
"Glass particles in footwear of members of the public in
south-eastern Australia-a survey." _Forensic Science
International_, *116*(2), 149-156.
doi:10.1016/S0379-0738(00)00355-8
<https://doi.org/10.1016/S0379-0738%2800%2900355-8>.Fit a zeta distribution with fit() and the
zetaModel() model descriptor:
New analyses should use this generic fit() interface;
older distribution-specific fitting functions are retained only for
backward compatibility.
The fitted object contains the estimated shape parameter \(\alpha\), its standard error, and the fitted probability terms. Printing the object gives a concise summary of these quantities.
The print method displays the first 10 fitted
probabilities from the model by default.
The estimated shape parameter is 4.9544
The standard error of shape parameter is 0.2366
The first 10 fitted values are:
P0 P1 P2 P3 P4
9.631547e-01 3.106447e-02 4.167082e-03 1.001917e-03 3.316637e-04
P5 P6 P7 P8 P9
1.344002e-04 6.262053e-05 3.231467e-05 1.802885e-05 1.069709e-05 For probability terms beyond those printed by default,
probfun() creates a function that evaluates the fitted
model at any requested term:
Calling the returned function with a value of \(k\) gives the corresponding fitted \(P_k\):
P5
0.0001344002 We can also fit a zero-inflated zeta model through the same generic
interface by supplying zizModel(). As before, we can choose
a variable name to store the results in.
The estimated mixing parameter, pi, is 0.8465
The estimated shape parameter is 2.8846
The first 10 fitted values are:
P0 P1 P2 P3 P4
0.9716490911 0.0169404164 0.0052597614 0.0022938450 0.0012050764
P5 P6 P7 P8 P9
0.0007122067 0.0004565511 0.0003106019 0.0002211302 0.0001631754 The same generic fitting interface is used for both zeta and ZIZ models.
In the example above we fit a zero-inflated model to Roux et al.’s data and print the resulting fit. We get estimates of the parameters and a default set of fitted values. The output is interesting because the value of \(\hat{\pi}\) shows that the zero part of the zero-inflated model is picking up about 85% of the zeros. It is interesting to compare the estimates from the raw frequencies, the zeta model, and the ZIZ model. The estimates are shown in Table \(\ref{tab:ex1}\).
We can see from Table \(\ref{tab:ex1}\) that we now have a non-zero estimate for \(P_5\), but this comes at the cost of smaller probabilities for the preceding terms \(P_0\)–\(P_4\), which is not necessarily a negative. The survey data is dominated by zeros. However, we think it likely that the raw sample estimates for \(P_0\)–\(P_4\) are overestimates. The model reduces the estimated value, which is in line with our thinking. Interestingly, the effect of including the zero-inflation factor is to increase nearly all of the probabilities, with the exception of \(P_1\). A natural question to ask is “Which model is correct?” The answer, unhelpfully, is “Neither”, because these are simply models. They can still help us without us having to believe that they are true.
The fitPS package provides a confint method
for the fitted value. The method returns both a Wald confidence interval
and a profile likelihood interval. The two intervals are returned as
elements of a list named wald and
prof, respectively.
2.5% 97.5%
4.490761 5.418099 2.5% 97.5%
4.520495 5.451277 For a two-parameter zero-inflated zeta fit,
plotUncertainty() shows joint parameter uncertainty. For an
MLE fit it draws profile-likelihood confidence regions. Use
confint() when numerical interval or region information is
required.
The same plotting call can be used after an ordinary bootstrap fit. In that case the plot shows a smoothed region based on the distribution of bootstrap parameter estimates.
> boot.fit = fit(
> roux,
> model = zizModel(),
> method = "bootstrap",
> B = 2000,
> seed = 123,
> silent = TRUE
> )
> plotUncertainty(boot.fit, level = c(0.80, 0.95))For new analyses, obtain the bootstrap distribution with
fit(..., method = "bootstrap") and use
plotUncertainty() to display the corresponding parameter
uncertainty.
We can use the methodology that has been demonstrated so far to compare surveys. One reason for comparing surveys is to explore the hypothesis that there is no difference in the underlying true value of \(\alpha\). If there is insufficient evidence to reject this hypothesis, then one may feel justified in combining data from two surveys. In the first instance we will take an ad hoc approach, and then treat this problem more formally. In our ad hoc approach we will compare confidence intervals for two surveys. If these confidence intervals overlap, then we might conclude that there is insufficient evidence in the data to suggest that the estimates of \(\alpha\) are different. We will illustrate this with the surveys conducted by Lau et al. (1997) and Jackson et al. (2013). Lau et al. (1997) surveyed the clothing of 213 Canadian high school students and observed two sets of clothing with one fragment on each. Similarly, Jackson et al. (2013) surveyed 232 “randomly” selected members of the population of New South Wales in Australia. We place “randomly” in quotes because this was not a true random sample, but rather a convenience sample. That being said, it is unlikely that using a truly random mechanism would have significantly changed the results.
Visual inspection of these surveys would suggest that they are fairly
similar. We can fit a zeta distribution to each survey, and then compute
a confidence interval for each survey. Again, these data sets are
included in the fitPS package.
> lau = Psurveys$lau
> jackson = Psurveys$jackson
> fit.lau = fit(lau, model = zetaModel())
> fit.jackson = fit(jackson, model = zetaModel())
> confint(fit.lau)$wald 2.5% 97.5%
4.948422 8.823417 2.5% 97.5%
4.443639 6.623473 We can see from the output that there is substantial overlap between these two Wald confidence intervals. The results using profile likelihood intervals lead to the same conclusion but are not shown. We can test this more formally. Specifically, we wish to test the null hypothesis that
\[ H_0: \alpha_1 = \alpha_2 \quad \mbox{or equivalently} \quad H_0: \alpha_1 - \alpha_2 = 0, \]
where \(\alpha_1\) is the true value of \(\alpha\) for the Lau et al. data, and \(\alpha_2\) is the true value of \(\alpha\) for the Jackson et al. data. We choose a two-tailed alternative, meaning we are not concerned about the sign of any difference, but simply the magnitude of the difference. That is,
\[ H_1: \alpha_1 \neq \alpha_2 \quad \mbox{or equivalently} \quad H_1: \alpha_1 - \alpha_2 \neq 0. \]
We test this hypothesis by constructing a test statistic and then computing a P-value under the assumption that the null hypothesis is true. We are interested in the difference between the two population values of \(\alpha\). We estimate this by computing the difference in the sample estimates. That is, our estimate of \(\alpha_1 - \alpha_2\) is given by \(\hat{\alpha}_1 - \hat{\alpha}_2\), where \(\hat{\alpha}_1\) and \(\hat{\alpha}_2\) are the maximum likelihood estimates based on the survey data. We scale this difference by the estimated standard deviation in the difference, that is, by the standard error of the difference, \(\mathrm{se}(\hat{\alpha}_1 - \hat{\alpha}_2)\). We estimate this–to keep the statistical theory to a minimum–as the square root of the sum of the two estimated variances, i.e.
\[ \mathrm{se}(\hat{\alpha}_1 - \hat{\alpha}_2) = \sqrt{\hat{V}(\hat{\alpha}_1) + \hat{V}(\hat{\alpha}_2)}. \]
Our test statistic is then
\[ Z_0 = \frac{\hat{\alpha}_1 - \hat{\alpha}_2}{\mathrm{se}(\hat{\alpha}_1 - \hat{\alpha}_2)}. \]
It can be shown that this test statistic follows an approximate normal distribution under the null hypothesis, which means our P-value can be computed by evaluating
\[ P = \Pr(Z > |Z_0|) = 2(1 - \Pr(Z < |Z_0|)). \]
All this theory has been integrated into a function called
compareSurveys.
Two-sided Wald test
data: lau and jackson
z = 1.1923, p-value = 0.2331
alternative hypothesis: true difference in shape parameters is not equal to 0
sample estimates:
Shape of lau Shape of jackson
6.885919 5.533556 The P-value is 0.23 (2 d.p.). This is significantly larger than
either 0.05 or 0.01. Based on this, we would conclude that there is
insufficient evidence to reject the null hypothesis of a common value of
\(\alpha\), and therefore it may be
sensible to combine data from these two surveys. We could have also used
the theory of likelihood ratio tests to test this hypothesis, but that
is beyond the scope of this article. We note, however, that the
fitPS package contains a function called
compareSurveysLRT which can compare two or more surveys
simultaneously using a likelihood ratio test.
Curran, J. M., Buzzini, P., and Trejos, T. (2024). Estimating probability terms for the background presence of glass when considering activity in forensic casework. Forensic Science International, 364, 112221. https://doi.org/10.1016/j.forsciint.2024.112221
Coulson, S. A., Buckleton, J. S., Gummer, A. B., and Triggs, C. M. (2001). Glass on clothing and shoes of members of the general population and people suspected of breaking crimes. Science & Justice, 41(1), 39-48. https://doi.org/10.1016/S1355-0306(01)71847-3
Evett, I. W., and Buckleton, J. S. (1990). The interpretation of glass evidence. A practical approach. Journal of the Forensic Science Society, 30(4), 215-223.
Jackson, F., Maynard, P., Cavanagh-Steer, K., Dusting, T., and Roux, C. (2013). A survey of glass found on the headwear and head hair of a random population vs. people working with glass. Forensic Science International, 226(1), 125-131.
Lau, L., Beveridge, A. D., Callowhill, B. C., Conners, N., Foster, K., Groves, R. J., Ohashi, K. N., Sumner, A. M., and Wong, H. (1997). The frequency of occurrence of paint and glass on the clothing of high school students. Canadian Society of Forensic Science Journal, 30(4), 233-240.
Roux, C., Kirk, R., Benson, S., Van Haren, T., and Petterd, C. I. (2001). Glass particles in footwear of members of the public in south-eastern Australia: a survey. Forensic Science International, 116(2), 149-156.