Using the eesim package

Generating exposure data

The first task of the package is generating exposure data. This can be done with the sim_exposure function. In this function, the user can specify whether he or she would like to generate exposure data that is binary or continuous (exposure_type). For continuous exposure data, the user must specify the mean (central) and standard deviation (sd) of the exposure data. For example, the following call simulates a dataframe of exposure data for an exposure that is normally distributed, with a mean value of 50, a standard deviation of 5, and no long-term or seasonal trends:

x_cont <- sim_exposure(n = 1000, central = 50, sd = 5, exposure_type = "continuous") 
x_cont %>% slice(1:5)

##         date        x
## 1 2001-01-01 49.26604
## 2 2001-01-02 42.91050
## 3 2001-01-03 52.53690
## 4 2001-01-04 51.53205
## 5 2001-01-05 52.42711

ggplot(x_cont, aes(x = date, y = x)) + geom_point(alpha = 0.2) + 
  theme_classic()

You can plot a calendar plot of this simulated exposure time series using the calendar_plot function that comes with the package. Within this function, the type of data (“continuous” or “discrete”) must be specified:

calendar_plot(x_cont, type = "continuous")

You can similarly use the sim_exposure function to simulate a binary exposure (e.g., occurence of an extreme event). For binary exposure data, the central argument of sim_exposure must also be expressed, but in this case it gives the probability of exposure on a study day:

x_bin <- sim_exposure(n = 1000, central = 0.05, exposure_type = "binary")
x_bin %>% slice(1:5)

##         date x
## 1 2001-01-01 0
## 2 2001-01-02 0
## 3 2001-01-03 1
## 4 2001-01-04 0
## 5 2001-01-05 0

Again, the calendar_plot function can be used to visualize the generated time series. In the case of binary exposure data, the labels to be used in the legend for each outcome level must also be specified using the labels argument:

calendar_plot(x_bin, type = "discrete", labels = c("Not exposed", "Exposed"))

So far, these sim_exposure calls have been used to simulate basic exposure data, without long-term or seasonal trends. However, for environmental epidemiology applications, exposure data often has a seasonal trend and / or long-term trend, and these temporal trends can serve as confounders in assessing the association between time-varying environmental exposures and health outcomes. The sim_exposure function therefore includes options to generate exposure data with long-term and seasonal trends relevant to environmental time series studies, through the trend argument.

The default for sim_exposure is to simulate the exposure data without a time trend (trend = "no trend"). However, we have also built in several time trends from which a user can to choose to simulate exposure data with a time trend, either seasonal or long-term or both, based on trend patterns used in a simulation study of case-crossover studies as a method of controlling for seasonal and long-term trends in environmental epidemiology studies (Bateson and Schwartz 1999). These trend patterns differ slightly depending on whether the user is simulating binary or continuous data. Below are plots of the built-in trends for continuous exposure data from which the user may choose.

You can use the amp argument to adjust the seasonal trend in any of the patterns with a seasonal trend. For example, here are plots of trends using tren = "cos1linear" for different values of amp:

The long-term trend in expected values can be changed in a similar way with the exposure_trend argument in eesim.

Here is an example of generating continuous exposure data with a “cos1linear” trend for an exposure with a mean value of 50 and a standard deviation of 10:

testexp <- sim_exposure(n = 365 * 3, central = 50, sd = 10, trend = "cos1linear",
                        exposure_type = "continuous")
a <- ggplot(testexp, aes(x = date, y = x)) +  
  geom_point(alpha = 0.5, size = 0.8) + 
  coord_cartesian(ylim = c(0,110)) + 
  labs(title = "Exposure with a 'cos1linear' trend", x = "Date", y="Exposure") + 
  theme_classic()
b <- calendar_plot(testexp, type = "continuous") + 
  ggtitle("Calendar plot of simulated exposure data") + 
  theme(legend.position = "bottom")
grid.arrange(a, b, ncol = 1)

Here is an example of changing the seasonal trend by changing the value for amp (the default value is 0.6) to simulate exposure data for an exposure with a smaller seasonal trend and with higher exposures typical in the summer than the winter:

small_amp <- sim_exposure(n = 365 * 3, central = 50, sd = 10, trend = "cos1linear",
                        amp = -0.3, exposure_type = "continuous")
a <- ggplot(small_amp, aes(x = date, y = x)) +  
  geom_point(alpha = 0.5, size = 0.8) + 
  coord_cartesian(ylim = c(0,110)) + 
  labs(title = "Exposure with a 'cos1linear' trend", x = "Date", y="Exposure") + 
  theme_classic()
b <- calendar_plot(small_amp, type = "continuous") + 
  ggtitle("Calendar plot of simulated exposure data") + 
  theme(legend.position = "bottom")
grid.arrange(a, b, ncol = 1)

The trend options are similar for binary exposure, but exclude “curvilinear” and “cos1linear”. Further, binary exposures can also be simulated using a “monthly” trend (trend = "monthly"), in which the probability of exposure can vary by month. When using this “monthly” trend option, the cental argument to sim_exposure should include a vector with 12 separate probabilities (the first is for January, the second for February, etc.) rather than a single probability. Here is an example of generating binary exposure data with a monthly trend, starting from June 1, 2002, with higher probability of the exposure in summer months than in winter months:

testbin <- sim_exposure(n=1000, central = c(.05, .05, .1, .2, .4, .4, .5, .7, .5, .2, .1, .05),
                        trend = "monthly", exposure_type = "binary", 
                        start.date = "2002-06-01")
a <- testbin %>% 
  mutate(x = factor(x, levels = c(0, 1), labels = c("Not exposed", "Exposed"))) %>% 
  ggplot(aes(x = date, y = x)) + 
  geom_jitter(alpha = 0.5, size = 0.7, fill = NA, width = 0, height = 0.1) + 
  theme_classic() + 
  labs(x = "Date", y = "Exposure")
b <- calendar_plot(testbin, type = "discrete", labels = c("Not exposed", "Exposed")) + 
  ggtitle("Calendar plot of simulated exposure data") + 
  theme(legend.position = "bottom")
grid.arrange(a, b, ncol = 1)

The sim_exposure function works by first calculating the expected exposure on any date in the simulated time series (figure below, left). This expected value is a mean for a continuous exposure and a probability for a binary exposure. The sim_exposure function then draws random values from the appropriate distribution (normal distribution for a continuous exposure, binomial distribution for a binary exposure) based on this day-specific expected exposure value and, in the case of continuous exposure, the standard deviation of the exposure (figure below, right). For continuous exposure data, the standard deviation specified in the call to eesim should measure the standard deviation of each point from its expected value (i.e., from the expected line shown on the left below), not the overall standard deviation of exposure values across all days in the simulated data.

Later in this vignette, we show how you can further customize this step of generating exposure data through the use of a user-created function, allowing extensive further flexibility in simulating exposure data.

Generating outcome data

Next, the sim_outcome function simulates outcome data. The health data can have an underlying seasonal and / or long term trend in its baseline value, and then that baseline is adjusted for the risk associated with exposure, based on the generated exposure data for that day. The baseline outcome count for a given day (\(B_t\)) are based on a user-specified trend and user-specified average outcome per day over the simulated time period. Further, the expected outcome count on a given day is adjusted for exposure-related risk through a user-specified relative rate per unit increase in exposure (\(RR\)) and the simulated exposure for that day (\(X_t\)). The eesim function then uses the following equation to calculate the expected outcome count (\(\lambda\)) on a given day in the simulated time series, based on the expected baseline rate and exposure-related risk for that day:

\[ \log(\lambda_t) = \log{(B_t)}+\log{(RR)}*X_t \]

For a binary outcome, the baseline count on a given day (\(B_t\)) is the expected outcome count for the day if there is not an event (e.g., in a heat wave study, a non-heat wave day). For a continuous exposure, the baseline count on a given day (\(B_t\)) is the expected outcome count for the day if exposure is at its mean value.

Once the expected count (\(\lambda_t\)) on each day of the simulated time series is calculated using this equation, the simulated count on each day is drawn as a random variable from a Poisson distribution with mean \(\lambda_t\). Later we describe how customization can be used to simulate output counts in other ways.

Here is an example of generating health outcome data with an upward linear trend using exposure data with a “cos1” trend. In this case, there is a steady increase in the baseline outcome count over time, as well as a seasonal trend linked to the risk associated with the seasonally-varying exposure:

testexp2 <- sim_exposure(n = 1000, central = 100, sd = 10, trend = "cos1",
                         exposure_type = "continuous")
testout <- sim_outcome(exposure = testexp2, average_outcome = 22,
                       trend = "linear", rr = 1.01)

Here are plots of the resulting output:

As with the exposure simulation step, this step can also be extensively customized by using a user-created function. This customization will be demonstrated in a later section of the vignette.

Fitting models

Next, the eesim package uses this process to generate many simulated data sets and then to fit statistical models to these generated datasets. This step allows tests of model performance. You must create an R function that fits the model you’d like to fit to the simulated dataset. This function needs to follow certain input / output rules to work correctly in the eesim framework. First, it must input the simulated dataframe with the argument df. When writing the function, you should assume that this simulated dataframe has at least the columns date (in a Date format), x (numeric class, this gives a daily value for exposure, with 0 for unexposed and 1 for exposed in the case of binary exposure), and outcome (non-negative integer, this gives the simulated outcome count each day). Other arguments can also be passed to this function if desired. The function should fit a desired model to the simulated dataframe and then should return a numeric vector of length 6 with values, in this order, for the log relative risk point estimate from the model (Estimate), the standard error for this point estimate (Std. Error), the t-statistic for a hypothesis test with the null hypothesis that this estimate is zero (t value), a p-value for that test (Pr(>|t|)), and the lower and upper 95% confidence intervals for the point estimate (2.5 % and 97.5 %).

The spline_mod function that comes with the package is an example of such a function. In this case, the function fits a GLM to the simulated data, with a natural cubic spline used to control for long-term and seasonal trends in mortality. The function, in addition to inputing the dataframe of simulated data (df), also allows an argument to use to set the smoothness of the time spline (df_year). Since the function is included in the package, you can see its code by running the bare function name at the console:

spline_mod

## function (df, df_year = 7) 
## {
##     dgrs_free <- df_year * as.numeric(diff(df[c(1, nrow(df)), 
##         "date"]))/365.4
##     df$time <- scale(df$date, center = TRUE, scale = FALSE)
##     mod <- stats::glm(outcome ~ x + splines::ns(time, round(dgrs_free)), 
##         data = df, family = stats::quasipoisson(link = "log"))
##     out_1 <- summary(mod)$coef[2, ]
##     out_2 <- stats::confint.default(mod)[2, ]
##     out <- c(out_1, out_2)
##     return(out)
## }
## <bytecode: 0x7f974761c7b8>
## <environment: namespace:eesim>

Here are examples of applying this function to a simulated dataframe:

# Create simulated data
sims <- create_sims(n_reps = 10, n = 100, central = 100, sd = 10,
             exposure_type="continuous", exposure_trend = "cos1",
             exposure_amp = .6, average_outcome = 22,
             outcome_trend = "no trend", outcome_amp = .6, rr = 1.01)
head(sims[[1]])

##         date        x outcome
## 1 2000-01-01 116.5873      22
## 2 2000-01-02 106.9395      29
## 3 2000-01-03 126.3188      21
## 4 2000-01-04 124.0723      21
## 5 2000-01-05 104.3558      22
## 6 2000-01-06 123.9815      32

# Apply `spline_mod` to the data
spline_mod(df = sims[[1]])

##    Estimate  Std. Error     t value    Pr(>|t|)       2.5 %      97.5 % 
## 0.007501871 0.002222082 3.376055719 0.001062752 0.003146671 0.011857071

spline_mod(df = sims[[1]], df_year = 6)

##    Estimate  Std. Error     t value    Pr(>|t|)       2.5 %      97.5 % 
## 0.007501871 0.002222082 3.376055719 0.001062752 0.003146671 0.011857071

The format_function function can be used within the modeling function to get the output in the correct format if running a GLM or similar model.

Once you’ve created the function, you can input it in a call to eesim using the custom_model argument. You can pass any additional arguments (df_year in our example) through to the function using the custom_model_args argument. This argument takes a list with the argument name and value for each argument you wish to pass to the modeling function. For example, the following call passes the spline_mod function shown above as the function to use for modeling the simulated data as well as a value for its df_year argument:

ex_sim2 <- eesim(n_reps = 100, n = 365 * 5, 
                 central = c(0, 0, 0, 0, 0, 0.05, 0.12, 0.02, 0, 0, 0, 0),
                 exposure_type = "binary", exposure_trend = "monthly",
                 exposure_amp = -.6, average_outcome = 50,
                 outcome_trend = "cos1", outcome_amp = 0.2, 
                 rr = 1.05, start.date = "1996-01-01",
                 custom_model = spline_mod, custom_model_args = list(df_year = 7))

The eesim function does this by applying the modeling function across all simulated datasets using a function called fit_mods. If you’d like, you can run that function independently. The fit_mods function outputs a data frame with estimates of the log relative risk, p-values, and upper and lower 95% confidence bounds for each simulated data set.

Here is an example of fitting the spline model coded in the spline_mod function, with 7 degrees of freedom per year used to model long-term and seasonal trends (df_year = 7 passed in a list to the model with the custom_model_args argument):

fits <- fit_mods(data = sims, custom_model = spline_mod, 
                 custom_model_args = list(df_year = 7))
fits

##       Estimate   Std.Error  t.value      p.value    lower_ci   upper_ci
## 1  0.007501871 0.002222082 3.376056 1.062752e-03 0.003146671 0.01185707
## 2  0.011958167 0.002215918 5.396484 4.893943e-07 0.007615048 0.01630129
## 3  0.009081710 0.002624721 3.460067 8.073299e-04 0.003937352 0.01422607
## 4  0.007669370 0.002354554 3.257249 1.555672e-03 0.003054528 0.01228421
## 5  0.010571590 0.001877903 5.629464 1.791011e-07 0.006890967 0.01425221
## 6  0.009580041 0.002354242 4.069267 9.680293e-05 0.004965810 0.01419427
## 7  0.010950804 0.002332539 4.694800 8.870286e-06 0.006379111 0.01552250
## 8  0.007214138 0.002456299 2.936995 4.148696e-03 0.002399881 0.01202839
## 9  0.007301189 0.002295396 3.180797 1.978360e-03 0.002802295 0.01180008
## 10 0.010476794 0.002291069 4.572886 1.435180e-05 0.005986383 0.01496721

As a note, the output of the fit_mods function is the output given as the second element of the list returned by a call to eesim.

Evaluating the models

The final step of the eesim function is to evaluate model performance across all simulations with several different measures. Within the eesim function, the check_sims function takes as inputs the true relative risk as well as the results from fitting the modeling function to all the simulations using the fim_mods function. It returns values for the mean effect estimate (log relative risk) and relative risk estimates across all simulated data sets, variance of the estimates of beta, the mean of the variances of each estimated log relative risk, the relative bias of the mean of the log relative risks, the percent coverage confidence intervals of the true log relative risk, and the power of the test at the 5% significance level (see the table near the beginning of the vignette).

Here is an example of the use of the check_sims function:

check_sims(fits, true_rr = 1.01)

##      beta_hat      rr_hat var_across_betas mean_beta_var percent_bias
## 1 0.009230567 0.009230567     3.012105e-06  5.334198e-06     7.233563
##   coverage power
## 1        1     1

In a run of eesim, this output is given in the third element of the returned list.

Internally, the functions used for this model assessment are:

• beta_bias
• beta_var
• coverage_beta
• mean_beta
• power_beta

A few more details about how some of these assessments are measured are given below.

Variance across estimated log relative risk and mean variance of estimates

Two values are measured by the beta_var function. First, the variance across all estimates of log relative risk is measured across all the simulations, using the equation:

\[ \text{variance of estimates} = E\left[\left(\hat{\beta_i} - \frac{1}{n}\sum_{i = 1}^n{\hat{\beta_i}}\right)\right] \]

where \(\hat{\beta_i}\) is the estimated log relative risk for a single simulation out of \(n\) total simulations and \(E\) represents the expected value.

Second, the function measures the mean value of the variance estimated for \(\hat{\beta}\) for each simulation:

\[ \text{mean of estimate variances} = \frac{1}{n}\sum_{i = 1}^{n}{var(\hat{\beta_i})} \]

where \(var(\hat{\beta_i})\) is the estimated variance of the estimated log relative risk for a single simulation out of \(n\) total simulations.

Generating power curves

The other main functionality of the eesim package is to run through simulations under varying data generating scenarios to estimate expected power of an analysis under different scenarios. For example, you can explore how expected power varies for different expected effect sizes (relative risk of the outcome associated with a change in exposure) or for different average daily number of outcomes. This is run using the power_calc function in the package.

The power_calc function allows you to put in varying values for one of the following three specifications in the simulations:

• Relative risk (rr): How you strongly expect the exposure and outcome to be associated
• Number of days in the study (n): How long you expect the study to last (or how many days of historical data you expect to be able to collect)
• Average daily count of outcomes (average_outcome): For the outcome of interest, how common it is on average on days in the study (this will usually be strongly associated with the size of the population being studied)

For whichever of these you choose to vary, you can specify different values to test. The power_calc function then loops through those values and runs eesim for each of them. From this, it can estimate the power for each value of the varying parameter.

Here is an example of running a power calculation with varying number of days in the study (n). The values argument is used to specify different values of n we would like to test (here, it’s testing power for studies with daily data for between 1 and 21 years):

pow <- power_calc(varying = "n", values = floor(365.25 * seq(1, 21, by = 5)), n_reps = 20,
                  central = 100, sd = 10, rr = 1.001, exposure_type = "continuous",
                  exposure_trend = "cos1", exposure_amp = .6, average_outcome = 22,
                  outcome_trend = "no trend", outcome_amp = .6,
                  custom_model = spline_mod, plot = TRUE)

## This function may take a minute or two to run, especially with
## lots of replications (`n_reps`) or options for `values`.

This call returns a dataframe with the estimated power for each of the values of n tested:

pow

##   values power
## 1    365  0.05
## 2   2191  0.75
## 3   4017  0.95
## 4   5844  0.95
## 5   7670  1.00

Because the argument plot is set to TRUE, it also generates a power curve plot as a side effect, as shown above.

Here is another example, but this time we assume that the study will have 4,000 days of daily data, but we explore estimated power as the average daily outcome count varies:

pow2 <- power_calc(varying = "average_outcome", values = c(1, 5, 10, 20, 30, 40),
                   n_reps = 20,
                   central = 100, sd = 10, rr = 1.001, exposure_type = "continuous",
                   exposure_trend = "cos1", exposure_amp = .6, n = 4000,
                   outcome_trend = "no trend", outcome_amp = .6,
                   custom_model = spline_mod, plot = TRUE)

## This function may take a minute or two to run, especially with
## lots of replications (`n_reps`) or options for `values`.

pow2

##   values power
## 1      1  0.10
## 2      5  0.25
## 3     10  0.60
## 4     20  0.85
## 5     30  0.90
## 6     40  1.00

Because these power curves and calculations are based on simulated data, there will be some randomness to results. Curves will be smoother as more simulations are used for each run (n_reps), although this will also increase the time needed to run the simulation.

Customizing the exposure trend

To take advantage of any of the customization options, you need to write a function that follows certain input and output (i.e., interface) rules. First, you can use a custom function for the underlying trend in expected exposure. This function must take the inputs:

• n (the number of days to simulate)
• Either mean for a continuous exposure (the average value of the outcome) or prob for a binary exposure (the average probability of exposure)

The function can take any other additional inputs, as well, but any such extra arguments (as well as mean) will need to be input to the eesim function in a list for the cust_expdraw_args argument (example below). The value for n will pass through directly from the n value specified for the call to eesim. The function must output a numeric vector that gives the simulated exposure values for each day in the simulated data.

For example, the following function creates a custom exposure trend with a long-term and seasonal trend, similar to trends available through the default package options. However, this function specifies a minimum value that the exposure trend cannot fall below– if the base exposure value is every set below this minimum within the algorithm, the value is reset to the minimum before the final values are output. This function can be useful in cases where the exposure cannot fall below a certain value (for example, a pollution concentration could not be lower than 0). This custom exposure function can also be used to customize how values are simulated from the expected exposure on each day (based on the expected distribution of the exposure). In the case of the example ozone concentration data from Chicago shown earlier in this vignette, we may want to simulate exposure based on the assumption that the square root of exposure is normally distributed, which will prevent negative values and may also help to simulate occasional very high values.

above_min_trend <- function(n, mean, sd_of_sqrt, minimum = 0){
  day <- c(1:n)
  
  ## Calculate a baseline exposure for each day
  base <- mean + -10 * cos(2 * pi * (day / 365))
  base[base < minimum] <- minimum            # Reset any values below 0 to 0
  
  ## Simulate exposure values from the baseline
  sqrt_base <- sqrt(base)                   # Transform to square root
  sqrt_sim <- rnorm(n, mean = sqrt_base, sd = sd_of_sqrt)
  sqrt_sim ^ 2                              # Transform back
}

Here is an example of running this custom exposure simulation function over 5 years, with a smooth line added to the plot to help show the seasonal trend included:

above_min_trend(n = 365.25 * 5, mean = 20, minimum = 0, sd_of_sqrt = 0.9) %>% 
  tbl_df() %>% 
  mutate(day = 1:n()) %>% 
  ggplot(aes(x = day, y = value)) + 
  geom_point(alpha = 0.5, size = 0.8) + 
  theme_classic() + 
  geom_smooth(se = FALSE, span = 0.1, method = "loess", color = "red")

You can then pass this custom function into the eesim function using the cust_exp_func argument. The value for n input to the custom function will be the value you input to eesim for n. For any other arguments you want to pass to the function (in the function we just created, you’ll want to pass values for mean, minimum, and sd_of_sqrt), you can include specifications for these as a list for the cust_exp_args argument of eesim. For example, the following call would run a simulation using this custom function for exposure:

ex_sim2 <- eesim(n_reps = 1, n = round(365.25 * 5), 
                 exposure_type = "continuous",
                 cust_exp_func = above_min_trend,
                 cust_exp_args = list(mean = 20, minimum = 0, sd_of_sqrt = 0.9),
                 average_outcome = 50, rr = 1.01, 
                 custom_model = spline_mod, custom_model_args = list(df_year = 7))

Customizing the outcome simulation

There are three ways to customize the simulated outcome data: creating a custom baseline for outcome values, customizing the relationship between outcome and exposure, and, as with the exposure values, customizing the randomization of the outcome values.

The outcome baseline (\(B_t\)) is comprised of the values the user expects the outcomes to have on each day of the simulated dataset without risk associated with the exposure factored in. The user may write a function to specify the trend of the baseline, then use it as an input in sim_outcome or eesim. Here is an example of creating a custom baseline function and using it in the eesim function:

custombase <- function(n, slope, intercept){
  day <- c(1:n)
  baseline <- day * slope + intercept
  return(baseline)
}

#Example:
custombase(n=5, slope = .3, intercept = 55)

## [1] 55.3 55.6 55.9 56.2 56.5

ex_sim3 <- eesim(n_reps = 3, n = 1000, central = 100, sd = 10,
                exposure_type = "continuous", exposure_trend = "cos1",
                exposure_amp = .6, average_outcome = 22, rr = 1.01, 
                cust_base_func = custombase,
                cust_base_args = list(n=1000, slope = 5, intercept = 12),
                custom_model = spline_mod, custom_model_args = list(df_year = 2))

## This function may take a minute or two to run, especially if you
## are creating lots of replications (`n_reps`).

ggplot(ex_sim3$simulated_datasets[[1]], aes(x=date, y=outcome))+ geom_point() + geom_point(alpha = 0.5, size = 0.8) + 
  theme_classic() + 
  geom_smooth(se = FALSE, span = 0.1, method = "loess", color = "red")

The second way of customizing the outcome simulation is to use a custom function to incorporate the added risk from the exposure when calculating the expected daily outcome count for a day, \(\lambda_t\), from the inputs of exposure (\(X_t\)) and outcome baseline (\(B_t\)) for the day. Here is an example of creating a custom lambda, meaning a custom function relating relative risk and exposure to outcomes, and using it in eesim with the custom baseline function created above. The custom lambda function must input arguments exposure, rr, and baseline and output a vector of lambda values.

customlambda <- function(exposure, rr, constant, baseline){
  log_lambda <- log(baseline) + log(rr) * exposure + constant
  lambda <- exp(log_lambda)
  return(lambda)
}

ex_sim4 <- eesim(n_reps = 3, n = 1000, central = 100, sd = 10,
                exposure_type = "continuous", exposure_trend = "cos1",
                exposure_amp = .6, average_outcome = 22, rr = 1.01, 
                cust_base_func = custombase,
                cust_base_args = list(n=1000, slope = .5, intercept = 12),
                cust_lambda_func = customlambda, cust_lambda_args = list(constant=10),
                custom_model = spline_mod, custom_model_args = list(df_year = 2))

## This function may take a minute or two to run, especially if you
## are creating lots of replications (`n_reps`).

The third way to customize the outcome simulation is to customize the randomization of the outcome values from the trend created by relating the baseline outcomes and the exposure (what we have called lambda). When the cust_outdraw argument is not specified in the eesim function, the function draws outcome values from a Poisson distribution with mean lambda. A custom function for outcome draws must input values called n and lambda, and any other arguments must be included in the cust_outdraw_args argument. Here is an example of using the custom functions to specify a negative binomial distribution for outcome randomization:

custnbinom <- function(n, lambda, prob){
  out <- rnbinom(n=n, size=lambda, prob=prob)
  return(out)
}

ex_sim5 <- eesim(n_reps = 3, n = 1000, central = 100, sd = 10,
                exposure_type = "continuous", exposure_trend = "cos1",
                exposure_amp = .6, average_outcome = 22, rr = 1.01, 
                cust_base_func = custombase,
                cust_base_args = list(n=1000, slope = .5, intercept = 12),
                cust_lambda_func = customlambda, cust_lambda_args = list(constant=10),
                cust_outdraw = custnbinom, cust_outdraw_args = list(prob=.3), 
                custom_model = spline_mod, custom_model_args = list(df_year = 2))

## This function may take a minute or two to run, especially if you
## are creating lots of replications (`n_reps`).