--- title: "Setting up different types of specifications" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Setting up different types of specifications} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} # Evaluate the code only if the suggested packages are installed knitr::opts_chunk$set(eval = requireNamespace("tidyverse", quietly = TRUE)) ``` The function `setup()` seems rather simple, but it is in fact very versatile and allows to incorporate (almost) any analytic choice that one might engage in during data analysis. In this vignette, I am going to exemplify how this function can be used to create several analytical choices. ```{r, message = F, warning = F} library(tidyverse) library(specr) ``` # Different independent and dependent variables One analytic choice may refer to choosing different independent or dependent variables (e.g., different scales that measure the same construct) or different ways in how these variables are compute (e.g., mean vs. latent measurement, log-transformation, etc). ## Simply providing different variable as vectors In most cases, this simply means adding such variables to the data set and provide them as alternative choices in setup. ```{r} # Setup specs specs <- setup(data = example_data, x = c("x1", "x2", "x3", "x4"), y = c("y1", "y2", "y3"), model = "lm") # Summary of specifications summary(specs, rows = 12) # Run analysis and plot results results <- specr(specs) head(results$data) ``` ## Specific variable/model combinations However, in some cases, we may have different type of variables (e.g., a mean vs. a dichotomous variable). In some cases, for example when this refers to a dependent variable, such a choice of variables requires also different model estimation functions. So in our setup, we want to acknowledge for specific variable/model combinations. In the following, let's imagine that next to our standard dependent variables (`y1`-`y3`), we also have a dichotomous item (`y_dich`; here only recoded from `y1`). In this case, we want to specify and add a custom function - here a logistic model (`glm(formula, data, family = binomial())`) - to the setup. ```{r} # Dichotomous dependent variable data <- example_data %>% mutate(y_dich = ifelse(y1 > mean(y1), 1, 0)) # Specific function log_glm <- function(formula, data) { glm(formula, data, family = binomial()) } # Setup specs specs <- setup(data = data, x = c("x1", "x2"), y = c("y1", "y2", "y3", "y_dich"), model = c("lm", "log_glm")) # Check specs %>% as_tibble ``` As we can see, the `setup()` function produces all combinations between the variables and the model functions. This is not meaningful as the logistic model should only be used if the dependent variable is `y_dich`. In comparison to specr version 0.2.1, this is where the strength of version 1.0.0 comes in. Because we setup all specifications beforehand, we can simply remove those that we are not interested in or that we deem non-meaningful. Note: You won't be able to filter or subset the `specr.setup` object (essentially a list). Instead, you want to filter/subset the first position of the list: the data frame that includes the specifications (here `specs$specs` or `specs[[1]]`). ```{r} # Filter out models that are not meaningful (here only keep log_glm, when y == "y_dich") specs$specs <- specs$specs %>% filter(!(model == "log_glm" & y != "y_dich")) %>% filter(!(model == "lm" & y == "y_dich")) # Check results (only meaningful specifications remain) summary(specs, rows = 8) # Run analysis and plot results results <- specr(specs) head(results$data) ``` ## Adding composite and latent variables In some cases, we may want to add mean indices or even latent measurement models of our variables of interest to explore how latent vs. manifest measurement affects our results. This requires a bit more upfront work as we need to work with [lavaan](https://lavaan.ugent.be/) for the latent measurement. Furthermore, we need to adjust the extract functions (`fun1` and `fun2` in `setup()`) so that the extraction of parameters aligns across different type of models. The beauty of `broom::tidy()` and `broom::glance()` (which are used by default) is that they adapt to the different model types that are passed to `setup`. Although the same parameters are extracted for most models, they sometimes differ for specific models (e.g., structural equation models resulting from `lavaan::sem()`). This means we need to adjust them to extract the same parameters across different models. ```{r, warning = F, message = F, eval = requireNamespace("tidyverse", quietly = TRUE) && requireNamespace("lavaan", quietly = TRUE)} # Add mean (one choice) data <- data %>% mutate(x_mean = rowMeans(across(c(x1, x2, x3, x4)))) # Add custom function with latent measurement models to pass to "models" (another choice) custom_sem <- function(formula, data) { # Add latent measurement as list latent <- list(latent_x12 = "latent_x12 =~ x1 + x2") # Remove +1 from formula as lavaan doesn't know how to process it semformula <- str_remove_all(formula, "\\+ 1") # remove non-used latent measurement models from list by checking the formula valid <- purrr::keep(names(latent), ~ stringr::str_detect(formula, .x)) # Create new formula that includes latent measurement models formula <- paste(semformula, "\n", paste(latent[valid], collapse = " \n ")) # Pass formula to `sem()` lavaan::sem(formula, data) } # Create custom tidy function that extracts the same parameters from different models! tidy_new <- function(x) { if(inherits(x, "lavaan")) { broom::tidy(x, conf.int = TRUE) %>% select(term, estimate, conf.low, conf.high) %>% # select parameters you want to keep filter(grepl(" ~ ", term)) %>% # term needs to be adjusted separate(term, c("dv", "term"), sep = " ~ ") %>% # extract independent variable select(-dv) # remove dependent variable } else { broom::tidy(x, conf.int = TRUE) %>% select(term, estimate, conf.low, conf.high) # same parameters as above } } # Setup specs with new custom function specs <- setup(data = data, x = c("x1", "x2", "x3", "x4", "x_mean", "latent_x12"), y = c("y1", "y2"), model = c("lm", "custom_sem"), fun1 = tidy_new, # We pass the new extract function fun2 = NULL) # switch off "glance" as it produces different fit indices and wouldn't work # Quick check (still includes non-meaningful specifications) summary(specs, rows = 12) # Filter out non-meaningful specifications specs$specs <- specs$specs %>% filter(!(model == "custom_sem" & !grepl("latent", x))) %>% filter(!(model == "lm" & grepl("latent", x))) # Check again summary(specs, rows = 12) # Run analysis and plot results results <- specr(specs) plot(results, choices = c("x", "y")) ``` As we can see, we can now investigate the manifest and latent measures next to the individual items. # Control variables Another analytical choice may refer to the inclusion of covariates. Here, the `setup()` function is again very versatile and allows to include different combinations of covariates. ## Adding a bunch of control variables By providing just a vector of covariates, `setup()` produces all combinations of these covariates plus a specification without any covariates. If we add the argument `simplify = TRUE`, not all combinations between covariates are created. Instead, only no covariates, each individually, and all together are included. ```{r} # Setup specification that include all combinations of covariates specs1 <- setup(data = example_data, x = c("x1", "x2"), y = c("y1", "y2"), model = "lm", controls = c("c1", "c2", "c3", "c4")) # simply providing a vector of control variables # Setup specifications that include only no covariates, each individually, and all together specs2 <- setup(data = example_data, x = c("x1", "x2"), y = c("y1", "y2"), model = "lm", controls = c("c1", "c2", "c3", "c4"), simplify = TRUE) # Difference to specs1! # Check distinct(specs1$specs, controls) distinct(specs2$specs, controls) ``` ## Sets of covariates We can also add groups of covariates. The way setup creates combinations between such groups remains the same. ```{r} # Add groups of covariates specs3 <- setup(data = example_data, x = c("x1", "x2"), y = c("y1", "y2"), model = "lm", controls = c("c1 + c2", "c3 + c4")) # Check distinct(specs3$specs, controls) ``` The `setup()` function further has an argument `add_to_formula`, which allows to add covariates (or any other formula-relevant aspects, e.g., random effect structures) to all model functions. ```{r} # Add some control variables to all models specs4 <- setup(data = example_data, x = c("x1", "x2"), y = c("y1", "y2"), model = "lm", controls = c("c1", "c2"), add_to_formula = "c3") # Check (see how `c3` is added to each formula, but is not part of controls) specs4$specs[1:6,] ``` ## Dealing with duplicates in specifications Finally, we may sometimes decide to add one of the independent or dependent variables as covariates as well. In this case, it would not make sense to add e.g., `x1` as covariate if it is already the independent variable. In this new version of specr, the `setup()` function automatically detects such duplications and deletes such specifications. ```{r} # Adding a covariate that is also an independent or dependent variable specs5 <- setup(data = example_data, x = c("x1", "x2"), y = c("y1", "y2"), model = "lm", controls = c("x1", "y1")) # Check (see how only 9 specifications are kept) specs5$specs ``` # Subset analyses One of the most powerful arguments within the `setup()` framework is `subsets`, which can be used for so-called "subset analyses". With this argument, a variety of analytical choices can be included. ## Simple subset analysis based on grouping variables In simple cases, we may want to investigate whether a relationship differs across certain subgroups. In the example data, we can for example investigate whether the relationship between `x` and `y` differs across age groups (`group1`) and gender (`group2`). We can add such subset analyses by specifying a list in two alternative ways. The easiest way refers to simply providing named vectors that refer to the respective variables in the data set and their unique values: ``` list(group1 = c("young", "middle", "old"), group2 = c("female", "male")) ``` Alternatively, we can achieve the same thing by using the function `unique()`, which extracts the unique values from each variable in the data set: ``` list(group1 = unique(example_data$group1), group2 = unique(example_data$group2)) ``` These lists are then passed as combinatorial factors to the expand function within `setup()`. If two or more grouping variables are provided in this way (you can provide as many as you want, but it quickly becomes messy), `setup()` automatically adds them as "choice" columns to the specification setup data frame. Note how it automatically adds the value `NA` to make sure that also the simple subsets and not only the combinations between the two groups are included in the specifications. It further creates a new column "subsets", which represents the combination of the subsetting factors. ```{r} # Setup specifications specs <- setup(data = example_data, x = c("x1", "x2"), y = c("y1", "y2"), model = c("lm"), controls = "c1", subsets = list(group1 = unique(example_data$group1), group2 = unique(example_data$group2))) # Summary of specifications summary(specs) # Check subsets (in this case, 12 different types of subset analysis, # including using "all" subjects) distinct(specs$specs, subsets) # Run analysis and plot results results <- specr(specs) plot(results, choices = c("x", "y", "subsets")) ``` ## More complex analysis using the subsets functionality The true potential of this type of subsetting can only be grasped with another example. Let's imagine we want to include the removal of outliers as an analytical choice in our specification curve analysis. Perhaps we are unsure about the cut-off value and want to investigate different (arbitrary) ones. This can again be seen as a subset analysis. We first create variables that denote who is an outlier and who is not according to the rule we want to use. We then pass these variables as subsets to the `setup()` function. This again produces non-meaningful specifications (e.g., a subset that includes outliers, a subset that includes people who align with different outlier rules, etc.). So we want to keep only the meaningful specifications. In this case, these are just 7 different types of specifications. ```{r} # Create variables that denote outliers (here with a range of arbitrary thresholds) data <- data %>% mutate(outlier1 = ifelse(y1 < mean(y1) - 2*sd(y1) | y1 > mean(y1) + 2*sd(y1), "outlier", "2.0*SD"), outlier2 = ifelse(y1 < mean(y1) - 2.1*sd(y1) | y1 > mean(y1) + 2.1*sd(y1), "outlier", "2.1*SD"), outlier3 = ifelse(y1 < mean(y1) - 2.2*sd(y1) | y1 > mean(y1) + 2.2*sd(y1), "outlier", "2.2*SD"), outlier4 = ifelse(y1 < mean(y1) - 2.3*sd(y1) | y1 > mean(y1) + 2.3*sd(y1), "outlier", "2.3*SD"), outlier5 = ifelse(y1 < mean(y1) - 2.4*sd(y1) | y1 > mean(y1) + 2.4*sd(y1), "outlier", "2.4*SD"), outlier6 = ifelse(y1 < mean(y1) - 2.5*sd(y1) | y1 > mean(y1) + 2.5*sd(y1), "outlier", "2.5*SD")) # Setup specs specs <- setup(data = data, x = c("x1", "x2"), y = c("y1", "y2"), model = "lm", controls = c("c1", "c2"), subsets = list(outlier1 = c("2.0*SD"), outlier2 = c("2.1*SD"), outlier3 = c("2.2*SD"), outlier4 = c("2.3*SD"), outlier5 = c("2.4*SD"), outlier6 = c("2.5*SD"))) # Remove unnecessary combinations specs$specs <- specs$specs %>% filter(subsets == "2.0*SD" | subsets == "2.1*SD" | subsets == "2.2*SD" | subsets == "2.3*SD" | subsets == "2.4*SD" | subsets == "2.5*SD" | subsets == "all") # Check specifications (see how it contains only meaningful subsets?) summary(specs, rows = 7) # Run analysis and plot results results <- specr(specs) plot(results, choices = c("x", "y", "subsets")) ``` # Other types of models The argument `model` in `setup()` accepts any function that has the arguments `formula` and `data` (e.g., `lm`, `glm`, or `lme4::lmer`). Other models can be included by writing a small custom function (a "wrapper") that passes these two arguments on to the respective model function and sets everything else that is needed. Custom functions can be provided by name (`model = "my_function"`) or, as in the following examples, as a named list of functions. The names in the list are then used as labels for the models. As `specr()` extracts the estimate whose term equals the independent variable, the extraction function (`fun1`) has to return a term that is named like `x`. The independent variables should hence be numeric (factors result in terms such as "xB"; recode them as dummy variables instead). ## Weighted regression Survey data, for example, often require weights. If the weights are a variable in the data set, a custom function can pass them on to `lm()`. Here, we simulate a weight variable and compare unweighted and weighted models: ```{r} # Simulate weights set.seed(42) data_weights <- example_data %>% mutate(weight = runif(n(), min = .5, max = 1.5)) # Custom model function lm_weighted <- function(formula, data) { lm(formula, data = data, weights = weight) } # Setup specs specs <- setup(data = data_weights, x = c("x1", "x2"), y = "y1", model = list(unweighted = lm, weighted = lm_weighted), controls = c("c1", "c2")) # Run analysis and plot results results <- specr(specs) plot(results, choices = c("x", "model", "controls")) ``` For complex survey designs, wrappers around e.g. `survey::svyglm()` can be used in the same way. ## Survival models Survival outcomes consist of two variables (the survival time and the status). They can be used as dependent variable by writing them as `Surv(time, status)` (from the package `survival`), and Cox regressions can be fitted with a wrapper around `survival::coxph()`. With `exponentiate = TRUE`, the extraction function returns hazard ratios. As no effect corresponds to a hazard ratio of 1, we pass `null = 1` to `plot()` to highlight the significant specifications correctly. ```{r, eval = requireNamespace("tidyverse", quietly = TRUE) && requireNamespace("survival", quietly = TRUE)} library(survival) # Deaths in the data set `colon` from the package survival colon_death <- colon %>% filter(etype == 2) # Custom model function cox_model <- function(formula, data) { coxph(as.formula(formula), data = data) } # Setup specs specs <- setup(data = colon_death, x = c("age", "nodes"), y = "Surv(time, status)", model = list(cox = cox_model), controls = c("sex", "obstruct", "perfor"), fun1 = function(x) broom::tidy(x, exponentiate = TRUE, conf.int = TRUE)) # Run analysis and plot results (hazard ratios) results <- specr(specs) plot(results, null = 1, choices = c("x", "controls")) ``` ## Fixed effects models With `fixest::feols()`, fixed effects are added after a `|` in the formula. A wrapper can add the fixed effects to the formula that `setup()` creates and also set the type of standard errors: ```{r, eval = requireNamespace("fixest", quietly = TRUE)} # Custom model function (fixed effects for group3 and heteroskedasticity-robust # standard errors; with real panel data, they can be clustered e.g. via `cluster = ~id`) feols_fe <- function(formula, data) { fixest::feols(as.formula(paste(formula, "| group3")), vcov = "hetero", data = data) } # Setup specs specs <- setup(data = example_data, x = c("x1", "x2"), y = c("y1", "y2"), model = list(ols = lm, fixed_effects = feols_fe), controls = c("c1", "c2"), fun2 = NULL) # different fit indices for lm and feols # Run analysis and plot results results <- specr(specs) plot(results, choices = c("x", "y", "model", "controls")) ``` ## Multilevel models with p-values Models fitted with `lme4::lmer()` do not report p-values. Using `lmerTest::lmer()` instead, the estimates are reported with p-values (based on Satterthwaite's degrees of freedom), which are extracted with `broom.mixed::tidy()`. For more details on multilevel models, see `vignette("multilevel-models")`. ```{r, eval = FALSE} specs <- setup(data = example_data, x = c("x1", "x2"), y = c("y1", "y2"), model = list(lmer = lmerTest::lmer), controls = "c1", add_to_formula = "(1|group3)", fun1 = function(x) broom.mixed::tidy(x, conf.int = TRUE), fun2 = NULL) results <- specr(specs) plot(results) ``` There are of course many more types of analytical decisions. For some more examples, see some of the other vignettes. If you have a specific type of analytical decision that you don't know how to include, feel free to open an issue on github.