---
title: "Modeling variance components"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Modeling variance components}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)
```


# Introduction

This vignette provides a detailed guide to specific functions in the `BGmisc` package that aid in the identification and fitting of variance component models common in behavior genetics. We will explore key functions such as `identifyComponentModel`, providing practical examples and theoretical background. Identification ensures a unique set of parameters that define the model-implied covariance matrix, preventing free parameters from trading off one another.

## Loading Required Libraries

```{r setup, include=FALSE}
library(BGmisc)

if (!requireNamespace("OpenMx", quietly = TRUE)) {
  message("Please install OpenMx to run the model fitting examples.")
} else {
  library(OpenMx)
}
if (!requireNamespace("EasyMx", quietly = TRUE)) {
  message("Please install EasyMx to run the model fitting examples.")
} else {
  library(EasyMx)
}
```

Ensure that the `BGmisc` package is installed and loaded.

Ensure that the following dependencies are installed before proceeding as they provide us with behavior genetic data and models:

- `EasyMx`

- `OpenMx`


```r
library(BGmisc)
library(EasyMx)
library(OpenMx)
```

Note: If any of the libraries are not installed, you can install them using install.packages("`package_name`").

# Working with Variance Component Models

In this section, we will demonstrate core functions related to the identification and fitting of variance component models.

## Using `comp2vech` Function

The `comp2vech` function is used to vectorize a components model. The function is often used in conjunction with the identification process. In this example, we apply it to a list of matrices:


```{r}
comp2vech(list(
  matrix(c(1, .5, .5, 1), 2, 2),
  matrix(1, 2, 2)
))
```

The result showcases how the matrices have been transformed, reflecting their role in subsequent variance component analysis.

## Using `identifyComponentModel` Function

The `identifyComponentModel` function helps determine if a variance components model is identified. It accepts relatedness component matrices and returns information about identified and non-identified parameters. 

Here's an example using the classical twin model *with only MZ twins*:
```{r}
identifyComponentModel(
  A = list(matrix(1, 2, 2)),
  C = list(matrix(1, 2, 2)),
  E = diag(1, 2)
)
```

As you can see, the model is not identified. We need to add an additional group so that we have sufficient information. Let us add the rest of the classical twin model, in this case DZ twins.

```{r}
identifyComponentModel(
  A = list(matrix(c(1, .5, .5, 1), 2, 2), matrix(1, 2, 2)),
  C = list(matrix(1, 2, 2), matrix(1, 2, 2)),
  E = diag(1, 4)
)
```

As you can see the model is identified, now that we've added another group. Let us confirm by fitting a model. First we prepare the data.

```{r, include=FALSE}
library(dplyr)
if (!requireNamespace("OpenMx", quietly = TRUE)) {
  # if OpenMx isn't available
  n_subjects <- 500
  require(dplyr)

  df_summary_data <- data.frame(
    age_mean = 34.45,
    age_sd = 10.12345,
    ht1_mean = 1.662,
    ht1_sd = 0.0896277,
    ht2_mean = 1.694,
    ht2_sd = 0.09903498
  )
  set.seed(12345)
  twinData <- data.frame(
    ht1 = rnorm(n_subjects, mean = df_summary_data$ht1_mean, sd = df_summary_data$ht1_sd),
    zyg = rep(c(1, 3), each = n_subjects / 2)
  )
  twinData$ht2 <- twinData$ht1 * ifelse(twinData$zyg == 1, 1, 0.5) +
    rnorm(n_subjects,
      mean = 0, sd = df_summary_data$ht2_sd
    ) +
    .1 * rnorm(n_subjects,
      mean = df_summary_data$ht2_mean,
      sd = df_summary_data$ht2_sd
    )

  twinData$ht2[twinData$zyg == 3] <- twinData$ht2[twinData$zyg == 3] + .5 * rnorm(sum(twinData$zyg == 3), mean = df_summary_data$ht2_mean, sd = df_summary_data$ht2_sd)
} else {
  data(twinData, package = "OpenMx")
}
```

```{r}
library(dplyr)


selVars <- c("ht1", "ht2")

mzdzData <- subset(
  twinData, zyg %in% c(1, 3),
  c(selVars, "zyg")
)

mzdzData$RCoef <- c(1, NA, .5)[mzdzData$zyg]


mzData <- mzdzData %>% filter(zyg == 1)
```


Let us fit the data with MZ twins by themselves. 

```{r}
if (!requireNamespace("EasyMx", quietly = TRUE)) {
  print("Please install EasyMx to run the model fitting examples.")
} else {
  library(EasyMx)
  run1 <- emxTwinModel(
    model = "Cholesky",
    relatedness = "RCoef",
    data = mzData,
    use = selVars,
    run = TRUE, name = "TwCh"
  )

  summary(run1)
}
```

As you can see the model was unsuccessful because it was not identified. But when we add another group, so that the model is identified, the model now fits.

```{r}
if (!requireNamespace("EasyMx", quietly = TRUE)) {
  print("Please install EasyMx to run the model fitting examples.")
} else {
  library(EasyMx)
  run2 <- emxTwinModel(
    model = "Cholesky",
    relatedness = "RCoef",
    data = mzdzData,
    use = selVars,
    run = TRUE, name = "TwCh"
  )

  summary(run2)
}
```

## Extended Pedigrees

The variance component framework used above for twin data extends naturally to
more complex family structures, such as multi-generational pedigrees. In this
setting, BGmisc can be used to derive the appropriate relatedness coefficients
(e.g., additive genetic, shared environmental, and dominance relationships)
from a pedigree object, and these coefficients can then be supplied to a
structural equation modeling package (such as EasyMx/OpenMx) to fit the model.

In practice, the workflow mirrors the twin example: (1) prepare a pedigree
dataset with individual IDs and parental links, (2) use BGmisc functions to
compute the relevant relatedness matrices for the pedigree, (3) specify a
variance component model (e.g., ACE, ADE, or custom structures) that uses
these matrices, and (4) estimate the model parameters via maximum likelihood.
Please see the vignette on "Modeling with pedigrees" for a detailed example of this process.
