Data Cleaning

For this example, we will load the following packages.

# For easy data manipulation
library(dplyr)
# For kinship linkages
library(NlsyLinks)
# For discordant-kinship regression
library(discord)
# To clean data frame names
library(janitor)
# tidy up output
library(broom)
# pipe
library(magrittr)

data(data_flu_ses)

After preprocessing, we obtain a data frame that contain subject identifiers, demographic information such as race and sex, and behavioral variables such as flu vaccination totals and SES at age 40. A random slice of this dataset is shown below.

Using kinship data from the NlsyLinks package, we restructure the dataset to better lend itself to discordant kinship analysis. For each kin pair, the function CreatePairLinksSingleEntered() takes a data set like the one above, [a specification of the NLSY database and the kin’s relatedness], and the variables of interest. It returns a data frame where each row represents a kin pair, and each variable appears twice—once for each sibling—using suffixes to distinguish individuals.

In this example, we examine the relationship between total flu vaccinations (received between 2006 and 2016) and SES at age 40, focusing on full siblings. The variable names used in this linkage are drawn from the preprocessed data previewed above.

# Get kinship links for individuals with the following variables:
link_vars <- c(
  "FLU_total", "FLU_2008", "FLU_2010",
  "FLU_2012", "FLU_2014", "FLU_2016",
  "S00_H40", "RACE", "SEX"
)

We now link the subjects by the specified variables using CreatePairLinksSingleEntered(), from the {NlsyLinks} package.

# Specify NLSY database and kin relatedness
link_pairs <- Links79PairExpanded %>%
  filter(RelationshipPath == "Gen1Housemates" & RFull == 0.5)

df_link <- CreatePairLinksSingleEntered(
  outcomeDataset = data_flu_ses,
  linksPairDataset = link_pairs,
  outcomeNames = link_vars
)

We have saved this data frame as df_link. A random subset of this data is shown below::¹

This data is almost ready for analysis, but we want to ensure that the data are representative of actual trends. The FLU_total column is simply a sum of the biennial survey responses. So for a given sibling-pair, one or both individuals may not have responded to the survey indicating their vaccination status. If that’s the case, we want to exclude those siblings to reduce [non-response bias], by examining the biennial responses and removing any rows with missingness.

# Take the linked data, group by the sibling pairs and
# count the number of responses for flu each year. If there is an NA,
# then data is missing for one of the years, and we omit it.
consistent_kin <- df_link %>%
  group_by(SubjectTag_S1, SubjectTag_S2) %>%
  count(
    FLU_2008_S1, FLU_2010_S1,
    FLU_2012_S1, FLU_2014_S1,
    FLU_2016_S1, FLU_2008_S2,
    FLU_2010_S2, FLU_2012_S2,
    FLU_2014_S2, FLU_2016_S2
  ) %>%
  na.omit()

# Create the flu_modeling_data object with only consistent responders.
# Clean the column names with the {janitor} package.
flu_modeling_data <- semi_join(df_link,
  consistent_kin,
  by = c(
    "SubjectTag_S1",
    "SubjectTag_S2"
  )
) %>%
  clean_names()

To avoid violating assumptions of independence, in our analysis we specify that the kin-pairs should be from unique households (i.e. we randomly select one sibling pair per household).

flu_modeling_data <- flu_modeling_data %>%
  group_by(extended_id) %>%
  slice_sample() %>%
  ungroup()

The data we will use for modeling now contains additional information for each member of the kin pair, including sex and race of each individual, flu vaccination status for the biennial survey between 2006-2016, and a total flu vaccination count for that period. The total vaccination count ranges from 0 - 5, where 0 indicates that the individual did not get a vaccine in any year between 2006-2016 and 5 indicates that an individual got at least 5 vaccines between 2006-2016. Although our data set has individual years, we focused on the aggregate as we felt that was a measure of general tendency. A subset of the data to use in this regression looks like:

Modeling and Interpretation

To perform the regression using the {discord} package, we supply the data frame and specify the outcome and predictors. It also requires a kinship pair id, extended_id in our case, as well as pair identifiers – the column name suffixes that identify to which kin a column’s values correspond (“_s1” and “_s2” in our case).² Optional, although recommended, are columns containing sex and race information to control for as additional covariates. In our case, these columns are prefixed “race” and “sex”. Per the pre-processing script, these columns contain dummy variables where the reference group for race is “non-Black, non-Hispanic” and the reference group for sex is female.

By entering this information into the discord_regression() function, we can run the model as such:

# Setting a seed for reproducibility
set.seed(18)
flu_model_output <- discord_regression(
  data = flu_modeling_data,
  outcome = "flu_total",
  predictors = "s00_h40",
  id = "extended_id",
  sex = "sex",
  race = "race",
  pair_identifiers = c("_s1", "_s2")
)

The default output of discord_regression() is an lm object. The metrics for our regression can be summarized as follows:

Looking at this output, the intercept can be thought of as the average difference in outcomes between siblings, controlling for all other variables. That is, it looks like the average difference for two sisters of a non-minority ethnic background (the reference groups for sex and race) is approximately 1.5. The term flu_total_mean is essentially an extra component of the intercept that captures some non-linear trends and allows the difference score to change as a function of the average predictors. Here, this is the mean socioeconomic status for the siblings, s00_h40_mean. We also accounted for sex and race, neither of which have a statistically significant effect on the differences in flu vaccine shots between siblings (different families) or within a sibling pair (same family).

The most important metric from the output, though, is the difference score, s00_h40_diff. Here, it is statistically significant. An interpretation of this result might be, “the difference in socioeconomic status between siblings at age 40 is positively associated with the difference in the number of flu vaccinations received between 2006-2016.” This finding means that a sibling with 10% higher SES is expected to have 0.0553768 average in flu shots.

The goal of performing a discordant-kinship regression is to see whether there is a significant difference in some behavioral measure while controlling for as much gene-and-environmental variance as possible. In this section, we walked through an analysis showing a statistically significant difference in the number of flu shots a sibling received and their socioeconomic status. From this finding, we could not claim the relationship is causal. However, we cannot eliminate causality because there are statistically significant within- and between-family differences in our predictors and outcomes.