Unified Regression Calibration Wrapper (External Reliability Study)
Source:R/RC_EX.R
RC_ExReliab.RdA single formula interface for regression calibration in external reliability studies. The user simply specifies `link = "linear"`, `"logistic"`, or `"log"`, and the wrapper selects the appropriate model: * `"linear"` → Gaussian (identity link) * `"logistic"` → Binomial (logit link) * `"log"` → Poisson (log link)
Usage
RC_ExReliab(
formula,
main_data,
rep_data,
link = c("linear", "logistic", "log"),
return_details = FALSE
)Arguments
- formula
A formula or character string such as `Y ~ sbp(sbp2, sbp3) + chol(chol2, chol3) + age + weight`. Terms of the form `var(rep1, rep2, ...)` are treated as error-prone exposures with replicates in `rep_data`; other terms are treated as covariates W.
- main_data
Data frame holding the outcome, error-prone exposures, and covariates.
- rep_data
Data frame holding replicate columns referenced in `formula`.
- link
Character; one of `"linear"`, `"logistic"`, or `"log"`.
- return_details
Logical; if `TRUE`, return parsed, prepared, and RC internals.
Value
A list with: * `uncorrected`: naive regression estimates * `corrected` : sandwich-corrected regression calibration estimates * optional `details` if `return_details = TRUE`
Examples
library(mgcv)
#> Loading required package: nlme
#> This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
set.seed(123)
add_err <- function(v, sd = sqrt(0.4)) v + rnorm(length(v), 0, sd)
## --- Example 1: External 1Z 0W ---
x <- rnorm(3000)
z.main <- x[1:1500] + rnorm(1500, 0, sqrt(0.4))
z_rep <- rbind(
cbind(add_err(x[1501:2000]), add_err(x[1501:2000]), NA, NA),
cbind(add_err(x[2001:2400]), add_err(x[2001:2400]), add_err(x[2001:2400]), NA),
cbind(add_err(x[2401:3000]), add_err(x[2401:3000]),
add_err(x[2401:3000]), add_err(x[2401:3000]))
)
colnames(z_rep) <- paste0("z_", 1:4)
Y <- rbinom(1500, 1, plogis(-2.3 + log(1.5) * x[1:1500]))
main_data <- data.frame(Y = Y, z = z.main)
rep_data <- data.frame(z_rep, check.names = FALSE)
res1 <- RC_ExReliab(Y ~ z(z_1, z_2, z_3, z_4), main_data, rep_data, link = "logistic")
res1$corrected
#> Estimate Std. Error z value Pr(>|z|) OR CI.low
#> (Intercept) -2.439767 0.09612526 -25.381117 2.000000000 0.0871812 0.07221033
#> z 0.283015 0.10867516 2.604229 0.009208109 1.3271251 1.07252077
#> CI.high
#> (Intercept) 0.1052559
#> z 1.6421697
## --- Example 2: External 1Z 1W ---
x <- rnorm(3000)
W_main <- rnorm(1500)
W_rep <- rnorm(1500)
z.main <- x[1:1500] + rnorm(1500, 0, sqrt(0.4))
z_rep <- rbind(
cbind(add_err(x[1501:2000]), add_err(x[1501:2000]), NA, NA),
cbind(add_err(x[2001:2400]), add_err(x[2001:2400]), add_err(x[2001:2400]), NA),
cbind(add_err(x[2401:3000]), add_err(x[2401:3000]),
add_err(x[2401:3000]), add_err(x[2401:3000]))
)
colnames(z_rep) <- paste0("z_", 1:4)
Y <- rbinom(1500, 1, plogis(-2.3 + log(1.5) * x[1:1500] + 0.5 * W_main))
main_data <- data.frame(Y = Y, z = z.main, W = W_main)
rep_data <- data.frame(z_rep, W = W_rep, check.names = FALSE)
res2 <- RC_ExReliab(Y ~ z(z_1, z_2, z_3, z_4) + W, main_data, rep_data, link = "logistic")
res2$corrected
#> Estimate Std. Error z value Pr(>|z|) OR CI.low
#> (Intercept) -2.1748410 0.08906444 -24.418735 2.000000e+00 0.1136262 0.09542571
#> z 0.2759814 0.09232740 2.989160 2.797452e-03 1.3178234 1.09968052
#> W 0.4279343 0.08398317 5.095477 3.478635e-07 1.5340854 1.30125259
#> CI.high
#> (Intercept) 0.1352981
#> z 1.5792391
#> W 1.8085788
## --- Example 3: External 2Z 0W ---
x <- mgcv::rmvn(3000, c(0, 0), matrix(c(1, 0.3, 0.3, 1), 2))
z.main <- x[1:1500, ] + matrix(rnorm(1500 * 2, 0, sqrt(0.4)), 1500, 2)
colnames(z.main) <- c("z1", "z2")
z1_rep <- rbind(
cbind(add_err(x[1501:2000, 1]), add_err(x[1501:2000, 1]), NA, NA),
cbind(add_err(x[2001:2400, 1]), add_err(x[2001:2400, 1]), add_err(x[2001:2400, 1]), NA),
cbind(add_err(x[2401:3000, 1]), add_err(x[2401:3000, 1]),
add_err(x[2401:3000, 1]), add_err(x[2401:3000, 1]))
)
colnames(z1_rep) <- paste0("z1_", 1:4)
z2_rep <- rbind(
cbind(add_err(x[1501:2000, 2]), add_err(x[1501:2000, 2]), NA, NA),
cbind(add_err(x[2001:2400, 2]), add_err(x[2001:2400, 2]), add_err(x[2001:2400, 2]), NA),
cbind(add_err(x[2401:3000, 2]), add_err(x[2401:3000, 2]),
add_err(x[2401:3000, 2]), add_err(x[2401:3000, 2]))
)
colnames(z2_rep) <- paste0("z2_", 1:4)
Y <- rbinom(1500, 1, plogis(-2.3 + log(1.5) * rowSums(x[1:1500, ])))
main_data <- data.frame(Y = Y, z1 = z.main[, 1], z2 = z.main[, 2])
rep_data <- data.frame(z1_rep, z2_rep, check.names = FALSE)
res3 <- RC_ExReliab(
Y ~ z1(z1_1, z1_2, z1_3, z1_4) + z2(z2_1, z2_2, z2_3, z2_4),
main_data, rep_data, link = "logistic"
)
res3$corrected
#> Estimate Std. Error z value Pr(>|z|) OR CI.low
#> (Intercept) -2.2593965 0.09444947 -23.921749 2.000000e+00 0.1044135 0.086768
#> z1 0.2657353 0.11124472 2.388745 1.690601e-02 1.3043898 1.048851
#> z2 0.5326281 0.10950870 4.863797 1.151548e-06 1.7034032 1.374364
#> CI.high
#> (Intercept) 0.1256474
#> z1 1.6221866
#> z2 2.1112183