Skip to contents

Calculate Simple (Stratum-Specific) Effects from a Joint Variable Model

Source: R/inteffects.R
inteffects.Rd

Replicates the output of an interaction-term model (e.g., simple effects from `publish()` or `emmeans`) by calculating all stratum-specific comparisons from a joint-variable model. Allows output on log or ratio scale with consistent column structure, including scale-appropriate standard errors.

Usage

inteffects(
  joint_model,
  joint_var_name,
  factor1_name,
  factor2_name,
  factor1_levels,
  factor2_levels,
  level_separator = ".",
  scale = c("ratio", "log"),
  digits = 2,
  conf.level = 0.95
)

Arguments

joint_model

A fitted model object (e.g., `svycoxph`, `svyglm`) that used a combined joint variable (e.g., `~ A_B`).

joint_var_name

Character string. The name of the joint variable as used in the model formula (e.g., `"Race1_ObeseStatus"`).

factor1_name

Character string. The name of the *first* logical factor used to create the joint variable (e.g., `"Race1"`). Used for naming output rows.

factor2_name

Character string. The name of the *second* logical factor used to create the joint variable (e.g., `"ObeseStatus"`). Used for naming output rows.

factor1_levels

Character vector. All levels of the *first* logical factor, with the reference level *first*.

factor2_levels

Character vector. All levels of the *second* logical factor, with the reference level *first*.

level_separator

Character string. The separator used to create the joint level names (e.g., `"_"` or `"."`). Default is `"."` (the R default for `interaction()`).

scale

Character string. The scale for the output estimates, CIs, and SE. Options are `"ratio"` (default, e.g., OR/HR) or `"log"` (log-OR/log-HR).

digits

Integer. Number of decimal places for rounding estimates, SE, and CIs when `scale = "ratio"`. Default is 2. (Set to NULL for no rounding).

conf.level

Confidence level for the interval (default 0.95).

Value

A `tibble` (data frame) with columns:

  • `Comparison`: A description of the simple effect, styled similarly to the `Publish` package (e.g., "ObeseStatus(Obese vs Not Obese): Race1(White)").

  • `Estimate`: The point estimate on the specified `scale`.

  • `SE`: The standard error on the specified `scale`. Calculated using the delta method for `scale = "ratio"`.

  • `CI.low`: The lower bound of the confidence interval (on the specified `scale`).

  • `CI.upp`: The upper bound of the confidence interval (on the specified `scale`).

  • `p_value`: The p-value for the simple effect (tests H0: logEstimate=0 or Estimate=1).

Examples

# \donttest{
# --- Load required libraries for the example ---
if (requireNamespace("survey", quietly = TRUE) &&
    requireNamespace("NHANES", quietly = TRUE) &&
    requireNamespace("dplyr", quietly = TRUE) &&
    requireNamespace("tidyr", quietly = TRUE) &&
    requireNamespace("msm", quietly = TRUE)) { # msm needed for SE calculation

  library(survey)
  library(NHANES)
  library(dplyr)
  library(tidyr)
  library(msm) # Load msm

  # --- 1. Data Preparation (NHANES Example) ---
  data(NHANESraw)

  vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130",
                   "SDMVPSU", "SDMVSTRA", "WTMEC2YR")

  nhanes_adults_processed <- NHANESraw %>%
    filter(Age >= 20) %>%
    mutate(
      ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"),
                           levels = c("Not Obese", "Obese")),
      Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"),
                                levels = c("No", "Yes")),
      Race1 = relevel(as.factor(Race1), ref = "White")
    ) %>%
    select(all_of(vars_needed)) %>%
    drop_na()

  adult_design_binary <- svydesign(
    id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR,
    nest = TRUE, data = nhanes_adults_processed
  )

  # --- 2. Create Joint Variable and Fit Joint Model ---
  adult_design_binary <- update(adult_design_binary,
    Race1_ObeseStatus = interaction(Race1, ObeseStatus, sep = "_", drop = TRUE)
  )
  adult_design_binary <- update(adult_design_binary,
    Race1_ObeseStatus = relevel(Race1_ObeseStatus, ref = "White_Not Obese")
  )

  joint_model_logit <- svyglm(
    Hypertension_130 ~ Race1_ObeseStatus + Age,
    design = adult_design_binary, family = quasibinomial()
  )

  # --- 3. Run the function ---
  f1_levels <- levels(adult_design_binary$variables$Race1)
  f2_levels <- levels(adult_design_binary$variables$ObeseStatus)

  # --- Example 1: Output on Ratio scale ---
  simple_effects_ratio <- inteffects(
    joint_model = joint_model_logit,
    joint_var_name = "Race1_ObeseStatus",
    factor1_name = "Race1",
    factor2_name = "ObeseStatus",
    factor1_levels = f1_levels,
    factor2_levels = f2_levels,
    level_separator = "_",
    scale = "ratio"
  )
  print("--- Output on Ratio Scale ---")
  print(simple_effects_ratio, n = 50)

  # --- Example 2: Output on Log scale ---
  simple_effects_log <- inteffects(
    joint_model = joint_model_logit,
    joint_var_name = "Race1_ObeseStatus",
    factor1_name = "Race1",
    factor2_name = "ObeseStatus",
    factor1_levels = f1_levels,
    factor2_levels = f2_levels,
    level_separator = "_",
    scale = "log"
  )
  print("--- Output on Log Scale ---")
  print(simple_effects_log, n = 50)
}
#> [1] "--- Output on Ratio Scale ---"
#> # A tibble: 13 × 6
#>    Comparison                             Estimate    SE CI.low CI.upp `p-value`
#>    <chr>                                     <dbl> <dbl>  <dbl>  <dbl>     <dbl>
#>  1 ObeseStatus(Obese vs Not Obese): Race…     1.4   0.1    1.21   1.61  4.88e- 6
#>  2 ObeseStatus(Obese vs Not Obese): Race…     1.2   0.14   0.95   1.52  1.30e- 1
#>  3 ObeseStatus(Obese vs Not Obese): Race…     1.66  0.24   1.25   2.2   4.01e- 4
#>  4 ObeseStatus(Obese vs Not Obese): Race…     1.65  0.2    1.3    2.1   4.15e- 5
#>  5 ObeseStatus(Obese vs Not Obese): Race…     1.09  0.27   0.68   1.77  7.15e- 1
#>  6 Race1(Black vs White): ObeseStatus(No…     2.11  0.21   1.73   2.58  1.42e-13
#>  7 Race1(Hispanic vs White): ObeseStatus…     0.99  0.12   0.77   1.27  9.37e- 1
#>  8 Race1(Mexican vs White): ObeseStatus(…     1.17  0.12   0.95   1.43  1.37e- 1
#>  9 Race1(Other vs White): ObeseStatus(No…     1.05  0.1    0.87   1.28  6.00e- 1
#> 10 Race1(Black vs White): ObeseStatus(Ob…     1.81  0.23   1.41   2.33  3.75e- 6
#> 11 Race1(Hispanic vs White): ObeseStatus…     1.18  0.15   0.92   1.51  1.96e- 1
#> 12 Race1(Mexican vs White): ObeseStatus(…     1.38  0.2    1.03   1.84  2.96e- 2
#> 13 Race1(Other vs White): ObeseStatus(Ob…     0.82  0.19   0.52   1.31  4.14e- 1
#> [1] "--- Output on Log Scale ---"
#> # A tibble: 13 × 6
#>    Comparison                           Estimate     SE  CI.low CI.upp `p-value`
#>    <chr>                                   <dbl>  <dbl>   <dbl>  <dbl>     <dbl>
#>  1 ObeseStatus(Obese vs Not Obese): Ra…  0.335   0.0732  0.191   0.478  4.88e- 6
#>  2 ObeseStatus(Obese vs Not Obese): Ra…  0.182   0.120  -0.0538  0.418  1.30e- 1
#>  3 ObeseStatus(Obese vs Not Obese): Ra…  0.508   0.144   0.227   0.789  4.01e- 4
#>  4 ObeseStatus(Obese vs Not Obese): Ra…  0.500   0.122   0.261   0.740  4.15e- 5
#>  5 ObeseStatus(Obese vs Not Obese): Ra…  0.0896  0.246  -0.392   0.571  7.15e- 1
#>  6 Race1(Black vs White): ObeseStatus(…  0.748   0.101   0.550   0.946  1.42e-13
#>  7 Race1(Hispanic vs White): ObeseStat… -0.00992 0.125  -0.255   0.235  9.37e- 1
#>  8 Race1(Mexican vs White): ObeseStatu…  0.154   0.104  -0.0493  0.358  1.37e- 1
#>  9 Race1(Other vs White): ObeseStatus(…  0.0516  0.0983 -0.141   0.244  6.00e- 1
#> 10 Race1(Black vs White): ObeseStatus(…  0.595   0.129   0.343   0.848  3.75e- 6
#> 11 Race1(Hispanic vs White): ObeseStat…  0.164   0.127  -0.0845  0.412  1.96e- 1
#> 12 Race1(Mexican vs White): ObeseStatu…  0.320   0.147   0.0317  0.608  2.96e- 2
#> 13 Race1(Other vs White): ObeseStatus(… -0.193   0.237  -0.657   0.270  4.14e- 1
# }