Simulating Data for a Poisson Regression

In a Poisson regression, the outcome \(Y_i\) represents a count, such as the number of units of a medication given intra-op for \(i^{th}\) subject.

We assume: \(Y_i\)∼ Poisson (\(\lambda_{i}\)), where \(\lambda_{i}\) is the expected event count (or event rate, depending on how exposure time, is handled).

For a simple Poisson regression with a binary exposure \(𝑋_𝑖\) , the model is: \(log(\lambda_{i})=\beta_0+\beta_1X_1\))

The model can equivalently be written as:\(\lambda_{i}=exp(\beta_0+\beta_1X_1)\)

Baseline Rate

When \(𝑋_i=0\), the subject belongs to the control group. Therefore: \(\lambda_0=exp(\beta_0)\)

We refer to \(\lambda_0\) as the baseline rate. For example, suppose the baseline rate is 2 events per subject-year. Then:\(\beta_0=log(2)\)=0.693

and therefore:\(e^{\beta_{0}} = 2\)

Thus, subjects in the control group have an expected event rate of 2 events per person-year.

Incidence Rate Ratio

The incidence rate ratio (IRR) compares the expected event rates between two exposure groups.

For a subject in exposed/treated group \(X_i=1\), \(\lambda_{1} =e^{\beta_0+\beta_1(1)}\)

The IRR is therefore:

\(IRR=\frac{\lambda_1}{\lambda_0}=e^{\beta_{1}}\)

and equivalently:\(\beta_1=log(IRR)\)

For example, if the target IRR is 1.5: \(\beta_1=log(1.5) =0.405\) 0.405

The expected event rate in the exposed/treated group is then:

\(\lambda_{1}=2*1.5=3\)

Thus, in this example:

Baseline rate = 2 events per person-year

Target IRR = 1.5

Expected rate among exposed subjects = 3 events per person-year

The data-generating model can therefore be written as:

\(\lambda_i=Baseline Rate*IRR^{X_i}\)

This equation will be used to generate the simulated outcome.

Simulating a Poisson Dataset

The exposure variable is generated using a Bernoulli distribution:

Each subject has a 50% probability of being exposed:

exposure = 0: unexposed

exposure = 1: exposed

The expected event rate for each subject is then calculated as: \(\lambda_i=Baseline Rate*IRR^{X_i}\)

For an unexposed/control subject:

\(\lambda_{ic}=2×1.5^{0} =2\)

For an exposed/treated subject:

\(\lambda{ie} =2×1.5^{1}=3\) Finally, the observed count is randomly generated from a Poisson distribution:

set.seed(123)

# Number of subjects
n <- 100

# Baseline event rate
baseline_rate <- 2

# Target incidence rate ratio
target_irr <- 1.5

# Binary exposure
exposure <- rbinom(
  n = n,
  size = 1,
  prob = 0.50
)

# Expected event rate for each subject
lambda <- baseline_rate * target_irr^exposure

# Simulate count outcome
count <- rpois(
  n = n,
  lambda = lambda
)

# Combine into a data frame
dat <- data.frame(
  count = count,
  exposure = exposure
)

head(dat)
##   count exposure
## 1     2        0
## 2     2        1
## 3     2        0
## 4     6        1
## 5     3        1
## 6     4        0

Examine the Simulated Data

The observed means will generally not be exactly 2 and 3 because of random variation. However, with a sufficiently large sample, they should be reasonably close to these values.

We can also examine the number of subjects in each group

Because exposure was generated with a probability of 0.50, we expect approximately half of the subjects to be exposed and half to be unexposed.

aggregate(
  count ~ exposure,
  data = dat,
  FUN = mean
)
##   exposure    count
## 1        0 2.094340
## 2        1 2.978723
table(dat$exposure)
## 
##  0  1 
## 53 47

Simulation-Based Sample Size Calculation

We can now use the simulated data-generating process to determine the sample size required to achieve 80% power. For illustration, we assume:

  1. Baseline event rate = 2 events per person-year. Note in the code below we are simulating: \(Y_i\)~ Poisson\((\lambda_i)\), where \(\lambda_i\) is the expected count per subject. If we later introduce person-time and an offset, then baseline_rate becomes appropriate, so more approprite wording at this point is baseline mean not baselien rate.

  2. Target IRR = 1.5

  3. 50% of subjects are exposed

  4. Two-sided significance level = 0.05

  5. Target power = 80%

The goal is to find the smallest sample size N for which the probability of detecting the target IRR is approximately 80%.

set.seed(123)

# Set Simulation parameters

baseline_rate <- 2 
target_irr <- 1.5
exposure_prob <- 0.50

alpha <- 0.05
target_power <- 0.80

# Number of simulated studies
nsimul <- 500 # eventually at east 2000

# Function to Simulate one Study

#The function returns TRUE if the simulated study detects a statistically #significant association and FALSE otherwise

simulate_study <- function(n,
                           baseline_rate,
                           target_irr,
                           exposure_prob,
                           alpha = 0.05) {

  # Generate exposure
  exposure <- rbinom(
    n = n,
    size = 1,
    prob = exposure_prob
  )

  # Expected event rate
  lambda <- baseline_rate *
    target_irr^exposure

  # Generate Poisson outcome
  count <- rpois(
    n = n,
    lambda = lambda
  )

  # Fit Poisson regression
  fit <- glm(
    count ~ exposure,
    family = poisson(link = "log")
  )

  # Extract p-value for exposure
  p_value <- summary(fit)$coefficients[
    "exposure",
    "Pr(>|z|)"
  ]

  # Determine whether the result is statistically significant
  p_value < alpha
}

## Estimate Power of one fixed Sample Size, by repeating the simulation many times for a particular sample size.
estimate_power <- function(n,
                           baseline_rate,
                           target_irr,
                           exposure_prob,
                           nsim = nsimul,
                           alpha = 0.05) {

  results <- replicate(
    nsim,
    simulate_study(
      n = n,
      baseline_rate = baseline_rate,
      target_irr = target_irr,
      exposure_prob = exposure_prob,
      alpha = alpha
    )
  )

  mean(results)
}

## Evaluate a range of Sam-ple Sizes
sample_sizes <- seq(
  from = 50,
  to = 500,
  by = 5
)

power_results <- sapply(
  sample_sizes,
  estimate_power,
  baseline_rate = baseline_rate,
  target_irr = target_irr,
  exposure_prob = exposure_prob,
  nsim = nsimul,
  alpha = alpha
)

power_results
##  [1] 0.582 0.614 0.678 0.722 0.776 0.748 0.822 0.852 0.856 0.858 0.870 0.906
## [13] 0.932 0.918 0.940 0.916 0.942 0.958 0.976 0.978 0.976 0.976 0.978 0.988
## [25] 0.986 0.978 0.992 0.994 0.992 0.988 0.992 0.994 0.996 0.996 0.998 0.998
## [37] 0.996 1.000 1.000 0.998 0.996 1.000 0.998 1.000 1.000 1.000 1.000 1.000
## [49] 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## [61] 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## [73] 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## [85] 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Combine results in data frame
results <- data.frame(
  sample_size = sample_sizes,
  power = power_results
)
#require(KableExtra)
#kbl(results)
results
##    sample_size power
## 1           50 0.582
## 2           55 0.614
## 3           60 0.678
## 4           65 0.722
## 5           70 0.776
## 6           75 0.748
## 7           80 0.822
## 8           85 0.852
## 9           90 0.856
## 10          95 0.858
## 11         100 0.870
## 12         105 0.906
## 13         110 0.932
## 14         115 0.918
## 15         120 0.940
## 16         125 0.916
## 17         130 0.942
## 18         135 0.958
## 19         140 0.976
## 20         145 0.978
## 21         150 0.976
## 22         155 0.976
## 23         160 0.978
## 24         165 0.988
## 25         170 0.986
## 26         175 0.978
## 27         180 0.992
## 28         185 0.994
## 29         190 0.992
## 30         195 0.988
## 31         200 0.992
## 32         205 0.994
## 33         210 0.996
## 34         215 0.996
## 35         220 0.998
## 36         225 0.998
## 37         230 0.996
## 38         235 1.000
## 39         240 1.000
## 40         245 0.998
## 41         250 0.996
## 42         255 1.000
## 43         260 0.998
## 44         265 1.000
## 45         270 1.000
## 46         275 1.000
## 47         280 1.000
## 48         285 1.000
## 49         290 1.000
## 50         295 1.000
## 51         300 1.000
## 52         305 1.000
## 53         310 1.000
## 54         315 1.000
## 55         320 1.000
## 56         325 1.000
## 57         330 1.000
## 58         335 1.000
## 59         340 1.000
## 60         345 1.000
## 61         350 1.000
## 62         355 1.000
## 63         360 1.000
## 64         365 1.000
## 65         370 1.000
## 66         375 1.000
## 67         380 1.000
## 68         385 1.000
## 69         390 1.000
## 70         395 1.000
## 71         400 1.000
## 72         405 1.000
## 73         410 1.000
## 74         415 1.000
## 75         420 1.000
## 76         425 1.000
## 77         430 1.000
## 78         435 1.000
## 79         440 1.000
## 80         445 1.000
## 81         450 1.000
## 82         455 1.000
## 83         460 1.000
## 84         465 1.000
## 85         470 1.000
## 86         475 1.000
## 87         480 1.000
## 88         485 1.000
## 89         490 1.000
## 90         495 1.000
## 91         500 1.000

With an IRR of 1.5 and a baseline expected count of 2, the effect is relatively large, so high power can be achieved with a relatively small sample size.

Identify the Sample Size Achieving 80% Power

Below code finds first evaluated sample size for which estimated power is at least 80%

required_n <- results$sample_size[
  which(results$power >= target_power)[1]
]

required_power <- results$power[
  which(results$power >= target_power)[1]
]

required_n
## [1] 80
required_power
## [1] 0.822

Power Curve

plot(
  results$sample_size,
  results$power,
  type = "b",
  pch = 19,
  ylim = c(0, 1),
  xlab = "Sample size",
  ylab = "Estimated power",
  main = "Simulation-Based Power for Poisson Regression"
)

# Target power
abline(
  h = target_power,
  col = "red",
  lty = 2
)

# Required sample size
abline(
  v = required_n,
  col = "blue",
  lty = 2
)