EVOLUTION-MANAGER

Edit File: stepwiseqtl.html

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"><html xmlns="http://www.w3.org/1999/xhtml"><head><title>R: Stepwise selection for multiple QTL</title>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<link rel="stylesheet" type="text/css" href="R.css" />
</head><body>

<table width="100%" summary="page for stepwiseqtl {qtl}"><tr><td>stepwiseqtl {qtl}</td><td style="text-align: right;">R Documentation</td></tr></table>

<h2>Stepwise selection for multiple QTL</h2>

<h3>Description</h3>

Performs forward/backward selection to identify a multiple QTL model,
with model choice made via a penalized LOD score, with separate
penalties on main effects and interactions.

<h3>Usage</h3>

<pre>
stepwiseqtl(cross, chr, pheno.col=1, qtl, formula, max.qtl=10, covar=NULL,
 method=c("imp", "hk"), model=c("normal", "binary"),
 incl.markers=TRUE, refine.locations=TRUE,
 additive.only=FALSE, scan.pairs=FALSE, penalties,
 keeplodprofile=TRUE, keeptrace=FALSE, verbose=TRUE,
 tol=1e-4, maxit=1000, require.fullrank=FALSE)
</pre>

<h3>Arguments</h3>

<table summary="R argblock">
<tr valign="top"><td><code>cross</code></td>
<td>
An object of class <code>cross</code>. See
<code><a href="read.cross.html">read.cross</a></code> for details.
</td></tr>
<tr valign="top"><td><code>chr</code></td>
<td>
Optional vector indicating the chromosomes to consider in
search for QTL. 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 <code>-</code> to have all
chromosomes but those considered. A logical (TRUE/FALSE) vector may
also be used.
</td></tr>
<tr valign="top"><td><code>pheno.col</code></td>
<td>
Column number in the phenotype matrix which should be
used as the phenotype. One may also give character strings matching
the phenotype names. Finally, one may give a numeric vector of
phenotypes, in which case it must have the length equal to the number
of individuals in the cross, and there must be either non-integers or
values &lt; 1 or &gt; no. phenotypes; this last case may be useful for studying
transformations.
</td></tr>
<tr valign="top"><td><code>qtl</code></td>
<td>
Optional QTL object (of class <code>"qtl"</code>, as created by
<code><a href="makeqtl.html">makeqtl</a></code>) to use as a starting point.
</td></tr>
<tr valign="top"><td><code>formula</code></td>
<td>
Optional formula to define the QTL model to be used as a
starting point.
</td></tr>
<tr valign="top"><td><code>max.qtl</code></td>
<td>
Maximum number of QTL to which forward selection should
proceed.
</td></tr>
<tr valign="top"><td><code>covar</code></td>
<td>
Data frame of additive covariates.
</td></tr>
<tr valign="top"><td><code>method</code></td>
<td>
Indicates whether to use multiple imputation or
Haley-Knott regression.
</td></tr>
<tr valign="top"><td><code>model</code></td>
<td>
The phenotype model: the usual model or a model for binary
traits
</td></tr>
<tr valign="top"><td><code>incl.markers</code></td>
<td>
If FALSE, do calculations only at points on an
evenly spaced grid.
</td></tr>
<tr valign="top"><td><code>refine.locations</code></td>
<td>
If TRUE, use <code><a href="refineqtl.html">refineqtl</a></code> to
refine the QTL locations after each step of forward and backward
selection.
</td></tr>
<tr valign="top"><td><code>additive.only</code></td>
<td>
If TRUE, allow only additive QTL models; if FALSE,
consider also pairwise interactions among QTL.
</td></tr>
<tr valign="top"><td><code>scan.pairs</code></td>
<td>
If TRUE, perform a two-dimensional, two-QTL scan at
each step of forward selection.
</td></tr>
<tr valign="top"><td><code>penalties</code></td>
<td>
Vector of three (or six) values indicating the penalty
on the number of QTL terms. If three values, these are the penalties
on main effects and heavy and light penalties on interactions. If
six values, these include X-chr-specific penalties, and the values
are: main effect for autosomes, main effect for X chr, heavy penalty
on A:A interactions, light penalty on A:A interactions, penalty on
A:X interactions, and penalty on X:X interactions. See the Details
below. If missing, default values are used that are based on
simulations of backcrosses and intercrosses with genomes modeled
after that of the mouse.
</td></tr>
<tr valign="top"><td><code>keeplodprofile</code></td>
<td>
If TRUE, keep the LOD profiles from the last
iteration as attributes to the output.
</td></tr>
<tr valign="top"><td><code>keeptrace</code></td>
<td>
If TRUE, keep information on the sequence of models
visited through the course of forward and backward selection as an
attribute to the output.
</td></tr>
<tr valign="top"><td><code>verbose</code></td>
<td>
If TRUE, give feedback about progress. If
<code>verbose</code> is an integer &gt; 1, even more information is printed.
</td></tr>
<tr valign="top"><td><code>tol</code></td>
<td>
Tolerance for convergence for the binary trait model.
</td></tr>
<tr valign="top"><td><code>maxit</code></td>
<td>
Maximum number of iterations for fitting the binary trait
model.
</td></tr>
<tr valign="top"><td><code>require.fullrank</code></td>
<td>
If TRUE, give LOD=0 when covariate matrix in
the linear regression is not of full rank.
</td></tr>
</table>

<h3>Details</h3>

We seek to identify the model with maximal penalized LOD score. The
penalized LOD score, defined in Manichaikul et al. (2009),
is the LOD score for the model (the log10 likelihood
ratio comparing the model to the null model with no QTL) with
penalties on the number of QTL and QTL:QTL interactions.

We consider QTL models allowing pairwise interactions among QTL but
with an enforced hierarchy in which inclusion of a pairwise
interaction requires the inclusion of both of the corresponding main
effects. Additive covariates may be included, but currently we do not
explore QTL:covariate interactions. Also, the penalized LOD score
criterion is currently defined only for autosomal loci, and results
with the X chromosome should be considered with caution.

The penalized LOD score is of the form pLOD(g) = LOD(g)
 - Tm pm - Th ph - Tl pl where g denotes a model,
pm is the number of QTL in the model (&quot;main effects&quot;),
ph is the number of pairwise interactions that will be
given a heavy interaction penalty, pl is the number of pairwise
interactions that will be given a light interaction penalty,
Tm is the penalty on main effects, Th is the heavy
interaction penalty, and Tl is the light interaction
penalty. The <code>penalties</code> argument is the vector (Tm, Th, Tl). If Tl is missing (<code>penalties</code>
has a vector of length 2), we assume Tl = Th, and so
all pairwise interactions are assigned the same penalty.

The &quot;heavy&quot; and &quot;light&quot; interaction penalties can be a bit
confusing. Consider the clusters of QTL that are connected via one or
more pairwise interactions. To each such cluster, we assign at most
one &quot;light&quot; interaction penalty, and give all other pairwise
interactions the heavy interaction penalty. In other words, if
pi is the total number of pairwise interactions for a QTL
model, we let pl be the number of clusters of connected QTL
with at least one pairwise interaction, and then let ph = pi - pl.

Let us give an explicit example. Consider a model with 6 QTL, and
with interactions between QTL 2 and 3, QTL 4 and 5 and QTL 4 and 6
(so we have the model formula
<code>y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q2:Q3 + Q4:Q5 + Q4:Q6</code>).
There are three clusters of connected QTL: (1), (2,3) and (4,5,6). We
would assign 6 main effect penalties (Tm), 2 light
interaction penalties (Tl), and 1 heavy interaction penalty
(Th).

Manichaikul et al. (2009) described a system for deriving the
three penalties on the basis of permutation results from a
two-dimensional, two-QTL genome scan (as calculated with
<code><a href="scantwo.html">scantwo</a></code>). These may be calculated with the
function <code><a href="calc.penalties.html">calc.penalties</a></code>.

A forward/backward search method is used, with the aim to optimize the
penalized LOD score criterion. That is, we seek to identify the
model with maximal the penalized LOD score. The search algorithm was
based closely on an algorithm described by Zeng et al. (1999).

We use forward selection to a model of moderate size (say 10 QTL),
followed by backward elimination all the way to the null model. The
chosen model is that which optimizes the penalized LOD score
criterion, among all models visited. The detailed algorithm is as
follows. Note that if <code>additive.only=TRUE</code>, no pairwise
interactions are considered.

<ol>
<li> Start at the null model, and perform a single-QTL genome scan,
and choose the position giving the largest LOD score. If
<code>scan.pairs=TRUE</code>, start with a two-dimensional, two-QTL genome
scan instead. If an initial QTL model were defined through the
arguments <code>qtl</code> and <code>formula</code>, start with this model and
jump immediately to step 2.

</li>
<li> With a fixed QTL model in hand:

<ol>
<li> Scan for an additional additive QTL.

</li>
<li> For each QTL in the current model, scan for an additional
interacting QTL.

</li>
<li> If there are &gt;= 2 QTL in the current model,
consider adding one of the possible pairwise interactions.

</li>
<li> If <code>scan.pairs=TRUE</code> perform a two-dimensional, two-QTL
scan, seeking to add a pair of novel QTL, either additive or
interacting.

</li>
<li> Step to the model that gives the largest value for the model
comparison criterion, among those considered at the current step.

</li></ol>

</li>
<li> Refine the locations of the QTL in the current model (if
<code>refine.locations=TRUE</code>).

</li>
<li> Repeat steps 2 and 3 up to a model with some pre-determined
number of loci.

</li>
<li> Perform backward elimination, all the way back to the null
model. At each step, consider dropping one of the current main
effects or interactions; move to the model that maximizes the model
comparison criterion, among those considered at this step. Follow
this with a refinement of the locations of the QTL.

</li>
<li> Finally, choose the model having the largest model comparison
criterion, among all models visited.

</li></ol>

In this forward/backward algorithm, it is likely best to build up to
an overly large model and then prune it back. Note that there is no
&quot;stopping rule&quot;; the chosen model is that which optimizes the model
comparison criterion, among all models visited. The search can be
time consuming, particularly if a two-dimensional scan is performed at
each forward step. Such two-dimensional scans may be useful for
identifying QTL linked in repulsion (having effects of opposite sign)
or interacting QTL with limited marginal effects, but our limited
experience suggests that they are not necessary; important linked or
interacting QTL pairs can be picked up in the forward selection to a
large model, and will be retained in the backward elimination phase.

<h3>Value</h3>

The output is a representation of the best model, as measured by the
penalized LOD score (see Details), among all models visited.
This is QTL object (of class <code>"qtl"</code>, as produced by
<code><a href="makeqtl.html">makeqtl</a></code>), with attributes <code>"formula"</code>,
indicating the model formula, and <code>"pLOD"</code> indicating the
penalized LOD score.

If <code>keeplodprofile=TRUE</code>, LOD profiles from the last pass through
the refinement algorithm are retained as an attribute,
<code>"lodprofile"</code>, to the object. These may be plotted with
<code><a href="plotLodProfile.html">plotLodProfile</a></code>.

If <code>keeptrace=TRUE</code>, the output will contain an attribute
<code>"trace"</code> containing information on the best model at each step
of forward and backward elimination. This is a list of objects of
class <code>"compactqtl"</code>, which is similar to a QTL object (as
produced by <code><a href="makeqtl.html">makeqtl</a></code>) but containing just
a vector of chromosome IDs and positions for the QTL. Each will also
have attributes <code>"formula"</code> (containing the model formula) and
<code>"pLOD"</code> (containing the penalized LOD score.

<h3>Methods</h3>

<code>imp</code>: multiple imputation is used, as described by Sen
and Churchill (2001).

<code>hk</code>: Haley-Knott regression is used (regression of the
phenotypes on the multipoint QTL genotype probabilities), as described
by Haley and Knott (1992).

<h3>Author(s)</h3>

Karl W Broman, <a href="mailto:broman@wisc.edu">broman@wisc.edu</a>

<h3>References</h3>

Manichaikul, A., Moon, J. Y., Sen, &#346;, Yandell, B. S. and Broman,
K. W. (2009) A model selection approach for the identification of
quantitative trait loci in experimental crosses, allowing
epistasis. Genetics, 181, 1077&ndash;1086.

Broman, K. W. and Speed, T. P. (2002) A model selection approach for the
identification of quantitative trait loci in experimental crosses (with
discussion). J Roy Stat Soc B 64, 641&ndash;656, 731&ndash;775.

Haley, C. S. and Knott, S. A. (1992) A simple regression method for mapping
quantitative trait loci in line crosses using flanking markers.
Heredity 69, 315&ndash;324.

Sen, &#346;. and Churchill, G. A. (2001) A statistical framework for quantitative
trait mapping. Genetics 159, 371&ndash;387.

Zeng, Z.-B., Kao, C.-H. and Basten, C. J. (1999) Estimating the
genetic architecture of quantitative traits. Genetical
Research, 74, 279&ndash;289.

<code><a href="calc.penalties.html">calc.penalties</a></code>,
<code><a href="plotModel.html">plotModel</a></code>, <code><a href="makeqtl.html">makeqtl</a></code>,
<code><a href="fitqtl.html">fitqtl</a></code>, <code><a href="refineqtl.html">refineqtl</a></code>,
<code><a href="addqtl.html">addqtl</a></code>, <code><a href="addpair.html">addpair</a></code>

<h3>Examples</h3>

<pre>
data(fake.bc)

## Not run: fake.bc &lt;- calc.genoprob(fake.bc, step=2.5)

outsw &lt;- stepwiseqtl(fake.bc, max.qtl=3, method="hk", keeptrace=TRUE)

# best model
outsw
plotModel(outsw)

# path through model space
thetrace &lt;- attr(outsw, "trace")

# plot of these
par(mfrow=c(3,3))
for(i in seq(along=thetrace))
 plotModel(thetrace[[i]], main=paste("pLOD =",round(attr(thetrace[[i]],"pLOD"), 2)))
</pre>

<hr /><div style="text-align: center;">[Package qtl version 1.66 <a href="00Index.html">Index</a>]</div>
</body></html>