Skip to contents

Simulate from count model given covariates

Source: R/outcome_models.R
outcome_count.Rd

Simulate from count model with intensity $$\lambda = \text{exposure-time}\exp(\text{par}^\top X)$$ where \(X\) is the design matrix specified by the formula

Usage

outcome_count(
  data,
  mean = NULL,
  par = NULL,
  outcome.name = "y",
  exposure = 1,
  remove = c("id", "num"),
  zero.inflation = NULL,
  overdispersion = NULL,
  ...
)

Arguments

data

(data.table) Covariate data, usually the output of the covariate model of a Trial object.

mean

formula specifying design from 'data' or a function that maps x to the mean value. If NULL all main-effects of the covariates will be used

par

(numeric) Regression coefficients (default zero). Can be given as a named list corresponding to the column names of model.matrix

outcome.name

Name of outcome variable ("y")

exposure

Exposure times. Either a scalar, vector or function.

remove

variables that will be removed from input data (if formula is not specified)

zero.inflation

vector of probabilities or a function of the covariates 'x' including an extra column 'rate' with the rate parameter.

overdispersion

variance of gamma-frailty either given as a numeric vector or a function of the covariates 'x' with an extra column 'rate' holding the rate parameter 'rate'

...

Additional arguments passed to mean and exposure function

See also

outcome_binary outcome_continuous outcome_lp

Examples

covariates <- function(n) data.frame(a = rbinom(n, 1, 0.5), x = rnorm(n))
trial <- Trial$new(covariates = covariates, outcome = outcome_count)
trial$args_model(
  mean = ~ 1 + a + x,
  par = c(2.5, 0.65, 0),
  overdispersion = 1 / 2,
  exposure = 2 # identical exposure time for all subjects
)
est <- function(data) {
  glm(y ~ a + x + offset(log(exposure)), data, family = poisson())
}
trial$simulate(1e4) |> est()
#> 
#> Call:  glm(formula = y ~ a + x + offset(log(exposure)), family = poisson(), 
#>     data = data)
#> 
#> Coefficients:
#> (Intercept)            a            x  
#>    2.506904     0.650294    -0.004953  
#> 
#> Degrees of Freedom: 9999 Total (i.e. Null);  9997 Residual
#> Null Deviance:	    211800 
#> Residual Deviance: 175900 	AIC: 226600

# intercept + coef for x default to 0 and regression coef for a takes
# the provided value
trial$simulate(1e4, par = c(a = 0.65)) |> est()
#> 
#> Call:  glm(formula = y ~ a + x + offset(log(exposure)), family = poisson(), 
#>     data = data)
#> 
#> Coefficients:
#> (Intercept)            a            x  
#>   -0.023634     0.686864     0.005638  
#> 
#> Degrees of Freedom: 9999 Total (i.e. Null);  9997 Residual
#> Null Deviance:	    27500 
#> Residual Deviance: 24250 	AIC: 47990
# trial$simulate(1e4, mean = ~ 1 + a, par = c("(Intercept)" = 1))

# define mean model that directly works on whole covariate data, incl id and
# num columns
trial$simulate(1e4, mean = \(x) with(x, exp(1 + 0.5 * a))) |>
  est()
#> 
#> Call:  glm(formula = y ~ a + x + offset(log(exposure)), family = poisson(), 
#>     data = data)
#> 
#> Coefficients:
#> (Intercept)            a            x  
#>    1.016003     0.503870     0.006757  
#> 
#> Degrees of Freedom: 9999 Total (i.e. Null);  9997 Residual
#> Null Deviance:	    49180 
#> Residual Deviance: 44670 	AIC: 78870

# treatment-dependent exposure times
trial$simulate(1e4, exposure = function(dd) 1 - 0.5 * dd$a) |>
  head()
#>       id     a          x   num     y exposure
#>    <num> <int>      <num> <num> <num>    <num>
#> 1:     1     1 -1.4883673     0     2      0.5
#> 2:     2     1 -0.9731949     0    20      0.5
#> 3:     3     0  0.1568903     0    15      1.0
#> 4:     4     0  0.7359493     0     3      1.0
#> 5:     5     0  2.0626036     0    25      1.0
#> 6:     6     0 -1.5535437     0    21      1.0