Movatterモバイル変換

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 [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'sformula to calculate 95of the functions noted below are useful ancillaries tofiellerandfieller2, notablyfoldp,fpower,extractLT,andgetScaleCoef.

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

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

foldp: Calculates the ratio ofp+eps to1-code+eps

getRho Extracts estimates of the intra-classcorrelation from a glmmTMB model object with betabinomial error.See the vignette [timeMortality] for details of the parametizationused for thebetabinomial error.

getScaleCoef: Extracts the scale coefficients from a vectorthat has been scaled usingscale, as needed so that the scalingcan be undone.

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

scaleLocAdjust: Returns, forglmmTMB models with abetabinomial error, dispersion factors (i.e., multipliers for thebinomial variance) as functions of predicted values.

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

Author(s)

Maintainer: John Maindonaldjohn@statsresearch.co.nz (ORCID)

See Also

Useful links:

• https://github.com/jhmaindonald/qra

• Report bugs athttps://github.com/jhmaindonald/qra/issues

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 forHawCon in theColdData has furtherdetails.

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 forplot(x.lm, which=3),wherex.lm is a linear model or a generalized linear model.Plot the object returned to assess how successful the weights,determined using the functionscaleLocAdjust, have beenin adjusting for heterogenous variances.

Usage

checkDisp(x, span = 0.75)

Arguments

Value

A data frame, with:

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 weightsrange(royalFix[[2]])## Range of updated dispersion estimatesrange(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 mothto fumigation with methyl bromide, for the year 1988 only

Usage

data(codling1988)data(codling1989)

Format

A data frame with 77 observations (codling1988), and with 40observations (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 levelsBRAEBURNFUJIGRANNYRed Delicious andROYAL;and with 1989 levelsGala,Red DeliciousandSplendour

rep

replicate number, withinCultivar

cultRep

Cultivar/replicate combination

Details

The research that generated these data was in part funded by New Zealandpipfruit growers. The published analysis was funded by New Zealandpipfruit growers. See alsoDAAG::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 fitteddose-response lines for each of several levels of a factor,extractLT calls the functionfieller 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

Details

Fixed coefficients fromobj must be for intercepts andfor slopes. Starting the model formula with0+ will commonlydo what is required. The coefficientsfixef(obj)[a] are assumedto specify line intercepts, whilefixef(obj)[b] specify thecorresponding slopes. These replace the argumentsnEsts(subscripts for intercepts were1:nEsts) andslopeAdd(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)} elsemessage("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 confidenceinterval for a specified mortality proportion, commonly0.50, or 0.90, or 0.99. Here "dose" is a generic term forany measure of intensity of a treatment that is designedto 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

Details

See the internal code for details of the valueg.The calculation gives increasing wide confidence intervals asg approaches 1. Ifg>=1, there are no limits.The default value fordf.t is a rough guess at whatmight be reasonable. For models fitted usinglme4::lmer(),abilities in thelmerTest package can be used to determinea suitable degrees of freedom approximation — this does notextend to use withglmer() orglmmTMB.

Value

A vector, with elements

References

Joe Hirschberg & Jenny Lye (2010) A GeometricComparison 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 Statisticsof Biological Assay, and Some Applications. QuarterlyJournal 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.scaledfieller(0.99, b=coef(redDel.glm), vv=vv, link='cloglog')

Title Function to calculate ratio ofp+eps to1-p+eps.

Description

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

Usage

foldp(p, eps)

Arguments

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 becausethis does for power transformations what Tukey's folded logarithmdoes 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 positiveoffset that ensures that\frac{p+\epsilon}{1-p+\epsilon} isgreater than 0 and finite.

Usage

fpower(p, lambda, eps)

Arguments

Value

The transformed values offpower(p).

Examples

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

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

Description

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

Usage

getRho(obj, varMult = FALSE)

Arguments

Details

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

Value

ifvarMult==FALSE return (as a vector) the estimates\rho, else (varMult==TRUE) returnlist(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)} elsemessage("Example requires `glmmTMB` version >= 1.1.2: not available")

Extract scaling coefficients from vector returned byscale()

Description

The functionscale() replacesx by(x-a)/b,wherea ismean(x) andb issd(x).The quantitiesa andb are available as attributesof the object that is returned.

Usage

getScaleCoef(z)

Arguments

Details

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

Value

A vector, whose elements are the scaling coefficientsa andb, or ifscale=FALSE thena.

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 groupsmay be supplied.

Usage

gpsWithin(x, f)

Arguments

Value

Integer vector whose values, within each specifiedcategory, 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 ofeach combination of each of a several factors (e.g.,species, lifestages, temperatures), mortalities foreach of a number of values of "dose". See for examplethe dataset help pagecodling1989.

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

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 thefitODBOD 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 standardizeddeviance residuals. The inverses of values from the smooth, raisedto the power oflambda, are then used as prior weights toupdate the model. A value oflambda that is a little morethan 2.0 has often worked well.

Usage

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

Arguments

Details

This function is primarily for experimental use, in investigatingpossible ways to model a dispersion factor that varies with thefitted value.

Value

A list, with elements

Note

The dispersion estimates that correspond to the updatedmodel are obtained by dividing the dispersion value givenbysummary() for the updated model by the (prior) weightssupplied when the model was updated. The approach for obtainingvarying dispersion estimates is used because, empirically, ithas been found to work well for at least some sets of data. Inparticular, there seems no obvious theoretical basis for thechoice oflambda. In the example given, used because thedata 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 weightsrange(ROYALFix[[2]])## Range of updated dispersion estimatesrange(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 valuein a further analysis step, the inverse of the variance may beused as a weight. The y-ordinate is for the link functiontransformed value of a specified mortality proportion, commonly0.50, or 0.90, or 0.99

Usage

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

Arguments

Details

This function should only be used, in order to speed upcalculations that use the functionfieller(callfieller 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 toconfidence intervals that are a close approximation to thosegiven with the default argument (type="Fieller").

Value

A vector, with elements

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.scaledqra::varRatio(0.99, b=coef(redDel.glm), vv=vv, link="cloglog")