A short note on a bug that is easy to write, hard to see, and lives in a lot of hand-rolled survival and life-table code.
Fit a parametric mortality model (Gompertz, Gompertz–Makeham, Siler, …) to a life table by maximum likelihood and you eventually need the probability that an individual alive at the start of an age interval dies during it. If you have
| # Load necessary packages | |
| library(V.PhyloMaker2) | |
| library(Rcompadre) | |
| library(dplyr) | |
| library(ape) | |
| # Fetch COMPADRE data | |
| compadre <- cdb_fetch("compadre") | |
| # Extract unique species list. |
| # de Finetti diagram | |
| library(dplyr) | |
| library(ggplot2) | |
| library(ggtern) | |
| HWEfrequencies <- data.frame( | |
| p = seq(0, 1, 0.01) | |
| ) %>% | |
| mutate(q = 1 - p) %>% | |
| mutate(AA = p^2, Aa = 2 * p * q, aa = q^2) |
| # I first simulate a dataset to use for this example | |
| set.seed(123) | |
| years <- 2013:2015 | |
| n <- 100 | |
| df1 <- data.frame(year = rep(years, each = n), | |
| eggdate = c(rnorm(n, 40, 4), | |
| rnorm(n, 45, 5), | |
| rnorm(n, 67, 7))) | |
| # Make sure year is coded as a "factor" |
| # Packages | |
| library(ggplot2) | |
| # Read in some data for the example | |
| # The data are grade data and acceptance into grad school. | |
| # GRE (Graduate Record Exam scores), GPA (grade point average) | |
| mydata <- read.csv("https://stats.oarc.ucla.edu/stat/data/binary.csv") | |
| head(mydata) | |
| summary(mydata) | |
| mylogit <- glm(admit ~ gre + gpa + gre:gpa, data = mydata, family = "binomial") |
| # offsets | |
| # Simulate some data for this example | |
| set.seed(12) | |
| sig <- matrix(c(1.1, 1, 1, 1.1), 2, 2) | |
| df1 <- data.frame(MASS::mvrnorm(n = 100, c(30, 30), Sigma = sig)) | |
| names(df1) <- c("x", "y") | |
| # Plot the data, and a 1:1 line | |
| plot(df1$x, df1$y) | |
| abline(0, 1) |
| # How to subset the matrix database based on the matrices themselves. | |
| # This example counts NA values in the F (fecundity) matrix of each population | |
| # and uses that to subset. Uses the modern Rcompadre CompadreDB API. | |
| library(Rcompadre) | |
| compadre <- cdb_fetch("compadre") | |
| # matF() returns a list of the F matrices, one per row of the database. | |
| # Count the NA values in each, and store as a new metadata column. | |
| compadre$NAinFmat <- sapply(matF(compadre), function(m) sum(is.na(m))) |
| #bar heights | |
| x<- c(1,2,3,2,4) | |
| #standard error values | |
| sem <- c(.5,.4,.2,.4,.2) | |
| #bar labels | |
| lab <- c("A","B","C","D","E") |
| par(mar = c(1, 1, 1, 1)) | |
| plot(x = rep(1:5, 5), y = rep(1:5, each = 5), pch = 1:25, col = "black", bg = "yellow", | |
| axes = FALSE, xlab = "", ylab = "", | |
| cex = 2, xlim = c(0, 5.5), ylim = c(0, 5.5)) | |
| text(x = rep(1:5, 5), y = rep(1:5, each = 5), labels = 1:25, pos = 2, cex = .8) |
| # How to subset the COMADRE Animal Matrix Database (works with COMPADRE too) | |
| # | |
| # NOTE: this uses the legacy list-style COMADRE object (metadata / mat / matrixClass). | |
| # Modern workflows use Rcompadre: `db <- cdb_fetch("comadre")` then | |
| # `subset(db, Order %in% c(...))`. See https://github.com/jonesor/Rcompadre | |
| # Create a vector of indices to retain | |
| subsetID <- which(comadre$metadata$Order %in% c("Monotremata", "Didelphimorphia", | |
| "Paucituberculata", "Microbiotheria", | |
| "Dasyuromorphia", "Peramelemorphia", |