fitstahl: Fit Stahl interference model

In qtl: Tools for Analyzing QTL Experiments

Fit Stahl interference model

Description

Fit the Stahl model for crossover inference (or the chi-square model, which is a special case).

Usage

fitstahl(cross, chr, m, p, error.prob=0.0001, maxit=4000, tol=1e-4,
         maxm=15, verbose=TRUE)

Arguments

cross	An object of class cross. See read.cross for details.
chr	Optional vector indicating the chromosomes to consider. This should be a vector of character strings referring to chromosomes by name; numeric values are converted to strings. Refer to chromosomes with a preceding - to have all chromosomes but those considered. A logical (TRUE/FALSE) vector may also be used.
m	Interference parameter (a non-negative integer); if unspecified, this is estimated.
p	The proportion of chiasmata coming from the no interference mechanism in the Stahl model (0 <= p <= 1). p=0 gives the chi-square model. If unspecified, this is estimated.
error.prob	The genotyping error probability. If = NULL, it is estimated.
maxit	Maximum number of iterations to perform.
tol	Tolerance for determining convergence.
maxm	Maximum value of m to consider, if m is unspecified.
verbose	Logical; indicates whether to print tracing information.

Details

This function is currently only available for backcrosses and intercrosses.

The Stahl model of crossover interference (of which the chi-square model is a special case) is fit. In the chi-square model, points are tossed down onto the four-strand bundle according to a Poisson process, and every (m+1)st point is a chiasma. With the assumption of no chromatid interference, crossover locations on a random meiotic product are obtained by thinning the chiasma process. The parameter m (a non-negative integer) governs the strength of crossover interference, with m=0 corresponding to no interference.

In the Stahl model, chiasmata on the four-strand bundle are a superposition of chiasmata from two mechanisms, one following a chi-square model and one exhibiting no interference. An additional parameter, p, gives the proportion of chiasmata from the no interference mechanism.

If all of m, p, and error.prob are specified, any of them with length > 1 must all have the same length.

If m is unspecified, we do a grid search starting at 0 and stop when the likelihood decreases (thus assuming a single mode), or maxm is reached.

Value

A matrix with four columns: m, p, error.prob, and the log likelihood.

If specific values for m, p, error.prob are provided, the log likelihood for each set are given.

If some are left unspecified, the maximum likelihood estimates are provided in the results.

Author(s)

Karl W Broman, [e-mail]

References

Armstrong, N. J., McPeek, M. J. and Speed, T. P. (2006) Incorporating interference into linkage analysis for experimental crosses. Biostatistics 7, 374–386.

Zhao, H., Speed, T. P. and McPeek, M. S. (1995) Statistical analysis of crossover interference using the chi-square model. Genetics 139, 1045–1056.

See Also

est.map, sim.cross

Examples

# Simulate genetic map: one chromosome of length 200 cM with
# a 2 cM marker spacing
mymap <- sim.map(200, 51, anchor.tel=TRUE, include.x=FALSE,
                 sex.sp=FALSE, eq.spacing=TRUE)

# Simulate data under the chi-square model, no errors
mydata <- sim.cross(mymap, n.ind=250, type="bc",
                    error.prob=0, m=3, p=0)

# Fit the chi-square model for specified m's
## Not run: output <- fitstahl(mydata, m=1:5, p=0, error.prob=0)

plot(output$m, output$loglik, lwd=2, type="b")

# Find the MLE of m in the chi-square model
## Not run: mle <- fitstahl(mydata, p=0, error.prob=0)

## Not run: 
# Simulate data under the Stahl model, no errors
mydata <- sim.cross(mymap, n.ind=250, type="bc",
                    error.prob=0, m=3, p=0.1)

# Find MLE of m for the Stahl model with known p
mle.stahl <- fitstahl(mydata, p=0.1, error.prob=0)

# Fit the Stahl model with unknown p and m,
# get results for m=0, 1, 2, ..., 8
output <- fitstahl(mydata, m=0:8, error.prob=0)
plot(output$m, output$loglik, type="b", lwd=2)
## End(Not run)

qtl documentation built on Sept. 11, 2024, 5:43 p.m.