Changes to remedy some problems.

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: mortAAR
 Type: Package
 Title: Analysis of Archaeological Mortality Data
-Version: 1.1.6
+Version: 1.1.7
 Authors@R: c(
   person("Nils", "Mueller-Scheessel", email = "[email protected]", role = c("aut", "cre", "cph")), 
   person("Martin", "Hinz", email = "[email protected]", role = c("aut")),
@@ -19,7 +19,7 @@ Description: A collection of functions for the analysis of archaeological mortal
     <https://books.google.de/books?id=nG5FoO_becAC&lpg=PA27&ots=LG0b_xrx6O&dq=life%20table%20archaeology&pg=PA27#v=onepage&q&f=false>). 
     It takes demographic data in different formats and displays the result in a standard life table
     as well as plots the relevant indices (percentage of deaths, survivorship, probability of death, life expectancy, percentage of population).
-Date: 2024-03-06
+Date: 2024-11-28
 License: GPL-3 | file LICENSE
 Encoding: UTF-8
 LazyData: true
@@ -29,7 +29,9 @@ Imports:
     reshape2 (>= 1.4.2),
     methods (>= 3.3.3),
     tibble (>= 3.0.3),
-    rlang (>= 1.1.1)
+    rlang (>= 1.1.1),
+    flexsurv (>= 2.2.2),
+    stats (>= 4.3.1)
 RoxygenNote: 7.2.3
 Suggests:
     testthat,

NEWS.md

@@ -1,3 +1,8 @@
+# mortAAR 1.1.7
+- added function for Halley bands
+- added function for population simulation
+- added helper functions for the random generating and applying of age categories
+
 # mortAAR 1.1.6
 - fixed the error in plotting that occurred after the fix of "aes_string" of ggplot2
   and that still persisted

R/halley_band.R

@@ -5,7 +5,7 @@
 #' methods of sampling from an skeletal population. It basically involves
 #' sampling n times from the age-estimation of each individual and then
 #' only taking the 2.5th and 97.5th percentile into account. The space
-#' between they dubbed 'Halley band' but pointed out that it
+#' between, they dubbed 'Halley band' but pointed out that it
 #' is not to be confused with confidence intervals.
 #'
 #' @param x a data.frame with individuals and age estimations.
@@ -48,20 +48,26 @@
 #' sim_ranges <- random.cat()
 #'
 #' # apply random age categories to simulated ages
-#' sim_appl <- random.cat.apply(pop_sim$result, age = "age", age_ranges = sim_ranges, from = "from", to = "to")
+#' sim_appl <- random.cat.apply(pop_sim$result, age = "age",
+#' age_ranges = sim_ranges, from = "from", to = "to")
 #'
 #' # create halley bands
-#' demo <- halley.band(sim_appl, n = 1000, uncert = 0.95, agebeg = "from", ageend = "to", agerange = "excluded")
+#' demo <- halley.band(sim_appl, n = 1000, uncert = 0.95, agebeg = "from",
+#' ageend = "to", agerange = "excluded")
 #'
 #' # plot band with ggplot
 #' library(ggplot2)
-#' ggplot(demo) + geom_ribbon(aes(x = age, ymin = lower_dx, ymax = upper_dx), linetype = 0, fill = "grey")
-#' ggplot(demo) + geom_ribbon(aes(x = age, ymin = lower_lx, ymax = upper_lx), linetype = 0, fill = "grey")
-#' ggplot(demo) + geom_ribbon(aes(x = age, ymin = lower_qx, ymax = upper_qx), linetype = 0, fill = "grey")
+#' ggplot(demo) + geom_ribbon(aes(x = age, ymin = lower_dx, ymax = upper_dx),
+#' linetype = 0, fill = "grey")
+#' ggplot(demo) + geom_ribbon(aes(x = age, ymin = lower_lx, ymax = upper_lx),
+#' linetype = 0, fill = "grey")
+#' ggplot(demo) + geom_ribbon(aes(x = age, ymin = lower_qx, ymax = upper_qx),
+#' linetype = 0, fill = "grey")
 
 #' @rdname halley.band
 #' @export
-halley.band <- function(x, n = 1000, uncert = 0.95, agebeg, ageend, agerange = "excluded") {
+halley.band <- function(x, n = 1000, uncert = 0.95, agebeg, ageend,
+                        agerange = "excluded") {
   asd <- data.frame(x)
 
   # Change the names of agebeg and ageend for further processes to "beg" and "ende".
@@ -78,21 +84,60 @@ halley.band <- function(x, n = 1000, uncert = 0.95, agebeg, ageend, agerange = "
 
   demo_sim_list <- list()
   for (i in 1:n) {
-    demo_sim_sim <- data.frame(ind = 1:nrow(asd)) %>%
-      dplyr::mutate(age = (round(runif(dplyr::n(), min = asd$beg, asd$ende) ) ) )
-    necdf <- data.frame(a = 1, demo_sim_sim %>% dplyr::group_by(age) %>% dplyr::summarize(Dx = dplyr::n()))
+    # Create the demo_sim_sim data frame
+    demo_sim_sim <- data.frame(ind = 1:nrow(asd))
+    demo_sim_sim$age <- round(stats::runif(nrow(asd), min = asd$beg, max = asd$ende))
 
-    necdf['dx'] <- necdf['Dx'] / sum(necdf['Dx']) * 100
-    necdf['lx'] <- c(100, 100 - cumsum(necdf[, 'dx']))[1:nrow(necdf)]
-    necdf['qx'] <- necdf['dx'] / necdf['lx'] * 100
-    necdf['Ax'] <- necdf[, 'a'] / 2
-    necdf['Lx'] <- necdf['a']* necdf['lx'] - ((necdf['a'] - necdf['Ax']) * necdf['dx'])
+    # Create necdf
+    age_counts <- table(demo_sim_sim$age)
+    necdf <- data.frame(
+      a = 1,
+      age = as.numeric(names(age_counts)),
+      Dx = as.numeric(age_counts)
+    )
+
+    # Add computed columns to necdf
+    necdf$dx <- necdf$Dx / sum(necdf$Dx) * 100
+    necdf$lx <- c(100, 100 - cumsum(necdf$dx))[1:nrow(necdf)]
+    necdf$qx <- necdf$dx / necdf$lx * 100
+    necdf$Ax <- necdf$a / 2
+    necdf$Lx <- necdf$a * necdf$lx - ((necdf$a - necdf$Ax) * necdf$dx)
+
+    # Store necdf in the list
     demo_sim_list[[i]] <- necdf
   }
-  demo_sim_all <- dplyr::bind_rows(demo_sim_list)
-  output <- demo_sim_all %>% dplyr::group_by(age) %>% dplyr::summarize(
-    lower_dx = quantile(dx, probs = low_q), upper_dx = quantile(dx, probs = up_q) ,
-    lower_qx = quantile(qx, probs = low_q), upper_qx = quantile(qx, probs = up_q) ,
-    lower_lx = quantile(lx, probs = low_q), upper_lx = quantile(lx, probs = up_q))
+
+  # Combine all data frames in the list into one
+  demo_sim_all <- do.call(rbind, demo_sim_list)
+
+  # Prepare for summarization
+  output <- data.frame()
+
+  # Loop through unique ages to calculate quantiles
+  unique_ages <- sort(unique(demo_sim_all$age))
+  for (age in unique_ages) {
+    # Subset data for the current age
+    age_data <- demo_sim_all[demo_sim_all$age == age, ]
+
+    # Calculate quantiles
+    lower_dx <- stats::quantile(age_data$dx, probs = low_q)
+    upper_dx <- stats::quantile(age_data$dx, probs = up_q)
+    lower_qx <- stats::quantile(age_data$qx, probs = low_q)
+    upper_qx <- stats::quantile(age_data$qx, probs = up_q)
+    lower_lx <- stats::quantile(age_data$lx, probs = low_q)
+    upper_lx <- stats::quantile(age_data$lx, probs = up_q)
+
+    # Append the results to the output
+    output <- rbind(output, data.frame(
+      age = age,
+      lower_dx = lower_dx,
+      upper_dx = upper_dx,
+      lower_qx = lower_qx,
+      upper_qx = upper_qx,
+      lower_lx = lower_lx,
+      upper_lx = upper_lx
+    ))
+  }
+
   return(output)
 }

R/population_simulation.R

@@ -60,18 +60,18 @@ pop.sim.gomp <- function(n, b = NULL, M = NULL, start_age = 15, max_age = 100) {
     M_1 <- 0
     M_2 <- 0
     while ( M < M_1 | M > M_2 ) {
-      b <- runif(n = 1, min = 0.02, max = 0.1)
-      a <- exp(rnorm(1, (-66.77 * (b - 0.0718) - 7.119), 0.0823))
+      b <- stats::runif(n = 1, min = 0.02, max = 0.1)
+      a <- exp(stats::rnorm(1, (-66.77 * (b - 0.0718) - 7.119), 0.0823))
       M_ <- 1 / b * log (b/a) + start_age
       M_1 <- M_ - 0.001
       M_2 <- M_ + 0.001
     }
   } else if (length(b) > 0) {
-    a <- exp(rnorm(1, (-66.77 * (b - 0.0718) - 7.119), sqrt(0.0823) ) )
+    a <- exp(stats::rnorm(1, (-66.77 * (b - 0.0718) - 7.119), sqrt(0.0823) ) )
     M <- 1 / b * log (b/a) + start_age
   } else {
-    b <- runif(n = 1, min = 0.02, max = 0.05)
-    a <- exp(rnorm(1, (-66.77 * (b - 0.0718) - 7.119), sqrt(0.0823) ) )
+    b <- stats::runif(n = 1, min = 0.02, max = 0.05)
+    a <- exp(stats::rnorm(1, (-66.77 * (b - 0.0718) - 7.119), sqrt(0.0823) ) )
     M <- 1 / b * log (b/a) + start_age
   }
   count = 1

R/random_categories.R

@@ -3,7 +3,7 @@
 #' Helper function that generates random age categories of absolute ages.
 #' It is mainly used together with the functions \code{pop.sim.gomp}
 #' and \code{random.cat.apply}. It will run until the number of categories
-#' are reached \emp{and} there are no gaps in the sequence left.
+#' are reached \emph{and} there are no gaps in the sequence left.
 #'
 #' @param n_cat numeric. Number of categories, default: 20.
 #'
@@ -33,14 +33,14 @@ random.cat <- function(n_cat = 20, min_age = 15, max_cat_low = 60, max_age = 74)
   seq_test_result <- 1
   sim_ranges <- data.frame()
   while (n_sim_ranges < n_cat | length(seq_test_result) > 0){
-    range_from <- round(runif(1, min = min_age, max = max_age)/5) * 5
+    range_from <- round(stats::runif(1, min = min_age, max = max_age)/5) * 5
     if(range_from > max_cat_low) {range_from <- max_cat_low}
 
     #define probabilities
     beta_1 <- range_from/40
     if(beta_1 == 0){beta_1 <- 0.01}
     beta_2 <- 3 - beta_1
-    range_probs <- dbeta(seq(1, 8, 1)/8.1, beta_1, beta_2)
+    range_probs <- stats::dbeta(seq(1, 8, 1)/8.1, beta_1, beta_2)
 
     range_to <- range_from + sample(seq(1, 8, 1), 1 , prob = range_probs) *5 -1
     if(range_to > max_cat_low) {range_to <- max_age}
@@ -69,6 +69,8 @@ random.cat <- function(n_cat = 20, min_age = 15, max_cat_low = 60, max_age = 74)
 #'
 #' @param x a data.frame with individual absolute ages.
 #'
+#' @param age the column containing the individual absolute ages.
+#'
 #' @param age_ranges a data.frame with age ranges.
 #'
 #' @param from numeric. Column name for the begin of an age range.
@@ -90,7 +92,8 @@ random.cat <- function(n_cat = 20, min_age = 15, max_cat_low = 60, max_age = 74)
 #' sim_ranges <- random.cat()
 #'
 #' # apply random age categories to simulated ages
-#' sim_appl <- random.cat.apply(pop_sim$result, age = "age", age_ranges = sim_ranges, from = "from", to = "to")
+#' sim_appl <- random.cat.apply(pop_sim$result, age = "age",
+#' age_ranges = sim_ranges, from = "from", to = "to")
 
 #' @rdname random.cat.apply
 #' @export