Unified Regression Calibration Wrapper (Internal Reliability Study)
Source:R/RC_IN.R
RC_InReliab.RdA single formula interface for regression calibration in internal 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_InReliab(
formula,
main_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 `main_data`; other terms are treated as covariates W.
- main_data
Data frame holding the outcome, replicate error-prone exposures, and any covariates.
- 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
set.seed(123)
add_err <- function(v, sd = sqrt(0.4)) v + rnorm(length(v), 0, sd)
## --- Example 1: Internal 1Z 0W ---
x <- rnorm(3000)
z <- rbind(
cbind(add_err(x[1:1500]), NA, NA, NA),
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) <- paste0("z_", 1:4)
Y <- rbinom(3000, 1, plogis(-2.65 + log(1.5) * x))
main_data <- data.frame(Y, z)
res1 <- RC_InReliab(Y ~ myz(z_1, z_2, z_3, z_4),
main_data = main_data,
link = "logistic")
res1$corrected
#> Estimate Std. Error z value Pr(>|z|) OR CI.low
#> (Intercept) -2.7171752 0.07701872 -35.279412 2.000000e+00 0.0660611 0.05680495
#> myz 0.3411573 0.08265864 4.127304 3.670418e-05 1.4065745 1.19619582
#> CI.high
#> (Intercept) 0.07682551
#> myz 1.65395313
## --- Example 2: Internal 1Z 1W ---
x <- rnorm(3000)
W1 <- rnorm(3000)
z <- rbind(
cbind(add_err(x[1:1500]), NA, NA, NA),
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) <- paste0("z_", 1:4)
Y <- rbinom(3000, 1, plogis(-2.65 + log(1.5) * x + 0.5 * W1))
main_data <- data.frame(Y, z, W1)
res2 <- RC_InReliab(Y ~ myz(z_1, z_2, z_3, z_4) + W1,
main_data = main_data,
link = "logistic")
res2$corrected
#> Estimate Std. Error z value Pr(>|z|) OR CI.low
#> (Intercept) -2.5679464 0.07450904 -34.464900 2.000000e+00 0.07669288 0.06627224
#> myz 0.4930941 0.07141417 6.904709 5.030643e-12 1.63737454 1.42350522
#> W1 0.4994671 0.06588726 7.580633 3.438739e-14 1.64784292 1.44820967
#> CI.high
#> (Intercept) 0.08875205
#> myz 1.88337587
#> W1 1.87499527
## --- Example 3: Internal 2Z 0W ---
x <- mgcv::rmvn(3000, c(0,0), matrix(c(1,0.3,0.3,1), 2))
z1 <- rbind(
cbind(add_err(x[1:1500, 1]), NA, NA, NA),
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) <- paste0("z1_", 1:4)
z2 <- rbind(
cbind(add_err(x[1:1500, 2]), NA, NA, NA),
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) <- paste0("z2_", 1:4)
Y <- rbinom(3000, 1, plogis(-2.65 + log(1.5) * rowSums(x)))
main_data <- data.frame(Y, z1, z2)
res3 <- RC_InReliab(
Y ~ myz1(z1_1, z1_2, z1_3, z1_4) + myz2(z2_1, z2_2, z2_3, z2_4),
main_data = main_data,
link = "logistic")
res3$corrected
#> Estimate Std. Error z value Pr(>|z|) OR CI.low
#> (Intercept) -2.6495693 0.07952853 -33.315959 2.000000e+00 0.07068165 0.06047984
#> myz1 0.4285456 0.08305847 5.159566 2.475230e-07 1.53502342 1.30441027
#> myz2 0.4835185 0.09087986 5.320414 1.035313e-07 1.62177059 1.35715958
#> CI.high
#> (Intercept) 0.08260431
#> myz1 1.80640779
#> myz2 1.93797389