Overview

In this tutorial, we will guide you through the use of the indexS(), subSV(), summaryGLS(), and summaryGLSbands() functions found in the GenomicSEM R package and introduced in Ennis et al. (XXXX).

summaryGLS() uses a generalized least squares approach that allows the user to conduct a regression analysis to detect patterns of change among parameter estimates with correlated estimation errors. Used together, indexS() and subSV() subset an LDSC object to include only the values of interest for summaryGLS(). summaryGLSbands() allows the user to model one primary set of regression predictors in a scatterplot with error bands, with the option to include control variables and a quadratic term.

Each function is accompanied by examples illustrating the analyses conducted in Ennis et al. (XXXX), in which these functions are introduced.

indexS()

indexS() returns a matrix of index values that map the parameter estimates of interest from the S matrix to their shared estimation errors in the V matrix. The returned matrix can then be subset to the parameters of interest. For example, indexS(MATRIX = LDSCoutput$S)[1:12, 13] selects the indices for the relationships between a set of 12 traits and one external trait. These indices can then be supplied to subSV() to obtain the subsetted S and V matrices for subsequent GLS estimation.

Arguments

Required

  • LDSC_OBJECT or MATRIX: Supply either an LDSC object using LDSC_OBJECT, from which the S matrix will be extracted, or an S matrix directly using MATRIX.

Optional

  • R: By default, the function returns index values from the S matrix that map the parameters onto their shared estimation errors in the V matrix. The user may instead work with the directly estimated genetic correlations acquired from rgmodel() in GenomicSEM, which have offset index values, by specifying R = TRUE.

Output

A symmetric matrix of index values with the same dimensions as the supplied S matrix.

Examples

The user may obtain the index values for the parameters describing the relationships between a set of traits and one external trait by selecting either the row or column for the external trait of interest. The matrix of index values is symmetric around the diagonal. Below is an example pulling the index values for 12 parameters relating one trait to 12 others:

indexS(MATRIX = LDSCoutput$S)[1:12, 13]
##  [1]  13  46  78 109 139 168 196 223 249 274 298 321

Here is an example of acquiring the same index values for the R matrix:

indexS(MATRIX = LDSCoutput$S, R = TRUE)[1:12, 13]
##  [1]  12  44  75 105 134 162 189 215 240 264 287 309

The user may want to use summaryGLS() to model cross-group or cross-trait trends. In this case, the diagonal of the subsetted index-value matrix can be used. Below is an example in which index values for cross-group parameters are extracted. In this hypothetical scenario, the summary statistics for a set of traits in one group are found across indices 1:5 and those for a second group are found across indices 6:10:

diag(indexS(MATRIX = LDSCoutput$S)[1:5, 6:10])
## [1]   6  40  73 105 136

If the user would like to model systematic trends in heritabilities across a subset of traits in S, the diagonal of the relevant block of the matrix is specified as follows:

diag(indexS(MATRIX = LDSCoutput$S)[1:12, 1:12])
##  [1]   1  35  68 100 131 161 190 218 245 271 296 320

subSV()

subSV() subsets an S and V matrix for GLS analyses. It uses the index values from indexS() to identify the parameters of interest in S and their cross-parameter estimation errors in V, which is more structurally complex.

Arguments

Required

  • LDSC_OBJECT or SMATRIX and VMATRIX: Supply either an LDSC object, from which the relevant S and V matrices will be extracted, or supply the S and V matrices directly.
  • INDEXVALS: A vector of index values, typically produced by indexS(), that maps the parameters of interest in S onto their corresponding sampling errors and covariances in V.

Optional

  • TYPE: Specifies which matrices should be extracted from the LDSC object. Use "S" for the unstandardized S and V matrices, "S_Stand" for the standardized matrices, or "R" for the directly estimated genetic correlation matrix from rgmodel(). The default is "S".

Output

A list containing:

  • subS: The selected parameter estimates.
  • subV: The corresponding sampling covariance matrix.

Examples

Here, we illustrate the use of subSV() using the LDSC object from Ennis et al. (XXXX), which includes the genetic correlations estimated or modeled in the paper. For this example, we subset the LDSC object to conduct the primary GLS analysis testing whether there is a linear and/or quadratic trend in the genetic correlations between a suite of neuroticism GWASs across the severity spectrum and MDD. Within this LDSC object, the 12 neuroticism cutpoints are the first 12 traits, and MDD is the 13th trait.

Using indexS(), we obtain the indices in V that correspond to the parameter estimates in S. When these indices are supplied to subSV(), the function returns the subsetted S and V matrices shown below:

subSV(
  LDSC_OBJECT = LDSCoutput,
  INDEXVALS = indexS(MATRIX = LDSCoutput$S)[1:12, 13],
  TYPE = "S_Stand"
)
## $subV
##              [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
##  [1,] 0.004525147 0.003516321 0.003239704 0.002765318 0.002473485 0.002381545
##  [2,] 0.003516321 0.003941051 0.003509050 0.003067823 0.002737499 0.002648608
##  [3,] 0.003239704 0.003509050 0.003960186 0.003489059 0.003058101 0.002991574
##  [4,] 0.002765318 0.003067823 0.003489059 0.003667373 0.003234418 0.003152736
##  [5,] 0.002473485 0.002737499 0.003058101 0.003234418 0.003378778 0.003243490
##  [6,] 0.002381545 0.002648608 0.002991574 0.003152736 0.003243490 0.003680699
##  [7,] 0.002308378 0.002660188 0.002917895 0.003042383 0.003102778 0.003475608
##  [8,] 0.002256821 0.002568257 0.002836095 0.002963695 0.003008038 0.003343570
##  [9,] 0.002085878 0.002461875 0.002690090 0.002811801 0.002843278 0.003119754
## [10,] 0.001918799 0.002324614 0.002523519 0.002654755 0.002702368 0.002990591
## [11,] 0.001778394 0.002162108 0.002402823 0.002560167 0.002602974 0.002939110
## [12,] 0.001750141 0.002210388 0.002540617 0.002964549 0.002961483 0.003328889
##              [,7]        [,8]        [,9]       [,10]       [,11]       [,12]
##  [1,] 0.002308378 0.002256821 0.002085878 0.001918799 0.001778394 0.001750141
##  [2,] 0.002660188 0.002568257 0.002461875 0.002324614 0.002162108 0.002210388
##  [3,] 0.002917895 0.002836095 0.002690090 0.002523519 0.002402823 0.002540617
##  [4,] 0.003042383 0.002963695 0.002811801 0.002654755 0.002560167 0.002964549
##  [5,] 0.003102778 0.003008038 0.002843278 0.002702368 0.002602974 0.002961483
##  [6,] 0.003475608 0.003343570 0.003119754 0.002990591 0.002939110 0.003328889
##  [7,] 0.003816842 0.003628525 0.003368926 0.003235128 0.003203597 0.003520683
##  [8,] 0.003628525 0.004168432 0.003931123 0.003823730 0.003713977 0.003987714
##  [9,] 0.003368926 0.003931123 0.004497136 0.004318204 0.004064820 0.004448955
## [10,] 0.003235128 0.003823730 0.004318204 0.005299458 0.005093393 0.005616733
## [11,] 0.003203597 0.003713977 0.004064820 0.005093393 0.007174898 0.007480572
## [12,] 0.003520683 0.003987714 0.004448955 0.005616733 0.007480572 0.013784695
## 
## $subS
##  [1] 0.4853433 0.4579246 0.5138622 0.5260669 0.5347160 0.5897716 0.6418681
##  [8] 0.6077084 0.6323309 0.6816521 0.7704100 0.8971379

To extract heritabilities, specify the diagonal of the block of traits to be modeled and use the unstandardized S and V matrices, as shown below:

subSV(
  LDSC_OBJECT = LDSCoutput,
  INDEXVALS = diag(indexS(MATRIX = LDSCoutput$S)[1:12, 1:12]),
  TYPE = "S"
)
## $subV
##               [,1]         [,2]         [,3]         [,4]         [,5]
##  [1,] 2.704202e-05 1.996753e-05 1.847191e-05 1.824910e-05 1.605820e-05
##  [2,] 1.996753e-05 2.751434e-05 2.430993e-05 2.379823e-05 2.056270e-05
##  [3,] 1.847191e-05 2.430993e-05 2.861398e-05 2.751379e-05 2.402241e-05
##  [4,] 1.824910e-05 2.379823e-05 2.751379e-05 3.132162e-05 2.698468e-05
##  [5,] 1.605820e-05 2.056270e-05 2.402241e-05 2.698468e-05 2.814617e-05
##  [6,] 1.501062e-05 1.910792e-05 2.207178e-05 2.475072e-05 2.494764e-05
##  [7,] 1.675471e-05 2.085260e-05 2.382657e-05 2.639904e-05 2.580722e-05
##  [8,] 1.544710e-05 1.837677e-05 2.166737e-05 2.389226e-05 2.349730e-05
##  [9,] 1.355800e-05 1.642057e-05 1.953899e-05 2.200789e-05 2.180589e-05
## [10,] 1.399402e-05 1.714505e-05 2.012701e-05 2.225750e-05 2.190175e-05
## [11,] 1.360250e-05 1.545840e-05 1.852190e-05 2.113016e-05 2.022905e-05
## [12,] 9.452816e-06 1.058275e-05 1.511169e-05 1.742131e-05 1.796487e-05
##               [,6]         [,7]         [,8]         [,9]        [,10]
##  [1,] 1.501062e-05 1.675471e-05 1.544710e-05 1.355800e-05 1.399402e-05
##  [2,] 1.910792e-05 2.085260e-05 1.837677e-05 1.642057e-05 1.714505e-05
##  [3,] 2.207178e-05 2.382657e-05 2.166737e-05 1.953899e-05 2.012701e-05
##  [4,] 2.475072e-05 2.639904e-05 2.389226e-05 2.200789e-05 2.225750e-05
##  [5,] 2.494764e-05 2.580722e-05 2.349730e-05 2.180589e-05 2.190175e-05
##  [6,] 2.727019e-05 2.719179e-05 2.484555e-05 2.262286e-05 2.318024e-05
##  [7,] 2.719179e-05 3.290943e-05 2.958371e-05 2.654778e-05 2.625035e-05
##  [8,] 2.484555e-05 2.958371e-05 3.353115e-05 3.029157e-05 2.926736e-05
##  [9,] 2.262286e-05 2.654778e-05 3.029157e-05 3.608514e-05 3.366224e-05
## [10,] 2.318024e-05 2.625035e-05 2.926736e-05 3.366224e-05 4.622626e-05
## [11,] 2.125675e-05 2.363781e-05 2.567604e-05 2.967184e-05 4.119582e-05
## [12,] 1.918614e-05 1.780262e-05 1.889673e-05 2.241044e-05 3.273765e-05
##              [,11]        [,12]
##  [1,] 1.360250e-05 9.452816e-06
##  [2,] 1.545840e-05 1.058275e-05
##  [3,] 1.852190e-05 1.511169e-05
##  [4,] 2.113016e-05 1.742131e-05
##  [5,] 2.022905e-05 1.796487e-05
##  [6,] 2.125675e-05 1.918614e-05
##  [7,] 2.363781e-05 1.780262e-05
##  [8,] 2.567604e-05 1.889673e-05
##  [9,] 2.967184e-05 2.241044e-05
## [10,] 4.119582e-05 3.273765e-05
## [11,] 6.417260e-05 4.866309e-05
## [12,] 4.866309e-05 1.304099e-04
## 
## $subS
##  [1] 0.10084849 0.11775859 0.12433008 0.12499084 0.12644499 0.12906139
##  [7] 0.13321876 0.13043655 0.12742782 0.12637293 0.11278042 0.08507253

summaryGLS()

summaryGLS() implements a generalized least squares approach to estimating trends in data with shared estimation errors. By default, the model includes an intercept, which may be toggled off with INTERCEPT = FALSE.

The function requires a vector of dependent variables, such as a subset of an S matrix, and its associated sampling covariance matrix, V. These may be supplied either as the output of subSV() or directly. The function also requires a vector or matrix of predictor values, with the number of rows equal to the number of dependent-variable observations and the number of columns equal to the number of predictor variables. If the user wants to add higher-order polynomial or interaction terms to the GLS regression, these must be created manually as additional predictor columns.

Do not include a column of ones in PREDICTORS. Use the INTERCEPT argument to specify whether the model contains an intercept.

Arguments

Required

  • OBJECT or Y and V_Y: Supply either a list containing subS and subV, such as the output of subSV(), or supply the dependent-variable vector with Y and its sampling covariance matrix with V_Y.
  • PREDICTORS: A numeric vector or matrix of predictor values. Its rows must align with the entries in OBJECT$subS or Y.

Optional

  • INTERCEPT: Logical. If TRUE, the default, the function adds an intercept. If FALSE, the regression is fit without an intercept.

Output

summaryGLS() returns a matrix containing the beta estimate, p value, standard error, and z statistic for each estimated coefficient. The number following each beta corresponds to the column number for that set of predictors.

Example

Here we show the output of summaryGLS() for the primary GLS analysis from Ennis et al. (XXXX), testing trends in the genetic correlations between neuroticism and MDD across the severity spectrum of negative emotionality. We supply summaryGLS() with the subsetted S and V matrices obtained using subSV(). These contain the 12 genetic correlations between the 12 neuroticism cutpoints and MDD, which is the 13th trait in the LDSC object. The index values mapping the genetic correlations in S to their shared estimation errors in V are obtained using indexS().

# Case prevalences for each neuroticism cutpoint in the primary GLS analysis
PREVS <- c(
  0.85556773, 0.75867410, 0.64767354,
  0.53828939, 0.43558773, 0.34194621,
  0.25818629, 0.19038580, 0.13159727,
  0.08571119, 0.04947786, 0.02166147
)

PREDICTORS <- -qnorm(PREVS)

summaryGLS(
  OBJECT = subSV(
    LDSC_OBJECT = LDSCoutput,
    INDEXVALS = indexS(MATRIX = LDSCoutput$S, R = FALSE)[1:12, 13],
    TYPE = "S_Stand"
  ),
  PREDICTORS = PREDICTORS,
  INTERCEPT = TRUE
)
##        betas        pvals         SE         Z
## b0 0.5439339 8.405734e-24 0.05407569 10.058752
## b1 0.1071030 1.922670e-04 0.02872190  3.728966

summaryGLSbands()

summaryGLSbands() is an extension of summaryGLS() that, in addition to returning the regression parameters, produces a regression plot with error bands for one focal predictor. A scatterplot is generated with the parameters of interest mapped onto the y-axis and the supplied predictor values mapped onto the x-axis. The axis labels may be defined using XLAB and YLAB, which are empty by default. The x- and y-axis ranges may be modified using XCOORDS and YCOORDS, which default to c(-2, 2) and c(0, 1), respectively.

A set of control variables may be included using CONTROLVARS. The coefficients for these control variables are included in the GLS model but are not used to draw the focal regression line in the scatterplot.

Error bands are provided by default and may be toggled off by setting BANDS = FALSE. By default, the bands are set to plus or minus one standard error and may be changed using BAND_SIZE.

Arguments

Required

  • OBJECT or Y and V_Y: Supply either a list containing subS and subV, such as the output of subSV(), or supply the dependent-variable vector with Y and its sampling covariance matrix with V_Y.
  • PREDICTORS: The focal predictor, supplied as a numeric vector or one-column matrix. Its rows must align with the outcome.

Optional

  • INTERVALS: The number of evenly spaced x-values used to draw the fitted line and error bands. The default is 20.
  • CONTROLVARS: An optional numeric vector or matrix of control variables. Its rows must align with the outcome and focal predictor.
  • INTERCEPT: Logical. If TRUE, the default, the function adds an intercept. If FALSE, the regression is fit without an intercept.
  • QUAD: Logical. If TRUE, a squared term for the focal predictor is included. The default is FALSE.
  • BANDS: Logical. If TRUE, the default, error bands are displayed.
  • BAND_SIZE: The multiplier applied to the pointwise standard errors. The default value of 1 produces bands of plus or minus one standard error.
  • XLAB: Label for the x-axis. The default is an empty string.
  • YLAB: Label for the y-axis. The default is an empty string.
  • TITLELAB: Plot title. The default is an empty string.
  • XCOORDS: A two-value vector giving the displayed limits of the x-axis. The default is c(-2, 2).
  • YCOORDS: A two-value vector giving the displayed limits of the y-axis. The default is c(0, 1).

Do not include a column of ones in PREDICTORS or CONTROLVARS. Use the INTERCEPT argument to specify whether the model contains an intercept.

Output

The function prints the same GLS coefficient table returned by summaryGLS() and displays the regression plot. When BANDS = TRUE, the plot includes the requested pointwise error bands.

Example

The following example models the genetic correlations between the 12 neuroticism cut points and a MDD, which is the 15th trait in the LDSC object:

PREDICTORS <- -qnorm(c(
  0.85556773, 0.75867410, 0.64767354,
  0.53828939, 0.43558773, 0.34194621,
  0.25818629, 0.19038580, 0.13159727,
  0.08571119, 0.04947786, 0.02166147
))

summaryGLSbands(
  OBJECT = subSV(
    LDSC_OBJECT = LDSCoutput,
    INDEXVALS = indexS(MATRIX = LDSCoutput$S, R = FALSE)[1:12, 13],
    TYPE = "S_Stand"
  ),
  PREDICTORS = PREDICTORS,
  INTERVALS = 30,
  XLAB = "Neuroticism Severity (Z-Distribution)",
  YLAB = "Genetic Correlation w/ \n MDD",
  QUAD = TRUE,
  XCOORDS = c(-1.2, 2.1),
  YCOORDS = c(0, 1)
)
##         betas        pvals         SE        Z
## b0 0.51646528 3.187222e-20 0.05606158 9.212464
## b1 0.09112062 2.373223e-03 0.02998322 3.039054
## b2 0.03988152 6.326221e-02 0.02147234 1.857344

For more details on GLS regression, please see:

Aitken, A. C. (1935). On least squares and linear combination of observations. Proceedings of the Royal Society of Edinburgh, 55, 42–48. https://doi.org/10.1017/S0370164600014346

Savalei, V. (2014). Understanding Robust Corrections in Structural Equation Modeling. Structural Equation Modeling: A Multidisciplinary Journal, 21(1), 149–160. https://doi.org/10.1080/10705511.2013.824793

https://en.wikipedia.org/wiki/Generalized_least_squares