xoi (version 0.72)

fitGamma: Fit Gamma model

Description

Fit the gamma model for crossover interference to data on crossover locations.

Usage

fitGamma(
  d,
  censor = NULL,
  nu = NULL,
  lo = NULL,
  hi = NULL,
  se = FALSE,
  supint = FALSE,
  rescale = FALSE,
  drop = 1.5,
  tol = 0.00001,
  maxit = 1000,
  max.conv = 25,
  integr.tol = 0.00000001,
  max.subd = 1000,
  min.subd = 10,
  h = 0.1,
  hstep = 1.5
)

Value

If nu is specified, we return a data frame with two columns: nu and the corresponding log (base e) likelihood. If rescale=TRUE, the maximum log likelihood is subtracted off, so that its maximum is at 0.

If lo and hi is specified, the output contains a single row with the MLE of \(\nu\) and the corresponding log likelihood. If se=TRUE, we also include the estimated SE. If supint=TRUE, we include two additional rows with the lower and upper limits of the likelihood support interval.

Arguments

d

A vector of inter-crossover distances in cM. This should include distances from start of chromosome to first crossover, last crossover to end of chromosome, and chromosome length, if there are no crossovers.

Alternatively, this may be a matrix with the first column being the distances and second column being the censoring types (censor).

censor

A vector of the same length as d, indicating the censoring type for each distance. 0 = uncensored, 1 = right-censored, 2 = initial crossover on chromosome, 3 = whole chromosome.

nu

A vector of interference parameters (\(\nu\)) at which to calculate the log likelihood. If NULL, lo and hi must be specified.

lo

If nu is unspecified, lo indicates the lower value of the interval in which to search for the MLE. If supint=TRUE, this should be below the lower limit of the support interval.

hi

If nu is unspecified, hi indicates the upper value of the interval in which to search for the MLE. If supint=TRUE, this should be above the upper limit of the support interval.

se

If TRUE and nu was not specified, an estimated SE (based on the second derivative of the log likelihood) is estimated.

supint

If TRUE and nu was not specified, a likelihood support interval is calculated, with drop being the amount to drop in log (base 10).

rescale

If TRUE and nu was specified, re-scale the log likelihoods so that the maximum is at 0.

drop

If supint was specified, this indicates the amount to drop in log (base 10) for the likelihood support interval.

tol

Tolerance for converence to calculate the likelihood, SE, and likelihood support interval.

maxit

Maximum number of iterations in estimating the SE and likelihood support interval.

max.conv

Maximum limit for summation in the convolutions to get inter-crossover distance distribution from the inter-chiasma distance distributions. This should be greater than the maximum number of chiasmata on the 4-strand bundle.

integr.tol

Tolerance for convergence of numerical integration.

max.subd

Maximum number of subdivisions in numerical integration.

min.subd

Minimum number of subdivisions in numerical integration.

h

Step used in estimating the second derivative of the log likelihood.

hstep

factor by which h is decreased in each iteration of the estimation of the second derivative of the log likelihood.

Author

Karl W Broman,

Details

See Broman and Weber (2000) for details of the method.

We use R's stats::integrate() function for numerical integrals, stats::optimize() for optimizing the likelihood, and stats::uniroot() for identifying the endpoints of the likelihood support interval.

References

Broman, K. W. and Weber, J. L. (2000) Characterization of human crossover interference. Am. J. Hum. Genet. 66, 1911--1926.

Broman, K. W., Rowe, L. B., Churchill, G. A. and Paigen, K. (2002) Crossover interference in the mouse. Genetics 160, 1123--1131.

McPeek, M. S. and Speed, T. P. (1995) Modeling interference in genetic recombination. Genetics 139, 1031--1044.

See Also

qtl::fitstahl()

Examples

data(bssbsb)
bssbsb <- bssbsb[,1:50]

xodist <- convertxoloc(find.breaks(bssbsb, chr=1))

# plot a rough log likelihood curve
if (FALSE) out <- fitGamma(xodist, nu=seq(1, 19, by=2))
out <- fitGamma(xodist, nu=seq(1, 19, by=2), tol=0.001)
plot(out, type="l", lwd=2)

# get MLE
if (FALSE) mle <- fitGamma(xodist, lo=8, hi=12)
mle <- fitGamma(xodist, lo=8, hi=12, tol=0.001)
mle

abline(v=mle[1], h=mle[2], col="blue", lty=2)

# get MLE and SE
if (FALSE) mle <- fitGamma(xodist, lo=9.5, hi=10.5, se=TRUE)
mle <- fitGamma(xodist, lo=9.5, hi=10.5, se=TRUE, tol=0.001)
mle

# get MLE and 10^1.5 support interval
if (FALSE) int <- fitGamma(xodist, lo=1, hi=20, supint=TRUE)
int <- fitGamma(xodist, lo=1, hi=20, supint=TRUE, tol=0.001)
int
abline(v=mle[2:3,1], h=mle[2:3,2], col="red", lty=2)