---
title: "Modeling and Relatedness"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{modelingandrelatedness}
  %\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 within 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)
require(EasyMx)
require(OpenMx)
```

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}
require(dplyr)
# require(purrr)

data(twinData, package = "OpenMx")
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}
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}
run2 <- emxTwinModel(
  model = "Cholesky",
  relatedness = "RCoef",
  data = mzdzData,
  use = selVars,
  run = TRUE, name = "TwCh"
)

summary(run2)
```