polygenic_hglm: Estimation of polygenic model
In GenABEL: genome-wide SNP association analysis

Description Usage Arguments Details Author(s) References See Also Examples

Estimation of polygenic model using a hierarchical generalized linear model (HGLM; Lee and Nelder 1996. hglm package; Ronnegard et al. 2010).

1 2	polygenic_hglm(formula, kinship.matrix, data, family = gaussian(), conv = 1e-06, maxit = 100, ...)

`formula`	Formula describing fixed effects to be used in analysis, e.g. y ~ a + b means that outcome (y) depends on two covariates, a and b. If no covariates used in analysis, skip the right-hand 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 `gwaa.data-class` or a data frame with outcome and covariates.
`family`	a description of the error distribution and link function to be used in the mean part of the model (see `family` for details of family functions).
`conv`	'conv' parameter of `hglm` (stricter than default, for great precision, use 1e-8).
`maxit`	'maxit' parameter of `hglm` (stricter than default).
`...`	other parameters passed to `hglm` call

The algorithm gives extended quasi-likelihood (EQL) estimates for restricted maximum likelihood (REML) (Ronnegard et al. 2010). Implemented is the inter-connected GLM interpretation of HGLM (Lee and Nelder 2001). Testing on fixed and random effects estimates are directly done using inter-connected GLMs. Testing on dispersion parameters (variance components) can be done by extracting the profile likelihood function value of REML.

Xia Shen, Yurii Aulchenko

Ronnegard, L, Shen, X and Alam, M (2010). hglm: A Package For Fitting Hierarchical Generalized Linear Models. The R Journal, 2(2), 20-28.

Lee, Y and Nelder JA (2001). Hierarchical generalised linear models: A synthesis of generalised linear models, random-effect models and structured dispersions. Biometrika, 88(4), 987-1006.

Lee, Y and Nelder JA (1996). Hierarchical Generalized Linear Models. Journal of the Royal Statistical Society. Series B, 58(4), 619-678.

polygenic, mmscore, grammar

require(GenABEL.data)
data(ge03d2ex.clean)
set.seed(1)
df <- ge03d2ex.clean[,sample(autosomal(ge03d2ex.clean),1000)]
gkin <- ibs(df,w="freq")

# ----- for quantitative traits
h2ht <- polygenic_hglm(height ~ sex + age, kin = gkin, df)
# ----- estimate of heritability
h2ht$esth2
# ----- other parameters
h2ht$h2an

# ----- test the significance of fixed effects
# ----- (e.g. quantitative trait)
h2ht <- polygenic_hglm(height ~ sex + age, kin = gkin, df)
# ----- summary with standard errors and p-values
summary(h2ht$hglm)
# ----- output for the fixed effects part
# ...
# MEAN MODEL
# Summary of the fixed effects estimates
#              Estimate Std. Error t value Pr(>|t|)
# (Intercept) 169.53525    2.57624  65.807  < 2e-16 ***
# sex          14.10694    1.41163   9.993  4.8e-15 ***
# age          -0.15332    0.05208  -2.944  0.00441 **
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Note: P-values are based on 69 degrees of freedom
# ...
# ----- extract fixed effects estimates and standard errors
h2ht$h2an

# ----- test the significance of polygenic effect
nullht <- lm(height ~ sex + age, df)
l1 <- h2ht$ProfLogLik
l0 <- as.numeric(logLik(nullht))
# the likelihood ratio test (LRT) statistic
(S <- -2*(l0 - l1))
# 5% right-tailed significance cutoff
# for 50:50 mixture distribution between point mass 0 and \chi^{2}(1)
# Self, SG and Liang KY (1987) Journal of the American Statistical Association.
qchisq(((1 - .05) - .50)/.50, 1)
# p-value
pchisq(S, 1, lower.tail = FALSE)/2

# ----- for binary traits
h2dm <- polygenic_hglm(dm2 ~ sex + age, kin = gkin, df, family = binomial(link = 'logit'))
# ----- estimated parameters
h2dm$h2an

Loading required package: MASS
Loading required package: GenABEL.data
Loading required package: hglm
Loading required package: Matrix
Loading required package: hglm.data
Loading required package: sp

hglm: Hierarchical Generalized Linear Models
Version 2.2-1 (2019-04-04) installed
Authors: Moudud Alam, Lars Ronnegard, Xia Shen
Maintainer: Xia Shen <xia.shen@ki.se>

Use citation("hglm") to know how to cite our work.

Discussion: https://r-forge.r-project.org/forum/?group_id=558
BugReports: https://r-forge.r-project.org/tracker/?group_id=558
VideoTutorials: http://www.youtube.com/playlist?list=PLn1OmZECD-n15vnYzvJDy5GxjNpVV5Jr8
[1] 0.3963493
$estimate
(Intercept)         sex         age          h2    totalVar 
168.9831028  14.2472113  -0.1430366   0.3963493  56.7957659 

$se
(Intercept)         sex         age          h2    totalVar 
 2.61089707  1.43324414  0.05286632          NA          NA 

Call: 
hglm.default(X = desmat, y = y, Z = L, family = family, conv = conv, 
    maxit = maxit)

----------
MEAN MODEL
----------

Summary of the fixed effects estimates:

             Estimate Std. Error t-value Pr(>|t|)    
(Intercept) 168.98310    2.61090  64.722  < 2e-16 ***
sex          14.24721    1.43324   9.941  4.4e-15 ***
age          -0.14304    0.05287  -2.706  0.00853 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 71 degrees of freedom

Summary of the random effects estimates:

    Estimate Std. Error
Z.1  -2.3095     2.8234
Z.2   2.3019     2.8823
Z.3  -4.9848     2.9131
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).

----------------
DISPERSION MODEL
----------------

NOTE: h-likelihood estimates through EQL can be biased.

Dispersion parameter for the mean model:
[1] 34.2848

Model estimates for the dispersion term:

Link = log 

Effects:
  Estimate Std. Error 
    3.5347     0.1676 

Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).

Dispersion parameter for the random effects:
[1] 22.51

Dispersion model for the random effects:

Link = log

Effects:
.|Random1 
  Estimate Std. Error 
    3.1140     0.2213 

Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).

EQL estimation converged in 2 iterations.
$estimate
(Intercept)         sex         age          h2    totalVar 
168.9831028  14.2472113  -0.1430366   0.3963493  56.7957659 

$se
(Intercept)         sex         age          h2    totalVar 
 2.61089707  1.43324414  0.05286632          NA          NA 

[1] 207.3487
[1] 2.705543
[1] 2.601859e-47
$estimate
 (Intercept)          sex          age           h2     totalVar 
-0.683884516  0.567843025  0.017869168  0.003164765  1.328178296 

$se
(Intercept)         sex         age          h2    totalVar 
 0.86099867  0.45838409  0.01766143          NA          NA