---
title: "Kinship Links"
author: Mason Garrison
output: rmarkdown::html_vignette
bibliography: references.bib
vignette: >
  %\VignetteIndexEntry{Kinship Links}
  %\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 demonstrates how to use the discord package with non-NLSY data. Specifically, it shows how to construct kinship links from a pedigree, estimate additive genetic relatedness and shared environmental exposure, and fit a discordant-kinship regression model. We use data and tools from the {BGmisc} package to illustrate the process.


# Loading Packages and Data

We begin by loading the required packages and loading a built-in toy dataset from {BGmisc}.


```{r echo=TRUE, message=FALSE, warning=FALSE}
library(BGmisc)
library(tidyverse)
library(discord)

data(potter)
```


We rename the family ID column to avoid naming conflicts and generate a pedigree-encoded data frame. This creates a consistent family identifier for use in subsequent functions.

```{r}
df_potter <- potter

names(df_potter)[names(df_potter) == "famID"] <- "oldfam"

df_potter <- ped2fam(df_potter,
  famID = "famID",
  personID = "personID"
)
```

Next, we compute two component matrices based on the pedigree:

- add: the additive genetic relatedness matrix

- cn: the shared environment matrix, indicating whether kin were raised together (1) or apart (0)


```{r}
add <- ped2add(df_potter)
cn <- ped2cn(df_potter)
```

We then verify and, if necessary, repair sex coding to ensure compatibility with downstream pedigree-based operations.

```{r}
df_potter <- checkSex(df_potter,
  code_male = 1,
  code_female = 0,
  verbose = FALSE, repair = TRUE
)
```

```{r, fig.cap="Pedigree plot of the Potter dataset", fig.width=8, fig.height=6, echo=FALSE}
plotPedigree(df_potter) %>% suppressMessages()
```

# Create Kinship Links

The `com2links` function generates a data frame of kinship links based on the pedigree data. The `writetodisk` argument is set to FALSE to avoid writing the output to disk, and we filter out upper triangular values to focus on unique pairs.

We use the `com2links()` function to convert the additive genetic and shared environment component matrices into a long-form data frame of kin pairs. Each row in the resulting data frame represents a pair of individuals, along with their additive genetic relatedness and shared environmental status. Self-pairs and duplicate (upper-triangular) entries are removed.


```{r}
df_links <- com2links(
  writetodisk = FALSE,
  ad_ped_matrix = add, cn_ped_matrix = cn,
  drop_upper_triangular = TRUE
) %>%
  filter(ID1 != ID2)

df_links %>%
  slice(1:10) %>%
  knitr::kable()
```

From the full set of kin links, we extract two subsets based on expected genetic and environmental structure:

- Full siblings: additive relatedness = 0.5 and shared environment = 1

- Cousins: additive relatedness = 0.125 and shared environment = 0

```{r}
df_sim <- df_links %>%
  filter(addRel == .5) %>% # only full siblings %>%
  filter(cnuRel == 1) # only kin raised in the same home

df_cousin <- df_links %>%
  filter(addRel == .125) %>% # only cousins %>%
  filter(cnuRel == 0) # only kin raised in separate homes
```

# Simulate Phenotypes and Fit Model

To illustrate the modeling framework, we simulate outcome variables for the cousin pairs using `kinsim()`. The simulated data reflect a known variance structure: additive genetic correlation = 0.125, no shared environment, and residual (unique environment) variance = 0.4. Latent component scores are excluded from the final dataset, but they can be useful for debugging and understanding the underlying structure of the data.

```{r}
set.seed(1234)
syn_df <- discord::kinsim(
  mu_all = c(1, 1), cov_a = .4,
  cov_e = .4, c_all = 0,
  r_vector = df_cousin$addRel
) %>%
  select(-c(
    A1_1, A1_2, A2_1, A2_2,
    C1_1, C1_2, C2_1, C2_2,
    E1_1, E1_2, E2_1, E2_2,
    r
  ))
```

We bind the simulated outcome data to the cousin link data to prepare it for modeling.

```{r}
data_demo <- cbind(df_cousin, syn_df)


summary(data_demo)
```

We then use `discord_regression()` to fit a discordant-kinship model, predicting y1 from y2. Based on the structure of the data, we expect that there will be a significant association between the two outcome variables, as there is a known overlapping non-shared environment covariance.

The model is fit using the `discord_regression()` function, which takes the following arguments:

```{r}
model_output <- discord_regression(
  data = data_demo,
  outcome = "y1",
  predictors = "y2",
  id = "id",
  sex = NULL,
  race = NULL,
  pair_identifiers = c("_1", "_2")
)
summary(model_output)
```


The output of the model includes estimates of the regression coefficients, standard errors, and p-values for the association between the two outcome variables. 

This example demonstrates how the discord package can be used with custom kinship data derived from pedigrees. By extracting component matrices for additive genetic relatedness and shared environment using {BGmisc}, and converting them to linked pair data, we can apply discordant-kinship regression outside of the NLSY context. This enables flexible applications of the method to any dataset where kin structure is known or derivable.

# References