Package 'qra'

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

---

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)

---

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")

---