# Analysis Clinic Case 004: significant interaction with nonsignificant main effects
#
# Purpose: demonstrate that lower-order coefficients in an interaction model are
# conditional effects at the other predictor's zero value, not unconditional
# effects. Re-expressing the moderator changes that reference point without
# changing the interaction coefficient, fitted values, or substantive model.

set.seed(20261003)

n <- 240L
x <- as.numeric(scale(rnorm(n), center = TRUE, scale = FALSE))
z <- as.numeric(scale(rnorm(n), center = TRUE, scale = FALSE))
error <- rnorm(n, mean = 0, sd = 1)
y <- 10 + 1.25 * x * z + error
dat <- data.frame(y, x, z)

fit <- lm(y ~ x * z, data = dat)
coefficient_table <- coef(summary(fit))

# In y = b0 + b1*x + b2*z + b3*x*z, the conditional slope of x at z = m is
# b1 + b3*m. Its variance uses both coefficient variances and their covariance.
simple_slope <- function(fit, moderator_value) {
  b <- coef(fit)
  v <- vcov(fit)
  estimate <- unname(b["x"] + moderator_value * b["x:z"])
  variance <- unname(
    v["x", "x"] +
      moderator_value^2 * v["x:z", "x:z"] +
      2 * moderator_value * v["x", "x:z"]
  )
  standard_error <- sqrt(variance)
  t_value <- estimate / standard_error
  p_value <- 2 * pt(abs(t_value), df = df.residual(fit), lower.tail = FALSE)
  c(estimate = estimate, standard_error = standard_error, t_value = t_value, p_value = p_value)
}

z_sd <- sd(dat$z)
probe_values <- c(low = -z_sd, mean = 0, high = z_sd)
simple_slopes <- t(vapply(probe_values, function(value) simple_slope(fit, value), numeric(4)))

# Shift z so zero now means one SD above its original mean. This is the same
# fitted model expressed around a different reference point.
dat$z_high <- dat$z - z_sd
fit_recentered <- lm(y ~ x * z_high, data = dat)
recentered_table <- coef(summary(fit_recentered))
maximum_prediction_difference <- max(abs(fitted(fit) - fitted(fit_recentered)))

cat("R version:", R.version.string, "\n")
cat("observations:", nrow(dat), "\n")
cat("mean of x:", round(mean(dat$x), 12), "\n")
cat("mean of z:", round(mean(dat$z), 12), "\n")
cat("SD of z:", round(z_sd, 6), "\n\n")

cat("ORIGINAL INTERACTION MODEL: y ~ x * z\n")
for (term in rownames(coefficient_table)) {
  cat(
    term,
    "estimate =", format(coefficient_table[term, "Estimate"], digits = 7),
    "SE =", format(coefficient_table[term, "Std. Error"], digits = 7),
    "p =", format.pval(coefficient_table[term, "Pr(>|t|)"], digits = 4),
    "\n"
  )
}
cat("\n")

cat("CONDITIONAL SLOPE OF x\n")
for (label in rownames(simple_slopes)) {
  cat(
    "z at", label,
    "estimate =", format(simple_slopes[label, "estimate"], digits = 7),
    "SE =", format(simple_slopes[label, "standard_error"], digits = 7),
    "p =", format.pval(simple_slopes[label, "p_value"], digits = 4),
    "\n"
  )
}
cat("\n")

cat("RECENTERED MODEL: zero of z_high equals original z at +1 SD\n")
cat(
  "x estimate =", format(recentered_table["x", "Estimate"], digits = 7),
  "SE =", format(recentered_table["x", "Std. Error"], digits = 7),
  "p =", format.pval(recentered_table["x", "Pr(>|t|)"], digits = 4),
  "\n"
)
cat(
  "interaction estimate =", format(recentered_table["x:z_high", "Estimate"], digits = 7),
  "p =", format.pval(recentered_table["x:z_high", "Pr(>|t|)"], digits = 4),
  "\n"
)
cat("maximum fitted-value difference:", format(maximum_prediction_difference, scientific = TRUE), "\n\n")

cat("INTERPRETATION\n")
cat("The interaction is significant while both lower-order coefficients at the centered zero values are not.\n")
cat("The slope of x is negative at a low value of z, near zero at mean z, and positive at a high value of z.\n")
cat("Recentering z changes the reference-point coefficient for x but not the interaction or fitted values.\n")
cat("The interaction should be interpreted with conditional estimates and uncertainty, not by requiring significant lower-order terms.\n")

stopifnot(
  abs(mean(dat$x)) < 1e-12,
  abs(mean(dat$z)) < 1e-12,
  coefficient_table["x", "Pr(>|t|)"] > 0.10,
  coefficient_table["z", "Pr(>|t|)"] > 0.10,
  coefficient_table["x:z", "Pr(>|t|)"] < 1e-10,
  simple_slopes["low", "estimate"] < 0,
  simple_slopes["low", "p_value"] < 0.001,
  simple_slopes["mean", "p_value"] > 0.10,
  simple_slopes["high", "estimate"] > 0,
  simple_slopes["high", "p_value"] < 0.001,
  abs(recentered_table["x", "Estimate"] - simple_slopes["high", "estimate"]) < 1e-10,
  abs(recentered_table["x:z_high", "Estimate"] - coefficient_table["x:z", "Estimate"]) < 1e-10,
  maximum_prediction_difference < 1e-10
)