*! version 1.0.0 30apr2008
help haplologit

Title

haplologit -- Haplotype-effects logistic regression for case-control data

Syntax

haplologit depvar [indepvars] [if] [in] ,
snpvars:(varlist) [options]

Model
* snpvars(varlist) specify SNP variables

inheritance(haplologit##inheritance:inhmode) specify mode of
inheritance; default is inheritance(additive)

riskhap(haplologit##riskhap_spec:riskhap_spec) specify a
single risk haplotype

riskhap#(riskhap_spec) specify #th risk haplotype

hfthreshold(#) retain observations with initial haplotype
frequencies exceeding hfthreshold(); default is
max(2/N,0.001)

constraints(numlist) apply specified linear constraints on environmental factors indepvars

collinear keep collinear variables

noconstant suppress constant term

Reporting
level(#) set confidence level; default is
level(95)

or report odds ratios

happrefix(string) use string as a prefix when labeling haplotypes in the output; default is happrefix(hap_)

alldots show all iterations (except ml) as dots

nocoef suppress coefficients table

nofreq suppress haplotype-frequency table

noheader suppress output header

EM options
emsample(controls|cases|all) obtain initial haplotype frequencies from the specified sample; default is emsample(controls)

emiterate(#) number of EM iterations; default is 500

emtolerance(#) EM convergence tolerance; default is 1e-6

eminit(matname) specify matrix containing starting values of haplotype frequencies for EM estimation

sort sort haplotypes by frequencies in the EM haplotype-frequency table; default is to sort by haplotypes (in a binary order)

emlog show EM iteration log

emdots show EM iterations as dots

noemshow suppress output from EM estimation

noemtable suppress EM haplotype-frequency table

Max options
maximize_options control the maximization process; seldom
used

* snpvars(varlist) is required.

inhmode

additive additive mode of inheritance; the default

dominant dominant mode of inheritance

recessive recessive mode of inheritance

where riskhap_spec is

riskhap_str|# [, riskhap_suboptions]

riskhap_str specifies the binary representation of a risk haplotype
enclosed in quotes or # specifies a risk haplotype index (position
of a risk haplotype in the ordered sequence of 2^M possible haplotypes at
M SNP sites).

riskhap_suboptions

interaction(varlist) specify interaction variables

noconstant suppress constant term; seldom used

Description

haplologit estimates haplotype effects and haplotype-environment
interactions from case-control genetic (SNP-based) data for one of three types
of genetic (haplotype risk) models: additive, dominant, or recessive. It fits
haplotype-effects logistic regression using the retrospective
profile-likelihood method in a special case of a rare disease and a single
candidate gene in Hardy-Weinberg equilibrium, under the assumption of
gene-environment independence. haplologit handles phased, unphased, and
missing genotypes and allows specifying multiple risk haplotypes.

Options

Model

snpvars(varlist) is required; it specifies SNP variables (variables
recording subjects' SNP genotypes). The SNP variables must contain values of
0, 1, 2, or missing (.). A missing value (.) indicates missing information at
a SNP site, other values represent the number of copies of a mutant (minor)
allele at a SNP site in a subject's pair of homologous chromosomes.

inheritance(inhmode) specifies a mode of inheritance (a genetic model).
The default is the additive risk model in which having two copies of a
risk haplotype in a pair of homologous chromosomes results in a two-fold
effect of the risk haplotype on a disease. The dominant risk model
assumes that having one or two copies of a risk haplotype has the same effect
on a disease. The recessive model assumes that only having two copies
of a risk haplotype has an effect on a disease.

riskhap(riskhap_spec) requests to include effects of the specified risk
haplotype in a regression model. riskhap() is a synonym for
riskhap1().

riskhap#(riskhap_spec) requests to include effects
of the #th risk haplotype in a regression model. If interaction(varlist) is specified in riskhap#() the
respective interaction effects of the risk haplotype with the covariates
specified in varlist are included (in addition to main effect) in the
regression model. If noconstant is used, the main haplotype effect is
omitted and only haplotype-environment interaction effects are included in the
model (seldom used).

interaction(varlist) specifies variables to be interacted with the
specified risk haplotype.

noconstant requests that the constant term (the main effect of a risk
haplotype) is not included in the model (seldom used). This option requires
interaction().

hfthreshold(#) specifies to retain in the computation only subjects'
diplotypes with initial frequencies of constituent haplotypes exceeding
#. The default is max(2/N,0.001) where N is the total number of
cases and controls.

constraints(numlist),
collinear; see [R] estimation options.
constraints() may only be used to define linear constraints on
environmental covariates indepvars.

noconstant suppresses the constant (intercept) term; seldom used.

Reporting

level(#); see [R] estimation options.

or reports the estimated coefficients transformed to odds ratios, i.e.,
exp(b) rather than b. Standard errors and confidence intervals are similarly
transformed. This option affects how results are displayed, not how they are
estimated. or may be specified at estimation or when replaying
previously estimated results.

happrefix(string) uses the specified string as a prefix when
labeling haplotypes in the output (except for the EM output). The default
prefix is hap_.

alldots specifies that iterations from all (possibly time-consuming)
computations be shown as dots except for the ml iterations. alldots
implies emdots.

nocoef specifies that the coefficient table not be displayed.

nofreq specifies that the haplotype-frequency table not be displayed.

noheader suppresses the output header, either at estimation or upon
replay.

EM options

emsample(controls|cases|all) requests that the initial haplotype
frequencies are estimated from the control sample, case sample, or combined
case-control sample. The default is to use the control sample.

emiterate(#) specifies the number of EM iterations to perform. The
default is 500.

emtolerance(#) specifies the convergence tolerance for the EM algorithm.
The default is 1e-6. The EM algorithm terminates when the maximum relative
change in estimated haplotype frequencies between successive iterations is
less than #.

eminit(matname) specifies the 1xL matrix matname containing
starting values of haplotype frequencies for EM estimation. If M is the
number of SNP loci (SNP variables), then L=2^M-1 . By default, all
haplotypes are assumed to be equally likely, that is all haplotype frequencies
are set to 1/2^M.

sort requests that haplotypes are displayed in descending order of
frequencies in the EM haplotype-frequency table. By default, haplotypes are
displayed according to their binary ordering.

emlog specifies that the EM iteration log be shown. The EM iteration
log is, by default, not displayed.

emdots specifies that the EM iterations be shown as dots. This option
can be convenient when the EM algorithm requires many iterations to converge.

noemshow suppresses the output from the EM estimation.

noemtable suppresses the EM haplotype-frequency table.

Max options

maximize_options:
difficult,
iterate(#),
[no]log,
trace,
hessian,
gradient,
showstep,
tolerance(#),
ltolerance(#),
gtolerance(#),
nrtolerance(#),
nonrtolerance,
shownrtolerance;
see maximize.

By default convergence is declared when the nrtolerance() criterion and
either of the tolerance() or ltolerance() criterion has been met.
If nonrtolerance is specified, then convergence is declared when either
of the tolerance() or ltolerance() criterion has been met.

If gtolerance() is specified, then the gtolerance() criterion must
be met in addition to any other required criteria for convergence to be
declared. See Tolerance options for more
information.

Remarks

haplologit is designed for the analysis of haplotype-disease
associations in case-control data. It utilizes the retrospective profile
likelihood approach of Spinka et al. (2005) and Lin and Zeng (2006) to
estimate haplotype and haplotype-environment effects. This approach is more
efficient for the analysis of case-control data than the conventional
prospective methods. The gain in efficiency is obtained from utilizing the
available information about genetic distribution (Hardy-Weinberg equilibrium,
independence of environmental factors) in the construction of the
retrospective likelihood of the data. If only haplotype main effects are
specified (no environmental effects), no profiling of the likelihood is needed
and the retrospective likelihood approach of Epstein and Satten (2003) is
used.

Subjects' genetic information is described by a sequence of single nucleotide
polymorphisms (SNPs), known to be genetic markers, along the gene of interest.
Specifically, this information consists of pairs of subjects' SNP
haplotypes where a haplotype is a set of nearby SNPs on the same
chromosome. In practice, the true haplotype pairs (diplotypes) are not
directly observed. Instead, a SNP genotype (a combination of two
homologous SNP haplotypes) is observed. The observed SNP genotype data is
supplied to haplologit as subjects' genetic information via required
option snpvars(). The genotype data must be recorded in the so-called
SNP variables (one variable for each SNP locus), containing values of 0, 1, 2,
and missing (.) only. A missing value (.) indicates missing information at a
SNP site, other values record the number of copies of a mutant (minor) allele
at a SNP site in a subject's pair of homologous chromosomes. Data must be in
the wide form, that is a single observation per subject.

The effect of a single gene in HWE on the disease is considered. A risk
haplotype (or causal haplotype) is a target haplotype whose effect on a
disease is of interest. The effects of risk haplotypes can be modeled
according to one of three genetic models, specified in option
inheritance() as the mode of inheritance: additive (the default),
dominant, or recessive. Genetic covariates are viewed as functions of
subjects' SNP genotype data and risk haplotypes. Specifically, they depend on
the number of copies of a risk haplotype present in the subject's diplotype.
Their functional forms are determined by the selected genetic model. For
example, under the additive risk model, having two copies of a risk haplotype
in a subject's diplotype doubles the effect of this haplotype on a disease
compared to having only one copy. In contrast, under the dominant risk model
having one or two copies has the same effect on the disease. Under the
recessive model, only having two copies of a risk haplotype has an effect on
the disease. haplologit uses genetic covariates indirectly in the
computation via the supplied information about SNP genotypes (option
snpvars()), the genetic model (option inheritance()), and the risk
haplotypes.

A risk haplotype may be specified as a string of a sequence of zeros and ones
(binary representation) or as a haplotype index (position of a risk haplotype
in the ordered sequence of 2^M possible haplotypes at M SNP sites). Risk
haplotypes are specified in options riskhap1(), riskhap2(), and so
on. By default, if no risk haplotypes are specified, haplologit uses
the most frequent haplotype estimated from the control sample (or the sample
specified in emsample()) as the risk haplotype. Environmental
covariates may be specified as indepvars following the dependent variable
depvar in the above syntax. The interaction effects of environmental
factors and risk haplotypes may be included by using riskhap#()'s
suboption interaction().

The distributional assumptions on genetic covariates are Hardy-Weinberg
equilibrium and independence with environmental covariates. Environmental
covariates can be both continuous and discrete and their distribution is left
unspecified.

haplologit's estimation process consists of three stages: (1) data
management, (2) initial estimation of haplotype frequencies, and (3) the
estimation of haplotype and optionally environmental effects. During the data
management stage haplologit performs data manipulations necessary for
handling unphased and missing SNP genotypes in the computation. At the second
stage the initial haplotype frequencies are estimated from the sample
specified in emsample() by using the EM algorithm. Only the haplotypes
with the estimated initial frequencies exceeding a default threshold (or an
alternate threshold specified in hfthreshold()) are retained for further
estimation. This is necessary for numerical stability of the algorithm. At
the third stage, the coefficients for environmental covariates, risk
haplotypes, and their interactions are estimated simultaneously with the
haplotype frequencies by Newton-Raphson. The command displays information
from and optionally progress at each of the three steps.

The execution time of haplologit increases significantly with an
increased number of SNP loci and haplotype and environmental effects. The
presence of many subjects with missing and/or unphased genotypes increases the
execution time as well.

For more details and methodology see Marchenko et al. (2008).

Examples

The dataset cc_snp.dta contains fictional data on 90 individuals (45
cases and 45 controls) genotyped at 2 SNP sites. The environmental covariates
are age and gender. The disease indicator variable is
status. Subjects' SNP genotypes are recorded in variables snp1
and snp2.

1. Additive main effects of the most frequent haplotype. Use
age as the environmental covariate, status as the dependent
variable, and 2 SNP variables snp1 and snp2.

. use cc_snp.dta
. haplologit status age, snpvars(snp1 snp2)

From the output, the estimated frequencies of the four haplotypes "00", "01",
"10", and "11" from the control sample are 0.025, 0.386, 0.497, and 0.092,
respectively. All haplotype frequencies exceed the default threshold of 0.022
and thus all haplotypes are used in the estimation. The regression includes
main effects of age and the default risk haplotype "10" (the most
frequent haplotype in the control sample). Although the obtained results are
not statistically significant, the estimated effect (log odds ratio) of
age, 0.019, suggests an increase in the risk of the disease with age,
whereas the presence of a risk haplotype "10" (coefficient for hap_10 is
-0.552) decreases this risk.

The left part of the header of the coefficient table reports general
information about the fitted model. This includes the type of the genetic
model used in the computation (the default is additive), the assumed genotype
distribution of the data (Hardy-Weinberg equilibrium), and variable names
containing subjects' SNP genotypes (snp1 and snp2).

The right part of the header provides general information about the data.
Number of obs displays the number of observations (subjects) used in the
computation. Subjects are determined by the rows of the dataset. In this
example all observations in the dataset (90) are used in the computation.
haplologit also reports the number of subjects with phased, unphased,
and incomplete (missing in at least on SNP variable) genotypes. Among 90
subjects used in the computation, 49 have phased genotypes, 34 have unphased
genotypes, and 7 have incomplete genotypes.

Under the rare-disease assumption the intercept parameter of the logistic
model (b0) is unidentifiable (see, for example, Spinka et al. (2005)); it is
confounded with the unknown marginal probability of a disease for a population
Pr(D=1)=p1. haplologit reports the estimate of the "retrospective"
constant term in the output, _cons. If p1 is known, the estimate of the
intercept b0 may be retrieved from the formula
_cons=b0+ln(N1/N0)-ln(p1/p0) where N1 is the number of cases, N0
is the number of controls, and p0=1-p1 .

2a. Additive main effects of 01. Here we investigate the main
effect of haplotype 01. We specify the risk haplotype "01" in option
riskhap().

. haplologit status age, snpvars(snp1 snp2) riskhap("01")

Alternatively, we can specify the haplotype index 2 instead of "01" in
riskhap() and obtain the same results.

. haplologit status age, snpvars(snp1 snp2) riskhap(2)

2b. Reporting odds ratios. We can obtain the above results
displayed as odds ratios by using option or on replay.

. haplologit, or

3. Additive main and interaction effects of 01 with age. In
this example, following the syntax of
riskhap_spec, we specify inter(age)
in option riskhap() to include the interaction effect of risk haplotype
01 with age.

. haplologit status age, snpvars(snp1 snp2) riskhap("01", inter(age))

4. Dominant main and interaction effects of 01 with age.
Here we change the mode of inheritance from the default additive to dominant
by using option inheritance(dominant). In fact, we use its abbreviated
version inher(d).

. haplologit status age, snpvars(snp1 snp2) riskhap("01", inter(age)) inher(d)

5. Joint additive main effects of 01 and 10. We specify risk
haplotypes 01 and 10 in options riskhap1() and riskhap2(),
respectively.

. haplologit status age, snpvars(snp1 snp2) riskhap1("01") riskhap2("10")

6. Joint additive main and interaction effects of 01 and 10 with
covariates age and gender. We add variable gender
to the list of independent variables and specify inter(age gender) in
options riskhap1() and riskhap2() to include their interaction
effects with the risk haplotypes.

. haplologit status age gender, snp(snp1 snp2) riskhap1("01", inter(age gender)) riskhap2("10", inter(age gender))

Saved results

haplologit saves the following in e():

Scalars

e(N) number of observations (subjects used in the computation)

e(N_phased) number of subjects with phased genotypes

e(N_unphased) number of subjects with unphased genotypes

e(N_miss) number of subjects with incomplete genotypes (missing in at least one SNP variable)

e(ll) retrospective (profile) log-likelihood

e(converged) 1 if converged, 0 otherwise

e(df_m) model degrees of freedom

e(chi2) chi-squared

e(p) significance of model test

e(em_N) number of observations at EM stage

e(em_ll) EM log-likelihood

e(cutoff) haplotype-frequency threshold

e(rc) return code

Macros

e(cmd) haplologit

e(cmdline) command as typed

e(depvar) name of dependent variable

e(snpvars) names of SNP variables

e(inheritance) mode of inheritance

e(genepop) genetic distribution

e(emsample) a sample used to obtain initial haplotype frequencies

e(happrefix) prefix used to label haplotypes in the output

Matrices

e(b) coefficient vector

e(V) variance-covariance matrix of the estimators

e(em_freq) initial haplotype frequency vector

Functions

e(sample) marks estimation sample

For other results see [R] maximize.

References

Epstein, M. P., and G. A. Satten. 2003. Inference on haplotype effects in
case-control studies using unphased genotype data. American Journal
of Human Genetics 73: 1316-1329.

Lin, D. Y., and D. Zeng. 2006. Likelihood-based inference on haplotype effects
in genetic association studies (with discussion). Journal of
the American Statistical Association 101: 89-118.

Marchenko, Y. V., Carroll, R. J., Lin, D. Y., Amos, C. I. , and R. G. Gutierrez.
2008. Semiparametric analysis of case-control genetic data in the presence of
environmental factors. Stata Journal ??: ??-??.

Spinka, C., Carroll, R. J., and Chatterjee, N. 2005. Analysis of case-control
studies of genetic and environmental factors with missing genetic information
and haplotype-phase ambiguity. Genetic Epidemiology 29: 108-127.