In johow/kinlab: Kinship and Reproduction in Historical Populations

Funktionsumfang des kinlab package für R, welches aus einem Forschungsprojekt mit obigem Titel (und einer Förderung durch die Deutsche Forschungsgemeinschaft in den Jahren 2012 bis 2015) hervorgegangen ist.

Stammbaumstrukuren

Aufbauend auf Konzepten zur genetischen Reproduktion (identity by descent) lassen sich Verwandtschaftsunterschiede innerhalb eines Stammbaums effizient in einer Verwandtschaftsmatrix abbilden (d.h. als sparse matrix, Therneau, 2015). In den Daten für die beiden untersuchten Populationen der historischen Krummhörn (18.-19. Jahrhundert) bzw. der frühen Québec-Region (17./18. Jahrhundert) lassen sich auf diese Weise jeweils die große Mehrheit aller dokumentierten Individuen einer einzigen großen Genealogie zuordnen -- mit über 74 000 Mitgliedern im Falle der Krummhörn-Datenbank (\glqq{}KH\grqq{}) bzw. über 300 000 Mitgliedern im Falle des Bevölkerungsregisters von Québec. Für die KH lassen sich innerhalb dieses großen Ausgangsstammbaum für die Zeit zwischen den Jahren 1720 und 1874 über 2000 Mütter identifizieren, welche sog. ,,vollständig dokumentierte'' Familien gegründet haben, und dabei -- genau wie ihr(e) Ehepartner -- ebenfalls aus ,,vollständig dokumentierten'' Familien stammen. Das ist Voraussetzung, um die elterliche Verwandtschaft der gegründeten Indexfamilien hinreichend genau identifizieren zu können. Für die Demonstrationszwecke dieser Seite musste allerdings die Datenbasis auf eine Stichprobe von lediglich 32 randomisierten Müttern beschränkt werden, um nicht gleichzeitig auch die gesamte KH-Datenbank mitveröffentlichen zu müssen.

Zeitgenossen in ,,begrenzten'' Stammbäumen

Für komplexe Analysen (sog. spatiotemporal event histories) sind diese riesigen Ausgangsstammbäume allerdings denkbar ungeeignet, weil die enorme Datenfülle ganz überwiegend Ereignisse von nicht- bzw. äußerst gering verwandter Individuen aus (früher oder später) angeheirateten Verzweigungen eines individuellen Stammbaums enthält. Im Gegensatz zu den lediglich sehr geringen Verwandtschaftsunterschiede im groben Populationsmittel, welche aufgrund der fehlenden Informationen für die Gründergeneration auch kaum als valide gelten können, spielen Unterschiede in der Verwandtschaft höchstwahrscheinlich eher in einem bestimmten \glqq{}genealogischen Nahfeld\grqq{}, d.h. für engere Familienangehörige von Mutter und Vater eine größere Rolle.

Daher ist es sinnvoll, Mitglieder eines bestimmten Verwandtennetzwerkes auf einen bestimmten Mindestwert in der Verwandtschaft zu einer Indexmutter bzw. ihre(n) Partner(n) zu beschränken. In diesem Beispiel wurden daher lediglich Verwandte bis zum Grad von Cousins oder Cousinen in der Analyse mit berücksichtigt, dabei dient $\phi$ als Maß für die statistische Wahrscheinlichkeit einer gemeinsamen Abstammung zwischen zwei (jeweils \glqq{}zufällig gezogenen\grqq{}) Allelen eines bestimmtes Gen zwischen zwei Individuen (also hier gilt $\phi >=0.0625$).

Dieser Wert kann durch den Benutzer frei gewählt werden und sollte sicherlich selbst Gegenstand weiterer Untersuchungen sein, allerdings steht die entsprechende Funktion (welche auch schnell zu rechenaufwändig für eine solche Web-Anwendung werden kann) hier aufgrund der benötigten Datenbasis aus lizenzrechtlichen Gründen bisher nicht zur Verfügung.

Verwandtschaftsunterschiede lassen sich z.B. als abgestufte Farbtöne zwischen grün (\glqq{}Verwandtschaft relativ hoch\grqq{}) und rot (\glqq{}Verwandtschaft relativ niedrig, d.h. i.d.R. 0\grqq{}) visualisieren, dabei kann auch der Überlebensstatus der Individuen z.B. durch gefüllte Symbole repräsentiert werden, und die Stammbaummitglieder können mit zusätzlichen Informationen (s.u.) beschriftet werden.

library(shiny)
inputPanel(
  selectInput("id_mother", label = "ID of mother:",
              choices = as.numeric(names(kh.data::kh_ped)), selected = as.numeric(names(kh.data::kh_ped)[1])),
  dateInput("evdat_adjust", label = "Date of interest:", value = "1800-01-01", min = "1720-01-01", max = "1874-12-31",
  format = "yyyy-mm-dd", startview = "decade", weekstart = 0,
  language = "en"), 
  shiny::checkboxInput("X", label = "show X chromosomal relatedness instead of autosomal relatedness", value = FALSE)
)

renderPlot({
  kinlab::plot_pedigree(input$id_mother, paste(as.Date(input$evdat_adjust, origin = "1970-01-01")),  kh.data::kh_ind, ped = kh.data::kh_ped[[paste(input$id_mother)]], evmat = kh.data::kh_mat, label = NULL, cex = 0.6, chrtype = ifelse(input$X==FALSE, "autosome", "X"))
  legend("top", pch = c(22,16, 16, 16, 16, 1), col = 
         c("gray", "gray", "green", "red", "gray", "gray"), title = paste0("Sex, relatedness, and survival status"),legend = c("male", "female", 
                  "related", "unrelated", 
                  "alive", "not alive"), horiz = TRUE)
})

 renderText({
  paste0("Actual spouse of ID ",input$id_mother, " on ", as.Date(input$evdat_adjust, origin = "1970-01-01"), " is ID ", kinlab::actual_spouse(input$id_mother, paste(as.Date(input$evdat_adjust, origin = "1970-01-01")), kh.data::kh_ind, kh.data::kh_fam), ".\n")
})



renderPlot({
  kinlab::plot_pedigree(
    kinlab::actual_spouse(input$id_mother, paste(as.Date(input$evdat_adjust, origin = "1970-01-01")), kh.data::kh_ind, kh.data::kh_fam), evdat =  paste(as.Date(input$evdat_adjust, origin = "1970-01-01")),  df_ind = kh.data::kh_ind, ped = kh.data::kh_ped[[paste(input$id_mother)]], evmat = kh.data::kh_mat, label = NULL, cex = 0.6, chrtype = ifelse(input$X==FALSE, "autosome", "X"))
  legend("top", pch = c(22,16, 16, 16, 16, 1), col = 
         c("gray", "gray", "green", "red", "gray", "gray"), title = paste0("Sex, relatedness, and survival status"),legend = c("male", "female", 
                  "related", "unrelated", 
                  "alive", "not alive"), horiz = TRUE)
})

library(shiny)
library(ggplot2)
library(ggmap)
inputPanel(
  selectInput("id_mother1", label = "ID of mother:",
              choices = as.numeric(names(kh.data::kh_ped)), selected = as.numeric(names(kh.data::kh_ped)[1])),
  selectInput("event_choice", label = "Choose date of event", c("birthday", "first_child", "death"))
)

Geographische Distanzen

Für eine Vielzahl von Stammbaummitgliedern gestatten die Daten zu deren individuellen Ereignissen eine zuverlässige Schätzung ihres Aufenthaltortes zu bestimmten Zeitpunkten. Unter Verwendung von web services lassen sich leicht die exakten Geopositionen für alle rezenten (bzw. je nach Anbieter wie Google Maps API etc. verfügbaren) Ortschaften innerhalb von beiden Populationen abrufen und zur Distanzberechnung zwischen solchen Individuen zu einem bestimmten Zeitpunkt nutzen.

renderPlot({ 
kinlab::plot_kinmap(id=as.numeric(input$id_mother1),
 evdat=kinlab::as_date(kh.data::kh_mat[paste(input$id_mother1), ifelse(input$event_choice=="birthday", 1, 
                                                                   ifelse(input$event_choice=="death",
                                                                          dim(kh.data::kh_mat)[[2]], 2)),1]),
                      list_kin=kh.data::kh_kin,
                      list_geo=kh.data::kh_geo,
                      my_map=kh.data::kh_geo[[2]][[1]],
                      spit_results=FALSE,
                      throw_plots=TRUE)
})

Tabellen zu Verwandtenetzwerken an spezifischen Zeitpunkten

In die obige Darstellung des Stammbaums sind bereits zusätzliche Angaben zu den lebenden Mitgliedern zu einem bestimmten Zeitpunkt mit eingeflossen, nämlich das individuelle Alter (in Jahren) und der (geschätzte) geographische Aufenthaltsort (als Kürzel für das betreffende Kirchspiel). Diese individuellen Angaben zu den jeweils lebenden Verwandten lassen sich prinzipiell zu jedem beliebigen Zeitpunkt für alle Mütter innerhalb der Stichprobe z.B. in Tabellenform für weitere Analysen bereitstellen:

renderTable({ rbind(head(subset(kinlab::get_kinset(x=as.numeric(input$id_mother1),
 evdat=kinlab::as_date(kh.data::kh_mat[paste(input$id_mother1), ifelse(input$event_choice=="birthday", 1, 
                                                                   ifelse(input$event_choice=="death",
                                                                          dim(kh.data::kh_mat)[[2]], 2)),1]),
 ped = kh.data::kh_ped[[paste(input$id_mother1)]], df_ind = kh.data::kh_ind, df_fam = kh.data::kh_fam , evmat=kh.data::kh_mat, map_dist = kh.data::kh_geo[["dist"]]), relmom >= 0.0625 | (!is.na(reldad) & reldad >= 0.0625))), c(rep("...", ncol(head(kinlab::get_kinset(x=as.numeric(input$id_mother1),
 evdat=kinlab::as_date(kh.data::kh_mat[paste(input$id_mother1), ifelse(input$event_choice=="birthday", 1, 
                                                                   ifelse(input$event_choice=="death",
dim(kh.data::kh_mat)[[2]], 2)),1]),
 ped = kh.data::kh_ped[[paste(input$id_mother1)]], df_ind = kh.data::kh_ind, df_fam = kh.data::kh_fam , evmat=kh.data::kh_mat, map_dist = kh.data::kh_geo[["dist"]]))))), row.names = NULL)
})

 renderText({
  paste0("Kin network of ID ",input$id_mother1, " on ",kinlab::as_date(kh.data::kh_mat[paste(input$id_mother1), ifelse(input$event_choice=="birthday", 1, 
                                                                   ifelse(input$event_choice=="death",
                                                                          dim(kh.data::kh_mat)[[2]], 2)),1]), " (\"",input$event_choice,  "\") has ", nrow(subset(kinlab::get_kinset(x=as.numeric(input$id_mother1),
 evdat=kinlab::as_date(kh.data::kh_mat[paste(input$id_mother1), ifelse(input$event_choice=="birthday", 1, 
                                                                   ifelse(input$event_choice=="death",
                                                                          dim(kh.data::kh_mat)[[2]], 2)),1]),
 ped = kh.data::kh_ped[[paste(input$id_mother1)]], df_ind = kh.data::kh_ind, df_fam = kh.data::kh_fam , evmat=kh.data::kh_mat, map_dist = kh.data::kh_geo[["dist"]]), relmom >= 0.0625 | (!is.na(reldad) & reldad >= 0.0625))), " members.\n")
})

 # downloadHandler() takes two arguments, both functions.
# The content function is passed a filename as an argument, and
#   it should write out data to that filename.
output$downloadData <- downloadHandler(

  # This function returns a string which tells the client
  # browser what name to use when saving the file.
  filename = function() {
    paste(input$dataset, input$filetype, sep = ".")
  },

  # This function should write data to a file given to it by
  # the argument 'file'.
  content = function(file) {
    sep <- switch(input$filetype, "csv" = ",", "tsv" = "\t")

    # Write to a file specified by the 'file' argument
    write.table(datasetInput(), file, sep = sep,
                row.names = FALSE)
  }
)

[Die Ausgabe wurde hier aus aus Platzgründenauf die ersten sechs Zeilen beschränkt. Die Auswahl und das Format der angezeigten Spalten sind natürlich noch nicht ausgereift. Es kommt außerdem noch in den je nach eingegebenem Datum relativ seltenen Fällen, in welchen gar keine als ,,Verwandte'' definierte Mitglieder eines Stammbaums gelebt haben, aufgrund des genutzten Webservers zu Fehlermeldungen, welche aber an sich vermeidbar sind...]

Multistate Life History Models

Untersuchungen zu Einflüssen auf Parity Progression Ratios und mütterliches Überleben lassen sich auch gut anhand von sog. Multistate Event History-Modellen durchführen, bei denen die einzelnen Paritäten (Null bis Maximum) als transiente Stadien und der Tod (oder Abschluss der fruchtbaren Phase) als absorbierendes Stadium definiert sind:

wellPanel(
checkboxGroupInput("id_input",
                   "IDs of mothers in sample:",
                   choices =  names(kh.data::kh_ped),
                   selected = names(kh.data::kh_ped), 
                   inline=TRUE)
)

renderPlot({ 
test <- kinlab::built_ms(id_sample = as.numeric(input$id_input),
                         evmat = kh.data::kh_mat)


tmat <-  attr(test, "trans")
# new data, to check whether results are the same for transition 1 as
# those in appendix E.1 of Therneau & Grambsch (2000)
library(mstate)
library(colorspace)
statecols <- heat_hcl(dim(tmat)[[1]], c = c(80, 30), l = c(30, 90), power = c(1/5, 2))[c(dim(tmat)[[1]]:1)]
ord <- 1:dim(tmat)[[1]]

# probtrans
cx <- coxph(Surv(Tstart,Tstop,status)~1+strata(trans),
            data=test,method="breslow")
HvH <- msfit(cx,trans=tmat)
summary(cx)
pt <- probtrans(HvH,predt=0)
# predictions from state 1
pt[[1]]
plot(pt, ord = ord, xlab = "Age", xlim = c(15,50), las = 1, type = "filled", col = statecols[ord])
title("Maternal Life Histories in Sample")#
})

Vergleich von Restreproduktionswerten

Durch Berechnung von Überlebenswahrscheinlichkeiten und der Parity Progression Ratios lassen sich auch Unterschiede in den mittleren Restreproduktionswerten zwischen bestimmten Gruppen auswerten. Hier als einfaches Beispiel als Unterschied zwischen ,,frühen'' und ,,späten'' Familien im Hinblick auf ein bestimmtes Datum.

library(shiny)
inputPanel(dateInput("evdat_adjust2", label = "Date of interest (Before-After)", value = "1800-01-01", min = "1750-01-01", max = "1829-12-31",
  format = "yyyy-mm-dd", startview = "decade", weekstart = 0,
  language = "en")
)
renderPlot({ 

kinlab::plot_rrv(df = data.frame(momid = as.numeric(unique(
  kh.data::kh_ind$momid[kh.data::kh_ind$momid > 0])),
 time = ifelse(as.numeric(format(as.Date(unlist(
lapply(as.numeric(unique(kh.data::kh_ind$momid[
  kh.data::kh_ind$momid > 0])), kinlab::get_date, kh.data::kh_ind)),
origin="1970-01-01"), "%Y"))<as.numeric(format(input$evdat_adjust2, "%Y")), "early", "late"))
, var = "time", df_ped = kh.data::kh_ind, 
                evmat = kh.data::kh_mat, evmat_bak = kh.data::kh_mat)
})

Weitere Arbeitsplanung

Unbedingt erforderlich ist aus meiner Sicht eine Anbindung der Tabellen, welche das spezifische Verwandtennetzwerk zu einem bestimmten Zeitpunkt beschreiben an die gegen Ende beschriebenen Multistate-Modelle zu mütterlichen Biographien. Theoretisch sollten hier doch deutliche Zusammenhänge sichtbar werden, oder? Für alle Ideen, Vorschläge, Kritik wäre ich sehr dankbar...

johow/kinlab documentation built on July 5, 2019, 4:23 p.m.