Description Usage Arguments Details Value Note Author(s) References See Also Examples
This function maximises the likelihood of the data under polygenic model with covariates an reports twice negative maximum likelihood estimates and the inverse of the variancecovariance matrix at the point of ML.
1 2 3 4 5 6 7  polygenic(formula, kinship.matrix, data, fixh2,
starth2 = 0.3, trait.type = "gaussian",
opt.method = "nlm", scaleh2 = 1, quiet = FALSE,
steptol = 1e08, gradtol = 1e08, optimbou = 8,
fglschecks = TRUE, maxnfgls = 8, maxdiffgls = 1e04,
patchBasedOnFGLS = TRUE, llfun = "polylik_eigen",
eigenOfRel, ...)

formula 
Formula describing fixed effects to be used in the analysis, e.g. y ~ a + b means that outcome (y) depends on two covariates, a and b. If no covariates used in the analysis, skip the righthand side of the equation. 
kinship.matrix 
Kinship matrix, as provided by e.g. ibs(,weight="freq"), or estimated outside of GenABEL from pedigree data. 
data 
An (optional) object of

fixh2 
Optional value of heritability to be used, instead of maximisation. The uses of this option are twofold: (a) testing significance of heritability and (b) using a priori known heritability to derive the rest of MLEs and var.cov. matrix. 
starth2 
Starting value for h2 estimate 
trait.type 
"gaussian" or "binomial" 
opt.method 
"nlm" or "optim". These two use
different optimisation functions. We suggest using the
default 
scaleh2 
Only relevant when "nlm" optimisation function is used. "scaleh2" is the heritability scaling parameter, regulating how "big" are parameter changes in h2 with respect to changes in other parameters. As other parameters are estimated from previous regression, these are expected to change little from the initial estimate. The default value of 1000 proved to work rather well under a range of conditions. 
quiet 
If FALSE (default), details of optimisation process are reported 
steptol 
steptal parameter of "nlm" 
gradtol 
gradtol parameter of "nlm" 
optimbou 
fixed effects boundary scale parameter for 'optim' 
fglschecks 
additional check for convergence on/off (convergence between estimates obtained and that from FGLS) 
maxnfgls 
number of fgls checks to perform 
maxdiffgls 
max difference allowed in fgls checks 
patchBasedOnFGLS 
if FGLS checks not passed, 'patch' fixed effect estimates based on FGLS expectation 
llfun 
function to compute likelihood (default 'polylik_eigen', also available – but not recommended – 'polylik') 
eigenOfRel 
results of eigen(relationship matrix = 2*kinship.matrix). Passing this can decrease computational time substantially if multiple traits are analysed using the same kinship matrix. This option will not work if any NA's are found in the trait and/or covariates and if the dimensions of the 'eigen'object, trait, covariates, kinship do not match. 
... 
Optional arguments to be passed to

One of the major uses of this function is to estimate
residuals of the trait and the inverse of the
variancecovariance matrix for further use in analysis
with mmscore
and grammar
.
Also, it can be used for a variant of GRAMMAR analysis,
which allows for permutations for GW significance by use
of environmental residuals as an analysis trait with
qtscore
.
"Environmental residuals" (not to be mistaken with just "residuals") are the residual where both the effect of covariates AND the estimated polygenic effect (breeding values) are factored out. This thus provides an estimate of the trait value contributed by environment (or, turning this other way around, the part of the trait not explained by covariates and by the polygene). Polygenic residuals are estimated as
σ^2 V^{1} (Y  (\hat{μ} + \hat{β} C_1 + ...))
where sigma^2 is the residual variance, V^{1} is the InvSigma (inverse of the varcov matrix at the maximum of polygenic model) and (Y  (\hat{μ} + \hat{β} C_1 + ...)) is the trait values adjusted for covariates (also at at the maximum of polygenic model likelihood).
It can also be used for heritability analysis. If you want to test significance of heritability, estimate the model and write down the function minimum reported at the "h2an" element of the output (this is twice the negative MaxLikelihood). Then do a next round of estimation, but set fixh2=0. The difference between your function minima gives a test distributed as chisquared with 1 d.f.
The way to compute the likelihood is partly based on the paper of Thompson (see refs), namely instead of taking the inverse of the varcov matrix every time, eigenvectors of the inverse of G (taken only once) are used.
A list with values
h2an 
A list supplied by the

esth2 
Estimate (or fixed value) of heritability 
residualY 
Residuals from analysis, based on covariate effects only; NOTE: these are NOT grammar "environmental residuals"! 
pgresidualY 
Environmental residuals from analysis, based on covariate effects and predicted breeding value. 
grresidualY 
GRAMMAR+ transformed trait residuals 
grammarGamma 
list with GRAMMARgamma correction factors 
InvSigma 
Inverse of the
variancecovariance matrix, computed at the MLEs – these
are used in 
call 
The details of call 
measuredIDs 
Logical values for IDs who were used in analysis (traits and all covariates measured) == TRUE 
convFGLS 
was convergence achieved according to FGLS criterionE 
Presence of twins may complicate your analysis. Check the
kinship matrix for singularities, or rather use
check.marker
for identification of twin
samples. Take special care in interpretation.
If a trait (no covariates) is used, make sure that the order of IDs in the kinship.matrix is exactly the same as in the outcome
Please note that there is alternative to 'polygenic',
polygenic_hglm
, which is faster than
polygenic() with the llfun='polylik' option, but slightly
slower than the default polygenic().
Yurii Aulchenko, Gulnara Svischeva
Thompson EA, Shaw RG (1990) Pedigree analysis for quantitative traits: variance components without matrix inversion. Biometrics 46, 399413.
for original GRAMMAR
Aulchenko YS, de Koning DJ, Haley C. Genomewide rapid association using mixed model and regression: a fast and simple method for genomewide pedigreebased quantitative trait loci association analysis. Genetics. 2007 177(1):57785.
for GRAMMARGC
Amin N, van Duijn CM, Aulchenko YS. A genomic background based method for association analysis in related individuals. PLoS ONE. 2007 Dec 5;2(12):e1274.
for GRAMMARGamma
Svischeva G, Axenovich TI, Belonogova NM, van Duijn CM, Aulchenko YS. Rapid variance componentsbased method for wholegenome association analysis. Nature Genetics. 2012 44:11661170. doi:10.1038/ng.2410
for GRAMMAR+ transformation
Belonogova NM, Svishcheva GR, van Duijn CM, Aulchenko YS, Axenovich TI (2013) RegionBased Association Analysis of Human Quantitative Traits in Related Individuals. PLoS ONE 8(6): e65395. doi:10.1371/journal.pone.0065395
polygenic_hglm
, mmscore
,
grammar
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18  # note that procedure runs on CLEAN data
require(GenABEL.data)
data(ge03d2ex.clean)
gkin < ibs(ge03d2ex.clean,w="freq")
h2ht < polygenic(height ~ sex + age, kin=gkin, ge03d2ex.clean)
# estimate of heritability
h2ht$esth2
# other parameters
h2ht$h2an
# the minimum twice negative loglikelihood
h2ht$h2an$minimum
# twice maximum loglikelihood
h2ht$h2an$minimum
# for binary trait (experimental)
h2dm < polygenic(dm2 ~ sex + age, kin=gkin, ge03d2ex.clean, trait="binomial")
# estimated parameters
h2dm$h2an

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.