diff --git a/.Rbuildignore b/.Rbuildignore
index cf44746..7b0fa46 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -2,3 +2,5 @@
^\.Rproj\.user$
^README\.Rmd$
^LICENSE\.md$
+^cran-comments\.md$
+^CRAN-SUBMISSION$
diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION
new file mode 100644
index 0000000..04b6fae
--- /dev/null
+++ b/CRAN-SUBMISSION
@@ -0,0 +1,3 @@
+Version: 1.0.3
+Date: 2023-06-04 13:27:38 UTC
+SHA: f54cd023b2e6bb2de1e74d1d0c89c13828149c44
diff --git a/DESCRIPTION b/DESCRIPTION
index f14db16..62a4398 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,26 +1,29 @@
Package: DemoKin
-Title: Demokin
-Description: Estimate population kin counts and its distribution by type and age
-Version: 1.0.0
+Title: Estimate Population Kin Distribution
+Description: Estimate population kin counts and its distribution by type, age and sex.
+ The package implements one-sex and two-sex framework for studying living-death availability,
+ with time varying rates or not, and multi-stage model.
+Version: 1.0.3
Authors@R: c(
person("Iván", "Williams", email = "act.ivanwilliams@gmail.com", role = "cre"),
person("Diego", "Alburez-Gutierrez", email = "alburezgutierrez@demogr.mpg.de", role = "aut"),
- person("Xi", "Song", email = "xisong@sas.upenn.edu", role = "ctb"))
+ person("Xi", "Song", email = "xisong@sas.upenn.edu", role = "ctb"),
+ person("Caswell", "Hal", email = "caswell@demogr.mpg.de", role = "ctb"))
License: MIT + file LICENSE
Encoding: UTF-8
LazyData: true
Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
Suggests:
knitr,
rmarkdown,
- testthat (>= 3.0.0)
+ testthat (>= 3.0.0),
+ ggplot2
VignetteBuilder: knitr
Imports:
dplyr,
tidyr,
purrr,
- HMDHFDplus,
progress,
matrixcalc,
Matrix,
@@ -28,7 +31,9 @@ Imports:
stats,
igraph,
magrittr,
+ data.table,
lifecycle
+URL: https://github.com/IvanWilli/DemoKin
BugReports: https://github.com/IvanWilli/DemoKin/issues
Depends:
R (>= 2.10)
diff --git a/NAMESPACE b/NAMESPACE
index f5c7de1..9a087d8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,12 +1,19 @@
# Generated by roxygen2: do not edit by hand
export("%>%")
-export(demokin_codes)
-export(get_HMDHFD)
export(kin)
+export(kin2sex)
export(kin_multi_stage)
export(kin_time_invariant)
+export(kin_time_invariant_2sex)
+export(kin_time_invariant_2sex_cod)
export(kin_time_variant)
+export(kin_time_variant_2sex)
+export(kin_time_variant_2sex_cod)
+export(output_period_cohort_combination)
export(plot_diagram)
export(rename_kin)
+export(timevarying_kin)
+export(timevarying_kin_2sex)
+export(timevarying_kin_2sex_cod)
importFrom(magrittr,"%>%")
diff --git a/NEWS.md b/NEWS.md
index 1db9e6f..cfc0438 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,16 @@
+# DemoKin 1.0.3
+
+# DemoKin 1.0.2
+
+# DemoKin 1.0.1
+
# DemoKin 1.0.0
* Added a `NEWS.md` file to track changes to the package.
* Change stable/non-stable references to time varying/non-varying rates.
* Add multi-state process.
+# DemoKin 1.0.1
+* Submitted to CRAN
+* Death counts are placed in the age where Focal experience the death.
+* Aggregated kin types are allowed (`s` for older and younger sisters, for example).
diff --git a/R/aux_funs.R b/R/aux_funs.R
index e8c745b..a8e7bfd 100644
--- a/R/aux_funs.R
+++ b/R/aux_funs.R
@@ -1,49 +1,17 @@
-
-#' print kin codes
-#' @description Print kin codes and labels
-#' @export
-demokin_codes <- function(){
- codes <- c("coa", "cya", "d", "gd", "ggd", "ggm", "gm", "m", "nos", "nys", "oa", "ya", "os", "ys")
- caswell_codes <- c("t", "v", "a", "b", "c", "h", "g", "d", "p", "q", "r", "s", "m", "n")
- labels <- c("Cousins from older aunt", "Cousins from younger aunt", "Daughter", "Grand-daughter", "Great-grand-daughter", "Great-grandmother", "Grandmother", "Mother", "Nieces from older sister", "Nieces from younger sister", "Aunt older than mother", "Aunt younger than mother", "Older sister", "Younger sister")
- data.frame(DemoKin = codes, Caswell = caswell_codes, Label = labels, row.names = NULL)
-}
-
#' rename kin
-#' @description Rename kin labels depending consolidate some types
-#' @export
-rename_kin <- function(df, consolidate_column = "no"){
-
- stopifnot("Argument 'consolidate_column' should be 'no' or a valid column name" = consolidate_column %in% c("no", colnames(df)))
-
- if(consolidate_column == "no"){
-
- relatives <- c("Cousins from older aunt", "Cousins from younger aunt", "Daughter", "Grand-daughter", "Great-grand-daughter", "Great-grandmother", "Grandmother", "Mother", "Nieces from older sister", "Nieces from younger sister", "Aunt older than mother", "Aunt younger than mother", "Older sister", "Younger sister")
- names(relatives) <- c("coa", "cya", "d", "gd", "ggd", "ggm", "gm", "m", "nos", "nys", "oa", "ya", "os", "ys")
- } else {
-
- # Combine kin types irrespective of whether they come from older
- # or younger sibling lines
- consolidate_vec <- c("c", "c", "d", "gd", "ggd", "ggm", "gm", "m", "n", "n", "a", "a", "s", "s")
- names(consolidate_vec) <- c("coa", "cya", "d", "gd", "ggd", "ggm", "gm", "m", "nos", "nys", "oa", "ya", "os", "ys")
-
- # Rename kin types from codes to actual words
- relatives <- c("Cousins", "Daughter", "Grand-daughter", "Great-grand-daughter", "Great-grandmother", "Grandmother", "Mother", "Nieces", "Aunt", "Sister")
- names(relatives) <- unique(consolidate_vec)
-
- df <- as.data.frame(df)
- df$count <- df[ , consolidate_column]
-
- df <-
- df %>%
- dplyr::mutate(kin = consolidate_vec[kin]) %>%
- dplyr::group_by(age_focal, kin) %>%
- dplyr::summarise(count = sum(count)) %>%
- dplyr::ungroup()
-
-
- }
- df$kin <- relatives[df$kin]
- df
+#' @description Add kin labels depending the sex
+#' @details See table `demokin_codes` to know label options.
+#' @param df data.frame. A data frame with variable `kin` with `DemoKin` codes to be labelled.
+#' @param sex character. "f" for female, "m" for male or "2sex" for both sex naming.
+#' @return Add a column with kin labels in the input data frame.
+#' @export
+rename_kin <- function(df, sex = "f"){
+ if(!"kin" %in% names(df)) stop("Input df needs a column named kin.")
+ if(sex == "f") demokin_codes_sex <- DemoKin::demokin_codes[,c("DemoKin", "Labels_female")]
+ if(sex == "m") demokin_codes_sex <- DemoKin::demokin_codes[,c("DemoKin", "Labels_male")]
+ if(sex == "2sex") demokin_codes_sex <- DemoKin::demokin_codes[,c("DemoKin", "Labels_2sex")]
+ colnames(demokin_codes_sex) <- c("kin", "kin_label")
+ df %>%
+ dplyr::left_join(demokin_codes_sex)
}
diff --git a/R/data.R b/R/data.R
index 1bacb90..43be162 100644
--- a/R/data.R
+++ b/R/data.R
@@ -127,3 +127,34 @@
#' @source
#' Caswell (2021)
"kin_svk1990_caswell2020"
+
+#' Fertility for France (2012) by sex in Caswell (2022).
+#'
+#' Fertility for France (2012) by sex in Caswell (2022).
+#' @docType data
+#' @format
+#' A data.frame with age specific fertility rates by age and sex.
+#'
+#' @source
+#' Caswell (2022)
+"fra_asfr_sex"
+
+#' Survival probability for France (2012) by sex in Caswell (2022).
+#'
+#' Survival probability for France (2012) by sex in Caswell (2022).
+#' @docType data
+#' @format
+#' A data.frame with survival probabilities by age and sex.
+#'
+#' @source
+#' Caswell (2022)
+"fra_surv_sex"
+
+#' DemoKin codes, Caswell (2020) codes, and useful labels.
+#'
+#' DemoKin codes, Caswell (2020) codes, and useful labels.
+#' @docType data
+#' @format
+#' A data.frame with codes and labels for distinction between kin types.
+
+"demokin_codes"
diff --git a/R/get_HMDHFD.R b/R/get_HMDHFD.R
deleted file mode 100644
index b4ee78e..0000000
--- a/R/get_HMDHFD.R
+++ /dev/null
@@ -1,125 +0,0 @@
-#' Get time serie matrix data from HMD/HFD
-
-#' @description Wrapper function to get data of female survival, fertlity and population
-#' of selected country on selected period.
-
-#' @param country numeric. Country code from rom HMD/HFD.
-#' @param max_year numeric. Latest year to get data.
-#' @param min_year integer. Older year to get data.
-#' @param user_HMD character. From HMD.
-#' @param user_HFD character. From HFD.
-#' @param pass_HMD character. From HMD.
-#' @param pass_HFD character. From HFD.
-#' @param OAG numeric. Open age group to standarize output.
-#' @return A list wiith female survival probability, survival function, fertility and poopulation age specific matrixes, with calendar year as colnames.
-#' @export
-
-get_HMDHFD <- function(country = "SWE",
- min_year = 1900,
- max_year = 2018,
- user_HMD = NULL,
- pass_HMD = NULL,
- user_HFD = NULL,
- pass_HFD = NULL,
- OAG = 100){
-
- if(any(c(is.null(user_HMD), is.null(user_HFD), is.null(pass_HMD), is.null(pass_HFD)))){
- stop("The function needs HMD and HMF access.")
- }
-
- # source HMD HFD -----------------------------------------------------------------
- pop <- HMDHFDplus::readHMDweb(CNTRY = country, "Population", user_HMD, pass_HMD, fixup = TRUE) %>%
- dplyr::select(Year, Age, N = Female1)%>%
- dplyr::filter(Year >= min_year, Year <= max_year)
- lt <- HMDHFDplus::readHMDweb(country, "fltper_1x1", user_HMD, pass_HMD, fixup = TRUE) %>%
- dplyr::filter(Year >= min_year, Year <= max_year)
- asfr <- HMDHFDplus::readHFDweb(country, "asfrRR", user_HFD, pass_HFD, fixup = TRUE)%>%
- dplyr::filter(Year >= min_year, Year <= max_year)
-
- # list of yearly Leslie matrix ---------------------------------------------------
-
- age = 0:OAG
- ages = length(age)
- w = last(age)
- last_year = max(lt$Year)
- years = min_year:last_year
-
- # survival probability
- px <- lt %>%
- dplyr::filter(Age<=OAG) %>%
- dplyr::mutate(px = 1 - qx,
- px = ifelse(Age==OAG, 0, px)) %>%
- dplyr::select(Year, Age, px) %>%
- tidyr::pivot_wider(names_from = "Year", values_from = "px") %>%
- dplyr::select(-Age) %>%
- as.matrix()
- rownames(px) = 0:OAG
-
- # survival function
- Lx <- lt %>%
- dplyr::filter(Age<=OAG) %>%
- dplyr::mutate(Lx = ifelse(Age==OAG, Tx, Lx)) %>%
- dplyr::select(Year, Age, Lx) %>%
- tidyr::pivot_wider(names_from = "Year", values_from = "Lx") %>%
- dplyr::select(-Age) %>%
- as.matrix()
-
- Sx <- rbind(Lx[c(-1,-ages),]/Lx[-c(w:ages),],
- Lx[ages,]/(Lx[w,]+Lx[ages,]),
- Lx[ages,]/(Lx[w,]+Lx[ages,]))
- rownames(Sx) = 0:w
-
- # fertility
- fx <- asfr %>%
- dplyr::filter(Year >= min_year) %>%
- dplyr::select(-OpenInterval) %>%
- rbind(
- expand.grid(Year = years,
- Age = c(0:(min(asfr$Age)-1),(max(asfr$Age)+1):OAG),
- ASFR = 0)) %>%
- dplyr::arrange(Year, Age) %>%
- tidyr::spread(Year, ASFR) %>%
- dplyr::select(-Age) %>%
- as.matrix()
- rownames(fx) = 0:OAG
-
- # population
- Nx <- pop %>%
- dplyr::mutate(Age = ifelse(Age>OAG, OAG, Age)) %>%
- dplyr::group_by(Year, Age) %>% summarise(N = sum(N)) %>%
- dplyr::filter(Age<=OAG, Year >= min_year) %>%
- dplyr::arrange(Year, Age) %>%
- tidyr::spread(Year, N) %>%
- dplyr::select(-Age) %>%
- as.matrix()
- rownames(Nx) = 0:OAG
-
- # only return data with values
- if(any(is.na(colSums(Sx)))){
- warning("Asked for data out of HMDHFD range")
- Sx <- Sx[,!is.na(colSums(Sx))]
- }
- if(any(is.na(colSums(fx)))){
- warning("Asked for data out of HMDHFD range")
- fx <- fx[,!is.na(colSums(fx))]
- }
- if(any(is.na(colSums(Nx)))){
- warning("Asked for data out of HMDHFD range")
- Nx <- Nx[,!is.na(colSums(Nx))]
- }
-
- return(list(px=px,
- Sx=Sx,
- fx=fx,
- Nx=Nx))
-}
-
-# save data
- # swe_px <- swe_data$px
- # swe_Sx <- swe_data$Sx
- # swe_asfr <-swe_data$fx
- # swe_pop <- swe_data$Nx
- # save(swe_px, file = "data/swe_px.rda")
- # save(swe_Sx, file = "data/swe_Sx.rda")
- # save(swe_asfr, file = "data/swe_asfr.rda")
- # save(swe_pop, file = "data/swe_pop.rda")
diff --git a/R/kin.R b/R/kin.R
index e9a5669..ac782c2 100644
--- a/R/kin.R
+++ b/R/kin.R
@@ -1,20 +1,27 @@
-#' Estimate kin counts
+#' Estimate kin counts in a one-sex framework.
-#' @description Implementation of Goodman-Keyfitz-Pullum equations in a matrix framework.
+#' @description Implementation of Goodman-Keyfitz-Pullum equations in a matrix framework. This produce a matrilineal (or patrilineal)
+#' kin count distribution by kin and age.
#' @details See Caswell (2019) and Caswell (2021) for details on formulas. One sex only (female by default).
-#' @param U numeric. A vector (atomic) or matrix with probabilities (or survival ratios, or transition between age class in a more general perspective) with rows as ages (and columns as years in case of matrix, being the name of each col the year).
-#' @param f numeric. Same as U but for fertility rates.
+#' @param p numeric. A vector (atomic) or matrix with probabilities (or survival ratios, or transition between age class
+#' in a more general perspective) with rows as ages (and columns as years in case of matrix, being the name of each col the year).
+#' @param f numeric. Same as `p` but for fertility rates.
#' @param time_invariant logical. Constant assumption for a given `year` rates. Default `TRUE`.
-#' @param N numeric. Same as U but for population distribution (counts or `%`). Optional.
-#' @param pi numeric. Same as U but for childbearing distribution (sum to 1). Optional.
+#' @param n numeric. Only for `time_invariant = FALSE`. Same as `p` but for population distribution (counts or `%`). Optional.
+#' @param pi numeric. Same as `U` but for childbearing distribution (sum to 1). Optional.
#' @param output_cohort integer. Vector of year cohorts for returning results. Should be within input data years range.
#' @param output_period integer. Vector of period years for returning results. Should be within input data years range.
#' @param output_kin character. kin types to return: "m" for mother, "d" for daughter,...
-#' @param birth_female numeric. Female portion at birth. This multiplies `f` argument. If `f` is already for female offspring, this needs to be set as 1.
+#' @param output_age_focal integer. Vector of ages to select (and make faster the run).
+#' @param birth_female numeric. Female portion at birth. This multiplies `f` argument. If `f` is already for female offspring,
+#' @param summary_kin logical. Whether or not include `kin_summary` table (see output details). Default `TRUE`.
+#' this needs to be set as 1.
#' @return A list with:
#' \itemize{
-#' \item{kin_full}{ a data frame with year, cohort, Focal´s age, related ages and type of kin (for example `d` is daughter, `oa` is older aunts, etc.), including living and dead kin at that age.}
-#' \item{kin_summary}{ a data frame with Focal´s age, related ages and type of kin, with indicators obtained processing `kin_full`, grouping by cohort or period (depending on the given arguments):}
+#' \item{kin_full}{ a data frame with year, cohort, Focal´s age, related ages and type of kin (for example `d` is daughter,
+#' `oa` is older aunts, etc.), including living and dead kin at that age.}
+#' \item{kin_summary}{ a data frame with Focal´s age, related ages and type of kin, with indicators obtained processing `kin_full`,
+#' grouping by cohort or period (depending on the given arguments):}
#' {\itemize{
#' \item{`count_living`}{: count of living kin at actual age of Focal}
#' \item{`mean_age`}{: mean age of each type of living kin.}
@@ -25,58 +32,76 @@
#' }
#' }
#' }
-
#' @export
-#'
-# get kin ----------------------------------------------------------------
-kin <- function(U = NULL, f = NULL,
- time_invariant = TRUE,
- N = NULL, pi = NULL,
- output_cohort = NULL, output_period = NULL, output_kin=NULL,
- birth_female = 1/2.04,
- stable = lifecycle::deprecated())
- {
+#' @examples
+#' # Kin expected matrilineal count for a Swedish female based on 2015 rates.
+#' swe_surv_2015 <- swe_px[,"2015"]
+#' swe_asfr_2015 <- swe_asfr[,"2015"]
+#' # Run kinship models
+#' swe_2015 <- kin(p = swe_surv_2015, f = swe_asfr_2015)
+#' head(swe_2015$kin_summary)
- age <- as.integer(rownames(U))
- years_data <- as.integer(colnames(U))
+kin <- function(p = NULL, f = NULL,
+ time_invariant = TRUE,
+ pi = NULL, n = NULL,
+ output_cohort = NULL, output_period = NULL, output_kin=NULL, output_age_focal = NULL,
+ birth_female = 1/2.04,
+ summary_kin = TRUE)
+ {
- if (lifecycle::is_present(stable)) {
- lifecycle::deprecate_warn("0.0.0.9000", "kin(stable)", details = "Used time_invariant")
- time_invariant <- stable
- }
+ # global vars
+ living<-dead<-age_kin<-age_focal<-cohort<-year<-total<-mean_age<-count_living<-sd_age<-count_dead<-mean_age_lost<-indicator<-value<-NULL
# kin to return
all_possible_kin <- c("coa", "cya", "d", "gd", "ggd", "ggm", "gm", "m", "nos", "nys", "oa", "ya", "os", "ys")
+ output_kin_asked <- output_kin
if(is.null(output_kin)){
output_kin <- all_possible_kin
}else{
+ if("s" %in% output_kin) output_kin <- c(output_kin, "os", "ys")
+ if("c" %in% output_kin) output_kin <- c(output_kin, "coa", "cya")
+ if("a" %in% output_kin) output_kin <- c(output_kin, "oa", "ya")
+ if("n" %in% output_kin) output_kin <- c(output_kin, "nos", "nys")
+ output_kin <- output_kin[!output_kin %in% c("s", "c", "a", "n")]
output_kin <- match.arg(tolower(output_kin), all_possible_kin, several.ok = TRUE)
}
- # if time dependent or not
+ # if is time dependent or not
+ age <- as.integer(rownames(p))
+ years_data <- as.integer(colnames(p))
if(time_invariant){
- if(!is.vector(U)) {
+ if(!is.vector(p)) {
output_period <- min(years_data)
- U <- U[,as.character(output_period)]
+ p <- p[,as.character(output_period)]
f <- f[,as.character(output_period)]
}
- kin_full <- kin_time_invariant(U = U, f = f,
+ kin_full <- kin_time_invariant(p = p, f = f, pi = pi,
output_kin = output_kin, birth_female = birth_female) %>%
dplyr::mutate(cohort = NA, year = NA)
}else{
if(!is.null(output_cohort) & !is.null(output_period)) stop("sorry, you can not select cohort and period. Choose one please")
- kin_full <- kin_time_variant(U = U, f = f, N = N, pi = pi,
+ kin_full <- kin_time_variant(p = p, f = f, pi = pi, n = n,
output_cohort = output_cohort, output_period = output_period,
output_kin = output_kin,
birth_female = birth_female)
message(paste0("Assuming stable population before ", min(years_data), "."))
}
- # reorder
- kin_full <- kin_full %>% dplyr::select(year, cohort, age_focal, kin, age_kin, living, dead)
+ # re-group if grouped type is asked
+ if(!is.null(output_kin_asked) & length(output_kin_asked)!=length(output_kin)){
+ if("s" %in% output_kin_asked) kin_full$kin[kin_full$kin %in% c("os", "ys")] <- "s"
+ if("c" %in% output_kin_asked) kin_full$kin[kin_full$kin %in% c("coa", "cya")] <- "c"
+ if("a" %in% output_kin_asked) kin_full$kin[kin_full$kin %in% c("oa", "ya")] <- "a"
+ if("n" %in% output_kin_asked) kin_full$kin[kin_full$kin %in% c("nos", "nys")] <- "n"
+ kin_full <- kin_full %>%
+ dplyr::summarise(living = sum(living), dead = sum(dead),
+ .by = c(kin, age_kin, age_focal, cohort, year))
+ }
- # summary
- # select period/cohort
+ # select period/cohort/age
+ if(!is.null(output_age_focal) & all(output_age_focal %in% 1:120)){
+ kin_full <- kin_full %>% dplyr::filter(age_focal %in% output_age_focal)
+ }
if(!is.null(output_cohort)){
agrupar <- "cohort"
} else if(!is.null(output_period)){
@@ -87,28 +112,35 @@ kin <- function(U = NULL, f = NULL,
agrupar_no_age_focal <- c("kin", agrupar)
agrupar <- c("age_focal", "kin", agrupar)
- kin_summary <- dplyr::bind_rows(
- kin_full %>%
- dplyr::rename(total=living) %>%
- dplyr::group_by(dplyr::across(dplyr::all_of(agrupar))) %>%
- dplyr::summarise(count_living = sum(total),
- mean_age = sum(total*age_kin)/sum(total),
- sd_age = (sum(total*age_kin^2)/sum(total)-mean_age^2)^.5) %>%
- tidyr::pivot_longer(count_living:sd_age, names_to = "indicator", "value"),
- kin_full %>%
- dplyr::rename(total=dead) %>%
- dplyr::group_by(dplyr::across(dplyr::all_of(agrupar))) %>%
- dplyr::summarise(count_dead = sum(total)) %>%
- dplyr::ungroup() %>%
- dplyr::group_by(dplyr::across(dplyr::all_of(agrupar_no_age_focal))) %>%
- dplyr::mutate(count_cum_dead = cumsum(count_dead),
- mean_age_lost = cumsum(count_dead * age_focal)/cumsum(count_dead)) %>%
- dplyr::ungroup() %>%
- tidyr::pivot_longer(count_dead:mean_age_lost, names_to = "indicator", "value")) %>%
+ # get summary indicators based on group variables. If it is asked
+ if(summary_kin){
+ kin_summary <- dplyr::bind_rows(
+ kin_full %>%
+ dplyr::rename(total=living) %>%
+ dplyr::group_by(dplyr::across(dplyr::all_of(agrupar))) %>%
+ dplyr::summarise(count_living = sum(total),
+ mean_age = sum(total*age_kin)/sum(total),
+ sd_age = (sum(total*age_kin^2)/sum(total)-mean_age^2)^.5) %>%
+ tidyr::pivot_longer(count_living:sd_age, names_to = "indicator", values_to = "value"),
+ kin_full %>%
+ dplyr::rename(total=dead) %>%
+ dplyr::group_by(dplyr::across(dplyr::all_of(agrupar))) %>%
+ dplyr::summarise(count_dead = sum(total)) %>%
+ dplyr::ungroup() %>%
+ dplyr::group_by(dplyr::across(dplyr::all_of(agrupar_no_age_focal))) %>%
+ dplyr::mutate(count_cum_dead = cumsum(count_dead),
+ mean_age_lost = cumsum(count_dead * age_focal)/cumsum(count_dead)) %>%
+ dplyr::ungroup() %>%
+ tidyr::pivot_longer(count_dead:mean_age_lost, names_to = "indicator", values_to = "value")) %>%
dplyr::ungroup() %>%
tidyr::pivot_wider(names_from = indicator, values_from = value)
# return
kin_out <- list(kin_full = kin_full, kin_summary = kin_summary)
+ }else{
+ # return
+ kin_out <- kin_full
+ }
+
return(kin_out)
}
diff --git a/R/kin2sex.R b/R/kin2sex.R
new file mode 100644
index 0000000..07dfa11
--- /dev/null
+++ b/R/kin2sex.R
@@ -0,0 +1,194 @@
+#' Estimate kin counts in a two-sex framework
+
+#' @description Implementation of two-sex matrix kinship model. This produces kin counts grouped by kin, age and sex of
+#' each relatives at each Focal´s age. For example, male cousins from aunts and uncles from different sibling's parents
+#' are grouped in one male count of cousins. Note that the output labels relative following female notation: the label `m`
+#' refers to either mothers or fathers, and column `sex_kin` determine the sex of the relative.
+#' @details See Caswell (2022) for details on formulas.
+#' @param pf numeric. A vector (atomic) or matrix with female probabilities (or survival ratios, or transition between age class in a more general perspective) with rows as ages (and columns as years in case of matrix, being the name of each col the year).
+#' @param pm numeric. A vector (atomic) or matrix with male probabilities (or survival ratios, or transition between age class in a more general perspective) with rows as ages (and columns as years in case of matrix, being the name of each col the year).
+#' @param ff numeric. Same as `pf` but for fertility rates.
+#' @param fm numeric. Same as `pm` but for fertility rates.
+#' @param time_invariant logical. Constant assumption for a given `year` rates. Default `TRUE`.
+#' @param sex_focal character. "f" for female or "m" for male.
+#' @param pif numeric. For using some specific age distribution of childbearing for mothers (same length as ages). Default `NULL`.
+#' @param pim numeric. For using some specific age distribution of childbearing for fathers (same length as ages). Default `NULL`.
+#' @param Hf numeric. A list where each list element (being the name of each list element the year) contains a matrix with cause-specific hazards for females with rows as causes and columns as ages, being the name of each col the age.
+#' @param Hm numeric. A list where each list element (being the name of each list element the year) contains a matrix with cause-specific hazards for males with rows as causes and columns as ages, being the name of each col the age.
+#' @param nf numeric. Only for `time_invariant = FALSE`. Same as `pf` but for population distribution (counts or `%`). Optional.
+#' @param nm numeric. Only for `time_invariant = FALSE`. Same as `pm` but for population distribution (counts or `%`). Optional.
+#' @param output_cohort integer. Vector of year cohorts for returning results. Should be within input data years range.
+#' @param output_period integer. Vector of period years for returning results. Should be within input data years range.
+#' @param output_kin character. kin types to return: "m" for mother, "d" for daughter,...
+#' @param output_age_focal integer. Vector of ages to select (and make faster the run).
+#' @param birth_female numeric. Female portion at birth. This multiplies `f` argument. If `f` is already for female offspring, this needs to be set as 1.
+#' @param summary_kin logical. Whether or not include `kin_summary` table (see output details). Default `TRUE`.
+#' @return A list with:
+#' \itemize{
+#' \item{kin_full}{ a data frame with year, cohort, Focal´s age, related ages and type of kin (for example `d` could be daughter or son depending `sex_kin`,
+#' `oa` is older aunts or uncles also depending `sex_kin` value, etc.), including living and dead kin at that age.}
+#' \item{kin_summary}{ a data frame with Focal´s age, related ages, sex and type of kin, with indicators obtained processing `kin_full`, grouping by cohort or period (depending on the given arguments):}
+#' {\itemize{
+#' \item{`count_living`}{: count of living kin at actual age of Focal}
+#' \item{`mean_age`}{: mean age of each type of living kin.}
+#' \item{`sd_age`}{: standard deviation of age of each type of living kin.}
+#' \item{`count_death`}{: count of dead kin at specific age of Focal.}
+#' \item{`count_cum_death`}{: cumulated count of dead kin until specific age of Focal.}
+#' \item{`mean_age_lost`}{: mean age where Focal lost her relative.}
+#' }
+#' }
+#' }
+#' @export
+#' @examples
+#' # Kin expected count by relative sex for a French female based on 2012 rates.
+#' fra_fert_f <- fra_asfr_sex[,"ff"]
+#' fra_fert_m <- fra_asfr_sex[,"fm"]
+#' fra_surv_f <- fra_surv_sex[,"pf"]
+#' fra_surv_m <- fra_surv_sex[,"pm"]
+#' fra_2012 <- kin2sex(fra_surv_f, fra_surv_m, fra_fert_f, fra_fert_m)
+#' head(fra_2012$kin_summary)
+#'
+# get kin ----------------------------------------------------------------
+kin2sex <- function(pf = NULL, pm = NULL, ff = NULL, fm = NULL,
+ time_invariant = TRUE,
+ sex_focal = "f",
+ birth_female = 1/2.04,
+ pif = NULL, pim = NULL,
+ nf = NULL, nm = NULL,
+ Hf = NULL, Hm = NULL,
+ output_cohort = NULL, output_period = NULL, output_kin=NULL,output_age_focal = NULL,
+ summary_kin = TRUE)
+ {
+
+ # global vars
+ living<-age_focal<-cohort<-year<-total<-mean_age<-count_living<-sd_age<-count_dead<-mean_age_lost<-indicator<-value<-sex_kin<-age_kin<-dead<-NULL
+ age <- as.integer(rownames(pf))
+ years_data <- as.integer(colnames(pf))
+
+ # kin to return
+ all_possible_kin <- c("coa", "cya", "d", "gd", "ggd", "ggm", "gm", "m", "nos", "nys", "oa", "ya", "os", "ys")
+ output_kin_asked <- output_kin
+ if(is.null(output_kin)){
+ output_kin <- all_possible_kin
+ }else{
+ if("s" %in% output_kin) output_kin <- c(output_kin, "os", "ys")
+ if("c" %in% output_kin) output_kin <- c(output_kin, "coa", "cya")
+ if("a" %in% output_kin) output_kin <- c(output_kin, "oa", "ya")
+ if("n" %in% output_kin) output_kin <- c(output_kin, "nos", "nys")
+ output_kin <- output_kin[!output_kin %in% c("s", "c", "a", "n")]
+ output_kin <- match.arg(tolower(output_kin), all_possible_kin, several.ok = TRUE)
+ }
+
+ # is cause of death specific or not
+ is_cod <- !is.null(Hf) & !is.null(Hm)
+
+ # if time dependent or not
+ if(time_invariant){
+ if(!is.vector(pf)) {
+ output_period <- min(years_data)
+ pf <- pf[,as.character(output_period)]
+ pm <- pm[,as.character(output_period)]
+ ff <- ff[,as.character(output_period)]
+ fm <- fm[,as.character(output_period)]
+ }
+ if(is_cod){
+ kin_full <- kin_time_invariant_2sex_cod(pf, pm, ff, fm,
+ sex_focal = sex_focal,
+ birth_female = birth_female,
+ pif = pif, pim = pim,
+ Hf = Hf, Hm = Hm,
+ output_kin = output_kin) %>%
+ dplyr::mutate(cohort = NA, year = NA)
+ }else{
+ kin_full <- kin_time_invariant_2sex(pf, pm, ff, fm,
+ sex_focal = sex_focal,
+ birth_female = birth_female,
+ pif = pif, pim = pim,
+ output_kin = output_kin) %>%
+ dplyr::mutate(cohort = NA, year = NA)
+ }
+
+ }else{
+ if(!is.null(output_cohort) & !is.null(output_period)) stop("sorry, you can not select cohort and period. Choose one please")
+ if(is_cod){
+ kin_full <- kin_time_variant_2sex_cod(pf = pf, pm = pm,
+ ff = ff, fm = fm,
+ sex_focal = sex_focal,
+ birth_female = birth_female,
+ pif = pif, pim = pim,
+ nf = nf, nm = nm,
+ Hf = Hf, Hm = Hm,
+ output_cohort = output_cohort, output_period = output_period,
+ output_kin = output_kin)
+ }else{
+ kin_full <- kin_time_variant_2sex(pf = pf, pm = pm,
+ ff = ff, fm = fm,
+ sex_focal = sex_focal,
+ birth_female = birth_female,
+ pif = pif, pim = pim,
+ nf = nf, nm = nm,
+ output_cohort = output_cohort, output_period = output_period,
+ output_kin = output_kin)
+ }
+ message(paste0("Assuming stable population before ", min(years_data), "."))
+ }
+
+ # reorder
+ kin_full <- kin_full %>% dplyr::select(year, cohort, age_focal, sex_kin, kin, age_kin, living, starts_with("dea"))
+
+ # re-group if grouped type is asked
+ if(!is.null(output_kin_asked) & length(output_kin_asked)!=length(output_kin)){
+ if("s" %in% output_kin_asked) kin_full$kin[kin_full$kin %in% c("os", "ys")] <- "s"
+ if("c" %in% output_kin_asked) kin_full$kin[kin_full$kin %in% c("coa", "cya")] <- "c"
+ if("a" %in% output_kin_asked) kin_full$kin[kin_full$kin %in% c("oa", "ya")] <- "a"
+ if("n" %in% output_kin_asked) kin_full$kin[kin_full$kin %in% c("nos", "nys")] <- "n"
+ kin_full <- kin_full %>%
+ dplyr::group_by(kin, age_kin, age_focal, sex_kin, cohort, year) %>%
+ dplyr::summarise_at(vars(c("living", dplyr::starts_with("dea"))), funs(sum)) %>%
+ dplyr::ungroup()
+ }
+
+ # summary
+ # select period/cohort/ge
+ if(!is.null(output_age_focal) & all(output_age_focal %in% 1:120)){
+ kin_full <- kin_full %>% dplyr::filter(age_focal %in% output_age_focal)
+ }
+ if(!is.null(output_cohort)){
+ agrupar <- "cohort"
+ } else if(!is.null(output_period)){
+ agrupar <- "year"
+ } else{
+ agrupar <- c("year", "cohort")
+ }
+ agrupar_no_age_focal <- c("kin", "sex_kin", agrupar)
+ agrupar <- c("age_focal", "kin", "sex_kin", agrupar)
+
+ # only return summary if is asked and is not cod
+ if(summary_kin & !is_cod){
+ kin_summary <- dplyr::bind_rows(
+ as.data.frame(kin_full) %>%
+ dplyr::rename(total=living) %>%
+ dplyr::group_by(dplyr::across(dplyr::all_of(agrupar))) %>%
+ dplyr::summarise(count_living = sum(total),
+ mean_age = sum(total*age_kin)/sum(total),
+ sd_age = (sum(total*age_kin^2)/sum(total)-mean_age^2)^.5) %>%
+ tidyr::pivot_longer(count_living:sd_age, names_to = "indicator", values_to = "value"),
+ as.data.frame(kin_full) %>%
+ dplyr::rename(total=dead) %>%
+ dplyr::group_by(dplyr::across(dplyr::all_of(agrupar))) %>%
+ dplyr::summarise(count_dead = sum(total)) %>%
+ dplyr::ungroup() %>%
+ dplyr::group_by(dplyr::across(dplyr::all_of(agrupar_no_age_focal))) %>%
+ dplyr::mutate(count_cum_dead = cumsum(count_dead),
+ mean_age_lost = cumsum(count_dead * age_focal)/cumsum(count_dead)) %>%
+ dplyr::ungroup() %>%
+ tidyr::pivot_longer(count_dead:mean_age_lost, names_to = "indicator", values_to = "value")) %>%
+ dplyr::ungroup() %>%
+ tidyr::pivot_wider(names_from = indicator, values_from = value)
+ kin_out <- list(kin_full = kin_full, kin_summary = kin_summary)
+ }else{
+ kin_out <- kin_full
+ }
+
+ return(kin_out)
+}
diff --git a/R/kin_multi_stage.R b/R/kin_multi_stage.R
index 7129134..829b232 100644
--- a/R/kin_multi_stage.R
+++ b/R/kin_multi_stage.R
@@ -2,110 +2,126 @@
#' @description Implementation of age-stage kin estimates (multi-state) by Caswell (2020). Stages are implied in length of input lists.
-#' @param U list. age elemnts with column-stochastic transition matrix with dimension for the state space, conditional on survival.
-#' @param f matrix. state-specific fertility (age in rows and states in columns).
-#' @param D matrix. survival probabilities by state (age in rows and states in columns)
-#' @param H matrix. assigns the offspring of individuals in some stage to the appropriate age class with 1 (age in rows and states in columns).
+#' @param U list. age elements with column-stochastic transition matrix with dimension for the state space, conditional on survival.
+#' @param f matrix. state-specific fertility (age in rows and states in columns). Is accepted also a list with for each age-class.
+#' @param D matrix. survival probabilities by state (age in rows and states in columns). Is accepted also a list for each state with survival matrices.
+#' @param H matrix. assigns the offspring of individuals in some stage to the appropriate age class (age in rows and states in columns). Is accepted also a list with a matrix for each state.
#' @param output_kin character. kin to return. For example "m" for mother, "d" for daughter. See the `vignette` for all kin types.
#' @param birth_female numeric. Female portion at birth.
+#' @param parity logical. parity states imply age distribution of mothers re-scaled to not have parity 0 when Focal born. Default `TRUE`.
#' @param list_output logical. Results as a list. Default `FALSE`.
#' @return A data frame with focal´s age, related ages and type of kin
#' (for example `d` is daughter, `oa` is older aunts, etc.), living and death kin counts, and specific stage. If `list_output = TRUE` then this is a list with elements as kin types.
#' @export
-#'
kin_multi_stage <- function(U = NULL, f = NULL, D = NULL, H = NULL,
- birth_female = 1/2.04,
- output_kin = NULL,
- list_output = FALSE){
+ birth_female = 1/2.04,
+ output_kin = NULL,
+ parity = FALSE,
+ list_output = FALSE){
+ # global vars
+ .<-age_kin<-stage_kin<-alive<-age_focal<-count<-NULL
+
+ # mandatory U as a list
if(!is.list(U)) stop("U must be a list with age length of elements, and stage transitiotn matrix for each one.")
- # stages and ages
+ # stages and age-classes
s <- ncol(U[[1]])
ages <- length(U)
age <- (1:ages)-1
- # build matrix structure from data.frame input
- H <- purrr::map(colnames(D), function(Y){
- Ht = matrix(0, nrow=ages, ncol=ages)
- Ht[1,] <- 1
- Ht
- })
- D <- purrr::map(colnames(D), function(Y){
+ # build H if it is not already a list
+ if(!is.list(H)){
+ H <- purrr::map(1:s, function(Y){
+ Ht = matrix(0, nrow=ages, ncol=ages)
+ Ht[1,] <- 1
+ Ht
+ })
+ }
+
+ # build D if it is not already a list
+ if(!is.list(D)){
+ D <- purrr::map(1:s, function(Y){
X <- D[,Y]
Dt = matrix(0, nrow=ages, ncol=ages)
Dt[row(Dt)-1 == col(Dt)] <- X[-ages]
Dt[ages, ages] = X[ages]
Dt
})
- f <- purrr::map(1:ages, function(Y){
- X <- f[Y,]
- ft = matrix(0, nrow=s, ncol=s)
- ft[1,] <- X
- ft
- })
-
- # build block matrix
+ }
+
+ # build f if it is not already a list
+ if(!is.list(f)){
+ f <- purrr::map(1:ages, function(Y){
+ X <- f[Y,]
+ ft = matrix(0, nrow=s, ncol=s)
+ ft[1,] <- X
+ ft
+ })
+ }
+
+ # build block matrices
bbU <- Matrix::bdiag(U)
bbF <- Matrix::bdiag(f) * birth_female
bbD <- Matrix::bdiag(D)
bbH <- Matrix::bdiag(H)
- # rearrange with conmutation matrix
+ # order transitions: first state within age, then age given state
K <- matrixcalc::commutation.matrix(s, ages)
Ut <- t(K) %*% bbD %*% K %*% bbU
ft <- t(K) %*% bbH %*% K %*% bbF
+
+ # focal transition but conditioned to survive
Gt <- Ut%*% MASS::ginv(diag(colSums(as.matrix(Ut))))
- # stable distribution mothers: age x stage
+ # stable distribution of mothers
At <- Ut + ft
A_decomp <- eigen(At)
lambda <- as.double(A_decomp$values[1])
wt <- as.double(A_decomp$vectors[,1])/sum(as.double(A_decomp$vectors[,1]))
pi <- wt*At[1,]/sum(wt*At[1,])
- # marginal mothers age
- Iom <- diag(1,ages, ages);
- ones <- t(rep(1,s))
+ # useful vectors and matrices
+ ones <- t(rep(1,s))
onesom <- t(rep(1,s*ages))
- piage <- kronecker(Iom,ones) %*% pi
-
- # momarray is an array with pit in each column
- momarray <- pi %*% matrix(1,1,ages)
- Iom = diag(1, ages)
- Is = diag(1, s)
- Isom = diag(1, s*ages)
- zsom = matrix(0, s*ages, s*ages)
- Z=Is;
- Z[1,1]=0;
- for(i in 1:ages){
- # imom = 1
- E <- Iom[,i] %*% t(Iom[i,]); # al cuadrado?
- momarray[,i] <- kronecker(E,Z) %*% momarray[,i]
- }
- # re-scale
- momarray <- momarray %*% MASS::ginv(diag(colSums(momarray)))
-
- # considering deaths
- phi = d = gd = ggd = m = gm = ggm = os = ys = nos = nys = oa = ya = coa = cya = matrix(0,ages*s*2,ages)
- phi[1,1] = 1
+ onesa <- t(rep(1,ages))
+ Iom <- diag(1, ages)
+ Is <- diag(1, s)
+ Isom <- diag(1, s*ages)
+ zsom <- matrix(0, s*ages, s*ages)
+
+ # momarray is an array with pi in each column
+ piage <- kronecker(Iom,ones) %*% pi
+ momarray <- pi %*% onesa
+
+ # considering deaths (no cumulated): reacreate block struct matrices
+ phi = d = gd = ggd = m = gm = ggm = os = ys = nos = nys = oa = ya = coa = cya = matrix(0,ages * s * 2, ages)
Mtt <- diag(as.numeric(onesom - onesom %*% Ut))
Utt <- rbind(cbind(Ut,zsom), cbind(Mtt,Isom)) %>% as.matrix()
ftt <- rbind(cbind(ft,zsom), cbind(zsom,zsom)) %>% as.matrix()
Gtt <- rbind(cbind(Gt,zsom), cbind(zsom,zsom)) %>% as.matrix()
sages <- 1:(ages*s)
- # no considering deaths
- # phi = d = gd = ggd = m = gm = ggm = os = ys = nos = nys = oa = ya = coa = cya = matrix(0,ages*s,ages)
- # phi[1,1] = 1
- # Utt = Ut %>% as.matrix()
- # ftt = ft %>% as.matrix()
- # Gtt = Gt %>% as.matrix()
+ # if parity: restriction to no initial mothers with 0 parity
+ if(parity){
+ Z=Is
+ Z[1,1]=0
+ for(i in 1:ages){
+ E <- Iom[,i] %*% t(Iom[i,])
+ momarray[,i] <- kronecker(E,Z) %*% momarray[,i]
+ }
+ # re-scale
+ momarray <- momarray %*% MASS::ginv(diag(colSums(momarray)))
+ # no 0 parity mothers: (momarray %*% piage)[seq(1,600,6)]
+ m[sages,1] = momarray %*% piage
+ }else{
+ m[sages,1] = pi
+ }
# focal´s trip
- m[sages,1] = momarray %*% piage;
+ phi[1,1] = 1
for(i in 1:(ages-1)){
phi[,i+1] = Gtt %*% phi[,i]
d[,i+1] = Utt %*% d[,i] + ftt %*% phi[,i]
@@ -145,7 +161,8 @@ kin_multi_stage <- function(U = NULL, f = NULL, D = NULL, H = NULL,
}
# get results
- kin_list <- list(d=d,gd=gd,ggd=ggd,m=m,gm=gm,ggm=ggm,os=os,ys=ys,
+ kin_list <- list(focal = phi,
+ d=d,gd=gd,ggd=ggd,m=m,gm=gm,ggm=ggm,os=os,ys=ys,
nos=nos,nys=nys,oa=oa,ya=ya,coa=coa,cya=cya)
# only selected kin
@@ -153,34 +170,32 @@ kin_multi_stage <- function(U = NULL, f = NULL, D = NULL, H = NULL,
kin_list <- kin_list %>% purrr::keep(names(.) %in% output_kin)
}
- # as data.frame
- kin <- purrr::map2(kin_list, names(kin_list),
- function(x,y){
- out <- as.data.frame(x)
- colnames(out) <- age
- out %>%
- dplyr::mutate(kin = y,
- age_kin = rep(sort(rep(age,s)),2),
- stage_kin = rep(rep(1:s,ages),2),
- alive = c(rep("living",s*ages),rep("dead",s*ages))
- # age_kin = sort(rep(age,s)),
- # stage_kin = rep(1:s,ages),
- # alive = c(rep("yes",s*ages))
- ) %>%
- tidyr::pivot_longer(c(-age_kin, -stage_kin, -kin, -alive), names_to = "age_focal", values_to = "count") %>%
- dplyr::mutate(age_focal = as.integer(age_focal)) %>%
- tidyr::pivot_wider(names_from = alive, values_from = count)
- }) %>%
+ # kin_full as data.frame
+ kin_full <- purrr::map2(kin_list, names(kin_list),
+ function(x,y){
+ # reassign deaths to Focal experienced age
+ x[(ages*s+1):(ages*s*2),1:(ages-1)] <- x[(ages*s+1):(ages*s*2),2:ages]
+ x[(ages*s+1):(ages*s*2),ages] <- 0
+ out <- as.data.frame(x)
+ colnames(out) <- age
+ out %>%
+ dplyr::mutate(kin = y,
+ age_kin = rep(sort(rep(age,s)),2),
+ stage_kin = rep(rep(1:s,ages),2),
+ alive = c(rep("living",s*ages),rep("dead",s*ages))) %>%
+ tidyr::pivot_longer(c(-age_kin, -stage_kin, -kin, -alive), names_to = "age_focal", values_to = "count") %>%
+ dplyr::mutate(age_focal = as.integer(age_focal)) %>%
+ tidyr::pivot_wider(names_from = alive, values_from = count)
+ }) %>%
purrr::reduce(rbind)
# results as list?
if(list_output) {
out <- kin_list
}else{
- out <- kin
+ out <- kin_full
}
+ # end
return(out)
}
-
-
diff --git a/R/kin_time_invariant.R b/R/kin_time_invariant.R
index d843c4e..06058b6 100644
--- a/R/kin_time_invariant.R
+++ b/R/kin_time_invariant.R
@@ -1,37 +1,40 @@
-#' Estimate kin counts in a time invariant framework
+#' Estimate kin counts in a time invariant framework for one-sex model (matrilineal/patrilineal)
-#' @description Implementation of Goodman-Keyfitz-Pullum equations adapted by Caswell (2019).
+#' @description Mtrix implementation of Goodman-Keyfitz-Pullum equations adapted by Caswell (2019).
-#' @param U numeric. A vector of survival probabilities with same length as ages.
+#' @param p numeric. A vector of survival probabilities with same length as ages.
#' @param f numeric. A vector of age-specific fertility rates with same length as ages.
#' @param birth_female numeric. Female portion at birth.
#' @param pi numeric. For using some specific non-stable age distribution of childbearing (same length as ages). Default `NULL`.
-#' @param output_kin character. kin to return. For example "m" for mother, "d" for daughter. See the `vignette` for all kin types.
+#' @param output_kin character. kin to return. For example "m" for mother, "d" for daughter. See `vignette` for all kin types.
#' @param list_output logical. Results as a list with `output_kin` elements, with focal´s age in columns and kin ages in rows (2 * ages, last chunk of ages for death experience). Default `FALSE`
#'
#' @return A data frame with focal´s age, related ages and type of kin
#' (for example `d` is daughter, `oa` is older aunts, etc.), alive and death. If `list_output = TRUE` then this is a list.
#' @export
-kin_time_invariant <- function(U = NULL, f = NULL,
+kin_time_invariant <- function(p = NULL, f = NULL,
birth_female = 1/2.04,
pi = NULL,
output_kin = NULL,
list_output = FALSE){
+ # global vars
+ .<-alive<-age_kin<-alive<-age_focal<-count<-NULL
+
# make matrix transition from vectors
- age = 0:(length(U)-1)
+ age = 0:(length(p)-1)
ages = length(age)
- Ut = Mt = zeros = Dcum = matrix(0, nrow=ages, ncol=ages)
- Ut[row(Ut)-1 == col(Ut)] <- U[-ages]
- Ut[ages, ages] = U[ages]
- diag(Mt) = 1 - U
+ Ut = Mt = zeros = matrix(0, nrow=ages, ncol=ages)
+ Ut[row(Ut)-1 == col(Ut)] <- p[-ages]
+ Ut[ages, ages] = p[ages]
+ diag(Mt) = 1 - p
Ut = rbind(cbind(Ut,zeros),
- cbind(Mt,Dcum))
+ cbind(Mt,zeros))
ft = matrix(0, nrow=ages*2, ncol=ages*2)
ft[1,1:ages] = f * birth_female
- # stable age distr
+ # stable age distribution in case no pi is given
if(is.null(pi)){
A = Ut[1:ages,1:ages] + ft[1:ages,1:ages]
A_decomp = eigen(A)
@@ -57,24 +60,20 @@ kin_time_invariant <- function(U = NULL, f = NULL,
ys[,i+1] = Ut %*% ys[,i] + ft %*% m[,i]
nys[,i+1] = Ut %*% nys[,i] + ft %*% ys[,i]
}
-
gm[1:ages,1] = m[1:ages,] %*% pi
for(i in 1:(ages-1)){
gm[,i+1] = Ut %*% gm[,i]
}
-
ggm[1:ages,1] = gm[1:ages,] %*% pi
for(i in 1:(ages-1)){
ggm[,i+1] = Ut %*% ggm[,i]
}
-
os[1:ages,1] = d[1:ages,] %*% pi
nos[1:ages,1] = gd[1:ages,] %*% pi
for(i in 1:(ages-1)){
os[,i+1] = Ut %*% os[,i]
nos[,i+1] = Ut %*% nos[,i] + ft %*% os[,i]
}
-
oa[1:ages,1] = os[1:ages,] %*% pi
ya[1:ages,1] = ys[1:ages,] %*% pi
coa[1:ages,1] = nos[1:ages,] %*% pi
@@ -95,9 +94,12 @@ kin_time_invariant <- function(U = NULL, f = NULL,
kin_list <- kin_list %>% purrr::keep(names(.) %in% output_kin)
}
- # as data.frame
+ # reshape as data.frame
kin <- purrr::map2(kin_list, names(kin_list),
function(x,y){
+ # reassign deaths to Focal experienced age
+ x[(ages+1):(ages*2),1:(ages-1)] <- x[(ages+1):(ages*2),2:ages]
+ x[(ages+1):(ages*2),ages] <- 0
out <- as.data.frame(x)
colnames(out) <- age
out %>%
@@ -111,13 +113,11 @@ kin_time_invariant <- function(U = NULL, f = NULL,
) %>%
purrr::reduce(rbind)
-
# results as list?
if(list_output) {
out <- kin_list
}else{
out <- kin
}
-
return(out)
}
diff --git a/R/kin_time_invariant_2sex.R b/R/kin_time_invariant_2sex.R
new file mode 100644
index 0000000..e9d7822
--- /dev/null
+++ b/R/kin_time_invariant_2sex.R
@@ -0,0 +1,165 @@
+#' Estimate kin counts in a time invariant framework for two-sex model.
+
+#' @description Two-sex matrix framework for kin count estimates.This produces kin counts grouped by kin, age and sex of
+#' each relatives at each Focal´s age. For example, male cousins from aunts and uncles from different sibling's parents
+#' are grouped in one male count of cousins.
+#' @details See Caswell (2022) for details on formulas.
+#' @param pf numeric. A vector of survival probabilities for females with same length as ages.
+#' @param ff numeric. A vector of age-specific fertility rates for females with same length as ages.
+#' @param pm numeric. A vector of survival probabilities for males with same length as ages.
+#' @param fm numeric. A vector of age-specific fertility rates for males with same length as ages.
+#' @param sex_focal character. "f" for female or "m" for male.
+#' @param birth_female numeric. Female portion at birth.
+#' @param pif numeric. For using some specific non-stable age distribution of childbearing for mothers (same length as ages). Default `NULL`.
+#' @param pim numeric. For using some specific non-stable age distribution of childbearing for fathers (same length as ages). Default `NULL`.
+#' @param output_kin character. kin to return, considering matrilineal names. For example "m" for parents, "d" for children, etc. See the `vignette` for all kin types.
+#' @param list_output logical. Results as a list with `output_kin` elements, with focal´s age in columns and kin ages in rows (2 * ages, last chunk of ages for death experience). Default `FALSE`
+#'
+#' @return A data frame with focal´s age, related ages and type of kin
+#' (for example `d` is children, `oa` is older aunts/uncles, etc.), sex, alive and death. If `list_output = TRUE` then this is a list.
+#' @export
+
+kin_time_invariant_2sex <- function(pf = NULL, pm = NULL,
+ ff = NULL, fm = NULL,
+ sex_focal = "f",
+ birth_female = 1/2.04,
+ pif = NULL, pim = NULL,
+ output_kin = NULL,
+ list_output = FALSE){
+
+ # global vars
+ .<-sex_kin<-alive<-count<-living<-dead<-age_kin<-age_focal<-cohort<-year<-total<-mean_age<-count_living<-sd_age<-count_dead<-mean_age_lost<-indicator<-value<-NULL
+
+ # same input length
+ if(!all(length(pf)==length(pm), length(pf)==length(ff), length(pf)==length(fm))) stop("Lengths of p's and f's should be the same")
+
+ # make matrix transition from vectors. Include death counts with matrix M
+ age = 0:(length(pf)-1)
+ ages = length(age)
+ agess = ages * 2
+ Uf = Um = Ff = Fm = Gt = zeros = matrix(0, nrow=ages, ncol=ages)
+ Uf[row(Uf)-1 == col(Uf)] <- pf[-ages]
+ Uf[ages, ages] = Uf[ages]
+ Um[row(Um)-1 == col(Um)] <- pm[-ages]
+ Um[ages, ages] = Um[ages]
+ Mm <- diag(1-pm)
+ Mf <- diag(1-pf)
+ Ut <- as.matrix(rbind(
+ cbind(Matrix::bdiag(Uf, Um), Matrix::bdiag(zeros, zeros)),
+ cbind(Matrix::bdiag(Mf, Mm), Matrix::bdiag(zeros, zeros))))
+ Ff[1,] = ff
+ Fm[1,] = fm
+ Ft <- Ft_star <- matrix(0, agess*2, agess*2)
+ Ft[1:agess,1:agess] <- rbind(cbind(birth_female * Ff, birth_female * Fm),
+ cbind((1-birth_female) * Ff, (1-birth_female) * Fm))
+
+ # mother and father do not reproduce independently to produce focal´s siblings. Assign to mother
+ Ft_star[1:agess,1:ages] <- rbind(birth_female * Ff, (1-birth_female) * Ff)
+
+ # parents age distribution under stable assumption in case no input
+ if(is.null(pim) | is.null(pif)){
+ A = Matrix::bdiag(Uf, Um) + Ft_star[1:agess,1:agess]
+ A_decomp = eigen(A)
+ lambda = as.double(A_decomp$values[1])
+ w = as.double(A_decomp$vectors[,1])/sum(as.double(A_decomp$vectors[,1]))
+ wf = w[1:ages]
+ wm = w[(ages+1):(2*ages)]
+ pif = wf * ff / sum(wf * ff)
+ pim = wm * fm / sum(wm * fm)
+ }
+
+ # initial count matrix (kin ages in rows and focal age in column)
+ phi = d = gd = ggd = m = gm = ggm = os = ys = nos = nys = oa = ya = coa = cya = matrix(0, agess*2, ages)
+
+ # locate focal at age 0 depending sex
+ sex_index <- ifelse(sex_focal == "f", 1, ages+1)
+ phi[sex_index, 1] <- 1
+
+ # G matrix moves focal by age
+ G <- matrix(0, nrow=ages, ncol=ages)
+ G[row(G)-1 == col(G)] <- 1
+ Gt <- matrix(0, agess*2, agess*2)
+ Gt[1:(agess), 1:(agess)] <- as.matrix(Matrix::bdiag(G, G))
+
+ # focal´s trip
+ # names of matrix count by kin refers to matrilineal as general reference
+ m[1:(agess),1] = c(pif, pim)
+ for(i in 1:(ages-1)){
+ # i = 1
+ phi[,i+1] = Gt %*% phi[,i]
+ d[,i+1] = Ut %*% d[,i] + Ft %*% phi[,i]
+ gd[,i+1] = Ut %*% gd[,i] + Ft %*% d[,i]
+ ggd[,i+1] = Ut %*% ggd[,i] + Ft %*% gd[,i]
+ m[,i+1] = Ut %*% m[,i]
+ ys[,i+1] = Ut %*% ys[,i] + Ft_star %*% m[,i]
+ nys[,i+1] = Ut %*% nys[,i] + Ft %*% ys[,i]
+ }
+
+ gm[1:(agess),1] = m[1:(agess),] %*% (pif + pim)
+ for(i in 1:(ages-1)){
+ gm[,i+1] = Ut %*% gm[,i]
+ }
+
+ ggm[1:(agess),1] = gm[1:(agess),] %*% (pif + pim)
+ for(i in 1:(ages-1)){
+ ggm[,i+1] = Ut %*% ggm[,i]
+ }
+
+ # atribuible to focal sex
+ pios = if(sex_focal == "f") pif else pim
+ os[1:(agess),1] = d[1:(agess),] %*% pif
+ nos[1:(agess),1] = gd[1:(agess),] %*% pif
+ for(i in 1:(ages-1)){
+ os[,i+1] = Ut %*% os[,i]
+ nos[,i+1] = Ut %*% nos[,i] + Ft %*% os[,i]
+ }
+
+ oa[1:(agess),1] = os[1:(agess),] %*% (pif + pim)
+ ya[1:(agess),1] = ys[1:(agess),] %*% (pif + pim)
+ coa[1:(agess),1] = nos[1:(agess),] %*% (pif + pim)
+ cya[1:(agess),1] = nys[1:(agess),] %*% (pif + pim)
+ for(i in 1:(ages-1)){
+ oa[,i+1] = Ut %*% oa[,i]
+ ya[,i+1] = Ut %*% ya[,i] + Ft_star %*% gm[,i]
+ coa[,i+1] = Ut %*% coa[,i] + Ft %*% oa[,i]
+ cya[,i+1] = Ut %*% cya[,i] + Ft %*% ya[,i]
+ }
+
+ # get results
+ kin_list <- list(d=d,gd=gd,ggd=ggd,m=m,gm=gm,ggm=ggm,os=os,ys=ys,
+ nos=nos,nys=nys,oa=oa,ya=ya,coa=coa,cya=cya)
+
+ # only selected kin
+ if(!is.null(output_kin)){
+ kin_list <- kin_list %>% purrr::keep(names(.) %in% output_kin)
+ }
+
+ # as data.frame
+ kin <- purrr::map2(kin_list, names(kin_list),
+ function(x,y){
+ # reassign deaths to Focal experienced age
+ x[(agess+1):(agess*2),1:(ages-1)] <- x[(agess+1):(agess*2),2:ages]
+ x[(agess+1):(agess*2),ages] <- 0
+ out <- as.data.frame(x)
+ colnames(out) <- age
+ out %>%
+ dplyr::mutate(kin = y,
+ age_kin = rep(age,4),
+ sex_kin = rep(c(rep("f",ages), rep("m",ages)),2),
+ alive = c(rep("living",2*ages), rep("dead",2*ages))) %>%
+ tidyr::pivot_longer(c(-age_kin, -kin, -sex_kin, -alive), names_to = "age_focal", values_to = "count") %>%
+ dplyr::mutate(age_focal = as.integer(age_focal)) %>%
+ tidyr::pivot_wider(names_from = alive, values_from = count)
+ }
+ ) %>%
+ purrr::reduce(rbind)
+
+ # results as list?
+ if(list_output) {
+ out <- kin_list
+ }else{
+ out <- kin
+ }
+
+ return(out)
+}
diff --git a/R/kin_time_invariant_2sex_cod.R b/R/kin_time_invariant_2sex_cod.R
new file mode 100644
index 0000000..e0d17d6
--- /dev/null
+++ b/R/kin_time_invariant_2sex_cod.R
@@ -0,0 +1,255 @@
+#' Estimate kin counts in a time invariant framework for two-sex model.
+
+#' @description Two-sex matrix framework for kin count and death estimates.This produces kin counts grouped by kin, age and sex of
+#' each relatives at each Focal´s age. For example, male cousins from aunts and uncles from different sibling's parents
+#' are grouped in one male count of cousins. This also produces kin deaths grouped by kin, age, sex of
+#' each relatives at each Focal´s age, and cause of death.
+#' @details See Caswell (2022) for details on formulas.
+#' @param pf numeric. A vector of survival probabilities for females with same length as ages.
+#' @param Hf numeric. A matrix with cause-specific hazards for females with rows as causes and columns as ages, being the name of each col the age.
+#' @param ff numeric. A vector of age-specific fertility rates for females with same length as ages.
+#' @param pm numeric. A vector of survival probabilities for males with same length as ages.
+#' @param fm numeric. A vector of age-specific fertility rates for males with same length as ages.
+#' @param Hm numeric. A matrix with cause-specific hazards for males with rows as causes and columns as ages, being the name of each col the age.
+#' @param sex_focal character. "f" for female or "m" for male.
+#' @param birth_female numeric. Female portion at birth.
+#' @param pif numeric. For using some specific non-stable age distribution of childbearing for mothers (same length as ages). Default `NULL`.
+#' @param pim numeric. For using some specific non-stable age distribution of childbearing for fathers (same length as ages). Default `NULL`.
+#' @param output_kin character. kin to return, considering matrilineal names. For example "m" for parents, "d" for children, etc. See the `vignette` for all kin types.
+#' @param list_output logical. Results as a list with `output_kin` elements, with focal´s age in columns and kin ages in rows (2 * ages, last chunk of ages for death experience). Default `FALSE`
+#'
+#' @return A data frame with focal´s age, related ages and type of kin
+#' (for example `d` is children, `oa` is older aunts/uncles, etc.), sex, alive and death. If `list_output = TRUE` then this is a list.
+#' @export
+
+# BEN: Added hazard matrices as inputs.
+# Assume that input of cause-specific mortality will be in terms of
+# matrices of cause-specific hazards for the two sexes (causes * ages).
+# Alternative: a matrix (causes * ages) containing the ratio mxi/mx.
+kin_time_invariant_2sex_cod <- function(pf = NULL,
+ pm = NULL,
+ ff = NULL,
+ fm = NULL,
+ Hf = NULL,
+ Hm = NULL,
+ sex_focal = "f",
+ birth_female = 1 / 2.04,
+ pif = NULL,
+ pim = NULL,
+ output_kin = NULL,
+ list_output = FALSE) {
+
+
+ # global vars
+ .<-sex_kin<-alive<-count<-living<-dead<-age_kin<-age_focal<-cohort<-year<-total<-mean_age<-count_living<-sd_age<-count_dead<-mean_age_lost<-indicator<-value<-NULL
+
+ # same input length
+
+ # BEN: Now we should also check the dimensions of the cause-specific hazard
+ # matrices.
+ if(!all(length(pf)==length(pm), length(pf)==length(ff), length(pf)==length(fm),
+ nrow(Hf)==nrow(Hm), ncol(Hf)==ncol(Hm), ncol(Hf)==length(pf))) stop("Number of age groups of p's, h's, and f's should match")
+
+ # make matrix transition from vectors. Include death counts with matrix M
+ age = 0:(length(pf)-1)
+ ages = length(age)
+ agess = ages * 2
+ Uf = Um = Ff = Fm = Gt = matrix(0, nrow=ages, ncol=ages)
+
+ # BEN: The zero matrix was deleted from line above and has
+ # to be made specific according to living/dead kin
+ # part of the block matrix Ut.
+ causes <- nrow(Hf) # number of causes of death
+ zeros_l <- matrix(0, nrow = ages, ncol = (causes*ages)) # zero matrix for living kin part
+ zeros_d = matrix(0, nrow = (causes*ages), ncol = (causes*ages)) # zero matrix for death kin part
+
+ Uf[row(Uf)-1 == col(Uf)] <- pf[-ages]
+
+ # BEN: What is the purpose of the following line? By default it is zero due to
+ # how the matrix is created
+ Uf[ages, ages] = Uf[ages]
+
+ Um[row(Um)-1 == col(Um)] <- pm[-ages]
+ Um[ages, ages] = Um[ages]
+
+ # BEN: Building of M, matrix of cause-specific prob. of dying.
+ # Hence, M = H D(h_tilde)^{-1} D(q)
+ # where h_tilde are the summed hazards for each age, and
+ # q = 1 - p
+ sum_hf <- t(rep(1, causes)) %*% Hf # h_tilde female
+ sum_hm <- t(rep(1, causes)) %*% Hm # h_tilde male
+ Mf <- Hf %*% solve(diag(c(sum_hf))) %*% diag(1-pf)
+ Mm <- Hm %*% solve(diag(c(sum_hm))) %*% diag(1-pm)
+ # Mm <- diag(1-pm)
+ # Mf <- diag(1-pf)
+
+ # BEN: In order to classify kin death by both cause and age at death,
+ # we need a mortality matrices M_hat of dimension
+ # ((causes*ages) * ages). See eq.12 in Caswell et al. (2024).
+ # Store columns of M as a list of vectors
+ Mf.cols <- lapply(1:ncol(Mf), function(j) return(Mf[,j]))
+ Mm.cols <- lapply(1:ncol(Mm), function(j) return(Mm[,j]))
+ # Create M_hat using the vectors as elements of the block diagonal
+ Ut <- as.matrix(rbind(
+ cbind(Matrix::bdiag(Uf, Um), Matrix::bdiag(zeros_l, zeros_l)),
+ cbind(Matrix::bdiag(Matrix::bdiag(Mf.cols), Matrix::bdiag(Mm.cols)), Matrix::bdiag(zeros_d, zeros_d))))
+
+ Ff[1,] = ff
+ Fm[1,] = fm
+
+ # BEN: Accounting for causes of death leads to have different dimensions
+ # in Ft and Ft_star.
+ Ft <- Ft_star <- matrix(0, (agess + agess*causes), (agess + agess*causes))
+ Ft[1:agess,1:agess] <- rbind(cbind(birth_female * Ff, birth_female * Fm),
+ cbind((1-birth_female) * Ff, (1-birth_female) * Fm))
+
+ # mother and father do not reproduce independently to produce focal´s siblings. Assign to mother
+ Ft_star[1:agess,1:ages] <- rbind(birth_female * Ff, (1-birth_female) * Ff)
+
+ # parents age distribution under stable assumption in case no input
+ if(is.null(pim) | is.null(pif)){
+ A = Matrix::bdiag(Uf, Um) + Ft_star[1:agess,1:agess]
+ A_decomp = eigen(A)
+ lambda = as.double(A_decomp$values[1])
+ w = as.double(A_decomp$vectors[,1])/sum(as.double(A_decomp$vectors[,1]))
+ wf = w[1:ages]
+ wm = w[(ages+1):(2*ages)]
+ pif = wf * ff / sum(wf * ff)
+ pim = wm * fm / sum(wm * fm)
+ }
+
+ # initial count matrix (kin ages in rows and focal age in column)
+ # BEN: Changed dimensions of lower part (dead kin) to account for death from causes.
+ phi = d = gd = ggd = m = gm = ggm = os = ys = nos = nys = oa = ya = coa = cya = matrix(0, (agess + agess*causes), ages)
+
+ # locate focal at age 0 depending sex
+ sex_index <- ifelse(sex_focal == "f", 1, ages+1)
+ phi[sex_index, 1] <- 1
+
+ # G matrix moves focal by age
+ G <- matrix(0, nrow=ages, ncol=ages)
+ G[row(G)-1 == col(G)] <- 1
+
+ # BEN: Changed dimensions
+ Gt <- matrix(0, (agess + agess*causes), (agess + agess*causes))
+
+ Gt[1:(agess), 1:(agess)] <- as.matrix(Matrix::bdiag(G, G))
+
+ # focal´s trip
+ # names of matrix count by kin refers to matrilineal as general reference
+ m[1:(agess),1] = c(pif, pim)
+ for(i in 1:(ages-1)){
+ # i = 1
+ phi[,i+1] = Gt %*% phi[,i]
+ d[,i+1] = Ut %*% d[,i] + Ft %*% phi[,i]
+ gd[,i+1] = Ut %*% gd[,i] + Ft %*% d[,i]
+ ggd[,i+1] = Ut %*% ggd[,i] + Ft %*% gd[,i]
+ m[,i+1] = Ut %*% m[,i]
+ ys[,i+1] = Ut %*% ys[,i] + Ft_star %*% m[,i]
+ nys[,i+1] = Ut %*% nys[,i] + Ft %*% ys[,i]
+ }
+
+ gm[1:(agess),1] = m[1:(agess),] %*% (pif + pim)
+ for(i in 1:(ages-1)){
+ gm[,i+1] = Ut %*% gm[,i]
+ }
+
+ ggm[1:(agess),1] = gm[1:(agess),] %*% (pif + pim)
+ for(i in 1:(ages-1)){
+ ggm[,i+1] = Ut %*% ggm[,i]
+ }
+
+ os[1:(agess),1] = d[1:(agess),] %*% pif
+ nos[1:(agess),1] = gd[1:(agess),] %*% pif
+ for(i in 1:(ages-1)){
+ os[,i+1] = Ut %*% os[,i]
+ nos[,i+1] = Ut %*% nos[,i] + Ft %*% os[,i]
+ }
+
+ oa[1:(agess),1] = os[1:(agess),] %*% (pif + pim)
+ ya[1:(agess),1] = ys[1:(agess),] %*% (pif + pim)
+ coa[1:(agess),1] = nos[1:(agess),] %*% (pif + pim)
+ cya[1:(agess),1] = nys[1:(agess),] %*% (pif + pim)
+ for(i in 1:(ages-1)){
+ oa[,i+1] = Ut %*% oa[,i]
+ ya[,i+1] = Ut %*% ya[,i] + Ft_star %*% gm[,i]
+ coa[,i+1] = Ut %*% coa[,i] + Ft %*% oa[,i]
+ cya[,i+1] = Ut %*% cya[,i] + Ft %*% ya[,i]
+ }
+
+ # get results
+ kin_list <- list(d=d,gd=gd,ggd=ggd,m=m,gm=gm,ggm=ggm,os=os,ys=ys,
+ nos=nos,nys=nys,oa=oa,ya=ya,coa=coa,cya=cya)
+
+ # only selected kin
+ if(!is.null(output_kin)){
+ kin_list <- kin_list %>% purrr::keep(names(.) %in% output_kin)
+ }
+
+ # as data.frame
+ kin <- purrr::map2(kin_list, names(kin_list),
+ function(x,y){
+
+ # BEN: Death take place in the same year and age!
+ # I adapted the code
+ # below such that it works with the new dimensions.
+
+ # reassign deaths to Focal experienced age
+ x[(agess+1):(agess + agess*causes),1:(ages-1)] <- x[(agess+1):(agess + agess*causes),2:ages]
+ x[(agess+1):(agess + agess*causes),ages] <- 0
+ out <- as.data.frame(x)
+ colnames(out) <- age
+ out %>%
+ # BEN: the matrices have different dimensions when
+ # we accounf for causes of death so what follows
+ # has been substantially changed.
+ dplyr::mutate(kin = y,
+ age_kin = c(rep(age,2), rep(rep(age,each=causes),2)),
+ sex_kin = c(rep(c("f", "m"),each=ages), rep(c("f", "m"),each=ages*causes)),
+ alive = c(rep("living",2*ages), rep(paste0("deadcause",1:causes),2*ages))) %>%
+ tidyr::pivot_longer(c(-age_kin, -kin, -sex_kin, -alive), names_to = "age_focal", values_to = "count") %>%
+ dplyr::mutate(age_focal = as.integer(age_focal)) %>%
+ tidyr::pivot_wider(names_from = alive, values_from = count)
+ }
+ ) %>%
+ purrr::reduce(rbind)
+
+ # results as list?
+ if(list_output) {
+ out <- kin_list
+ }else{
+ out <- kin
+ }
+
+ return(out)
+}
+
+## BEN: ========================================================================
+
+# Checks
+
+# No dead parent at birth: deadcausei=0 when age_focal==0
+# ff # fertility starts at age 13
+# kin |> filter(kin == "m", age_focal ==0, age_kin >= 12)
+#
+# # pi when age_focal==0 and age_kin when fx>0:
+# kin |> filter(kin == "m", age_kin >= 13, age_focal ==0)
+# pif[14:101]
+#
+# # mother dying from cause i at age x when focal is age==1 comes from nber of
+# # living mother age x when focal is age==1 multiplied by (1-pf[x])*(1/3)
+# kin |> filter(kin == "m", age_kin == 14, age_focal ==1)
+# 0.000246 * ((1-pf[15])*(1/3)) # mother
+# 0.0000486 * ((1-pm[15])*(1/3)) # father
+#
+# # Store to compare with kin_time_invariant_2sex.R
+# saveRDS(
+# kin,
+# here(
+# "checks",
+# "output_time_invariant_2sex.rds"
+# )
+# )
+
+
+## =============================================================================
diff --git a/R/kin_time_variant.R b/R/kin_time_variant.R
index bc962fc..71a5dfa 100644
--- a/R/kin_time_variant.R
+++ b/R/kin_time_variant.R
@@ -1,10 +1,10 @@
-#' Estimate kin counts in a time variant framework
+#' Estimate kin counts in a time variant framework (dynamic rates) for one-sex model (matrilineal/patrilineal)
-#' @description Implementation of time variant Goodman-Keyfitz-Pullum equations based on Caswell (2021).
-#'
-#' @param U numeric. A matrix of survival ratios with rows as ages and columns as years. Column names must be equal interval.
+#' @description Matrix implementation of time variant Goodman-Keyfitz-Pullum equations in a matrix framework.
+#' @details See Caswell (2021) for details on formulas.
+#' @param p numeric. A matrix of survival ratios with rows as ages and columns as years. Column names must be equal interval.
#' @param f numeric. A matrix of age-specific fertility rates with rows as ages and columns as years. Coincident with `U`.
-#' @param N numeric. A matrix of population with rows as ages and columns as years. Coincident with `U`.
+#' @param n numeric. A matrix of population with rows as ages and columns as years. Coincident with `U`.
#' @param pi numeric. A matrix with distribution of childbearing with rows as ages and columns as years. Coincident with `U`.
#' @param output_cohort integer. Year of birth of focal to return as output. Could be a vector. Should be within input data years range.
#' @param output_period integer. Year for which to return kinship structure. Could be a vector. Should be within input data years range.
@@ -12,82 +12,79 @@
#' @param birth_female numeric. Female portion at birth.
#' @param list_output logical. Results as a list with years elements (as a result of `output_cohort` and `output_period` combination), with a second list of `output_kin` elements, with focal´s age in columns and kin ages in rows (2 * ages, last chunk of ages for death experience). Default `FALSE`
-#' @return A data frame of population kinship structure, with focal's cohort, focal´s age, period year, type of relatives
+#' @return A data frame of population kinship structure, with Focal's cohort, focal´s age, period year, type of relatives
#' (for example `d` is daughter, `oa` is older aunts, etc.), living and death kin counts, and age of (living or time deceased) relatives. If `list_output = TRUE` then this is a list.
#' @export
-kin_time_variant <- function(U = NULL, f = NULL, N = NULL, pi = NULL,
+kin_time_variant <- function(p = NULL, f = NULL, pi = NULL, n = NULL,
output_cohort = NULL, output_period = NULL, output_kin = NULL,
birth_female = 1/2.04, list_output = FALSE){
+ # global vars
+ .<-living<-dead<-age_kin<-age_focal<-cohort<-year<-total<-mean_age<-count_living<-sd_age<-count_dead<-mean_age_lost<-indicator<-value<-NULL
+
# check input
- if(is.null(U) | is.null(f)) stop("You need values on U and/or f.")
+ if(is.null(p) | is.null(f)) stop("You need values on p and f.")
# diff years
- if(!any(as.integer(colnames(U)) == as.integer(colnames(f)))) stop("Data should be from same years.")
+ if(!any(as.integer(colnames(p)) == as.integer(colnames(f)))) stop("Make sure that p and f are matrices and have the same column names.")
# data should be from same interval years
- years_data <- as.integer(colnames(U))
- if(var(diff(years_data))!=0) stop("Data should be for same interval length years. Fill the gaps and run again")
+ years_data <- as.integer(colnames(p))
+ if(stats::var(diff(years_data))!=0) stop("The years given as column names in the p and f matrices must be equally spaced.")
# utils
- age <- 0:(nrow(U)-1)
+ age <- 0:(nrow(p)-1)
n_years_data <- length(years_data)
ages <- length(age)
om <- max(age)
zeros <- matrix(0, nrow=ages, ncol=ages)
# age distribution at childborn
+ pi_N_null_flag <- FALSE
if(is.null(pi)){
- if(is.null(N)){
+ if(is.null(n)){
# create pi and fill it during the loop
message("Stable assumption was made for calculating pi on each year because no input data.")
+ pi_N_null_flag <- TRUE
pi <- matrix(0, nrow=ages, ncol=n_years_data)
}else{
- pi <- rbind(t(t(N * f)/colSums(N * f)), matrix(0,ages,length(years_data)))
+ pi_N_null_flag <- FALSE
+ pi <- rbind(t(t(n * f)/colSums(n * f)), matrix(0,ages,length(years_data)))
}
}
- # get lists of matrix
- Ul = fl = list()
- for(t in 1:n_years_data){
- Ut = Mt = Dcum = matrix(0, nrow=ages, ncol=ages)
- Ut[row(Ut)-1 == col(Ut)] <- U[-ages,t]
- Ut[ages, ages]=U[ages,t]
- diag(Mt) = 1 - U[,t]
- Ul[[as.character(years_data[t])]] <- rbind(cbind(Ut,zeros),cbind(Mt,Dcum))
- ft = matrix(0, nrow=ages*2, ncol=ages*2)
- ft[1,1:ages] = f[,t] * birth_female
- fl[[as.character(years_data[t])]] <- ft
- }
- U <- Ul
- f <- fl
-
# loop over years (more performance here)
kin_all <- list()
pb <- progress::progress_bar$new(
format = "Running over input years [:bar] :percent",
- total = n_years_data, clear = FALSE, width = 60)
- for (iyear in 1:n_years_data){
- # print(iyear)
- Ut <- as.matrix(U[[iyear]])
- ft <- as.matrix(f[[iyear]])
- if(is.null(pi)){
+ total = n_years_data + 1, clear = FALSE, width = 50)
+ for (t in 1:n_years_data){
+ # build matrix
+ Ut = Mt = matrix(0, nrow=ages, ncol=ages)
+ Ut[row(Ut)-1 == col(Ut)] <- p[-ages,t]
+ Ut[ages, ages] = p[ages,t]
+ diag(Mt) = 1 - p[,t]
+ Ut = rbind(cbind(Ut,zeros),cbind(Mt,zeros))
+ ft = matrix(0, nrow=ages*2, ncol=ages*2)
+ ft[1,1:ages] = f[,t] * birth_female
+ if(pi_N_null_flag){
A <- Ut[1:ages,1:ages] + ft[1:ages,1:ages]
A_decomp = eigen(A)
w <- as.double(A_decomp$vectors[,1])/sum(as.double(A_decomp$vectors[,1]))
- pit <- pi[,iyear] <- w*A[1,]/sum(w*A[1,])
+ pit <- pi[,t] <- w*A[1,]/sum(w*A[1,])
}else{
- pit <- pi[,iyear]
+ pit <- pi[,t]
}
- if (iyear==1){
- U1 <- c(diag(Ut[-1,])[1:om],Ut[om,om])
+ # proj
+ if (t==1){
+ p1 <- c(diag(Ut[-1,])[1:om],Ut[om,om])
f1 <- ft[1,][1:ages]
pi1 <- pit[1:ages]
- kin_all[[1]] <- kin_time_invariant(U = U1, f = f1/birth_female, pi = pi1, birth_female = birth_female,
+ kin_all[[1]] <- kin_time_invariant(p = p1, f = f1/birth_female, pi = pi1, birth_female = birth_female,
list_output = TRUE)
}
- kin_all[[iyear+1]] <- timevarying_kin(Ut=Ut,ft=ft,pit=pit,ages,pkin=kin_all[[iyear]])
+ kin_all[[t+1]] <- timevarying_kin(Ut=Ut,ft=ft,pit=pit,ages,pkin=kin_all[[t]])
pb$tick()
}
@@ -110,23 +107,29 @@ kin_time_variant <- function(U = NULL, f = NULL, N = NULL, pi = NULL,
purrr::map(~ .[selected_kin_position])
# long format
- kin <- lapply(names(kin_list), function(Y){
+ message("Preparing output...")
+ kin <- lapply(names(kin_list), FUN = function(Y){
X <- kin_list[[Y]]
- X <- purrr::map2(X, names(X), function(x,y) as.data.frame(x) %>%
- dplyr::mutate(year = Y,
- kin=y,
- age_kin = rep(age,2),
- alive = c(rep("living",ages), rep("dead",ages)),
- .before=everything())) %>%
- dplyr::bind_rows() %>%
- stats::setNames(c("year","kin","age_kin","alive",as.character(age))) %>%
- tidyr::gather(age_focal, count,-age_kin, -kin, -year, -alive) %>%
- dplyr::mutate(age_focal = as.integer(age_focal),
- year = as.integer(year),
- cohort = year - age_focal) %>%
- dplyr::filter(age_focal %in% out_selected$age[out_selected$year==as.integer(Y)]) %>%
- tidyr::pivot_wider(names_from = alive, values_from = count)}) %>%
- dplyr::bind_rows()
+ X <- purrr::map2(X, names(X), function(x,y){
+ # reassign deaths to Focal experienced age
+ x[(ages+1):(ages*2),1:(ages-1)] <- x[(ages+1):(ages*2),2:ages]
+ x[(ages+1):(ages*2),ages] <- 0
+ x <- as.data.frame(x)
+ x$year <- Y
+ x$kin <- y
+ x$age_kin <- rep(age,2)
+ x$alive <- c(rep("living",ages), rep("dead",ages))
+ return(x)
+ }) %>%
+ data.table::rbindlist() %>%
+ stats::setNames(c(as.character(age), "year","kin","age_kin","alive")) %>%
+ data.table::melt(id.vars = c("year","kin","age_kin","alive"), variable.name = "age_focal", value.name = "count")
+ X$age_focal = as.integer(as.character(X$age_focal))
+ X$year = as.integer(X$year)
+ X$cohort = X$year - X$age_focal
+ X[X$age_focal %in% out_selected$age[out_selected$year==as.integer(Y)],] %>%
+ data.table::dcast(year + kin + age_kin + age_focal + cohort ~ alive, value.var = "count")
+ }) %>% data.table::rbindlist()
# results as list?
if(list_output) {
@@ -148,7 +151,8 @@ kin_time_variant <- function(U = NULL, f = NULL, N = NULL, pi = NULL,
#' @param pit numeric. A matrix with distribution of childbearing.
#' @param ages numeric.
#' @param pkin numeric. A list with kin count distribution in previous year.
-#
+#' @return A list of 14 types of kin matrices (kin age by Focal age) projected one time interval.
+#' @export
timevarying_kin<- function(Ut, ft, pit, ages, pkin){
# frequently used zero vector for initial condition
@@ -172,7 +176,8 @@ timevarying_kin<- function(Ut, ft, pit, ages, pkin){
coa[1:ages,1]= pkin[["nos"]][1:ages,] %*% pit[1:ages]
cya[1:ages,1]= pkin[["nys"]][1:ages,] %*% pit[1:ages]
- for (ix in 1:om){
+ # vers1
+ for(ix in 1:om){
d[,ix+1] = Ut %*% pkin[["d"]][,ix] + ft %*% I[,ix]
gd[,ix+1] = Ut %*% pkin[["gd"]][,ix] + ft %*% pkin[["d"]][,ix]
ggd[,ix+1] = Ut %*% pkin[["ggd"]][,ix] + ft %*% pkin[["gd"]][,ix]
@@ -195,10 +200,17 @@ timevarying_kin<- function(Ut, ft, pit, ages, pkin){
return(kin_list)
}
-#' defince apc combination to return
+#' APC combination to return
-#' @description defince apc to return.
-#'
+#' @description define APC combination to return in `kin` and `kin2sex`.
+#' @details Because returning all period and cohort data from a huge time-series would be hard memory consuming,
+#' this function is an auxiliary one to deal with selection from inputs `output_cohort` and `output_period`.
+#' @param output_cohort integer. A vector with selected calendar years.
+#' @param output_period integer. A vector with selected cohort years.
+#' @param age integer. A vector with ages from the kinship network to be filtered.
+#' @param years_data integer. A vector with years from the time-varying kinship network to be filtered.
+#' @return data.frame with years and ages to filter in `kin` and `kin_2sex` functions.
+#' @export
output_period_cohort_combination <- function(output_cohort = NULL, output_period = NULL, age = NULL, years_data = NULL){
# no specific
@@ -214,10 +226,11 @@ output_period_cohort_combination <- function(output_cohort = NULL, output_period
unlist(use.names = F))
}else{selected_cohorts_year_age <- c()}
- # period year combination
+ # period combination
if(!is.null(output_period)){selected_years_age <- expand.grid(age, output_period) %>% dplyr::rename(age=1,year=2)
}else{selected_years_age <- c()}
# end
return(dplyr::bind_rows(selected_years_age,selected_cohorts_year_age) %>% dplyr::distinct())
}
+
diff --git a/R/kin_time_variant_2sex.R b/R/kin_time_variant_2sex.R
new file mode 100644
index 0000000..e67db6f
--- /dev/null
+++ b/R/kin_time_variant_2sex.R
@@ -0,0 +1,254 @@
+#' Estimate kin counts in a time variant framework (dynamic rates) in a two-sex framework (Caswell, 2022)
+
+#' @description Two-sex matrix framework for kin count estimates with varying rates.
+#' This produces kin counts grouped by kin, age and sex of each relatives at each Focal´s age.
+#' For example, male cousins from aunts and uncles from different sibling's parents are grouped in one male count of cousins.
+#' @details See Caswell (2022) for details on formulas.
+#' @param pf numeric. A vector (atomic) or matrix with probabilities (or survival ratios, or transition between age class in a more general perspective) with rows as ages (and columns as years in case of matrix, being the name of each col the year).
+#' @param pm numeric. A vector (atomic) or matrix with probabilities (or survival ratios, or transition between age class in a more general perspective) with rows as ages (and columns as years in case of matrix, being the name of each col the year).
+#' @param ff numeric. Same as pf but for fertility rates.
+#' @param fm numeric. Same as pm but for fertility rates.
+#' @param sex_focal character. "f" for female or "m" for male.
+#' @param pif numeric. For using some specific age distribution of childbearing for mothers (same length as ages). Default `NULL`.
+#' @param pim numeric. For using some specific age distribution of childbearing for fathers (same length as ages). Default `NULL`.
+#' @param nf numeric. Same as pf but for population distribution (counts or `%`). Optional.
+#' @param nm numeric. Same as pm but for population distribution (counts or `%`). Optional.
+#' @param output_cohort integer. Vector of year cohorts for returning results. Should be within input data years range.
+#' @param output_period integer. Vector of period years for returning results. Should be within input data years range.
+#' @param output_kin character. kin types to return: "m" for mother, "d" for daughter,...
+#' @param birth_female numeric. Female portion at birth. This multiplies `f` argument. If `f` is already for female offspring, this needs to be set as 1.
+#' @param list_output logical. Results as a list with years elements (as a result of `output_cohort` and `output_period` combination), with a second list of `output_kin` elements, with focal´s age in columns and kin ages in rows (2 * ages, last chunk of ages for death experience). Default `FALSE`
+#' @return A data.frame with year, cohort, Focal´s age, related ages, sex and type of kin (for example `d` is daughter, `oa` is older aunts, etc.), including living and dead kin at that age and sex.
+#' @export
+
+kin_time_variant_2sex <- function(pf = NULL, pm = NULL,
+ ff = NULL, fm = NULL,
+ sex_focal = "f",
+ birth_female = 1/2.04,
+ pif = NULL, pim = NULL,
+ nf = NULL, nm = NULL,
+ output_cohort = NULL, output_period = NULL, output_kin = NULL,
+ list_output = FALSE){
+
+ # global vars
+ .<-living<-dead<-age_kin<-age_focal<-cohort<-year<-total<-mean_age<-count_living<-sd_age<-count_dead<-mean_age_lost<-indicator<-value<-NULL
+
+ # same input length
+ if(!all(dim(pf) == dim(pm), dim(pf) == dim(ff), dim(pf) == dim(fm))) stop("Dimension of P's and F's should be the same")
+
+ # data should be from same interval years
+ years_data <- as.integer(colnames(pf))
+ if(stats::var(diff(years_data))!=0) stop("Data should be for same interval length years. Fill the gaps and run again")
+
+ # utils
+ age <- 0:(nrow(pf)-1)
+ n_years_data <- length(years_data)
+ ages <- length(age)
+ agess <- ages*2
+ om <- max(age)
+ zeros <- matrix(0, nrow=ages, ncol=ages)
+
+ # age distribution at child born
+ Pif <- pif; no_Pif <- FALSE
+ Pim <- pim; no_Pim <- FALSE
+ if(is.null(pif)){
+ if(!is.null(nf)){
+ Pif <- t(t(nf * ff)/colSums(nf * ff))
+ }else{
+ Pif <- matrix(0, nrow=ages, ncol=n_years_data)
+ no_Pif <- TRUE
+ }
+ }
+ if(is.null(pim)){
+ if(!is.null(nm)){
+ Pim <- t(t(nm * fm)/colSums(nm * fm))
+ }else{
+ Pim <- matrix(0, nrow=ages, ncol=n_years_data)
+ no_Pim <- TRUE
+ }
+ }
+
+ # get lists of matrix
+ Ul = Fl = Fl_star = list()
+ kin_all <- list()
+ pb <- progress::progress_bar$new(
+ format = "Running over input years [:bar] :percent",
+ total = n_years_data + 1, clear = FALSE, width = 60)
+ for(t in 1:n_years_data){
+ # t = 1
+ Uf = Um = Fft = Fmt = Mm = Mf = Gt = zeros = matrix(0, nrow=ages, ncol=ages)
+ Uf[row(Uf)-1 == col(Uf)] <- pf[-ages,t]
+ Uf[ages, ages] = pf[ages,t]
+ Um[row(Um)-1 == col(Um)] <- pm[-ages,t]
+ Um[ages, ages] = pm[ages,t]
+ Mm <- diag(1-pm[,t])
+ Mf <- diag(1-pf[,t])
+ Ut <- as.matrix(rbind(
+ cbind(Matrix::bdiag(Uf, Um), Matrix::bdiag(zeros, zeros)),
+ cbind(Matrix::bdiag(Mf, Mm), Matrix::bdiag(zeros, zeros))))
+ Ul[[as.character(years_data[t])]] <- Ut
+ Fft[1,] = ff[,t]
+ Fmt[1,] = fm[,t]
+ Ft <- Ft_star <- matrix(0, agess*2, agess*2)
+ Ft[1:agess,1:agess] <- rbind(cbind(birth_female * Fft, birth_female * Fmt),
+ cbind((1-birth_female) * Fft, (1-birth_female) * Fmt))
+ Ft_star[1:agess,1:ages] <- rbind(birth_female * Fft, (1-birth_female) * Fft)
+ Fl[[as.character(years_data[t])]] <- Ft
+ Fl_star[[as.character(years_data[t])]] <- Ft_star
+ # parents age distribution under stable assumption in case no input
+ if(no_Pim | no_Pif){
+ A = Matrix::bdiag(Uf, Um) + Ft_star[1:agess,1:agess]
+ A_decomp = eigen(A)
+ lambda = as.double(A_decomp$values[1])
+ w = as.double(A_decomp$vectors[,1])/sum(as.double(A_decomp$vectors[,1]))
+ wf = w[1:ages]
+ wm = w[(ages+1):(2*ages)]
+ Pif[,t] = wf * ff[,t] / sum(wf * ff[,t])
+ Pim[,t] = wm * fm[,t] / sum(wm * fm[,t])
+ }
+
+ # project
+ Ut <- as.matrix(Ul[[t]])
+ Ft <- as.matrix(Fl[[t]])
+ Ft_star <- as.matrix(Fl_star[[t]])
+ pitf <- Pif[,t]
+ pitm <- Pim[,t]
+ pit <- c(pitf, pitm)
+ if (t==1){
+ p1f <- pf[,1]
+ p1m <- pm[,1]
+ f1f <- ff[,1]
+ f1m <- fm[,1]
+ pif1 <- Pif[,1]
+ pim1 <- Pim[,1]
+ kin_all[[1]] <- kin_time_invariant_2sex(pf = p1f, pm = p1m,
+ ff = f1f, fm = f1m,
+ sex_focal = sex_focal,
+ pif = pif1, pim = pim1,
+ birth_female = birth_female, list_output = TRUE)
+ }
+ kin_all[[t+1]] <- timevarying_kin_2sex(Ut=Ut, Ft=Ft, Ft_star=Ft_star, pit=pit, sex_focal, ages, pkin=kin_all[[t]])
+ pb$tick()
+ }
+
+ # filter years and kin that were selected
+ names(kin_all) <- as.character(years_data)
+
+ # combinations to return
+ out_selected <- output_period_cohort_combination(output_cohort, output_period, age = age, years_data = years_data)
+
+ possible_kin <- c("d","gd","ggd","m","gm","ggm","os","ys","nos","nys","oa","ya","coa","cya")
+ if(is.null(output_kin)){
+ selected_kin_position <- 1:length(possible_kin)
+ }else{
+ selected_kin_position <- which(possible_kin %in% output_kin)
+ }
+
+ # first filter
+ kin_list <- kin_all %>%
+ purrr::keep(names(.) %in% as.character(unique(out_selected$year))) %>%
+ purrr::map(~ .[selected_kin_position])
+ # long format
+ message(" Preparing output...")
+ kin <- lapply(names(kin_list), FUN = function(Y){
+ X <- kin_list[[Y]]
+ X <- purrr::map2(X, names(X), function(x,y){
+ # reassign deaths to Focal experienced age
+ x[(agess+1):(agess*2),1:(ages-1)] <- x[(agess+1):(agess*2),2:ages]
+ x[(agess+1):(agess*2),ages] <- 0
+ x <- data.table::as.data.table(x)
+ x$year <- Y
+ x$kin <- y
+ x$sex_kin <- rep(c(rep("f",ages), rep("m",ages)),2)
+ x$age_kin <- rep(age, 4)
+ x$alive <- c(rep("living",agess), rep("dead",agess))
+ return(x)
+ }) %>%
+ data.table::rbindlist() %>%
+ stats::setNames(c(as.character(age), "year","kin","sex_kin","age_kin","alive")) %>%
+ data.table::melt(id.vars = c("year","kin","sex_kin","age_kin","alive"), variable.name = "age_focal", value.name = "count")
+ X$age_focal = as.integer(as.character(X$age_focal))
+ X$year = as.integer(X$year)
+ X$cohort = X$year - X$age_focal
+ X <- X[X$age_focal %in% out_selected$age[out_selected$year==as.integer(Y)],]
+ X <- data.table::dcast(X, year + kin + sex_kin + age_kin + age_focal + cohort ~ alive, value.var = "count", fun.aggregate = sum)
+ }) %>% data.table::rbindlist()
+
+ # results as list?
+ if(list_output) {
+ out <- kin_list
+ }else{
+ out <- kin
+ }
+ return(out)
+}
+
+#' one time projection kin
+
+#' @description one time projection kin. internal function.
+#'
+#' @param Ut numeric. A matrix of survival probabilities (or ratios).
+#' @param Ft numeric. A matrix of age-specific fertility rates.
+#' @param Ft_star numeric. Ft but for female fertility.
+#' @param pit numeric. A matrix with distribution of childbearing.
+#' @param sex_focal character. "f" for female or "m" for male.
+#' @param ages numeric.
+#' @param pkin numeric. A list with kin count distribution in previous year.
+#' @return A list of 14 types of kin matrices (kin age by Focal age, blocked for two sex) projected one time interval.
+#' @export
+timevarying_kin_2sex<- function(Ut, Ft, Ft_star, pit, sex_focal, ages, pkin){
+
+ agess <- ages*2
+ om <- ages-1
+ pif <- pit[1:ages]
+ pim <- pit[(ages+1):agess]
+ phi = d = gd = ggd = m = gm = ggm = os = ys = nos = nys = oa = ya = coa = cya = matrix(0,agess*2,ages)
+ kin_list <- list(d=d,gd=gd,ggd=ggd,m=m,gm=gm,ggm=ggm,os=os,ys=ys,
+ nos=nos,nys=nys,oa=oa,ya=ya,coa=coa,cya=cya)
+
+ # G matrix moves focal by age
+ G <- matrix(0, nrow=ages, ncol=ages)
+ G[row(G)-1 == col(G)] <- 1
+ Gt <- matrix(0, agess*2, agess*2)
+ Gt[1:(agess), 1:(agess)] <- as.matrix(Matrix::bdiag(G, G))
+
+ # locate focal at age 0 depending sex
+ sex_index <- ifelse(sex_focal == "f", 1, ages+1)
+ phi[sex_index, 1] <- 1
+
+ # initial distribution
+ m[1:agess,1] = pit
+ gm[1:agess,1] = pkin[["m"]][1:agess,] %*% (pif + pim)
+ ggm[1:agess,1] = pkin[["gm"]][1:agess,] %*% (pif + pim)
+ oa[1:agess,1] = pkin[["os"]][1:agess,] %*% (pif + pim)
+ ya[1:agess,1] = pkin[["ys"]][1:agess,] %*% (pif + pim)
+ coa[1:agess,1] = pkin[["nos"]][1:agess,] %*% (pif + pim)
+ cya[1:agess,1] = pkin[["nys"]][1:agess,] %*% (pif + pim)
+ # atribuible to focal sex
+ pios = if(sex_focal == "f") pif else pim
+ os[1:agess,1] = pkin[["d"]][1:agess,] %*% pios
+ nos[1:agess,1] = pkin[["gd"]][1:ages,] %*% pios
+
+ for (ix in 1:om){
+ phi[,ix+1] = Gt %*% phi[, ix]
+ d[,ix+1] = Ut %*% pkin[["d"]][,ix] + Ft %*% phi[,ix]
+ gd[,ix+1] = Ut %*% pkin[["gd"]][,ix] + Ft %*% pkin[["d"]][,ix]
+ ggd[,ix+1] = Ut %*% pkin[["ggd"]][,ix] + Ft %*% pkin[["gd"]][,ix]
+ m[,ix+1] = Ut %*% pkin[["m"]][,ix]
+ gm[,ix+1] = Ut %*% pkin[["gm"]][,ix]
+ ggm[,ix+1] = Ut %*% pkin[["ggm"]][,ix]
+ os[,ix+1] = Ut %*% pkin[["os"]][,ix]
+ ys[,ix+1] = Ut %*% pkin[["ys"]][,ix] + Ft_star %*% pkin[["m"]][,ix]
+ nos[,ix+1] = Ut %*% pkin[["nos"]][,ix] + Ft %*% pkin[["os"]][,ix]
+ nys[,ix+1] = Ut %*% pkin[["nys"]][,ix] + Ft %*% pkin[["ys"]][,ix]
+ oa[,ix+1] = Ut %*% pkin[["oa"]][,ix]
+ ya[,ix+1] = Ut %*% pkin[["ya"]][,ix] + Ft_star %*% pkin[["gm"]][,ix]
+ coa[,ix+1] = Ut %*% pkin[["coa"]][,ix] + Ft %*% pkin[["oa"]][,ix]
+ cya[,ix+1] = Ut %*% pkin[["cya"]][,ix] + Ft %*% pkin[["ya"]][,ix]
+ }
+
+ kin_list <- list(d=d,gd=gd,ggd=ggd,m=m,gm=gm,ggm=ggm,os=os,ys=ys,
+ nos=nos,nys=nys,oa=oa,ya=ya,coa=coa,cya=cya)
+
+ return(kin_list)
+}
diff --git a/R/kin_time_variant_2sex_cod.R b/R/kin_time_variant_2sex_cod.R
new file mode 100644
index 0000000..562be20
--- /dev/null
+++ b/R/kin_time_variant_2sex_cod.R
@@ -0,0 +1,316 @@
+#' Estimate kin counts in a time variant framework (dynamic rates) in a two-sex framework (Caswell, 2022)
+
+#' @description Two-sex matrix framework for kin count estimates with varying rates.
+#' This produces kin counts grouped by kin, age and sex of each relatives at each Focal´s age.
+#' For example, male cousins from aunts and uncles from different sibling's parents are grouped in one male count of cousins. This also produces kin deaths grouped by kin, age, sex of
+#' each relatives at each Focal´s age, and cause of death.
+#' @details See Caswell (2022) for details on formulas.
+#' @param pf numeric. A vector (atomic) or matrix with probabilities (or survival ratios, or transition between age class in a more general perspective) with rows as ages (and columns as years in case of matrix, being the name of each col the year).
+#' @param pm numeric. A vector (atomic) or matrix with probabilities (or survival ratios, or transition between age class in a more general perspective) with rows as ages (and columns as years in case of matrix, being the name of each col the year).
+#' @param ff numeric. Same as pf but for fertility rates.
+#' @param fm numeric. Same as pm but for fertility rates.
+#' @param Hf numeric. A list where each list element (being the name of each list element the year) contains a matrix with cause-specific hazards for females with rows as causes and columns as ages, being the name of each col the age.
+#' @param Hm numeric. A list where each list element (being the name of each list element the year) contains a matrix with cause-specific hazards for males with rows as causes and columns as ages, being the name of each col the age.
+#' @param sex_focal character. "f" for female or "m" for male.
+#' @param pif numeric. For using some specific age distribution of childbearing for mothers (same length as ages). Default `NULL`.
+#' @param pim numeric. For using some specific age distribution of childbearing for fathers (same length as ages). Default `NULL`.
+#' @param nf numeric. Same as pf but for population distribution (counts or `%`). Optional.
+#' @param nm numeric. Same as pm but for population distribution (counts or `%`). Optional.
+#' @param output_cohort integer. Vector of year cohorts for returning results. Should be within input data years range.
+#' @param output_period integer. Vector of period years for returning results. Should be within input data years range.
+#' @param output_kin character. kin types to return: "m" for mother, "d" for daughter,...
+#' @param birth_female numeric. Female portion at birth. This multiplies `f` argument. If `f` is already for female offspring, this needs to be set as 1.
+#' @param list_output logical. Results as a list with years elements (as a result of `output_cohort` and `output_period` combination), with a second list of `output_kin` elements, with focal´s age in columns and kin ages in rows (2 * ages, last chunk of ages for death experience). Default `FALSE`
+#' @return A data.frame with year, cohort, Focal´s age, related ages, sex and type of kin (for example `d` is daughter, `oa` is older aunts, etc.), including living and dead kin at that age and sex.
+#' @export
+
+# BEN: Added hazard matrices as inputs.
+# Assume that input of cause-specific mortality will be in terms of
+# matrices of cause-specific hazards for the two sexes (causes * ages).
+# Alternative: a matrix (causes * ages) containing the ratio mxi/mx.
+kin_time_variant_2sex_cod <- function(pf = NULL, pm = NULL,
+ ff = NULL, fm = NULL,
+ Hf = NULL, Hm = NULL,
+ sex_focal = "f",
+ birth_female = 1/2.04,
+ pif = NULL, pim = NULL,
+ nf = NULL, nm = NULL,
+ output_cohort = NULL, output_period = NULL, output_kin = NULL,
+ list_output = FALSE){
+
+ # global vars
+ .<-living<-dead<-age_kin<-age_focal<-cohort<-year<-total<-mean_age<-count_living<-sd_age<-count_dead<-mean_age_lost<-indicator<-value<-NULL
+
+ # same input length
+
+ # BEN: Now we should also check the dimensions of the cause-specific hazard
+ # matrices.
+ if(!all(dim(pf) == dim(pm), dim(pf) == dim(ff), dim(pf) == dim(fm),
+ nrow(Hf)==nrow(Hm), ncol(Hf)==ncol(Hm), ncol(Hf)==nrow(pf),
+ length(Hf)==length(Hm), length(Hm)==ncol(pf))) stop("Dimension of P's, F's, and H's should match")
+
+ # data should be from same interval years
+ years_data <- as.integer(colnames(pf))
+ if(stats::var(diff(years_data))!=0) stop("Data should be for same interval length years. Fill the gaps and run again")
+
+ # utils
+ age <- 0:(nrow(pf)-1)
+ n_years_data <- length(years_data)
+ ages <- length(age)
+ agess <- ages*2
+ om <- max(age)
+
+ # BEN: The zero matrix was deleted from line above and has
+ # to be made specific according to living/dead kin
+ # part of the block matrix Ut.
+ causes <- nrow(Hf[[1]]) # number of causes of death
+ zeros_l <- matrix(0, nrow = ages, ncol = (causes*ages)) # zero matrix for living kin part
+ zeros_d = matrix(0, nrow = (causes*ages), ncol = (causes*ages)) # zero matrix for death kin part
+
+ # age distribution at child born
+ Pif <- pif; no_Pif <- FALSE
+ Pim <- pim; no_Pim <- FALSE
+ if(is.null(pif)){
+ if(!is.null(nf)){
+ Pif <- t(t(nf * ff)/colSums(nf * ff))
+ }else{
+ Pif <- matrix(0, nrow=ages, ncol=n_years_data)
+ no_Pif <- TRUE
+ }
+ }
+ if(is.null(pim)){
+ if(!is.null(nm)){
+ Pim <- t(t(nm * fm)/colSums(nm * fm))
+ }else{
+ Pim <- matrix(0, nrow=ages, ncol=n_years_data)
+ no_Pim <- TRUE
+ }
+ }
+
+ # get lists of matrix
+ Ul = Fl = Fl_star = list()
+ kin_all <- list()
+ pb <- progress::progress_bar$new(
+ format = "Running over input years [:bar] :percent",
+ total = n_years_data + 1, clear = FALSE, width = 60)
+
+ # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ # BEN: First load function at the end of script
+ # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ for(t in 1:n_years_data){
+ # t = 1
+ Uf = Um = Fft = Fmt = Mm = Mf = Gt = matrix(0, nrow=ages, ncol=ages)
+ Uf[row(Uf)-1 == col(Uf)] <- pf[-ages,t]
+ Uf[ages, ages] = pf[ages,t]
+ Um[row(Um)-1 == col(Um)] <- pm[-ages,t]
+ Um[ages, ages] = pm[ages,t]
+
+ # BEN: Building of M, matrix of cause-specific prob. of dying.
+ # Hence, M = H D(h_tilde)^{-1} D(q)
+ # where h_tilde are the summed hazards for each age, and
+ # q = 1 - p
+ sum_hf <- t(rep(1, causes)) %*% Hf[[t]] # h_tilde female
+ sum_hm <- t(rep(1, causes)) %*% Hm[[t]] # h_tilde male
+ Mf <- Hf[[t]] %*% solve(diag(c(sum_hf))) %*% diag(1-pf[,t])
+ Mm <- Hm[[t]] %*% solve(diag(c(sum_hm))) %*% diag(1-pm[,t])
+ # Mm <- diag(1-pm[,t])
+ # Mf <- diag(1-pf[,t])
+
+ # BEN: In order to classify kin death by both cause and age at death,
+ # we need a mortality matrices M_hat of dimension
+ # ((causes*ages) * ages). See eq.12 in Caswell et al. (2024).
+ # Store columns of M as a list of vectors
+ Mf.cols <- lapply(1:ncol(Mf), function(j) return(Mf[,j]))
+ Mm.cols <- lapply(1:ncol(Mm), function(j) return(Mm[,j]))
+ # Create M_hat using the vectors as elements of the block diagonal
+ Ut <- as.matrix(rbind(
+ cbind(Matrix::bdiag(Uf, Um), Matrix::bdiag(zeros_l, zeros_l)),
+ cbind(Matrix::bdiag(Matrix::bdiag(Mf.cols), Matrix::bdiag(Mm.cols)), Matrix::bdiag(zeros_d, zeros_d))))
+
+ Ul[[as.character(years_data[t])]] <- Ut
+ Fft[1,] = ff[,t]
+ Fmt[1,] = fm[,t]
+
+ # BEN: Accounting for causes of death leads to have different dimensions
+ # in Ft and Ft_star.
+ Ft <- Ft_star <- matrix(0, (agess + agess*causes), (agess + agess*causes))
+
+ Ft[1:agess,1:agess] <- rbind(cbind(birth_female * Fft, birth_female * Fmt),
+ cbind((1-birth_female) * Fft, (1-birth_female) * Fmt))
+ Ft_star[1:agess,1:ages] <- rbind(birth_female * Fft, (1-birth_female) * Fft)
+ Fl[[as.character(years_data[t])]] <- Ft
+ Fl_star[[as.character(years_data[t])]] <- Ft_star
+ # parents age distribution under stable assumption in case no input
+ if(no_Pim | no_Pif){
+ A = Matrix::bdiag(Uf, Um) + Ft_star[1:agess,1:agess]
+ A_decomp = eigen(A)
+ lambda = as.double(A_decomp$values[1])
+ w = as.double(A_decomp$vectors[,1])/sum(as.double(A_decomp$vectors[,1]))
+ wf = w[1:ages]
+ wm = w[(ages+1):(2*ages)]
+ Pif[,t] = wf * ff[,t] / sum(wf * ff[,t])
+ Pim[,t] = wm * fm[,t] / sum(wm * fm[,t])
+ }
+
+ # project
+ Ut <- as.matrix(Ul[[t]])
+ Ft <- as.matrix(Fl[[t]])
+ Ft_star <- as.matrix(Fl_star[[t]])
+ pitf <- Pif[,t]
+ pitm <- Pim[,t]
+ pit <- c(pitf, pitm)
+ if (t==1){
+ p1f <- pf[,1]
+ p1m <- pm[,1]
+ f1f <- ff[,1]
+ f1m <- fm[,1]
+ pif1 <- Pif[,1]
+ pim1 <- Pim[,1]
+
+ # BEN: Add Hf and Hm
+ H1f <- Hf[[1]]
+ H1m <- Hm[[1]]
+
+ # BEN: cod version !!!
+ kin_all[[1]] <- kin_time_invariant_2sex_cod(pf = p1f, pm = p1m,
+ ff = f1f, fm = f1m,
+ pif = pif1, pim = pim1,
+ Hf = H1f, Hm = H1m,
+ birth_female = birth_female, list_output = TRUE)
+ }
+ kin_all[[t+1]] <- timevarying_kin_2sex_cod(Ut=Ut, Ft=Ft, Ft_star=Ft_star, pit=pit, sex_focal, ages, pkin=kin_all[[t]])
+ pb$tick()
+ }
+
+ # filter years and kin that were selected
+ names(kin_all) <- as.character(years_data)
+
+ # combinations to return
+ out_selected <- output_period_cohort_combination(output_cohort, output_period, age = age, years_data = years_data)
+
+ possible_kin <- c("d","gd","ggd","m","gm","ggm","os","ys","nos","nys","oa","ya","coa","cya")
+ if(is.null(output_kin)){
+ selected_kin_position <- 1:length(possible_kin)
+ }else{
+ selected_kin_position <- which(possible_kin %in% output_kin)
+ }
+
+ # first filter
+ kin_list <- kin_all %>%
+ purrr::keep(names(.) %in% as.character(unique(out_selected$year))) %>%
+ purrr::map(~ .[selected_kin_position])
+ # long format
+ message("Preparing output...")
+ kin <- lapply(names(kin_list), FUN = function(Y){
+ X <- kin_list[[Y]]
+ X <- purrr::map2(X, names(X), function(x,y){
+ # reassign deaths to Focal experienced age
+ x[(agess+1):(agess + agess*causes),1:(ages-1)] <- x[(agess+1):(agess + agess*causes),2:ages]
+ x[(agess+1):(agess + agess*causes),ages] <- 0
+ x <- data.table::as.data.table(x)
+ x$year <- Y
+ x$kin <- y
+ x$sex_kin <- c(rep(c("f", "m"),each=ages), rep(c("f", "m"),each=ages*causes))
+ x$age_kin <- c(rep(age,2), rep(rep(age,each=causes),2))
+ x$alive <- c(rep("living",2*ages), rep(paste0("deadcause",1:causes),2*ages))
+ return(x)
+ }) %>%
+ data.table::rbindlist() %>%
+ stats::setNames(c(as.character(age), "year","kin","sex_kin","age_kin","alive")) %>%
+ data.table::melt(id.vars = c("year","kin","sex_kin","age_kin","alive"), variable.name = "age_focal", value.name = "count")
+ X$age_focal = as.integer(as.character(X$age_focal))
+ X$year = as.integer(X$year)
+ X$cohort = X$year - X$age_focal
+ X <- X[X$age_focal %in% out_selected$age[out_selected$year==as.integer(Y)],]
+ X <- data.table::dcast(X, year + kin + sex_kin + age_kin + age_focal + cohort ~ alive, value.var = "count", fun.aggregate = sum)
+ }) %>% data.table::rbindlist()
+
+ # results as list?
+ if(list_output) {
+ out <- kin_list
+ }else{
+ out <- kin
+ }
+ return(out)
+}
+
+#' one time projection kin
+
+#' @description one time projection kin. internal function.
+#'
+#' @param Ut numeric. A matrix of survival probabilities (or ratios).
+#' @param Ft numeric. A matrix of age-specific fertility rates.
+#' @param Ft_star numeric. Ft but for female fertility.
+#' @param pit numeric. A matrix with distribution of childbearing.
+#' @param sex_focal character. "f" for female or "m" for male.
+#' @param ages numeric.
+#' @param pkin numeric. A list with kin count distribution in previous year.
+#' @return A list of 14 types of kin matrices (kin age by Focal age, blocked for two sex) projected one time interval.
+#' @export
+timevarying_kin_2sex_cod<- function(Ut, Ft, Ft_star, pit, sex_focal, ages, pkin){
+
+ agess <- ages*2
+ om <- ages-1
+ pif <- pit[1:ages]
+ pim <- pit[(ages+1):agess]
+
+ # BEN : Add the number of CoD
+ causes <- nrow(Hf[[1]])
+
+ # BEN: Changed dimensions of lower part (dead kin) to account for death from causes.
+ phi = d = gd = ggd = m = gm = ggm = os = ys = nos = nys = oa = ya = coa = cya = matrix(0, (agess + agess*causes), ages)
+
+ kin_list <- list(d=d,gd=gd,ggd=ggd,m=m,gm=gm,ggm=ggm,os=os,ys=ys,
+ nos=nos,nys=nys,oa=oa,ya=ya,coa=coa,cya=cya)
+
+ # G matrix moves focal by age
+ G <- matrix(0, nrow=ages, ncol=ages)
+ G[row(G)-1 == col(G)] <- 1
+
+ # BEN: Changed dimensions
+ Gt <- matrix(0, (agess + agess*causes), (agess + agess*causes))
+
+ Gt[1:(agess), 1:(agess)] <- as.matrix(Matrix::bdiag(G, G))
+
+ # locate focal at age 0 depending sex
+ sex_index <- ifelse(sex_focal == "f", 1, ages+1)
+ phi[sex_index, 1] <- 1
+
+ # BEN: NOT SURE ABOUT WHAT IS HAPPENING BELOW
+ # Rows are multiplied by the sum of the pi?
+
+ # initial distribution
+ m[1:agess,1] = pit
+ gm[1:agess,1] = pkin[["m"]][1:agess,] %*% (pif + pim)
+ ggm[1:agess,1] = pkin[["gm"]][1:agess,] %*% (pif + pim)
+ os[1:agess,1] = pkin[["d"]][1:agess,] %*% pif
+ nos[1:agess,1] = pkin[["gd"]][1:ages,] %*% pif
+ oa[1:agess,1] = pkin[["os"]][1:agess,] %*% (pif + pim)
+ ya[1:agess,1] = pkin[["ys"]][1:agess,] %*% (pif + pim)
+ coa[1:agess,1] = pkin[["nos"]][1:agess,] %*% (pif + pim)
+ cya[1:agess,1] = pkin[["nys"]][1:agess,] %*% (pif + pim)
+
+ for (ix in 1:om){
+ phi[,ix+1] = Gt %*% phi[, ix]
+ d[,ix+1] = Ut %*% pkin[["d"]][,ix] + Ft %*% phi[,ix]
+ gd[,ix+1] = Ut %*% pkin[["gd"]][,ix] + Ft %*% pkin[["d"]][,ix]
+ ggd[,ix+1] = Ut %*% pkin[["ggd"]][,ix] + Ft %*% pkin[["gd"]][,ix]
+ m[,ix+1] = Ut %*% pkin[["m"]][,ix]
+ gm[,ix+1] = Ut %*% pkin[["gm"]][,ix]
+ ggm[,ix+1] = Ut %*% pkin[["ggm"]][,ix]
+ os[,ix+1] = Ut %*% pkin[["os"]][,ix]
+ ys[,ix+1] = Ut %*% pkin[["ys"]][,ix] + Ft_star %*% pkin[["m"]][,ix]
+ nos[,ix+1] = Ut %*% pkin[["nos"]][,ix] + Ft %*% pkin[["os"]][,ix]
+ nys[,ix+1] = Ut %*% pkin[["nys"]][,ix] + Ft %*% pkin[["ys"]][,ix]
+ oa[,ix+1] = Ut %*% pkin[["oa"]][,ix]
+ ya[,ix+1] = Ut %*% pkin[["ya"]][,ix] + Ft_star %*% pkin[["gm"]][,ix]
+ coa[,ix+1] = Ut %*% pkin[["coa"]][,ix] + Ft %*% pkin[["oa"]][,ix]
+ cya[,ix+1] = Ut %*% pkin[["cya"]][,ix] + Ft %*% pkin[["ya"]][,ix]
+ }
+
+ kin_list <- list(d=d,gd=gd,ggd=ggd,m=m,gm=gm,ggm=ggm,os=os,ys=ys,
+ nos=nos,nys=nys,oa=oa,ya=ya,coa=coa,cya=cya)
+
+ return(kin_list)
+}
diff --git a/R/plot_diagramm.R b/R/plot_diagramm.R
index 2497186..6bffa8c 100644
--- a/R/plot_diagramm.R
+++ b/R/plot_diagramm.R
@@ -1,90 +1,77 @@
#' plot a Kin diagram (network)
-#' @description Given estimation of kin counts from `kins` function, draw a network diagramm.
-#' @param kin_total data.frame. With columns `kin` with type and `count` with some measeure.
-#' @param rounding numeric. Estimation could have a lot of decimals. Rounding will make looks more clear the diagramm.
-#' @return A plot
+#' @description Draws a Keyfitz-style kinship diagram given a kinship object created by the `kin` function. Displays expected kin counts for a Focal aged 'a'.
+#' @param kin_total data.frame. values in column `kin` define the relative type - see `demokin_codes()`. Values in column `count` are the expected number of relatives.
+#' @param rounding numeric. Number of decimals to show in diagram.
+#' @return A Keyfitz-style kinship plot.
#' @export
-plot_diagram <- function(kin_total, rounding = 3){
-
- vertices <- data.frame(
- nodes = c("ggd", "gd", "d", "Focal", "m", "gm", "ggm", "oa", "coa", "os", "nos", "ya", "cya", "ys", "nys")
- , x = c(1, 1, 1, 1, 1, 1, 1, 0, -1, 0, -1, 2, 3, 2, 3)
- , y = c(0, 1, 2, 3, 4, 5, 6, 4, 3, 3, 2, 4, 3, 3, 2)
- )
-
- d <- data.frame(
- from = c("ggd", "gd", "d", "Focal", "m", "gm", "gm", "oa", "m", "os", "gm", "ya", "m", "ys")
- , to = c("gd", "d", "Focal", "m", "gm", "ggm", "oa", "coa", "os", "nos", "ya", "cya", "ys", "nys")
- )
-
- # Add values
- lookup <- c(with(kin_total, paste0(kin, " \n", round(count, rounding))), "Focal")
- names(lookup) <- c(kin_total$kin, "Focal")
-
- vertices$nodes <- lookup[vertices$nodes]
- d$from <- lookup[d$from]
- d$to <- lookup[d$to]
-
- # Plot
-
- b <- igraph::graph_from_data_frame(vertices = vertices, d= d, directed = FALSE)
-
- plot(
- b
- , vertex.size = 30
- , curved = 1
- , vertex.color = "#FFF1E2"
- , vertex.shape = "circle"
- , vertex.label.cex = 0.8
- , vertex.label.color = "black"
- , label.degree = -pi/2
- , edge.width = 2
- , edge.color = "black"
- )
-
-}
-
-# old function
-
-# plot_diagram <- function(kin_total, rounding = 3){
-# # https://cran.r-project.org/web/packages/DiagrammeR/vignettes/graphviz-mermaid.html
-# # https://color.hailpixel.com/#D9E9BE,BF62CB,94C2DB,79D297,CDA76A,C8695B
-#
-# kin_total <- kin_total %>% mutate(count = round(count,digits = rounding))
-#
-# DiagrammeR::mermaid(
-# paste0("graph TD
-#
-# GGM(ggm:
", kin_total$count[kin_total$kin=="ggm"] ,")
-# GGM ==> GM(gm:
", kin_total$count[kin_total$kin=="gm"] ,")
-# GM --> AOM(oa:
", kin_total$count[kin_total$kin=="oa"] ,")
-# GM ==> M(m:
", kin_total$count[kin_total$kin=="m"] ,")
-# GM --> AYM(ya:
", kin_total$count[kin_total$kin=="ya"] ,")
-# AOM --> CAOM(coa:
", kin_total$count[kin_total$kin=="coa"] ,")
-# M --> OS(os:
", kin_total$count[kin_total$kin=="os"] ,")
-# M ==> E((Ego))
-# M --> YS(ys:
", kin_total$count[kin_total$kin=="ys"] ,")
-# AYM --> CAYM(cya:
", kin_total$count[kin_total$kin=="cya"] ,")
-# OS --> NOS(nos:
", kin_total$count[kin_total$kin=="nos"] ,")
-# E ==> D(d:
", kin_total$count[kin_total$kin=="d"] ,")
-# YS --> NYS(nys:
", kin_total$count[kin_total$kin=="nys"] ,")
-# D ==> GD(gd:
", kin_total$count[kin_total$kin=="gd"] ,")
-# style GGM fill:#a1f590, stroke:#333, stroke-width:2px;
-# style GM fill:#a1f590, stroke:#333, stroke-width:2px, text-align: center;
-# style M fill:#a1f590, stroke:#333, stroke-width:2px, text-align: center
-# style D fill:#a1f590, stroke:#333, stroke-width:2px, text-align: center
-# style YS fill:#a1f590, stroke:#333, stroke-width:2px, text-align: center
-# style OS fill:#a1f590, stroke:#333, stroke-width:2px, text-align: center
-# style CAOM fill:#f1f0f5, stroke:#333, stroke-width:2px, text-align: center
-# style AYM fill:#f1f0f5, stroke:#333, stroke-width:2px, text-align: center
-# style AOM fill:#f1f0f5, stroke:#333, stroke-width:2px, text-align: center
-# style CAYM fill:#f1f0f5, stroke:#333, stroke-width:2px, text-align: center
-# style NOS fill:#f1f0f5, stroke:#333, stroke-width:2px, text-align: center
-# style NYS fill:#f1f0f5, stroke:#333, stroke-width:2px, text-align: center
-# style E fill:#FFF, stroke:#333, stroke-width:4px, text-align: center
-# style D fill:#a1f590, stroke:#333, stroke-width:2px, text-align: center
-# style GD fill:#a1f590, stroke:#333, stroke-width:2px, text-align: center"))
-# }
-
+plot_diagram <-
+ function (kin_total, rounding = 3) {
+ rels <- c("ggd", "gd", "d", "Focal", "m", "gm", "ggm", "oa", "coa", "os", "nos", "ya", "cya", "ys", "nys")
+ # check all types are in
+ if(!any(unique(kin_total$kin) %in% rels) | any(c("s", "c", "a", "n") %in% unique(kin_total$kin))) stop("You need all specific types. If some are missed or grouped, for example old and younger sisters in 's', this will fail.")
+ vertices <- data.frame(
+ nodes = rels
+ , x = c(1, 1, 1, 1, 1, 1, 1, 0, -1, 0, -1, 2, 3, 2, 3)
+ , y = c(0, 1, 2, 3, 4, 5, 6, 4, 3, 3, 2, 4, 3, 3, 2)
+ )
+ d <- data.frame(from = c("ggd", "gd", "d", "Focal", "m",
+ "gm", "gm", "oa", "m", "os", "gm", "ya", "m", "ys"),
+ to = c("gd", "d", "Focal", "m", "gm", "ggm", "oa", "coa",
+ "os", "nos", "ya", "cya", "ys", "nys"))
+ lookup <- c(with(kin_total, paste0(kin, " \n", round(count, rounding))), "Focal")
+ names(lookup) <- c(kin_total$kin, "Focal")
+ vertices$nodes <- lookup[vertices$nodes]
+ d$from <- lookup[d$from]
+ d$to <- lookup[d$to]
+ # to show full relative names
+ relatives <- c("Cousins from older aunt", "Cousins from younger aunt",
+ "Daughter", "Grand-daughter", "Great-grand-daughter",
+ "Great-grandmother", "Grandmother", "Mother", "Nieces from older sister",
+ "Nieces from younger sister", "Aunt older than mother",
+ "Aunt younger than mother", "Older sister", "Younger sister", "")
+ names(relatives) <- c("coa", "cya", "d", "gd", "ggd",
+ "ggm", "gm", "m", "nos", "nys", "oa", "ya", "os",
+ "ys", "Focal")
+ labs <- relatives[rels]
+ # Plot
+ b <- igraph::graph_from_data_frame(vertices = vertices, d= d, directed = FALSE)
+ b_auto_layout <- igraph::layout.auto(b)
+ b_auto_layout_scaled <- igraph::norm_coords(b_auto_layout, ymin=-1, ymax=1, xmin=-1, xmax=1)
+ plot(
+ b
+ , vertex.size = 70
+ , curved = 1
+ , vertex.color = "#FFF1E2"
+ , vertex.shape = "circle"
+ , vertex.label.cex = 0.8
+ , vertex.label.color = "black"
+ , edge.width = 2
+ , layout = b_auto_layout_scaled * 3
+ , rescale = FALSE
+ , xlim = c(-3.3,3.3)
+ , ylim = c(-3.1,3.1)
+ )
+ # Add relative names
+ # Thanks to Egor Kotov for this tip!
+ plot(
+ b
+ , vertex.size = 70
+ , curved = 1
+ , vertex.color = NA
+ , vertex.shape = "none"
+ , vertex.label = labs
+ , vertex.label.dist = -6.5
+ , vertex.label.cex = 0.8
+ , vertex.label.color = "black"
+ , vertex.label.degree = -pi/2
+ , edge.width = 2
+ , edge.color = NA
+ , layout = b_auto_layout_scaled * 3
+ , rescale = FALSE
+ , xlim = c(-3.3,3.3)
+ , ylim = c(-3.1,3.1)
+ , add = T
+ )
+ }
diff --git a/README.Rmd b/README.Rmd
index 5eb0592..a3f462f 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -1,9 +1,8 @@
---
output: github_document
+bibliography: vignettes\\references.bib
---
-
-
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
@@ -23,19 +22,25 @@ library(knitr)
::: {.column width="60%"}
-`DemoKin` uses matrix demographic methods to compute expected (average) kin counts from demographic rates under a range of scenarios and assumptions. The package is an R-language implementation of Caswell (2019), Caswell (2020), and Caswell and Song (2021). It draws on previous theoretical development by Goodman, Keyfitz and Pullum (1974).
+`DemoKin` uses matrix demographic methods to compute expected (average) kin counts from demographic rates under a range of scenarios and assumptions. The package is an R-language implementation of Caswell [-@caswell_formal_2019; -@caswell_formal_2020; -@caswell_formal_2022], and Caswell and Song [-@caswell_formal_2021]. It draws on previous theoretical development by Goodman, Keyfitz and Pullum [-@goodman_family_1974].
:::
::: {.column width="40%"}
-
+
:::
::::::::::::::
## Installation
-You can install the development version from GitHub with:
+Download the stable version [from CRAN](https://cran.r-project.org/web/packages/DemoKin/):
+
+``` {r, eval=FALSE, include = T}
+install.packages("DemoKin")
+```
+
+Or you can install the development version from GitHub:
``` {r, eval=FALSE}
# install.packages("devtools")
@@ -44,22 +49,22 @@ devtools::install_github("IvanWilli/DemoKin")
## Usage
-Consider an average Swedish woman called 'Focal'. For this exercise, we assume a female closed population in which everyone experiences the Swedish 2015 mortality and fertility rates at each age throughout their life (the 'time-invariant' assumption in Caswell [2019]).
+Consider an average Swedish woman called 'Focal.' For this exercise, we assume a female closed population in which everyone experiences the Swedish 2015 mortality and fertility rates at each age throughout their life; i.e., the 'time-invariant' assumption in Caswell [-@caswell_formal_2019].
We then ask:
-> How many living relatives does Focal have at each age?
+> What is the expected number of relatives of Focal over her life course?
Let's explore this using the Swedish data already included with `DemoKin`.
-```{r, fig.height=6, fig.width=8}
+```{r, fig.height=6, fig.width=8, message=FALSE, warning=FALSE}
library(DemoKin)
swe_surv_2015 <- swe_px[,"2015"]
swe_asfr_2015 <- swe_asfr[,"2015"]
-swe_2015 <- kin(U = swe_surv_2015, f = swe_asfr_2015, time_invariant = TRUE)
+swe_2015 <- kin(p = swe_surv_2015, f = swe_asfr_2015, time_invariant = TRUE)
```
-*px* is the survival probability by age from a life table and *f* are the age specific fertility raties by age (see `?kin` for details).
+*p* is the survival probability by age from a life table and *f* are the age specific fertility ratios by age (see `?kin` for details).
Now, we can visualize the implied kin counts (i.e., the average number of living kin) of Focal at age 35 using a network or 'Keyfitz' kinship diagram with the function `plot_diagram`:
@@ -75,12 +80,13 @@ plot_diagram(kin_total, rounding = 2)
Relatives are identified by a unique code:
```{r, fig.height=6, fig.width=8, echo=FALSE}
-kable(DemoKin::demokin_codes()[-2])
+# kable(DemoKin::demokin_codes[,c("DemoKin", "Labels_2sex")])
+kable(DemoKin::demokin_codes[,-c(2)])
```
## Vignette
-For more details, including an extension to time varying-populations rates, deceased kin, and multi-state models, see `vignette("Reference", package = "DemoKin")`.
+For more details, including an extension to time-variant rates, deceased kin, and multi-state models in a one-sex framework, see the [Reference_OneSex](https://cran.r-project.org/web/packages/DemoKin/vignettes/Reference_OneSex.html) vignette; also accessible from DemoKin: `vignette("Reference_OneSex", package = "DemoKin")`. For two-sex models, see the [Reference_TwoSex](https://cran.r-project.org/web/packages/DemoKin/vignettes/Reference_TwoSex.html) vignette; also accessible from DemoKin: `vignette("Reference_TwoSex", package = "DemoKin")`.
If the vignette does not load, you may need to install the package as `devtools::install_github("IvanWilli/DemoKin", build_vignettes = T)`.
## Citation
@@ -98,18 +104,3 @@ If you're interested in contributing, please get in touch, create an issue, or s
We look forward to hearing from you!
## References
-
-Caswell, H. 2019. The formal demography of kinship: A matrix formulation. Demographic Research 41:679–712. doi:10.4054/DemRes.2019.41.24.
-
-Caswell, H. 2020. The formal demography of kinship II: Multistate models, parity, and sibship. Demographic Research 42: 1097-1144. doi:10.4054/DemRes.2020.42.38.
-
-Caswell, Hal and Xi Song. 2021. “The Formal Demography of Kinship. III. Kinship Dynamics with Time-Varying Demographic Rates.” Demographic Research 45: 517–46. doi:10.4054/DemRes.2021.45.16.
-
-Goodman, L.A., Keyfitz, N., and Pullum, T.W. (1974). Family formation and the frequency of various kinship relationships. Theoretical Population Biology 5(1):1–27. doi:10.1016/0040-5809(74)90049-5.
-
-
-
-
-
-
-
diff --git a/README.md b/README.md
index 0142d8b..9255f10 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,4 @@
-
-
# DemoKin