qtl (version 1.66)

fitstahl: 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)

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.

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.

Author

Karl W Broman,

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.

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

set.seed(12201969)
# 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
if (FALSE) output <- fitstahl(mydata, m=1:5, p=0, error.prob=0)
output <- fitstahl(mydata, m=1, p=0, error.prob=0,tol=0.001)
plot(output$m, output$loglik, lwd=2, type="b")

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

if (FALSE) {
# 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)}