03c_stats_fitting.Rmd

# Dose-effect curve fitting {#stats-fitting}

## Fitting method

The technique suggested for determining the best fit coefficients is that of maximum likelihood [@Papworth1975; @Merkle1983]. Using this method, the best fit value for each coefficient is achieved by assuming a Poisson distribution and maximising the likelihood of the observations by the method of iteratively reweighted least squares. For overdispersed (non-Poisson) distributions, as obtained after high-LET radiation, the weights must take into account the overdispersion. If the data show a statistically significant trend of $\hat{\sigma}^{2}/\bar{y}$ with dose, then that trend should be used. Otherwise, the Poisson weight on each data point should be divided by the average value of $\hat{\sigma}^{2}/\bar{y}$.

Table \@ref(tab:th-dic-count-data)  gives example data used to construct a dose-effect curve for low-LET $\gamma$-radiation [@Barquinero1995], while Table \@ref(tab:th-dic-coeffs) shows the estimated coefficients when the data from Table \@ref(tab:th-dic-count-data) are fitted. The $p$-values of the $t$-test indicate that each parameter is statistically significant.

```{r th-dic-count-data, echo=FALSE}
dic_count_data <- here::here("data/count-data-IAEA.csv") |>
  utils::read.csv() |>
  biodosetools::calculate_aberr_table(type = "count")

dic_count_data <- dic_count_data |>
  as.matrix() |>
  formatC(format = "f", digits = 3) |>
  as.data.frame() |>
  dplyr::mutate(
    dplyr::across(
      .cols = c("N", "X", dplyr::starts_with("C")),
      .fns = as.integer
    )
  ) |>
  biodosetools:::fix_count_data_names(type = "count", output = "kable")

colnames(dic_count_data)[1] <- "$D \\text{ (Gy)}$"

dic_count_data |>
  knitr::kable(
    caption = "Dicentric distribution within cells, dispersion index, and $u$-value of the dose-effect curve constructed from cytogenetic results obtained from blood samples irradiated with $\\gamma$-rays [@Barquinero1995].",
    align = "c"
  )
```

For each dose analysed, the total number of cells scored ($N$), the total number of dicentrics observed ($X$), the cell distribution of dicentrics ($C_{0}, \dots, C_{5}$), and the dispersion index ($\hat{\sigma}^{2}/\bar{y}$) and $u$-test statistic ($u$) are presented (a sample is considered overdispersed if $u > 1.96$).

```{r th-dic-coeffs, echo=FALSE}
dic_coeffs <- here::here("data/dicentrics-coeffs.csv") |>
  data.table::fread() |>
  dplyr::mutate(
    unit = c("", "Gy⁻¹", "Gy⁻²"),
    est_formatted = paste0(
      "(",
      Estimate,
      " ± ",
      round(`Std. Error`, 5),
      ") ", unit
    ),
    `p-value` = scales::pvalue(`p-value`)
  ) |>
  dplyr::select(
    Coefficient,
    Estimate = est_formatted,
    "$t$-statistic" = `t-statistic`,
    "$p$-value" = `p-value`
  )

dic_coeffs |>
  knitr::kable(
    caption = "Fitted values of the coefficients of the LQ dose-effect calibration curve \\@ref(eq:dose-response-lq) constructed from the dicentric distribution in Table \\@ref(tab:th-dic-count-data).",
    align = "c",
    escape = TRUE
  )
```

## Goodness of fit

The goodness of fit of the fitted curve and significance of estimated coefficients should then be tested, for instance using an appropriate form of the $F$-test, $z$-test or $t$-test. {biodosetools} implements the $t$-test.

Let $\hat{\theta}$ be an estimator of the parameter $\theta \in \{ \alpha, \beta, C\}$ in the fit model. Then the $t$-statistic for this parameter is defined as
\begin{align}
  t_{\hat{\theta}} = \frac{\hat{\theta}}{\widehat{se}(\hat{\theta})},
  (\#eq:t-test)
\end{align}
where $\widehat{se}(\hat{\theta})$ is the standard error of $\hat\theta$.

## Dose-effect calibration

Adequate curve fitting requires a sufficient number of evaluated dose points to minimise the error. A minimum of 7 doses should be evaluated, 5 of them at doses equal or lower than 1 Gy, including the 0 Gy dose. Usually, for low-LET radiation (X-rays and $\gamma$-rays) the dose range evaluated is 0--5 Gy, indeed, beyond this dose there is evidence of saturation of the aberration yield which will lead to a distortion of the $\beta$ coefficient in equation \@ref(eq:dose-response-lq) [@ISO19238:2014; @Lloyd1983]. For high-LET radiation, a range of 0--2 Gy is suggested [@IAEA2011].

At higher doses, scoring should aim to detect 100 dicentrics at each dose. However, at lower doses this is difficult to achieve and instead several thousand cells per point should be scored; a number between \num{3000} and \num{5000} is suggested. In all cases, the actual number of cells scored should be dependent on the number of dose points in the low dose region, with the focus on minimising the error on the fitted curve [@Higueras2020].

Opinions vary on how to treat the background level of aberrations in fitting dose-effect data. In general there are three approaches: (a) a dose point at 0 Gy is included in the curve fitting procedure, (b) the zero dose point is ignored, or else (c) the zero dose point is represented in every fitting procedure by a standard background value. Nowadays almost all laboratories include a 0 Gy dose point, and it should always be included for dicentrics, translocations and micronuclei.

If the measured yield at zero dose is used as one of the data points for the curve fitting (as used in the curve fitting presented in Table \@ref(tab:th-dic-coeffs)), the background becomes a variable parameter. However, since the yield in unirradiated cells is usually low, often none are observed so the measured yield at zero dose is zero. As discussed, at low doses, the statistical resolution of the data points is generally low. Thus, including the zero dose point in the curve fitting procedure can sometimes lead to negative estimates of the background value ($C$) and negative linear coefficients ($\alpha$), which obviously have no biological basis. Some investigators resolve this problem by ignoring zero dose data points and constraining the curve to pass through the origin. Another way to solve this is to fit the calibration curve model with constrained maximum likelihood which forces the intercept parameter to be non-negative [@Oliveira2016].

There are, however, sufficient data published from surveys of subjects exposed only to background radiation to show that there is a small positive background level of aberrations. An alternative method adopted by some experts is therefore to use a small positive background value as a data point and to ascribe a large percentage of uncertainty to it. Ideally a laboratory should generate its own background data, although this requires the analysis of many thousands of cells. A consensus has emerged that the background level of dicentrics is ~ 0.5--1.0 per 1000 cells [@Lloyd1980] whilst for translocations [@Sigurdson2008] and for micronuclei [@Fenech1993] the control values are higher.

## Protracted and fractionated exposure

In case of a protracted or fractionated exposure the resulting chromosomal aberration yield may be lower than when receiving the same dose acutely. No dose rate or fractionation effect would be expected for high-LET radiation, where the dose-effect relationship is close to linear. However, with low-LET radiation, dose protraction has the impact of reducing the dose squared coefficient, $\beta$, in the yield \@ref(eq:dose-response-lq). It is assumed that this term refers to those aberrations generated by the interaction of two tracks or more. Exchange type aberrations, such as dicentrics have a fast kinetics formation [@Darroudi1998; @Durante1996; @Durante1998; @Pujol2020], and hence the effect on the interaction of two tracks can be modified by repair mechanisms that have time to function over the course of a long exposure or in the intervals between intermittent acute exposures. The decrease in the frequency of aberrations appears to follow a single exponential function with a mean time of about 2 hours.

In the early times of radiobiology and evaluating the effect of IR on \emph{Tradescantia} @Lea1942 suggested a time-dependent factor called the $G$ function \@ref(eq:protraction-function) to allow change of the dose squared coefficient and so account for the impacts of dose protraction.
\begin{align}
  \lambda = C + \alpha D + \beta G(x) D^{2},
  (#eq:dose-response-lq-protracted)
\end{align}
where $t$ is the time over which the irradiation occurred, and $t_{0}$ is the mean lifetime of the breaks, which has been shown to be of the order of ~ 2 hours [@Bauchinger1979; @Lloyd1984].

Currently in biological dosimetry this approach is still accepted and that is why the LQ dose-effect equation \@ref(eq:dose-response-lq) may be modified as follows:
\begin{align}
  G(x) = \frac{2}{x^{2}} \left( x - 1 + e^{-x} \right) ,
  \quad \text{and} \quad
  x = \frac{t}{t_{0}}.
  (#eq:protraction-function)
\end{align}

Given a value of $t$, we can have three protraction cases:

- Acute exposure when $t$ is approximately 0. In this case $\lim\limits_{x\to 0} G(x) = 1$. Therefore, the dose-effect equation becomes \@ref(eq:dose-response-lq).
- Highly protracted exposure when $t$ is high. In this case $G(x)$ reduces virtually to zero. Therefore, even if a high dose (> 1.0 Gy) is involved, the dose-effect equation becomes linear \@ref(eq:dose-response-l).
- Protracted exposure for any case in between. In this case the effects of $G(x)$ need to be considered.

All dose estimation methods implemented in {biodosetools} consider protraction, with the extreme of $G(x) = 1$ or $G(x) \approx 0$ being just particular cases. For this reason, all methods described in Chapter \@ref(stats-estimation) are derived using equation \@ref(eq:dose-response-lq-protracted) instead of \@ref(eq:dose-response-lq).

## Genomic equivalence for translocations

Curve fitting for translocations follow the same procedure described for dicentrics. However, Fluorescence \emph{In Situ} Hybridisation (FISH) analysis only evaluates the translocation frequency in a specific set of chromosomes painted. The conversion of this frequency to full genome equivalence is a recommended procedure to use when different data must be combined when results using different combinations of whole chromosome paintings are compared, e.g. in the frame of interlaboratory comparisons.

When the DNA probes used to paint different chromosome pairs are labelled with different fluorochromes, the genomic conversion factor $F_{p}/F_{G}$ is usually calculated by using the formula for the painted fractions of the genome [@Lucas2000]:
\begin{align}
  \frac{F_{p}}{F_{G}} = \frac{2}{0.974} \left[ \sum_{i} f_{i} (1 - f_{i}) - \sum_{i < j} f_{i} f_{j} \right],
  (#eq:genome-factor)
\end{align}
where $F_{G}$ is the full genome aberration frequency, $F_{p}$ is the translocation frequency detected by FISH, and $f_{i}$ is the fraction of genome corresponding to each chromosome $i$ used in the hybridisation, taking into account the gender of the subjects [@Morton1991]. The total number of interchromosomal exchanges is 0.974, using the same assumption of DNA proportionality (see calculations in @Lucas1992).

In the case when the DNA probes used to paint different chromosomes are all labelled with the same fluorochrome, the genomic conversion factor \@ref(eq:genome-factor) reduces to:
\begin{align}
  \frac{F_{p}}{F_{G}} = \frac{2}{0.974} f_{p} (1 - f_{p}) = 2.05 f_{p} (1 - f_{p}),
  (#eq:genome-factor-simple)
\end{align}
where $F_{G}$ is the full genome aberration frequency, $F_{p}$ is the translocation frequency detected by FISH, and $f_{p}$ is the fraction of genome covered by the combination of whole chromosome probes used in the hybridisation, also taking into account the gender of the subjects.

Translocations have higher background levels than dicentrics, which is partly due to the former being a persistent type of aberration. When attempting retrospective biological dosimetry, it is critical to account for the translocation background, especially at low doses. Because a pre-exposure control blood sample from the unintentionally irradiated person or a population study group is not available, an estimated value based on generic survey data must be utilised. Ideally, a laboratory would create its own control database, but this is a huge undertaking because it would have to cover a lot of confounding factors and, more importantly, a wide range of age groups. The greatest international database, split down by age, sex, ethnicity, geographic region, and smoking behaviours, is currently available thanks to a comprehensive meta-analysis published by @Sigurdson2008.

It is important to account for the background and subtract from the total number of translocations observed in an individual's lymphocyte the expected translocation rate given a set of confounding factors \@ref(eq:sigurdson-corrected-FG), the most important of which is age:
\begin{align}
  F'_{G} = \frac{X - X_{c}}{N (F_{p} / F_{G})},
  (#eq:sigurdson-corrected-FG)
\end{align}
where $F'_{G}$ is the corrected full genome aberration frequency, $X$ is the total number of translocations observed, $X_{c}$ is the total number of expected aberrations, $N$ is the total number of cells, and $F_{p} / F_{G}$ is the genomic conversion factor given by \@ref(eq:genome-factor).