I have a simulated data set where each row is an individual bird, and I've written some functions that determine whether each individual lives or dies. The conditions for whether or not an individual survives is based off of its age class (AHY or HY) and its sex (M or F). I created a function for each age/sex combination, and use pmap_chr inside of mutate/case_when, which should fill in a column called status. In my code this provides a value of 'live' or 'die.' Here's an abbreviated version of my dataset:
library(tidyverse)
agents <- structure(list(id = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18), sex = c("F", "F", "F", "F", "F", "M",
"M", "M", "M", "M", "F", "F", "M", "M", "M", "M", "M", "M"),
class = c("AHY", "AHY", "AHY", "AHY", "AHY", "AHY", "AHY",
"AHY", "AHY", "AHY", "HY", "HY", "HY", "HY", "HY", "HY",
"HY", "HY"), hDateCtr = c(-0.84852029241304, 0.558881154137435,
-0.909711659654365, 1.21158907137824, -0.56296057862019,
-0.0938267631033649, -1.54202245448139, -0.216209497586015,
1.33397180586089, 1.06880921448181, -0.935414346693485, -0.935414346693485,
-0.935414346693485, -0.935414346693485, 0.935414346693485,
0.935414346693485, 0.935414346693485, 0.935414346693485),
aDateCtr = c(-1.13245629117638, 1.13245629117638, -0.490731059509763,
1.13245629117638, -0.641725231666613, 1.13245629117638, -1.13245629117638,
1.13245629117638, -0.490731059509763, -0.641725231666613,
NA, NA, NA, NA, NA, NA, NA, NA), selfOrig = c("imm", "imm",
"imm", "imm", "imm", "imm", "imm", "imm", "imm", "imm", "local",
"local", "local", "local", "local", "local", "local", "local"
), sameSexOrig = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
"imm", "imm", "imm", "imm", "imm", "imm", "imm", "imm"),
success = c(TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE,
TRUE, FALSE, FALSE, NA, NA, NA, NA, NA, NA, NA, NA), paired = c(TRUE,
TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE,
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE)), class = c("tbl_df",
"tbl", "data.frame"), row.names = c(NA, -18L))
# A tibble: 18 x 9
id sex class hDateCtr aDateCtr selfOrig sameSexOrig success paired
<dbl> <chr> <chr> <dbl> <dbl> <chr> <chr> <lgl> <lgl>
1 1 F AHY -0.849 -1.13 imm NA TRUE TRUE
2 2 F AHY 0.559 1.13 imm NA TRUE TRUE
3 3 F AHY -0.910 -0.491 imm NA FALSE TRUE
4 4 F AHY 1.21 1.13 imm NA FALSE TRUE
5 5 F AHY -0.563 -0.642 imm NA FALSE TRUE
6 6 M AHY -0.0938 1.13 imm NA FALSE TRUE
7 7 M AHY -1.54 -1.13 imm NA TRUE TRUE
8 8 M AHY -0.216 1.13 imm NA TRUE TRUE
9 9 M AHY 1.33 -0.491 imm NA FALSE TRUE
10 10 M AHY 1.07 -0.642 imm NA FALSE TRUE
11 11 F HY -0.935 NA local imm NA FALSE
12 12 F HY -0.935 NA local imm NA FALSE
13 13 M HY -0.935 NA local imm NA FALSE
14 14 M HY -0.935 NA local imm NA FALSE
15 15 M HY 0.935 NA local imm NA FALSE
16 16 M HY 0.935 NA local imm NA FALSE
17 17 M HY 0.935 NA local imm NA FALSE
18 18 M HY 0.935 NA local imm NA FALSE
Here's an example of the mortality functions I've written that go into pmap_chr. These all work fine if I run the code below on a data set that has a single age class or sex:
hDateEffect <- TRUE
winterTemp <- -3
# hatchling mortality -----------------------------------------------------
hatchMortInt <- -4.67
hatchMortIntSD <- 0.39
hatchMortBeta1 <- 0.6
hatchMortBeta1SD <- 0.27
hatchMortBeta2 <- 1.12
hatchMortBeta2SD <- 0.36
hatchMortBeta3 <- -0.3
hatchMortBeta3SD <- 0.16
hatchMortBeta4 <- -0.3
hatchMortBeta4SD <- 0.16
# male mortality ----------------------------------------------------------
maleMortInt <- -2.09
maleMortIntSD <- 0.32
maleMortBeta1 <- 0.81
maleMortBeta1SD <- 0.34
maleMortBeta2 <- -1.36
maleMortBeta2SD <- 0.84
maleMortBeta3 <- 1.67
maleMortBeta3SD <- 0.32
# female mortality --------------------------------------------------------
femMortInt <- -0.93
femMortIntSD <- 0.87
femMortBeta1 <- 1.59
femMortBeta1SD <- 0.35
femMortBeta2 <- -1.77
femMortBeta2SD <- 0.78
# hatch-year female
HY_female_mortality <- function(hDateCtr, sameSexOrig, ...) {
intercept <- rnorm(1, hatchMortInt, hatchMortIntSD)
beta2 <- rnorm(1, hatchMortBeta2, hatchMortBeta2SD)
beta4 <- rnorm(1, hatchMortBeta4, hatchMortBeta4SD)
if(hDateEffect == TRUE) {
beta3 <- rnorm(1, hatchMortBeta3, hatchMortBeta3SD)
} else {
beta3 <- 0
}
if (sameSexOrig == 'local') {
linSurv <- intercept + beta2 + (beta3 * hDateCtr) + (beta4 * winterTemp)
} else {
linSurv <- intercept + (beta3 * hDateCtr) + (beta4 * winterTemp)
}
probSurv <- plogis(linSurv)
randDraw <- runif(1, 0, 1)
if (randDraw > probSurv) {
val <- 'die'
return(val)
} else {
val <- 'live'
return(val)
}
}
# hatch-year male
HY_male_mortality <- function(hDateCtr, sameSexOrig, ...) {
intercept <- rnorm(1, hatchMortInt, hatchMortIntSD)
beta1 <- rnorm(1, hatchMortBeta1, hatchMortBeta1SD)
beta2 <- rnorm(1, hatchMortBeta2, hatchMortBeta2SD)
beta4 <- rnorm(1, hatchMortBeta4, hatchMortBeta4SD)
if(hDateEffect == TRUE) {
beta3 <- rnorm(1, hatchMortBeta3, hatchMortBeta3SD)
} else {
beta3 <- 0
}
if (sameSexOrig == 'local') {
linSurv <- intercept + beta1 + beta2 + (beta3 * hDateCtr) + (beta4 * winterTemp)
} else {
linSurv <- intercept + beta1 + (beta3 * hDateCtr) + (beta4 * winterTemp)
}
probSurv <- plogis(linSurv)
randDraw <- runif(1, 0, 1)
if (randDraw > probSurv) {
val <- 'die'
return(val)
} else {
val <- 'live'
return(val)
}
}
# after-hatch-year mortality functions
# after-hatch-year male
AHY_male_mortality <- function(aDateCtr, success, selfOrig, ...) {
intercept <- rnorm(1, maleMortInt, maleMortIntSD)
beta1 <- rnorm(1, maleMortBeta1, maleMortBeta1SD)
beta3 <- rnorm(1, maleMortBeta3, maleMortBeta3SD)
if(hDateEffect == TRUE) {
beta2 <- rnorm(1, hatchMortBeta3, hatchMortBeta3SD)
} else {
beta2 <- 0
}
if (success == TRUE) {
linSurv <- intercept + beta1 + (beta2 * aDateCtr)
} else {
linSurv <- intercept + (beta2 * aDateCtr)
}
if (selfOrig == 'local') {
linSurv <- linSurv + beta3
} else {
linSurv <- linSurv
}
probSurv <- plogis(linSurv)
randDraw <- runif(1, 0, 1)
if (randDraw > probSurv) {
val <- 'die'
return(val)
} else {
val <- 'live'
return(val)
}
}
# after-hatch-year female
AHY_female_mortality <- function(aDateCtr, success, ...) {
intercept <- rnorm(1, femMortInt, femMortIntSD)
beta1 <- rnorm(1, femMortBeta1, femMortBeta1SD)
beta2 <- rnorm(1, femMortBeta2, femMortBeta2SD)
if (success == TRUE) {
linSurv <- intercept + beta1 + (beta2 * aDateCtr)
} else {
linSurv <- intercept + (beta2 * aDateCtr)
}
probSurv <- plogis(linSurv)
randDraw <- runif(1, 0, 1)
if (randDraw > probSurv) {
val <- 'die'
} else {
val <- 'live'
}
return(val)
}
Here's the pmap_chr part, which doesn't work across all of the age and sex class combinations:
agents %>%
mutate(
status = case_when(
class == 'HY' & sex == 'F' ~ pmap_chr(., HY_female_mortality),
class == 'HY' & sex == 'M' ~ pmap_chr(., HY_male_mortality),
class == 'AHY' & sex == 'M' ~ pmap_chr(., AHY_male_mortality),
class == 'AHY' & sex == 'F' ~ pmap_chr(., AHY_female_mortality)
)
)
But if I do this same thing instead for another logical called 'success' (so if (success == TRUE)) which is the one I actually need the condition to be based, it produces an error:
Error in mutate_impl(.data, dots) :
Evaluation error: missing value where TRUE/FALSE needed.
I'm at a loss for why these functions work separately, but not on the entire dataset containing all age and sex classes. I have examples of different processes (reproduction, immigration) where I do a similar thing (take the data set, write functions that are used inside pmap, which in turn are inside case_when and mutate).
pmap_chr(agents, HY_female_mortality)
throws an error for me, which is ultimately what you're calling in yourcase_when
statement. I think it's due to anNA
value being to yoursameSexOrig
value and the design of your user-defined functions. – zacksameSexOrig
does not containNA
values for those observations. – zackNAs
forsameSexOrig
but don't rely on that variable won't work? – Jason