Every model in this section is fit on clinical_sample.csv data. This clinical_sample.csv data is purely additive, so every interaction tested on it is null by design.

1 Load data

d <- read.csv("../source_material/data/clinical_sample.csv")
str(d)
## 'data.frame':    280 obs. of  6 variables:
##  $ id         : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ self_esteem: num  19 28.2 37.6 34.3 28.7 23 30.5 21.9 23.2 29.4 ...
##  $ neg_affect : num  33.1 19.5 24.8 20.1 30.8 23.2 19.7 24.1 18.9 21.3 ...
##  $ treatment  : int  0 0 1 1 0 0 1 0 1 1 ...
##  $ gender     : int  1 1 0 0 0 0 0 1 1 0 ...
##  $ depression : num  16.8 14.2 0.2 6.1 7.1 13 10.3 9.8 6.4 1.9 ...
summary(d)
##        id         self_esteem     neg_affect     treatment         gender     
##  Min.   :  1.0   Min.   :10.0   Min.   :10.0   Min.   :0.000   Min.   :0.000  
##  1st Qu.: 70.8   1st Qu.:25.0   1st Qu.:17.2   1st Qu.:0.000   1st Qu.:0.000  
##  Median :140.5   Median :28.6   Median :21.6   Median :1.000   Median :0.000  
##  Mean   :140.5   Mean   :28.7   Mean   :22.2   Mean   :0.507   Mean   :0.468  
##  3rd Qu.:210.2   3rd Qu.:32.1   3rd Qu.:26.9   3rd Qu.:1.000   3rd Qu.:1.000  
##  Max.   :280.0   Max.   :40.0   Max.   :41.4   Max.   :1.000   Max.   :1.000  
##    depression   
##  Min.   : 0.00  
##  1st Qu.: 5.38  
##  Median : 9.00  
##  Mean   : 9.11  
##  3rd Qu.:12.82  
##  Max.   :23.80

2 From covariance to correlation

2.1 Covariance matrix

vars <- c("self_esteem", "neg_affect", "depression")
cov_mat <- cov(d[, vars])
round(cov_mat, 2)
##             self_esteem neg_affect depression
## self_esteem       26.83     -11.90      -8.37
## neg_affect       -11.90      48.99      16.81
## depression        -8.37      16.81      26.52

2.1.1 Manual covariance function (cf. 01_Correlation_and_Regression.Rmd)

We can write our own covariance function, checked against cov().

mycov <- function(x, y) sum((x - mean(x)) * (y - mean(y))) / (length(x) - 1)
mycov(d$self_esteem, d$depression)
## [1] -8.374
cov(d$self_esteem, d$depression)   # should match exactly
## [1] -8.374

2.2 Correlation matrix

cor_mat <- cor(d[, vars])
round(cor_mat, 3)
##             self_esteem neg_affect depression
## self_esteem       1.000     -0.328     -0.314
## neg_affect       -0.328      1.000      0.466
## depression       -0.314      0.466      1.000

2.3 Scatterplot matrix (sanity check on linearity/outliers)

pairs(d[, vars], pch = 16, col = adjustcolor(ntnublue, 0.4))

3 Simple linear regression

3.1 Continuous predictor: self-esteem -> depression

m1 <- lm(depression ~ self_esteem, data = d)
summary(m1)
## 
## Call:
## lm(formula = depression ~ self_esteem, data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -10.78  -3.52  -0.43   3.58  13.38 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.0614     1.6500   10.95  < 2e-16 ***
## self_esteem  -0.3121     0.0566   -5.51  8.1e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.9 on 278 degrees of freedom
## Multiple R-squared:  0.0985, Adjusted R-squared:  0.0953 
## F-statistic: 30.4 on 1 and 278 DF,  p-value: 8.09e-08
confint(m1)   # 95% CIs quoted on every results table in the slides
##               2.5 %  97.5 %
## (Intercept) 14.8134 21.3094
## self_esteem -0.4236 -0.2006
# standardized beta (should equal r in simple regression)
lm.beta(m1)
## 
## Call:
## lm(formula = depression ~ self_esteem, data = d)
## 
## Standardized Coefficients::
## (Intercept) self_esteem 
##          NA     -0.3139
cor(d$self_esteem, d$depression)
## [1] -0.3139
ggplot(d, aes(self_esteem, depression)) +
  geom_point(alpha = 0.55, color = ntnublue) +
  geom_smooth(method = "lm", se = FALSE, color = red, linewidth = 1.2) +
  labs(x = "Self-esteem (RSES, 10-40)", y = "Depression severity (0-27)")

3.2 Dummy predictor: treatment -> depression

m2 <- lm(depression ~ treatment, data = d)
summary(m2)
## 
## Call:
## lm(formula = depression ~ treatment, data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -10.68  -3.44   0.01   3.55  13.12 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   10.675      0.419   25.48  < 2e-16 ***
## treatment     -3.085      0.588   -5.24  3.1e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.92 on 278 degrees of freedom
## Multiple R-squared:  0.09,   Adjusted R-squared:  0.0867 
## F-statistic: 27.5 on 1 and 278 DF,  p-value: 3.12e-07
confint(m2)
##              2.5 % 97.5 %
## (Intercept)  9.851 11.500
## treatment   -4.243 -1.927
d %>% group_by(treatment) %>%
  summarise(n = n(), M = mean(depression), SD = sd(depression))
## # A tibble: 2 × 4
##   treatment     n     M    SD
##       <int> <int> <dbl> <dbl>
## 1         0   138 10.7   4.98
## 2         1   142  7.59  4.87

4 Multiple linear regression

4.1 Two continuous predictors

m3 <- lm(depression ~ self_esteem + neg_affect, data = d)
summary(m3)
## 
## Call:
## lm(formula = depression ~ self_esteem + neg_affect, data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.010 -3.090 -0.192  3.065 14.924 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.6087     2.0721    3.67  0.00029 ***
## self_esteem  -0.1792     0.0549   -3.26  0.00123 ** 
## neg_affect    0.2996     0.0406    7.37  1.9e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.49 on 277 degrees of freedom
## Multiple R-squared:  0.246,  Adjusted R-squared:  0.241 
## F-statistic: 45.3 on 2 and 277 DF,  p-value: <2e-16
confint(m3)
##               2.5 %   97.5 %
## (Intercept)  3.5296 11.68784
## self_esteem -0.2873 -0.07115
## neg_affect   0.2196  0.37955
lm.beta(m3)
## 
## Call:
## lm(formula = depression ~ self_esteem + neg_affect, data = d)
## 
## Standardized Coefficients::
## (Intercept) self_esteem  neg_affect 
##          NA     -0.1803      0.4071

4.2 Two categorical predictors (= two-way factorial ANOVA)

m4 <- lm(depression ~ treatment + gender, data = d)
summary(m4)
## 
## Call:
## lm(formula = depression ~ treatment + gender, data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.810 -3.473 -0.078  3.364 13.717 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.810      0.490   20.04  < 2e-16 ***
## treatment     -3.128      0.578   -5.41  1.4e-07 ***
## gender         1.895      0.580    3.27   0.0012 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.84 on 277 degrees of freedom
## Multiple R-squared:  0.124,  Adjusted R-squared:  0.117 
## F-statistic: 19.6 on 2 and 277 DF,  p-value: 1.12e-08
confint(m4)
##               2.5 % 97.5 %
## (Intercept)  8.8467 10.774
## treatment   -4.2664 -1.989
## gender       0.7538  3.036

4.3 Mixed predictors (continuous + categorical = ANCOVA)

m5 <- lm(depression ~ self_esteem + treatment, data = d)
summary(m5)
## 
## Call:
## lm(formula = depression ~ self_esteem + treatment, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.677  -2.937  -0.144   3.203  12.504 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.5410     1.5927   12.27  < 2e-16 ***
## self_esteem  -0.3096     0.0539   -5.75  2.4e-08 ***
## treatment    -3.0583     0.5571   -5.49  9.1e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.66 on 277 degrees of freedom
## Multiple R-squared:  0.187,  Adjusted R-squared:  0.181 
## F-statistic: 31.9 on 2 and 277 DF,  p-value: 3.54e-13
confint(m5)
##               2.5 %  97.5 %
## (Intercept) 16.4056 22.6763
## self_esteem -0.4157 -0.2036
## treatment   -4.1550 -1.9615

4.4 Full model (used for diagnostics later)

m6 <- lm(depression ~ self_esteem + neg_affect + treatment + gender, data = d)
summary(m6)
## 
## Call:
## lm(formula = depression ~ self_esteem + neg_affect + treatment + 
##     gender, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.188  -3.053   0.087   2.564  12.437 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.6906     1.9283    3.99  8.5e-05 ***
## self_esteem  -0.1709     0.0500   -3.42  0.00073 ***
## neg_affect    0.3109     0.0371    8.38  2.8e-15 ***
## treatment    -3.1159     0.4888   -6.37  7.7e-10 ***
## gender        2.1552     0.4909    4.39  1.6e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.09 on 275 degrees of freedom
## Multiple R-squared:  0.379,  Adjusted R-squared:  0.37 
## F-statistic: 41.9 on 4 and 275 DF,  p-value: <2e-16
confint(m6)
##               2.5 %   97.5 %
## (Intercept)  3.8945 11.48664
## self_esteem -0.2694 -0.07239
## neg_affect   0.2379  0.38394
## treatment   -4.0782 -2.15356
## gender       1.1889  3.12154

4.5 Orthogonal predictors: a concrete demonstration

treatment is randomized, so it is (by design) almost perfectly uncorrelated with the continuous predictors. That lets us show, very concretely, that when two predictors are orthogonal, the joint-model coefficient can be recovered exactly two different ways – jointly, or by residualizing first and regressing second (Kline, 2016, uses the same cor.test(residuals(m), predictor) orthogonality check

cor(d$treatment, d$self_esteem)   # ~0: near-orthogonal by randomization
## [1] 0.008419
cor(d$treatment, d$neg_affect)    # ~0
## [1] -0.0006069
# (1) Jointly
m_joint <- lm(depression ~ self_esteem + treatment, data = d)
coef(m_joint)["treatment"]
## treatment 
##    -3.058
# (2) Sequentially: residualize depression on self_esteem FIRST,
#     then regress those residuals on treatment ALONE
m_step1 <- lm(depression ~ self_esteem, data = d)
resid_dep_on_se <- resid(m_step1)
m_step2 <- lm(resid_dep_on_se ~ treatment, data = d)
coef(m_step2)["treatment"]   # identical to (1)
## treatment 
##    -3.058
# Orthogonality check on residuals (Kline, 2016 style)
cor.test(resid(m1), d$self_esteem)   # r should be ~0 by construction of OLS
## 
##  Pearson's product-moment correlation
## 
## data:  resid(m1) and d$self_esteem
## t = -3.4e-15, df = 278, p-value = 1
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1172  0.1172
## sample estimates:
##        cor 
## -2.011e-16

5 Interaction and Moderation Effects

Every model in this section (except the “back to our real data” check at the end) is fit on interaction_demo.csv, a second, purpose-built companion dataset – same variable names/scales as clinical_sample.csv, but with three true two-way interactions built in. clinical_sample.csv itself is purely additive, so every interaction tested on it is null by construction.

d2 <- read.csv("../source_material/data/interaction_demo.csv")
str(d2)
## 'data.frame':    280 obs. of  6 variables:
##  $ id         : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ self_esteem: num  25.6 23.4 33.6 28.8 34.6 21.8 33.3 28.5 19.3 30.3 ...
##  $ neg_affect : num  16.5 33 19.7 37.1 17 23.9 22.8 12.7 18.7 22.9 ...
##  $ treatment  : int  1 1 1 0 1 0 0 0 0 0 ...
##  $ gender     : int  1 0 1 1 0 0 0 0 0 1 ...
##  $ depression : num  9.9 19.9 9.7 20.5 8.5 15.6 13.4 5.6 8.5 16.6 ...

5.1 Continuous x continuous: self-esteem x negative affect

m7 <- lm(depression ~ self_esteem * neg_affect, data = d2)
summary(m7)
## 
## Call:
## lm(formula = depression ~ self_esteem * neg_affect, data = d2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.464 -2.838 -0.062  2.834 11.014 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             -2.1604     5.1780   -0.42   0.6768    
## self_esteem              0.2485     0.1718    1.45   0.1492    
## neg_affect               0.9199     0.2312    3.98  8.8e-05 ***
## self_esteem:neg_affect  -0.0221     0.0079   -2.80   0.0054 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.13 on 276 degrees of freedom
## Multiple R-squared:  0.258,  Adjusted R-squared:  0.25 
## F-statistic: 31.9 on 3 and 276 DF,  p-value: <2e-16
confint(m7)
##                            2.5 %    97.5 %
## (Intercept)            -12.35392  8.033054
## self_esteem             -0.08970  0.586776
## neg_affect               0.46477  1.374974
## self_esteem:neg_affect  -0.03768 -0.006588

Simple slopes of self-esteem at low/mean/high negative affect (Equation: \(B_1+B_3w\)), plus the region-of-significance (Johnson-Neyman) boundary:

b <- coef(m7); V <- vcov(m7)
na_mean <- mean(d2$neg_affect); na_sd <- sd(d2$neg_affect)
simple_slope <- function(w) {
  slope <- b["self_esteem"] + b["self_esteem:neg_affect"] * w
  se <- sqrt(V["self_esteem", "self_esteem"] + w^2 * V["self_esteem:neg_affect", "self_esteem:neg_affect"] +
             2 * w * V["self_esteem", "self_esteem:neg_affect"])
  t <- slope / se
  p <- 2 * (1 - pt(abs(t), df = m7$df.residual))
  c(neg_affect = w, slope = slope, se = se, t = t, p = p)
}
rbind(low = simple_slope(na_mean - na_sd), mean = simple_slope(na_mean), high = simple_slope(na_mean + na_sd))
##      neg_affect slope.self_esteem      se t.self_esteem p.self_esteem
## low       14.97          -0.08288 0.07132        -1.162     2.462e-01
## mean      21.56          -0.22873 0.05651        -4.047     6.736e-05
## high      28.15          -0.37458 0.08196        -4.570     7.363e-06

5.1.1 The Johnson-Neyman technique

Reproduces the “The Johnson-Neyman Technique” slide: instead of evaluating the simple slope at three arbitrary points, solve directly for the moderator value(s) where its \(t\)-statistic crosses \(\pm t_{crit}\).

tcrit <- qt(0.975, m7$df.residual)
b1 <- b["self_esteem"]; b3 <- b["self_esteem:neg_affect"]
Vb1 <- V["self_esteem", "self_esteem"]; Vb3 <- V["self_esteem:neg_affect", "self_esteem:neg_affect"]
Cb13 <- V["self_esteem", "self_esteem:neg_affect"]
A <- b3^2 - tcrit^2 * Vb3; B <- 2*b1*b3 - 2*tcrit^2*Cb13; C <- b1^2 - tcrit^2*Vb1
roots <- sort((-B + c(1, -1) * sqrt(B^2 - 4*A*C)) / (2*A))
na_range <- range(d2$neg_affect)

jn_boundary <- roots[roots >= na_range[1] & roots <= na_range[2]]
data.frame(root = roots,
           in_observed_range = roots >= na_range[1] & roots <= na_range[2])
##     root in_observed_range
## 1 -12.58             FALSE
## 2  16.85              TRUE
cat(sprintf("Observed neg_affect range: %.1f-%.1f\n", na_range[1], na_range[2]))
## Observed neg_affect range: 10.0-37.1
cat(sprintf("Johnson-Neyman boundary (in range): %.2f\n", jn_boundary))
## Johnson-Neyman boundary (in range): 16.85

Only the in-range root (16.85) is a genuine region-of-significance boundary; the other is outside the data and not interpretable. This matches the pick-a-point table directly: Low (\(w=15\)) sits just below the boundary and was non-significant; Mean and High (\(w=22, 28\)) sit above it and were both significant.

w <- seq(na_range[1], na_range[2], length.out = 400)
jn_df <- data.frame(
  w = w,
  slope = b1 + b3 * w,
  se = sqrt(Vb1 + w^2 * Vb3 + 2 * w * Cb13)
) %>%
  mutate(lo = slope - tcrit * se, hi = slope + tcrit * se,
         sig = w >= jn_boundary)

ggplot(jn_df, aes(w, slope)) +
  geom_ribbon(aes(ymin = lo, ymax = hi, fill = sig), alpha = 0.18) +
  geom_line(color = ntnublue, linewidth = 1) +
  geom_hline(yintercept = 0, linewidth = 0.4) +
  geom_vline(xintercept = jn_boundary, color = red, linetype = "dotted") +
  scale_fill_manual(values = c(`FALSE` = "grey60", `TRUE` = ntnublue),
                     labels = c("Non-significant (CI includes 0)", "Significant (CI excludes 0)"),
                     name = NULL) +
  labs(x = "Negative affect (moderator)",
       y = "Simple slope of self-esteem on depression") +
  theme(legend.position = "bottom")

5.2 Continuous x dummy: self-esteem x gender

m8 <- lm(depression ~ self_esteem * gender, data = d2)
summary(m8)
## 
## Call:
## lm(formula = depression ~ self_esteem * gender, data = d2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.480  -2.925  -0.128   2.714  10.765 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         24.5674     2.2210   11.06  < 2e-16 ***
## self_esteem         -0.5078     0.0765   -6.63  1.7e-10 ***
## gender              -7.1588     3.2670   -2.19   0.0293 *  
## self_esteem:gender   0.3407     0.1112    3.06   0.0024 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.27 on 276 degrees of freedom
## Multiple R-squared:  0.207,  Adjusted R-squared:  0.198 
## F-statistic:   24 on 3 and 276 DF,  p-value: 8.27e-14
confint(m8)
##                       2.5 %  97.5 %
## (Intercept)         20.1952 28.9396
## self_esteem         -0.6585 -0.3571
## gender             -13.5901 -0.7274
## self_esteem:gender   0.1218  0.5596
coef(m8)["self_esteem"]                              # men's slope
## self_esteem 
##     -0.5078
coef(m8)["self_esteem"] + coef(m8)["self_esteem:gender"]  # women's slope
## self_esteem 
##     -0.1671

5.3 Dummy x dummy: treatment x gender

m9 <- lm(depression ~ treatment * gender, data = d2)
summary(m9)
## 
## Call:
## lm(formula = depression ~ treatment * gender, data = d2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.165 -3.159 -0.105  2.456 11.495 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        10.805      0.494   21.87  < 2e-16 ***
## treatment          -1.641      0.716   -2.29   0.0228 *  
## gender              3.666      0.686    5.35  1.9e-07 ***
## treatment:gender   -3.174      1.032   -3.07   0.0023 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.28 on 276 degrees of freedom
## Multiple R-squared:  0.204,  Adjusted R-squared:  0.195 
## F-statistic: 23.6 on 3 and 276 DF,  p-value: 1.33e-13
confint(m9)
##                   2.5 %  97.5 %
## (Intercept)       9.833 11.7778
## treatment        -3.051 -0.2304
## gender            2.317  5.0159
## treatment:gender -5.206 -1.1413
# cell means
aggregate(depression ~ treatment + gender, data = d2, FUN = mean)
##   treatment gender depression
## 1         0      0     10.805
## 2         1      0      9.165
## 3         0      1     14.472
## 4         1      1      9.657
# simple effects of treatment, by gender
coef(m9)["treatment"]                                  # men
## treatment 
##    -1.641
coef(m9)["treatment"] + coef(m9)["treatment:gender"]    # women
## treatment 
##    -4.814

5.4 Centering: raw vs. mean-centered (self-esteem x negative affect)

d2$se_c <- d2$self_esteem - mean(d2$self_esteem)
d2$na_c <- d2$neg_affect - mean(d2$neg_affect)
m7_centered <- lm(depression ~ se_c * na_c, data = d2)
summary(m7_centered)
## 
## Call:
## lm(formula = depression ~ se_c * na_c, data = d2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.464 -2.838 -0.062  2.834 11.014 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.0427     0.2565   43.05  < 2e-16 ***
## se_c         -0.2287     0.0565   -4.05  6.7e-05 ***
## na_c          0.2782     0.0393    7.09  1.1e-11 ***
## se_c:na_c    -0.0221     0.0079   -2.80   0.0054 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.13 on 276 degrees of freedom
## Multiple R-squared:  0.258,  Adjusted R-squared:  0.25 
## F-statistic: 31.9 on 3 and 276 DF,  p-value: <2e-16
# Interaction coefficient, R^2, and RSE are identical to the raw model;
# only the intercept and main effects change meaning.
identical(round(unname(coef(m7)["self_esteem:neg_affect"]), 4),
          round(unname(coef(m7_centered)["se_c:na_c"]), 4))
## [1] TRUE

5.5 Back to clinical_sample.csv: a non-significant interaction

m_real_int <- lm(depression ~ self_esteem * neg_affect, data = d)
summary(m_real_int)   # interaction term is non-significant (p = .120)
## 
## Call:
## lm(formula = depression ~ self_esteem * neg_affect, data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.817 -3.126 -0.158  3.144 15.020 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)   
## (Intercept)            -0.12488    5.37482   -0.02   0.9815   
## self_esteem             0.08298    0.17691    0.47   0.6394   
## neg_affect              0.64028    0.22231    2.88   0.0043 **
## self_esteem:neg_affect -0.01176    0.00755   -1.56   0.1202   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.48 on 276 degrees of freedom
## Multiple R-squared:  0.253,  Adjusted R-squared:  0.245 
## F-statistic: 31.2 on 3 and 276 DF,  p-value: <2e-16
# Dropping it recovers Example 1 (m3) from the Multiple Linear Regression
# section above exactly:
identical(round(unname(coef(m3)), 4),
          round(unname(coef(lm(depression ~ self_esteem + neg_affect, data = d))), 4))
## [1] TRUE

6 Assessing model fit

# R2 / adjusted R2 by hand, verified against summary()
r2 <- summary(m3)$r.squared
n <- nrow(d); k <- 2
adj_r2 <- 1 - (1 - r2) * (n - 1) / (n - k - 1)
c(R2 = r2, adjR2_handcalc = adj_r2, adjR2_lm = summary(m3)$adj.r.squared)
##             R2 adjR2_handcalc       adjR2_lm 
##         0.2464         0.2410         0.2410
# overall F-test
summary(m3)$fstatistic
##  value  numdf  dendf 
##  45.29   2.00 277.00

7 Relative importance of predictors

7.1 Standardized betas (recap)

lm.beta(m3)
## 
## Call:
## lm(formula = depression ~ self_esteem + neg_affect, data = d)
## 
## Standardized Coefficients::
## (Intercept) self_esteem  neg_affect 
##          NA     -0.1803      0.4071

7.2 Partial and part (semipartial) correlations

pcor(d[, c("depression", "self_esteem", "neg_affect")])$estimate
##             depression self_esteem neg_affect
## depression      1.0000     -0.1925     0.4050
## self_esteem    -0.1925      1.0000    -0.2165
## neg_affect      0.4050     -0.2165     1.0000
spcor(d[, c("depression", "self_esteem", "neg_affect")])$estimate
##             depression self_esteem neg_affect
## depression      1.0000     -0.1703     0.3846
## self_esteem    -0.1818      1.0000    -0.2056
## neg_affect      0.3826     -0.1916     1.0000

7.2.1 Full worked derivation (not just the final numbers)

For pedagogical rigor, the slides shows every step of this calculation – reproduced here in full so the arithmetic can be checked line by line.

r_dep_se <- cor(d$depression, d$self_esteem)
r_dep_na <- cor(d$depression, d$neg_affect)
r_se_na  <- cor(d$self_esteem, d$neg_affect)
cat(sprintf("r(dep,se) = %.3f\nr(dep,na) = %.3f\nr(se,na)  = %.3f\n",
            r_dep_se, r_dep_na, r_se_na))
## r(dep,se) = -0.314
## r(dep,na) = 0.466
## r(se,na)  = -0.328
# Partial correlation of self_esteem with depression, controlling neg_affect
partial_se <- (r_dep_se - r_dep_na * r_se_na) /
  sqrt((1 - r_se_na^2) * (1 - r_dep_na^2))
cat(sprintf(
  "\npartial r = (%.3f - (%.3f)(%.3f)) / sqrt((1-%.3f^2)(1-%.3f^2))\n          = (%.3f - %.3f) / sqrt(%.3f * %.3f)\n          = %.3f / %.3f = %.3f\n",
  r_dep_se, r_dep_na, r_se_na, r_se_na, r_dep_na,
  r_dep_se, r_dep_na * r_se_na, 1 - r_se_na^2, 1 - r_dep_na^2,
  r_dep_se - r_dep_na * r_se_na,
  sqrt((1 - r_se_na^2) * (1 - r_dep_na^2)), partial_se))
## 
## partial r = (-0.314 - (0.466)(-0.328)) / sqrt((1--0.328^2)(1-0.466^2))
##           = (-0.314 - -0.153) / sqrt(0.892 * 0.783)
##           = -0.161 / 0.836 = -0.192
# Part (semipartial) correlation: same numerator, smaller (X-only) denominator
part_se <- (r_dep_se - r_dep_na * r_se_na) / sqrt(1 - r_se_na^2)
cat(sprintf(
  "\npart r = (%.3f - %.3f) / sqrt(%.3f) = %.3f\n",
  r_dep_se, r_dep_na * r_se_na, 1 - r_se_na^2, part_se))
## 
## part r = (-0.314 - -0.153) / sqrt(0.892) = -0.170
# Cross-check against ppcor
stopifnot(abs(partial_se - pcor(d[, c("depression","self_esteem","neg_affect")])$estimate["self_esteem","depression"]) < 1e-2)

7.3 Incremental R2

r2_full <- summary(m3)$r.squared
r2_se_only <- summary(lm(depression ~ self_esteem, d))$r.squared
r2_na_only <- summary(lm(depression ~ neg_affect, d))$r.squared
c(delta_R2_adding_NA_after_SE = r2_full - r2_se_only,
  delta_R2_adding_SE_after_NA = r2_full - r2_na_only)
## delta_R2_adding_NA_after_SE delta_R2_adding_SE_after_NA 
##                     0.14789                     0.02899

7.4 LMG relative importance (relaimpo)

calc.relimp(m3, type = c("lmg", "first", "last", "betasq", "pratt"))
## Response variable: depression 
## Total response variance: 26.52 
## Analysis based on 280 observations 
## 
## 2 Regressors: 
## self_esteem neg_affect 
## Proportion of variance explained by model: 24.64%
## Metrics are not normalized (rela=FALSE). 
## 
## Relative importance metrics: 
## 
##                 lmg    last   first  betasq   pratt
## self_esteem 0.06376 0.02899 0.09853 0.03249 0.05658
## neg_affect  0.18266 0.14789 0.21743 0.16575 0.18984
## 
## Average coefficients for different model sizes: 
## 
##                  1X     2Xs
## self_esteem -0.3121 -0.1792
## neg_affect   0.3431  0.2996

7.5 Hoffman’s product measure (and its limits)

std <- lm.beta(m3)$standardized.coefficients
r_se <- cor(d$self_esteem, d$depression)
r_na <- cor(d$neg_affect, d$depression)
hoffman <- c(self_esteem = std["self_esteem"] * r_se,
             neg_affect  = std["neg_affect"]  * r_na)
hoffman
## self_esteem.self_esteem   neg_affect.neg_affect 
##                 0.05658                 0.18984
sum(hoffman); r2_full
## [1] 0.2464
## [1] 0.2464

8 Standardization

8.1 Worked example: standardizing one score by hand

xbar <- mean(d$self_esteem); sx <- sd(d$self_esteem)
one_person <- 35
z <- (one_person - xbar) / sx
cat(sprintf("z = (%.0f - %.2f) / %.2f = %.2f\n", one_person, xbar, sx, z))
## z = (35 - 28.68) / 5.18 = 1.22
# Full-sample standardized model, for comparison with lm.beta(m3) above
d_std <- d %>% mutate(across(c(self_esteem, neg_affect, depression), scale))
m3_std <- lm(depression ~ self_esteem + neg_affect, data = d_std)
coef(m3_std)
## (Intercept) self_esteem  neg_affect 
##   1.358e-15  -1.803e-01   4.071e-01

9 Assumptions: LINE diagnostics on the real model

par(mfrow = c(2, 2))
plot(m6)

par(mfrow = c(1, 1))

shapiro.test(resid(m6))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m6)
## W = 0.99, p-value = 0.3

9.1 A fix, worked: log-transforming a right-skewed predictor

New illustrative variable (not the running clinical example): a right-skewed stress_events count whose true relationship with depression is with its log, not its raw value. This reproduces the before/after figure used in the slides.

set.seed(42)
n <- 280
stress_events <- pmin(pmax(round(rgamma(n, shape = 1.4, scale = 3.2)), 0), 30)
depression_sim <- pmin(pmax(3 + 3.1 * log1p(stress_events) + rnorm(n, 0, 2.2), 0), 27)
dt <- data.frame(stress_events, depression = depression_sim)

m_raw <- lm(depression ~ stress_events, data = dt)
m_log <- lm(depression ~ log1p(stress_events), data = dt)

summary(m_raw)$r.squared; summary(m_log)$r.squared
## [1] 0.4365
## [1] 0.4983
par(mfrow = c(2, 2))
plot(dt$stress_events, dt$depression, pch = 16, col = adjustcolor(ntnublue, 0.5),
     xlab = "Stressful life events (raw)", ylab = "Depression", main = "Raw")
abline(m_raw, col = red, lwd = 2)
plot(fitted(m_raw), resid(m_raw), pch = 16, col = adjustcolor(ntnublue, 0.5),
     xlab = "Fitted", ylab = "Residual", main = "Raw: residuals vs. fitted")
abline(h = 0, col = red)
plot(log1p(dt$stress_events), dt$depression, pch = 16, col = adjustcolor(ntnublue, 0.5),
     xlab = "log(stress+1)", ylab = "Depression", main = "Log-transformed")
abline(m_log, col = red, lwd = 2)
plot(fitted(m_log), resid(m_log), pch = 16, col = adjustcolor(ntnublue, 0.5),
     xlab = "Fitted", ylab = "Residual", main = "Transformed: residuals vs. fitted")
abline(h = 0, col = red)

par(mfrow = c(1, 1))

10 Suppression effects

suppression_report <- function(path, y, x, z) {
  dd <- read.csv(path)
  cat("\nCorrelation matrix:\n")
  print(round(cor(dd[, c(y, x, z)]), 3))
  m_simple <- lm(reformulate(x, y), dd)
  m_full <- lm(reformulate(c(x, z), y), dd)
  cat("\nSimple model:\n"); print(coef(summary(m_simple)))
  cat("\nFull model:\n"); print(coef(summary(m_full)))
  cat("\nStandardized betas (full model):\n"); print(lm.beta(m_full)$standardized.coefficients)
  cat("\nR2 simple:", summary(m_simple)$r.squared, " | R2 full:", summary(m_full)$r.squared, "\n")
}

10.1 Classical suppression

suppression_report("../source_material/data/suppression_classical.csv",
                    "depression", "neg_affect", "social_desirability")
## 
## Correlation matrix:
##                     depression neg_affect social_desirability
## depression               1.000      0.502              -0.001
## neg_affect               0.502      1.000               0.398
## social_desirability     -0.001      0.398               1.000
## 
## Simple model:
##             Estimate Std. Error t value  Pr(>|t|)
## (Intercept)   2.0639    0.90998   2.268 2.418e-02
## neg_affect    0.3607    0.03943   9.149 2.179e-17
## 
## Full model:
##                     Estimate Std. Error t value  Pr(>|t|)
## (Intercept)           4.5560    1.07050   4.256 2.958e-05
## neg_affect            0.4289    0.04167  10.293 6.586e-21
## social_desirability  -0.1996    0.04856  -4.109 5.401e-05
## 
## Standardized betas (full model):
##         (Intercept)          neg_affect social_desirability 
##                  NA              0.5972             -0.2384 
## 
## R2 simple: 0.2523  | R2 full: 0.3002

10.2 Net suppression

suppression_report("../source_material/data/suppression_net.csv",
                    "depression", "rumination", "neg_affect")
## 
## Correlation matrix:
##            depression rumination neg_affect
## depression      1.000      0.304      0.507
## rumination      0.304      1.000      0.697
## neg_affect      0.507      0.697      1.000
## 
## Simple model:
##             Estimate Std. Error t value  Pr(>|t|)
## (Intercept)   5.3728    0.97062   5.535 7.888e-08
## rumination    0.3856    0.07678   5.022 9.768e-07
## 
## Full model:
##             Estimate Std. Error t value  Pr(>|t|)
## (Intercept)   2.2831    0.96795   2.359 1.912e-02
## rumination   -0.1229    0.09675  -1.270 2.051e-01
## neg_affect    0.4178    0.05542   7.539 8.985e-13
## 
## Standardized betas (full model):
## (Intercept)  rumination  neg_affect 
##          NA    -0.09686     0.57477 
## 
## R2 simple: 0.09232  | R2 full: 0.2621

10.3 Cooperative suppression

suppression_report("../source_material/data/suppression_cooperative.csv",
                    "well_being", "self_esteem", "optimism")
## 
## Correlation matrix:
##             well_being self_esteem optimism
## well_being       1.000       0.496    0.511
## self_esteem      0.496       1.000   -0.395
## optimism         0.511      -0.395    1.000
## 
## Simple model:
##             Estimate Std. Error t value  Pr(>|t|)
## (Intercept)    8.399    1.30963   6.414 7.139e-10
## self_esteem    0.400    0.04451   8.988 6.607e-17
## 
## Full model:
##             Estimate Std. Error t value  Pr(>|t|)
## (Intercept) -18.9264    1.09641  -17.26 2.453e-44
## self_esteem   0.6668    0.02254   29.59 3.699e-83
## optimism      5.5967    0.18669   29.98 2.882e-84
## 
## Standardized betas (full model):
## (Intercept) self_esteem    optimism 
##          NA      0.8263      0.8373 
## 
## R2 simple: 0.2457  | R2 full: 0.8374

11 Bridge to lavaan

m_lm <- lm(depression ~ self_esteem, data = d)
m_lavaan <- sem('depression ~ self_esteem', data = d, meanstructure = TRUE)

summary(m_lm)$coefficients
##             Estimate Std. Error t value  Pr(>|t|)
## (Intercept)  18.0614    1.64996  10.947 2.001e-23
## self_esteem  -0.3121    0.05662  -5.512 8.088e-08
parTable(m_lavaan)
##   id         lhs op         rhs user block group free ustart exo label plabel
## 1  1  depression  ~ self_esteem    1     1     1    1     NA   0         .p1.
## 2  2  depression ~~  depression    0     1     1    2     NA   0         .p2.
## 3  3 self_esteem ~~ self_esteem    0     1     1    0     NA   1         .p3.
## 4  4  depression ~1                0     1     1    3     NA   0         .p4.
## 5  5 self_esteem ~1                0     1     1    0     NA   1         .p5.
##    start    est    se
## 1 -0.312 -0.312 0.056
## 2 23.826 23.826 2.014
## 3 26.734 26.734 0.000
## 4 18.061 18.061 1.644
## 5 28.678 28.678 0.000

12 Session info

sessionInfo()
## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS 27.0
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Oslo
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] kableExtra_1.4.0 knitr_1.51       broom_1.0.12     lavaan_0.6-21   
##  [5] relaimpo_2.2-7   mitools_2.4      survey_4.5       survival_3.8-6  
##  [9] Matrix_1.7-5     boot_1.3-32      ppcor_1.1        MASS_7.3-65     
## [13] lm.beta_1.7-3    lubridate_1.9.5  forcats_1.0.1    stringr_1.6.0   
## [17] dplyr_1.2.0      purrr_1.2.1      readr_2.2.0      tidyr_1.3.2     
## [21] tibble_3.3.1     ggplot2_4.0.2    tidyverse_2.0.0  pacman_0.5.1    
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       xfun_0.57          bslib_0.10.0       lattice_0.22-9    
##  [5] tzdb_0.5.0         quadprog_1.5-8     vctrs_0.7.2        tools_4.4.2       
##  [9] generics_0.1.4     parallel_4.4.2     stats4_4.4.2       pkgconfig_2.0.3   
## [13] RColorBrewer_1.1-3 S7_0.2.1           lifecycle_1.0.5    compiler_4.4.2    
## [17] farver_2.1.2       textshaping_1.0.5  mnormt_2.1.2       htmltools_0.5.9   
## [21] sass_0.4.10        yaml_2.3.12        pillar_1.11.1      jquerylib_0.1.4   
## [25] cachem_1.1.0       nlme_3.1-168       tidyselect_1.2.1   digest_0.6.39     
## [29] stringi_1.8.7      labeling_0.4.3     splines_4.4.2      fastmap_1.2.0     
## [33] cli_3.6.5          magrittr_2.0.4     pbivnorm_0.6.0     corpcor_1.6.10    
## [37] withr_3.0.2        scales_1.4.0       backports_1.5.0    timechange_0.4.0  
## [41] rmarkdown_2.30     otel_0.2.0         hms_1.1.4          evaluate_1.0.5    
## [45] viridisLite_0.4.3  mgcv_1.9-4         rlang_1.1.7        Rcpp_1.1.1        
## [49] xtable_1.8-8       glue_1.8.0         DBI_1.3.0          xml2_1.5.2        
## [53] svglite_2.2.2      rstudioapi_0.18.0  jsonlite_2.0.0     R6_2.6.1          
## [57] systemfonts_1.3.2