Simulating data for exploratory analyses

First, we simulate some data to be able to enable showcasing of the functionalities. For this we use the glm_data() function from the package, where the user can specify an expression alongside variables and a family of the response to then simulate a response from a GLM with linear predictor given by the expression provided.

n <- 1000
b0 <- 1
b1 <- 3
b2 <- 2

# Simulate data with a non-linear effect
dat_treat <- glm_data(
  Y ~ b0+b1*sin(W)^2+b2*A,
  W = runif(n, min = -2, max = 2),
  A = rbinom(n, 1, prob = 1/2),
  family = gaussian() # Default value
)

Fitting rctglm() without prognostic covariate adjustment

The rctglm() function estimates any specified estimand using plug-in estimation for randomised clinical trials and estimates the variance using the influence function of the marginal effect estimand.

The interface of rctglm() is similar to that of the stats::glm() function but with an added mandatory specification of

• The randomisation variable in data, usually being the (name of) the treatment variable
• The randomisation ratio - the probability of being allocated to group 1 (rather than 0)
  • As a default, a ratio of 1’s in data is used
• An estimand function
  • As a default, the function takes the average treatment effect (ATE) as the estimand

Thus, we can estimate the ATE by simply writing the below:

Note that as a default, verbose = 2, meaning that information about the algorithm is printed to the console. However, here we suppress this behavior. See more in vignette("model-fit").

ate <- rctglm(formula = Y ~ A * W,
              exposure_indicator = A,
              exposure_prob = 1/2,
              data = dat_treat,
              family = "gaussian") # Default value

This creates an rctglm object which prints as

ate
#> 
#> Object of class rctglm 
#> 
#> Call:  rctglm(formula = Y ~ A * W, exposure_indicator = A, exposure_prob = 1/2, 
#>     data = dat_treat, family = "gaussian")
#> 
#> Counterfactual control mean (psi_0=E[Y|X, A=0]) estimate: 2.776
#> Counterfactual control mean (psi_1=E[Y|X, A=1]) estimate: 4.867
#> Estimand function r: psi1 - psi0
#> Estimand (r(psi_1, psi_0)) estimate (SE): 2.091 (0.09209)

The structure of such an rctglm object is broken down in the Value section of the documentation in rctglm().

Methods available are estimand (or the shorthand est) which prints a data.frame with and estimate of the estimand and its standard error. A method for coef is also available to extract coefficients from the underlying glm fit.

est(ate)
#>   Estimate Std. Error
#> 1 2.091095 0.09209306

See more info in the documentation page rctglm_methods().

Using prognostic covariate adjustment

The rctglm_with_prognosticscore() function uses the fit_best_learner() function to fit a prognostic model to historical data and then uses the prognostic model to predict \[\begin{align} \mathbb{E}[Y|X,A=0] \end{align}\]

for all observations in the current data set. These prognostic scores are then used as a covariate in the GLM when running rctglm().

Allowing the use of complex non-linear models to create such a prognostic score allows utilising information from potentially many variables, “catching” non-linear relationships and then using all this information in the GLM model using a single covariate adjustment.

We simulate some historical data to showcase the use of this function as well:

dat_notreat <- glm_data(
  Y ~ b0+b1*sin(W)^2,
  W = runif(n, min = -2, max = 2),
  family = gaussian # Default value
)

The call to rctglm_with_prognosticscore() is the same as to rctglm() but with an added specification of

• (Historical) data to fit the prognostic model using fit_best_learner()
• A formula used when fitting the prognostic model
  • Default uses all covariates in the data.
• (Optionally) number folds in cross validation and a list of learners for fitting the best learner

Thus, a simple call which estimates the average treatment effect, adjusting for a prognostic score, is seen below:

ate_prog <- rctglm_with_prognosticscore(
  formula = Y ~ A * W,
  exposure_indicator = A,
  exposure_prob = 1/2,
  data = dat_treat,
  family = gaussian(link = "identity"), # Default value
  data_hist = dat_notreat)

Quick results of the fit can be seen by printing the object:

ate_prog
#> 
#> Object of class rctglm_prog 
#> 
#> Call:  rctglm_with_prognosticscore(formula = Y ~ A * W, exposure_indicator = A, 
#>     exposure_prob = 1/2, data = dat_treat, family = gaussian(link = "identity"), 
#>     data_hist = dat_notreat)
#> 
#> Counterfactual control mean (psi_0=E[Y|X, A=0]) estimate: 2.827
#> Counterfactual control mean (psi_1=E[Y|X, A=1]) estimate: 4.821
#> Estimand function r: psi1 - psi0
#> Estimand (r(psi_1, psi_0)) estimate (SE): 1.994 (0.06405)

It’s evident that in this case where there is a non-linear relationship between the covariate we observe and the response, adjusting for the prognostic score reduces the standard error of our estimand approximation by quite a bit.

For marginal effects

The method proposed in Powering RCTs for marginal effects with GLMs using prognostic score adjustment by Højbjerre-Frandsen et. al (2025), which can be used to estimate the power when estimating any marginal effect, is implemented in the function power_marginaleffect().

According to the conservative approach in the article, if wanting to conduct power analyses to figure out how many participants is needed for an upcoming trial, where you are planning to use prognostic covariate adjustment, predictions should be obtained from a discrete super learner identical to the one planned to use for generating prognostic scores when adjusting in the analysis when estimating the marginal effect.

Here we showcase a simple use of a glm(), but fx. fit_best_learner() could be used to fit a discrete super learner as the prediction model. Could also add steps to get out-of-sample (OOS) predictions (see examples).

pred_mod <- glm(Y ~ W + A, data = dat_treat)
preds <- predict(pred_mod, dat_treat)

power_marginaleffect(
  response = dat_treat$Y,
  predictions = preds,
  target_effect = 0.4,
  exposure_prob = 1/2
)
#> [1] 0.9056303

Specific to linear models

vanc <- variance_ancova(Y ~ A + W, data = dat_treat)
vanc
#> [1] 2.120278

power_gs(variance = vanc, n = 100, ate = 0.8)
#> [1] 0.7765394
power_nc(variance = vanc, n = 100, df = 97, ate = 0.8)
#> [1] 0.7763057
samplesize_gs(variance = vanc, ate = 0.8, power = 0.9)
#> [1] 141.1623