I have this function:
mweFitModelsGLMER <- function(longData,
strModelName = "ID Rand Model",
strFormula = "phMeta_Toxicity ~ 1 + (1 | phMeta_UID)",
...)
{
output = list() # Initialise output
Ptn <- unique(longData$Ptn)
Model <- longData %>%
lme4::glmer(formula = stats::as.formula(strFormula),
...)
#Do a few other things like prepare summary table etc.
output$Models <-
# list(
tibble(
Protein = Ptn,
Formula = strFormula,
`Model Name` = strModelName,
Model = list(Model),
)
#)
return(output)
}
Whose primary purpose is to fit a mixed model, do a few other things and then output a big tibble. This works well on its own with the data at the end of the post:
mweFitModelsGLMER(allData, family = "binomial", na.action = na.pass)
Now I try to create a helper function, which can take in a named list of formula specifications, pass it on to the function above and return the output:
mweMultipleModelsGLMER <- function(longData,
namLstFormula = list(Model1 = "phMeta_Toxicity ~ 1 + (1 | phMeta_UID)",
Model2 = "phMeta_Toxicity ~ value + (1 | phMeta_UID)"),
...)
{
purrr::imap(
namLstFormula,
~ mweFitModelsGLMER(
longData = longData,
strModelName = paste(.y),
strFormula = .x,
...
)
)
}
And run it like so:
mweMultipleModelsGLMER(allData,
namLstFormula =
list(Model1 = "phMeta_Toxicity ~ 1 + (1 | phMeta_UID)",
Model2 = "phMeta_Toxicity ~ value + (1 | phMeta_UID)"),
family = "binomial")
To which, R responds with
Error in
map2(): ℹ In index: 1. ℹ With name: Model1. Caused by error inlme4::glmer(): ! 'control' is not a list; use glmerControl()
Ok, then I be explicit in my call..
mweMultipleModelsGLMER(
allData,
namLstFormula = list(Model1 = "phMeta_Toxicity ~ 1 + (1 | phMeta_UID)",
Model2 = "phMeta_Toxicity ~ value + (1 | phMeta_UID)"),
family = "binomial",
control = lme4::glmerControl(optimizer = "bobyqa")
)
To which R comes back to me with:
Error in
map2(): ℹ In index: 1. ℹ With name: Model1. Caused by error inlme4::glmer(): ! 'control' is not a list; use glmerControl()
Backtrace:
▆
1. ├─global mweMultipleModelsGLMER(...)
2. │ └─purrr::imap(...)
3. │ └─purrr::map2(.x, vec_index(.x), .f, ...)
4. │ └─purrr:::map2_("list", .x, .y, .f, ..., .progress = .progress)
5. │ ├─purrr:::with_indexed_errors(...)
6. │ │ └─base::withCallingHandlers(...)
7. │ ├─purrr:::call_with_cleanup(...)
8. │ └─.f(.x[[i]], .y[[i]], ...)
9. │ └─global mweFitModelsGLMER(...)
10. │ └─longData %>% ...
11. └─lme4::glmer(., formula = stats::as.formula(strFormula), ...)
12. └─base::stop("'control' is not a list; use glmerControl()")
For reference, I did help("glmer")
glmer(formula, data = NULL, family = gaussian , control = glmerControl() , start = NULL , verbose = 0L , nAGQ = 1L , subset, weights, na.action, offset, contrasts = NULL , mustart, etastart , devFunOnly = FALSE)
What am I missing?
Data below:
library(tidyverse)
allData <- structure(list(phMeta_UID = c("Pat 1 BRS", "Pat 1 BRS", "Pat 1 BRS",
"Pat 1 BRS", "Pat 1 BRS", "Pat 1 BRS", "Pat 1 BRS", "Pat 1 BRS",
"Pat 1 BRS", "Pat 1 SHF", "Pat 1 SHF", "Pat 1 SHF", "Pat 1 SHF",
"Pat 1 SHF", "Pat 1 SHF", "Pat 1 SHF", "Pat 2 BRS", "Pat 2 BRS",
"Pat 2 BRS", "Pat 2 SHF", "Pat 2 SHF", "Pat 2 SHF", "Pat 2 SHF",
"Pat 2 SHF", "Pat 2 SHF", "Pat 3 SHF", "Pat 3 SHF", "Pat 3 SHF",
"Pat 3 SHF", "Pat 3 SHF", "Pat 3 SHF", "Pat 3 SHF", "Pat 3 SHF",
"Pat 1 SHF", "Pat 2 BRS", "Pat 2 BRS", "Pat 2 BRS", "Pat 2 BRS",
"Pat 2 BRS", "Pat 2 BRS", "Pat 2 SHF"), phMeta_Time = c(0, 0.5,
1, 2, 3, 4, 6, 8, 12, 0, 0.5, 2, 3, 4, 8, 12, 0, 0.5, 2, 0, 0.5,
2, 3, 4, 10, 3, 4, 8, 9, 12, 0, 0.5, 2, 1, 1, 3, 4, 6, 8, 12,
1), phMeta_Batch = c(1, 1, 2, 1, 1, 1, 2, 2, 2, 2, 2, 1, 2, 1,
1, 2, 1, 2, 2, 1, 1, 2, 2, 2, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1,
1, 2, 1, 1, 2, 2), phMeta_Site = c("BRS", "BRS", "BRS", "BRS",
"BRS", "BRS", "BRS", "BRS", "BRS", "SHF", "SHF", "SHF", "SHF",
"SHF", "SHF", "SHF", "BRS", "BRS", "BRS", "SHF", "SHF", "SHF",
"SHF", "SHF", "SHF", "SHF", "SHF", "SHF", "SHF", "SHF", "SHF",
"SHF", "SHF", "SHF", "BRS", "BRS", "BRS", "BRS", "BRS", "BRS",
"SHF"), phMeta_Toxicity = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L), levels = c("N", "Y"), class = "factor"), phMeta_Patient = c(4,
4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 2, 2, 2,
2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 1, 5, 5, 5, 5, 5, 5, 2), phMeta_SiteXPatient = c(1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 1, 2, 2, 2, 2, 2, 2, 2), phMeta_Phlebotomy = c(NA,
NA, NA, NA, NA, NA, NA, NA, NA, "Venous", "Capillary", "Venous",
NA, "Venous", "Capillary", "Venous", NA, NA, NA, "Venous", "Capillary",
"Capillary", "Capillary", "Capillary", "Venous", "Venous", "Venous",
"Venous", NA, "Venous", "Venous", "Capillary", "Venous", "Venous",
NA, NA, NA, NA, NA, NA, "Capillary"), value = structure(c(0.821944197243299,
-0.0198715543825022, 0.631293040769747, 0.0849008002934989, 0.0887812578740912,
1.05453893286552, -2.12824809442977, 0.276304235154362, 0.859985670512456,
-1.04724028808727, -0.00411159277206202, 1.67249503748148, -0.54509175851945,
-1.46818604842327, -0.498407201908304, 0.51962029081445, -1.08569307903582,
-0.270959349353233, -0.211864905388811, 0.158115347225517, 0.0227322978830837,
0.852349233070034, -0.401612245382643, -2.98257656282869, -0.191542564781942,
-1.45732401444245, 0.494605681417659, 0.464925863604591, 0.856056504259303,
1.63695935481179, -1.33574175565861, 0.42694933523218, -0.0213328145592944,
-0.966697791972374, -0.634006734239892, 0.128202810199108, 0.960861331383678,
-0.258051551124902, 1.25488311517846, 1.08015428255721, 1.18190128745978
), dim = c(41L, 1L)), Ptn = c("phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand",
"phOSI_ICOS ligand", "phOSI_ICOS ligand", "phOSI_ICOS ligand"
)), row.names = c(NA, -41L), class = c("tbl_df", "tbl", "data.frame"
))
The function works if I ignore the list names:
fFitMultipleModelsGLMER2 <- function(nestlistData,
strPredictorROC = "phMeta_Toxicity",
LstFormula = c("phMeta_Toxicity ~ 1 + (1 | phMeta_UID)",
"phMeta_Toxicity ~ value + (1 | phMeta_UID)"),
...)
{
purrr::map(
namLstFormula,
~fFitModelsGLMER(
longData = nestlistData,
strPredictor = strPredictorROC,
strFormula = .x,
...
)
)
}
EDIT2: Editing to improve discoverability
The issue is the use of a purrr-style anonymous function, i.e. using a formula, and .... As an unintended side-effect, instead of the optional arguments you passed via ... to the "outer" function, you are passing the .x and .y arguments to the "inner" mweFitModelsGLMER a second time via ....
This can be seen from the following example where I simply print the content of ... inside the imap call:
library(purrr)
foo <- function(longData,
namLstFormula = list(
Model1 = "phMeta_Toxicity ~ 1 + (1 | phMeta_UID)",
Model2 = "phMeta_Toxicity ~ value + (1 | phMeta_UID)"
),
...) {
purrr::imap(
namLstFormula,
~ print(list(...))
)
}
baz <- foo(allData,
namLstFormula =
list(
Model1 = "phMeta_Toxicity ~ 1 + (1 | phMeta_UID)",
Model2 = "phMeta_Toxicity ~ value + (1 | phMeta_UID)"
),
family = "binomial"
)
#> [[1]]
#> [1] "phMeta_Toxicity ~ 1 + (1 | phMeta_UID)"
#>
#> [[2]]
#> [1] "Model1"
#>
#> [[1]]
#> [1] "phMeta_Toxicity ~ value + (1 | phMeta_UID)"
#>
#> [[2]]
#> [1] "Model2"
To fix your issue use function(.x, .y) or \(.x, .y) instead of a formula:
mweMultipleModelsGLMER <- function(longData,
namLstFormula = list(
Model1 = "phMeta_Toxicity ~ 1 + (1 | phMeta_UID)",
Model2 = "phMeta_Toxicity ~ value + (1 | phMeta_UID)"
),
...) {
purrr::imap(
namLstFormula,
\(.x, .y)
mweFitModelsGLMER(
longData = longData,
strModelName = .y,
strFormula = .x,
...
)
)
}
mweMultipleModelsGLMER(allData,
namLstFormula =
list(
Model1 = "phMeta_Toxicity ~ 1 + (1 | phMeta_UID)",
Model2 = "phMeta_Toxicity ~ value + (1 | phMeta_UID)"
),
family = "binomial"
)
#> $Model1
#> $Model1$Models
#> # A tibble: 1 × 4
#> Protein Formula `Model Name` Model
#> <chr> <chr> <chr> <list>
#> 1 phOSI_ICOS ligand phMeta_Toxicity ~ 1 + (1 | phMeta_U… Model1 <glmerMod>
#>
#>
#> $Model2
#> $Model2$Models
#> # A tibble: 1 × 4
#> Protein Formula `Model Name` Model
#> <chr> <chr> <chr> <list>
#> 1 phOSI_ICOS ligand phMeta_Toxicity ~ value + (1 | phMe… Model2 <glmerMod>