rGLS: GWAS for Studies having Repeated Measurements on Related...
In RepeatABEL: GWAS for Multiple Observations on Related Individuals

Description Usage Arguments Details Author(s) Examples

It is used to perform genome-wide association studies on individuals that are both related and have repeated measurements. The function computes score statistic based p-values for a linear mixed model including random polygenic effects and a random effect for repeated measurements. A p-value is computed for each marker and the null hypothesis tested is a zero additive marker effect.

1 2	rGLS(formula.FixedEffects = y ~ 1, genabel.data, phenotype.data, id.name = "id", GRM = NULL, V = NULL, memory = 1e+08)

`formula.FixedEffects`	Formula including the response variable and cofactors as fixed effects.
`genabel.data`	An GenABEL object including marker information. This object has one observtion per individuals.
`phenotype.data`	A data frame including the repeated observations and IDs.
`id.name`	The column name of the IDs in phen.data
`GRM`	An optional genetic relationship matrix (GRM) can be included as input. Otherwise the GRM is computed within the function.
`V`	An optional (co)variance matrix can be included as input. Otherwise it is computed using the hglm function.
`memory`	Used to optimize computations. The maximum number of elements in a matrix that can be stored efficiently.

A generalized squares (GLS) is fitted for each marker given a (co)variance matrix V. The computations are made fast by transforming the GLS to an ordinary least-squares (OLS) problem using an eigen-decomposition of V. The OLS are computed using QR-factorization. If V is not specified then a model including random polygenic effects and permanent environmental effects is fitted (using the hglm package) to compute V. A GenABEL object (scan.gwaa class) is returned (including also the hglm results). Let e.g. GWAS1 be an object returned by the rGLS function. Then a Manhattan plot can be produced by calling plot(GWAS1) and the top SNPs using summary(GWAS1). Both of these functions are generic GenABEL functions.
The results from the fitted linear mixed model without any SNP effect included are produced by calling summary(GWAS1@call$hglm).

Lars Ronnegard

data(Phen.Data) #Phenotype data with repeated observations
 data(gen.data) #GenABEL object including IDs and marker genotypes
 GWAS1 <- rGLS(y ~ age + sex, genabel.data = gen.data, phenotype.data = Phen.Data)
 plot(GWAS1, main="")
 summary(GWAS1)
 #Summary for variance component estimation without SNP effects
 summary(GWAS1@call$hglm)

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

hglm: Hierarchical Generalized Linear Models
Version 2.1-1 (2015-08-28) 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
Loading required package: GenABEL
Loading required package: GenABEL.data
[1] "GRM ready"
[1] "Variance component estimation ready"
[1] "Rotation matrix ready"
[1] "Rotate LMM started"
[1] "Rotate LMM ready"
Summary for top 10 results, sorted by P1df 
          Chromosome Position Strand A1 A2       effB   se_effB chi2.1df
rs120315           1  3725352      +  T  C  2.1696316 0.4258979       NA
rs9670687          1  4779800      -  C  A -2.0412127 0.5351582       NA
rs891586           2  8024318      +  A  G -1.6669267 0.4478499       NA
rs1352451          2  8022651      -  A  C  1.6169429 0.4564078       NA
rs1252282          2  8167984      -  A  T -0.9982089 0.2825403       NA
rs7633966          1  3721952      +  C  T  1.8586037 0.5325691       NA
rs9922492          1  3729983      +  C  T -1.8586037 0.5325691       NA
rs4378234          1  4800348      -  G  T -1.9012073 0.5509388       NA
rs7499832          3 10242953      -  G  C -1.0897248 0.3182487       NA
rs6561272          1  3715538      -  T  G -1.9073305 0.5726170       NA
                  P1df Pc1df
rs120315  3.501204e-07    NA
rs9670687 1.366120e-04    NA
rs891586  1.975997e-04    NA
rs1352451 3.959643e-04    NA
rs1252282 4.109058e-04    NA
rs7633966 4.832322e-04    NA
rs9922492 4.832322e-04    NA
rs4378234 5.588239e-04    NA
rs7499832 6.167705e-04    NA
rs6561272 8.656541e-04    NA
Call: 
hglm.default(X = X, y = y, Z = Z, maxit = 200, RandC = c(ncol(Z.GRM), 
    ncol(Z.indx)))

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

Summary of the fixed effects estimates:

             Estimate Std. Error t-value Pr(>|t|)    
(Intercept) 102.05662    0.78674 129.721  < 2e-16 ***
age           0.06504    0.01451   4.483 1.06e-05 ***
sex           0.85160    0.38016   2.240   0.0258 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 290 degrees of freedom

Summary of the random effects estimates:

    Estimate Std. Error
Z.1  -2.1532     0.7458
Z.2   0.0267     0.7619
Z.3  -0.2035     0.7559
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).

Summary of the random effects estimates:

      Estimate Std. Error
Z.101   0.3941     0.9913
Z.102   0.0778     0.8626
Z.103   1.5205     0.8765
...
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] 3.090352

Model estimates for the dispersion term:

Link = log 

Effects:
  Estimate Std. Error 
    1.1283     0.0830 

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

Dispersion parameter for the random effects:
[1] 1.376 1.354

Dispersion model for the random effects:

Link = log

Effects:
.|Random1 
  Estimate Std. Error 
    0.3194     0.2495 

.|Random2 
  Estimate Std. Error 
    0.3032     0.2407 

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

EQL estimation converged in 2 iterations.