--- title: "Getting started with fitPS" output: rmarkdown::pdf_document vignette: > %\VignetteIndexEntry{Getting started with fitPS} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, echo=FALSE} options(width = 70) knitr::opts_chunk$set( message = FALSE, comment = "", highlight = TRUE, prompt = TRUE, tidy = TRUE, tidy.opts = list(arrow = FALSE), warning = FALSE ) suppressWarnings(suppressPackageStartupMessages(library(fitPS))) ``` ## Fitting a zeta distribution 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: ```{r load-roux} data("Psurveys") roux = Psurveys$roux ``` Printing the survey object gives a compact summary of the observed frequencies: ```{r print-roux} roux ``` Fit a zeta distribution with `fit()` and the `zetaModel()` model descriptor: ```{r fit-roux} fit = fit(roux, model = zetaModel()) ``` 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. ### Using the fitted distribution to estimate P terms The `print` method displays the first 10 fitted probabilities from the model by default. ```{r print-fit} fit ``` For probability terms beyond those printed by default, `probfun()` creates a function that evaluates the fitted model at any requested term: ```{r probability-function} P = probfun(fit) ``` Calling the returned function with a value of $k$ gives the corresponding fitted $P_k$: ```{r probability-example} P(5) ``` ## Fitting a zero-inflated zeta distribution 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. ```{r fit-zero-inflated} fit.zi = fit(roux, model = zizModel()) fit.zi ``` 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 `r round(100 * fit.zi$pi)`% 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}. ```{r ex1, echo=FALSE, results='asis'} raw = c(roux$data$rn / sum(roux$data$rn), 0) tableData = cbind( 0:5, raw, fitted(fit, 6), fitted(fit.zi, 6) ) colnames(tableData) = c( "$k$", "$P_k^{raw}$", "$P_k^{zeta}$", "$P_k^{ZIZ}$" ) probabilityTable = xtable::xtable( tableData, align = "ccccc", digits = c(0, 0, 4, 4, 4), caption = "Estimated probability that k groups of glass would be found in shoes of a random member of the population based on the data of Roux et al. (2001), the raw frequencies, and those produced from the zeta and ZIZ models respectively.", label = "tab:ex1" ) print( probabilityTable, type = "latex", include.rownames = FALSE, sanitize.text.function = identity, comment = FALSE ) ``` 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. ## Confidence intervals for the parameter estimates 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. ```{r confidence-intervals} ci = confint(fit) ci$wald ci$prof ``` ### Profile-likelihood and bootstrap confidence regions for the zero-inflated zeta 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. ```{r confidence-region, eval=FALSE} plotUncertainty(fit.zi, level = c(0.80, 0.95)) cr = confint(fit.zi, level = c(0.80, 0.95)) ``` 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. ```{r boot-confidence-region, eval=FALSE} 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. ## Comparing two surveys 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. ```{r compare-data-table, echo=FALSE, results='asis'} lau = Psurveys$lau jackson = Psurveys$jackson surveyTableData = cbind( 0:1, lau$data$rn, jackson$data$rn ) colnames(surveyTableData) = c( "$k$", "Lau et al. (1997)", "Jackson et al. (2013)" ) surveyTable = xtable::xtable( surveyTableData, align = "ccrr", digits = 0, caption = "Survey results from Lau et al. (1997) and Jackson et al. (2013).", label = "tab:survey-comparison" ) print( surveyTable, type = "latex", include.rownames = FALSE, sanitize.text.function = identity, comment = FALSE ) ``` 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. ```{r compare-fits} lau = Psurveys$lau jackson = Psurveys$jackson fit.lau = fit(lau, model = zetaModel()) fit.jackson = fit(jackson, model = zetaModel()) confint(fit.lau)$wald confint(fit.jackson)$wald ``` 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`. ```{r compare-surveys} compareSurveys(lau, jackson) ``` The P-value is `r round(compareSurveys(lau, jackson)$p.value, 2)` (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. ## References 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.