Type: Package
Title: Quantal Response Analysis for Dose-Mortality Data
Version: 0.2.8.1
Maintainer: John Maindonald <john@statsresearch.co.nz>
Description: Functions are provided that implement the use of the Fieller's formula methodology, for calculating a confidence interval for a ratio of (commonly, correlated) means. See Fieller (1954) <doi:10.1111/j.2517-6161.1954.tb00159.x>. Here, the application of primary interest is to studies of insect mortality response to increasing doses of a fumigant, or, e.g., to time in coolstorage. The formula is used to calculate a confidence interval for the dose or time required to achieve a specified mortality proportion, commonly 0.5 or 0.99. Vignettes demonstrate link functions that may be considered, checks on fitted models, and alternative choices of error family. Note in particular the betabinomial error family. See also Maindonald, Waddell, and Petry (2001) <doi:10.1016/S0925-5214(01)00082-5>.
Encoding: UTF-8
License: GPL-3
Depends: R (≥ 4.1.0), lattice, latticeExtra, knitr, rmarkdown
Imports: lme4, splines, ggplot2
Suggests: fitODBOD, VGAM, glmmTMB (≥ 1.1.2), gamlss, prettydoc, DHARMa, kableExtra (≥ 1.2), plotrix, dfoptim, optimx, bookdown
URL: https://github.com/jhmaindonald/qra
BugReports: https://github.com/jhmaindonald/qra/issues
VignetteBuilder: knitr, rmarkdown, bookdown, prettydoc
LazyData: TRUE
RoxygenNote: 7.2.3
NeedsCompilation: no
Packaged: 2025-05-22 01:38:44 UTC; johnm1
Author: John Maindonald ORCID iD [aut, cre]
Repository: CRAN
Date/Publication: 2025-05-22 04:20:02 UTC

qra: A package for calculations that relate to dose-mortality, or time-mortality, or other such models.

Description

The qra package provides the functions: checkDisp, getRho, extractLT, getScaleCoef, scaleLocAdjust, fieller, gpsWithin, varRatio, foldp, graphSum, fpower

Details

Vignettes provide examples of the use of the functions.

qra functions

fieller: Calculates lethal dose estimates, using Fieller's formula to calculate 95 of the functions noted below are useful ancillaries to fieller and fieller2, notably foldp, fpower, extractLT, and getScaleCoef.

fieller2: Use when a folded power link has been used. See fpower.

extractLT: Obtains complete set of LT or LD estimate, when it is convenient to get results from several models at the same time.

foldp: Calculates the ratio of p+eps to 1-code+eps

getRho Extracts estimates of the intra-class correlation from a glmmTMB model object with betabinomial error. See the vignette [timeMortality] for details of the parametization used for the betabinomial error.

getScaleCoef: Extracts the scale coefficients from a vector that has been scaled using scale, as needed so that the scaling can be undone.

gpsWithin Renumbers group identifiers so that they run from 1 to number of groups within for each level of the specified factor.

scaleLocAdjust: Returns, for glmmTMB models with a betabinomial error, dispersion factors (i.e., multipliers for the binomial variance) as functions of predicted values.

varRatio: Returns a first order approximation to the variance of the $y$-ordinate to slope ratio. This is used in the type="Delta" approximation, for calculation of LT and LD confidence intervals. Primarily, this is provided for purposes of comparison, to make it easy to show how poor the approximation can be, and to warn against its general dewvuse!

Author(s)

Maintainer: John Maindonald john@statsresearch.co.nz (ORCID)

See Also

Useful links:

Hawaiian Contemporary Cold Treatment Dataset

Description

The counts of live/dead were derived by injecting a known number of individuals of the target life stage into citrus fruits, subjecting them to treatment and then counting the number of individuals emerging.

Usage

data("HawCon")

Format

A data frame with 106 observations on the following 10 variables.

Species

Species of fruitfly

CN

Common name, in abbreviated form. MedFly is ‘Mediterranean Fruit Fly’. MelonFly is ‘Melon Fly’

LifestageTrt

Lifestage treated

RepNumber

Replicate number

PropDead

Fraction dead

TrtTime

Treatment time (days)

Dead

a numeric vector

Live

a numeric vector

Total

a numeric vector

Details

The help page for HawCon in the ColdData has further details.

Source

Dr Peter Follett

References

A paper is in the course of preparation.

Examples

data(HawCon)
str(HawCon)

Reproduce data for the linear model scale-location diagnostic plot

Description

The values returned are those used for plot(x.lm, which=3), where x.lm is a linear model or a generalized linear model. Plot the object returned to assess how successful the weights, determined using the function scaleLocAdjust, have been in adjusting for heterogenous variances.

Usage

checkDisp(x, span = 0.75)

Arguments

x

Model fitted using lm() or glm()

span

span parameter for use in smoothing the square root of standardized deviance residuals.

Value

A data frame, with:

linpred

Predicted values, on the scale of the linear predictor

absrSmooth

Smoothed values of the square roots of absolute values of standardised deviance residuals.

Examples

royal <- subset(qra::codling1988, Cultivar=="ROYAL")
royal.glm <- glm(cbind(dead,total-dead)~ct, data=royal,
                 family=quasibinomial(link='cloglog'))
royalFix <- qra::scaleLocAdjust(royal.glm, lambda=2)
## Check range of indicated prior weights
range(royalFix[[2]])
## Range of updated dispersion estimates
range(summary(royalFix[[1]])[['dispersion']]/royalFix[[2]])
xy <- qra::checkDisp(royalFix[[1]])
plot(xy)

Dose-mortality data, for fumigation of codling moth with methyl bromide

Description

Data are from trials that studied the mortality response of codling moth to fumigation with methyl bromide, for the year 1988 only

Usage

data(codling1988)
data(codling1989)

Format

A data frame with 77 observations (codling1988), and with 40 observations (codling1989), on the following 8 variables.

dose

Injected dose of methyl bromide, in gm per cubic meter

ct

Concentration-time sum

total

Number of insects in chamber

dead

Number of insects dying

PropDead

Proportion dying

Cultivar

a factor with 1988 levels BRAEBURN FUJI GRANNY Red Delicious and ROYAL; and with 1989 levels Gala, Red Delicious and Splendour

rep

replicate number, within Cultivar

cultRep

Cultivar/replicate combination

Details

The research that generated these data was in part funded by New Zealand pipfruit growers. The published analysis was funded by New Zealand pipfruit growers. See also DAAG::sorption.

Source

Maindonald, J.H.; Waddell, B.C.; Petry, R.J. 2001. Apple cultivar effects on codling moth (Lepidoptera: Tortricidae) egg mortality following fumigation with methyl bromide. Postharvest Biology and Technology 22: 99-110.

Obtain complete set of LT or LD estimates

Description

When supplied with a model object that has fitted dose-response lines for each of several levels of a factor, extractLT calls the function fieller to calculate lethal time

Usage

extractLT(
  obj,
  a = 1:3,
  b = 4:6,
  link = NULL,
  logscale = FALSE,
  p = 0.99,
  eps = 0,
  offset = 0,
  df.t = NULL
)

extractLTpwr(
  obj,
  a = 1:3,
  b = 1:3,
  link = "fpower",
  logscale = FALSE,
  p = 0.99,
  lambda = 0,
  eps = 0.015,
  offset = 0,
  df.t = NULL
)

Arguments

obj

merMod object, created using lmer() or glmerMod object, created using glmer().

a

Subscripts for intercepts.

b

Subscripts for corresponding slopes.

link

Link function, for use with objects where no link was specified in the function call, but it is required to back-transform a transformation that was performed prior to the function call. Otherwise leave as link=NULL, and the link function will be extracted as family(obj)[['link']]. For a folded power function, with extractLTpwr(), the only available link is fpower, and the exponent lambda must be specified.

logscale

Logical. Specify TRUE, if LT values are to be back-transformed from a logarithmic scale.

p

Target response proportion.

eps

Replace prob by prob+eps before transformation.

offset

Use to undo scaling of time or dose variable. This is passed to the fieller function that extractLT calls.

df.t

Degrees of freedom for a t-distribution approximation for 't' or 'z' statistics. If NULL, a conservative (low) value will be used. For linear (but not generalized linear) models and mixed models, approximations are implemented in the afex package. See vignette('introduction-mixed-models', package="afex"), page 19.

lambda

(extractLTpwr only) Power for power function.

Details

Fixed coefficients from obj must be for intercepts and for slopes. Starting the model formula with 0+ will commonly do what is required. The coefficients fixef(obj)[a] are assumed to specify line intercepts, while fixef(obj)[b] specify the corresponding slopes. These replace the arguments nEsts (subscripts for intercepts were 1:nEsts) and slopeAdd (subscripts for slopes were (nEsts+1):(nEsts+slopeAdd)).

Value

Matrix holding LD or LD estimates.

Examples

pcheck <- suppressWarnings(requireNamespace("glmmTMB", quietly = TRUE))
if(pcheck) pcheck & packageVersion("glmmTMB") >= "1.1.2"
if(pcheck){
form <- cbind(Dead,Live)~0+trtGp/TrtTime+(1|trtGpRep)
HawMed <- droplevels(subset(HawCon, CN=="MedFly"&LifestageTrt!="Egg"))
HawMed <- within(HawMed,
                 {trtGp <- factor(paste0(CN,LifestageTrt, sep=":"))
                 trtGpRep <- paste0(CN,LifestageTrt,":",RepNumber)
                 scTime <- scale(TrtTime) })
HawMedbb.cll <- glmmTMB::glmmTMB(form, dispformula=~trtGp+splines::ns(scTime,2),
                                 family=glmmTMB::betabinomial(link="cloglog"),
                                 data=HawMed)
round(qra::extractLT(p=0.99, obj=HawMedbb.cll, link="cloglog",
               a=1:3, b=4:6, eps=0, df.t=NULL)[,-2], 2)} else
message("Example requires `glmmTMB` version >= 1.1.2: not available")

Confidence Limits for Lethal Dose Estimate From Dose-response Line

Description

This uses Fieller's formula to calculate a confidence interval for a specified mortality proportion, commonly 0.50, or 0.90, or 0.99. Here "dose" is a generic term for any measure of intensity of a treatment that is designed to induce insect death.

Usage

fieller(
  phat,
  b,
  vv,
  df.t = Inf,
  offset = 0,
  logscale = FALSE,
  link = "logit",
  eps = 0,
  type = c("Fieller", "Delta"),
  maxg = 0.99
)

fieller2(
  phat,
  b,
  vv,
  df.t = Inf,
  offset = 0,
  logscale = FALSE,
  link = "fpower",
  lambda = 0,
  eps = 0,
  type = c("Fieller", "Delta"),
  maxg = 0.99
)

Arguments

phat

Mortality proportion

b

Length 2 vector of intercept and slope

vv

Variance-covariance matrix for intercept and slope

df.t

Degrees of freedom for variance-covariance matrix

offset

Offset to be added to intercept. This can be of length 2, in order to return values on the original scale, in the case where b and vv are for values that have been scaled by subtracting offset[1] and dividing by offset[2].

logscale

Should confidence limits be back transformed from log scale?

link

Link function that transforms expected mortalities to the scale of the linear predictor

eps

If eps>0 phat is replaced by \frac{p+\epsilon}{1+2*\epsilon} before applying the transformation.

type

The default is to use Fieller's formula. The Delta (type="Delta") method, which relies on a first order Taylor series approximation to the variance, is provided so that it can be used for comparative purposes. It can be reliably used only in cases where the interval has been shown to be essentially the same as given by type="Fieller"!

maxg

Maximum value of g for which a confidence interval will be calculated. Must be < 1.

lambda

The power \lambda, when using the link="fpower". (This applies to fieller2 only.)

Details

See the internal code for details of the value g. The calculation gives increasing wide confidence intervals as g approaches 1. If g>=1, there are no limits. The default value for df.t is a rough guess at what might be reasonable. For models fitted using lme4::lmer(), abilities in the lmerTest package can be used to determine a suitable degrees of freedom approximation — this does not extend to use with glmer() or glmmTMB.

Value

A vector, with elements

est

Estimate

var

Variance, calculated using the Delta method

lwr

Lower bound of confidence interval

upr

upper bound of confidence interval

g

If g is close to 0 (perhaps g < 0.05), confidence intervals will be similar to those calculated using the Delta method, and the variance can reasonably be used for normal theory inference.

References

Joe Hirschberg & Jenny Lye (2010) A Geometric Comparison of the Delta and Fieller Confidence Intervals, The American Statistician, 64:3, 234-241, DOI: 10.1198/ tast.2010.08130

E C Fieller (1944). A Fundamental Formula in the Statistics of Biological Assay, and Some Applications. Quarterly Journal of Pharmacy and Pharmacology, 17, 117-123.

David J Finney (1978). Statistical Method in Biological Assay (3rd ed.), London, Charles Griffin and Company.

See Also

varRatio

Examples

redDel <- subset(qra::codling1988, Cultivar=="Red Delicious")
redDel.glm <- glm(cbind(dead,total-dead)~ct, data=redDel,
                  family=quasibinomial(link='cloglog'))
vv <- summary(redDel.glm)$cov.scaled
fieller(0.99, b=coef(redDel.glm), vv=vv, link='cloglog')

Title Function to calculate ratio of p+eps to 1-p+eps.

Description

This is a convenience function that returns \frac{p+\epsilon}{1-p+\epsilon}. It calculates the argument that is supplied to the log function in Tukey's ‘flog’.

Usage

foldp(p, eps)

Arguments

p

Proportion

eps

Offset. The choice eps=0.01 has the same effect as replacing \frac{r}{n-r} by \frac{r+0.5}{n-r+0.5} when n=50, or by \frac{r+1}{n-r+1} when n=100

Value

(p+eps)/(1-p+eps)

Examples

foldp(c(0.2,0.75), 0)

Folded Power Transformation

Description

The name “folded Power Transformation” is used because this does for power transformations what Tukey's folded logarithm does for the logarithmic tranformation. The function calculates

f(p, \lambda, \epsilon) = \frac{p+\epsilon}{1-p+\epsilon}^\lambda

where \lambda is the power and \epsilon is a positive offset that ensures that \frac{p+\epsilon}{1-p+\epsilon} is greater than 0 and finite.

Usage

fpower(p, lambda, eps)

Arguments

p

Mortality proportion

lambda

Power lambda for the power transformation

eps

If eps>0 phat is replaced by \frac{p+\epsilon}{1+\epsilon} before applying the power transformation.

Value

The transformed values of fpower(p).

Examples

p <- (0:10)/10
ytrans <- fpower(p, lambda=0.25, eps=1/450)

Extract estimates of the intra-class correlation from a glmmTMB model object with beta-binomial error.

Description

The intra-class correlation is calculated as (1+exp(\theta))^{-1}, where \theta is the estimate given by the formula specified in the argument dispformula.

Usage

getRho(obj, varMult = FALSE)

Arguments

obj

glmmTMB model object with betabinomial error, and with a 'dispformula' argument supplied.

varMult

If TRUE return, in addition to rho, the factor mult by which the variance is inflated relative to the binomial.

Details

The variance for the betabinomial model is then obtained by multiplying the binomial variance by 1+(n-1)\rho, where $n$ is the binomial 'size'.

Value

if varMult==FALSE return (as a vector) the estimates \rho, else (varMult==TRUE) return list(rho, mult).

Examples

pcheck <- suppressWarnings(requireNamespace("glmmTMB", quietly = TRUE))
if(pcheck) pcheck & packageVersion("glmmTMB") >= "1.1.2"
if(pcheck){
form <- cbind(Dead,Live)~0+trtGp/TrtTime+(1|trtGpRep)
HawMed <- droplevels(subset(HawCon, CN=="MedFly"&LifestageTrt!="Egg"))
HawMed <- within(HawMed,
                 {trtGp <- factor(paste0(CN,LifestageTrt, sep=":"))
                 trtGpRep <- paste0(CN,LifestageTrt,":",RepNumber)
                 scTime <- scale(TrtTime) })
HawMedbb.TMB <- glmmTMB::glmmTMB(form, dispformula=~trtGp+splines::ns(scTime,2),
                                 family=glmmTMB::betabinomial(link="cloglog"),
                                 data=HawMed)
rho <- qra::getRho(HawMedbb.TMB)} else
message("Example requires `glmmTMB` version >= 1.1.2: not available")

Extract scaling coefficients from vector returned by scale()

Description

The function scale() replaces x by (x-a)/b, where a is mean(x) and b is sd(x). The quantities a and b are available as attributes of the object that is returned.

Usage

getScaleCoef(z)

Arguments

z

Object returned by scale()

Details

Use of a scaled explanatory variable can be helpful in getting a model to fit. The scaling coefficient(s) will then be needed when the fitted model is used with explanatory variable values on the original scale.

Value

A vector, whose elements are the scaling coefficients a and b, or if scale=FALSE then a.

Examples

z <- scale(1:10)
qra::getScaleCoef(z)

Use given vector to identify groups with specified categories

Description

Any one-dimensional object whose values distinguish groups may be supplied.

Usage

gpsWithin(x, f)

Arguments

x

One-dimensional object whose values distinguish groups

f

One-dimensional object or list of objects, the combinations of whose values specify categories within which groups are to be defined.

Value

Integer vector whose values, within each specified category, run from 1 to the number of groups

Examples

repnum <- with(qra::codling1988, gpsWithin(cultRep, Cultivar))
table(codling1988$Cultivar,repnum)

Draw graphs of insect mortality or other exposure-response data

Description

Datasets that are in mind hold, for each replicate of each combination of each of a several factors (e.g., species, lifestages, temperatures), mortalities for each of a number of values of "dose". See for example the dataset help page codling1989.

Usage

graphSum(
  df,
  subSet = NULL,
  link = "cloglog",
  logScale = FALSE,
  dead = "Dead",
  tot = "Tot",
  dosevar = "logCT",
  Rep = "Rep",
  fitRep = NULL,
  fitPanel = NULL,
  byFacet = ~Species,
  layout = NULL,
  maint = "Codling Moth, MeBr",
  ptSize = 2,
  xzeroOffsetFrac = 0.08,
  yzeroOneOffsets = c(-0.08, 0.08),
  yEps = 0.005,
  xlab = expression(bold("CT ") * "(gm.h." * m^{
     -3
 } * ")"),
  ylabel = NULL,
  ytiklab = c(0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)
)

Arguments

df

Data frame from which data will be taken

subSet

NULL, or an expression, such as for example expression(LifeStage=='Eggs')) that evaluates to a logical that specifies the required data subset. If not NULL then the subsetting information is pasted on after the main title

link

Link function. If character, obtain from make.link. Alternatively, a function may be supplied as argument.

logScale

Logical, indicating whether the dose ($x$-variable) is on a log scale.

dead

Character; name of column holding number dead

tot

Character; column holding total number

dosevar

Character; column holding "dose" values

Rep

Character; NULL, or column holding replicate number, within panel

fitRep

Character; NULL, or column holding replicate fitted values

fitPanel

Character; NULL, or column holding panel fitted values

byFacet

Graphics formula specifying factor combination that determines panel

layout

Graphics formula that can be supplied to grid_facet

maint

Main title

ptSize

Pointsize, by default 2

xzeroOffsetFrac

$x$-axis zero offset fraction, required when scale is logarithmic

yzeroOneOffsets

Length two vector, giving 0 100 mortalities, on the scale of the link function.

yEps

Fractional increase at bottom and top of $y$ user range to accommodate points for mortalities of 0 and 1.

xlab

Expression specifying x-axis label

ylabel

If not NULL, $y$-axis label

ytiklab

Place $y$ axis tiks and labels at these mortalities

Value

No return value, called for side effects

Kerrich Coin Toss Trial Outcomes

Description

A data set containing 2,000 trials of coin flips from statistician John Edmund Kerrich's 1940s experiments while imprisoned by the Nazis during World War Two.

Usage

data("kerrich")

Format

The format is: List of 1 $ : chr [1:2000] "0" "0" "0" "1" ...

Source

https://en.wikipedia.org/wiki/John_Edmund_Kerrich

References

Kerrich, J. E. (1950). An experimental introduction to the theory of probability. Belgisk Import Company.

Examples

data(kerrich)

Number of males among first 12 in families of 13 children

Description

The number of male children among the first 12 children of family size 13 in 6115 families taken from the hospital records in the nineteenth century Saxony (Lindsey (1995), p.59). The thirteenth child is ignored to assuage the effect of families non-randomly stopping when a desired gender is reached.

Usage

data("malesINfirst12")

Format

A data frame with 13 observations on the following 2 variables.

No_of_Males

a numeric vector

freq

a numeric vector

Details

Data are available in the fitODBOD package.

Source

fitODBOD package

References

Edwards, A. W. F. (1958). An analysis of Geissler's data on the human sex ratio. Annals of human genetics, 23(1), 6-15.

Geissler, A. (1889) Beiträge zur Frage des Geschlechtsverhältnisses der Geborenen. Z. Köngl. Sächs. Statist. Bur., 35, 1±24.

Lindsey, J. K., & Altham, P. M. E. (1998). Analysis of the human sex ratio by using overdispersion models. Journal of the Royal Statistical Society: Series C (Applied Statistics), 47(1), 149-157.

Examples

data(malesINfirst12)
boxplot(freq ~ No_of_Males, data=malesINfirst12)

Incidence of ray blight disease of pyrethrum

Description

An assessment of the incidence of ray blight disease of pyrethrum in 62 sampling units, containing 6 plants each.

Usage

data("rayBlight")

Format

The format is: int [1:62] 4 6 6 6 6 6 6 6 4 6 ...

Source

epiphy package.

References

Pethybridge SJ, Esker P, Hay F, Wilson C, Nutter FW. 2005. Spatiotemporal description of epidemics caused by Phoma ligulicola in Tasmanian pyrethrum fields. Phytopathology 95, 648-658.

Examples

data(rayBlight)
barplot(table(rayBlight))

Estimate dispersion as a function of predicted values

Description

A loess smooth is applied to the square roots of the standardized deviance residuals. The inverses of values from the smooth, raised to the power of lambda, are then used as prior weights to update the model. A value of lambda that is a little more than 2.0 has often worked well.

Usage

scaleLocAdjust(x, lambda = 2, span = 0.75)

Arguments

x

Model fitted using glm or, possibly lm

lambda

Power of smooth of square roots of absolute values of residuals, to try for values whose inverses will be used as weights

span

span parameter for use in smoothing the square root of standardized deviance residuals.

Details

This function is primarily for experimental use, in investigating possible ways to model a dispersion factor that varies with the fitted value.

Value

A list, with elements

model

Model updated to use the newly calculated weights

estDisp

Estimated dispersions

Note

The dispersion estimates that correspond to the updated model are obtained by dividing the dispersion value given by summary() for the updated model by the (prior) weights supplied when the model was updated. The approach for obtaining varying dispersion estimates is used because, empirically, it has been found to work well for at least some sets of data. In particular, there seems no obvious theoretical basis for the choice of lambda. In the example given, used because the data is publicly available, the method has limited success.

See Also

checkDisp

Examples

ROYAL <- subset(qra::codling1988, Cultivar=="ROYAL")
ROYAL.glm <- glm(cbind(dead,total-dead)~ct, data=ROYAL,
                  family=quasibinomial(link='cloglog'))
ROYALFix <- qra::scaleLocAdjust(ROYAL.glm)
## Check range of indicated prior weights
range(ROYALFix[[2]])
## Range of updated dispersion estimates
range(summary(ROYALFix[[1]])[['dispersion']]/ROYALFix[[2]])

First order approximation to variance of y-ordinate to slope ratio

Description

In contexts where an LD99 estimate will be used as a data value in a further analysis step, the inverse of the variance may be used as a weight. The y-ordinate is for the link function transformed value of a specified mortality proportion, commonly 0.50, or 0.90, or 0.99

Usage

varRatio(phat = 0.99, b, vv, link = "cloglog")

Arguments

phat

Mortality proportion

b

Length 2 vector of intercept and slope

vv

Variance-covariance matrix for intercept and slope

link

Link function that transforms expected mortalities to the scale of the linear predictor

Details

This function should only be used, in order to speed up calculations that use the function fieller (call fieller with (type="Delta")), in a context where it is to be used many times, and where a check has been made that its use leads to confidence intervals that are a close approximation to those given with the default argument (type="Fieller").

Value

A vector, with elements

xhat

Estimate

var

Variance, calculated using the Delta method, See the help page for fieller for further details and references.

Examples

redDel <- subset(qra::codling1988, Cultivar=="Red Delicious")
redDel.glm <- glm(cbind(dead,total-dead)~ct, data=redDel,
                  family=quasibinomial(link='cloglog'))
vv <- summary(redDel.glm)$cov.scaled
qra::varRatio(0.99, b=coef(redDel.glm), vv=vv, link="cloglog")