Use beginning salary (salbegin) and age
(age) to predict current salary (salary). This
document includes the regression model, analysis, coefficient
interpretations, model-fit interpretation, a multicollinearity check,
and an interaction test.
Uncomment one import command below and replace the
example path with the location of your data. The imported data frame
must be named dat.
library(haven)
dat <- read_sav(
"/Users/carlos/Library/CloudStorage/OneDrive-TexasTechUniversity/QM 3/Lecture Notes & Assignments/Week 2/LN2-1 data_Multiple Regression Lab_bankdata.sav"
)
The following code confirms that the required variables exist. It
then uses listwise deletion, meaning that a participant is excluded if
salary, salbegin, or age is
missing.
required_variables <- c("salary", "salbegin", "age")
if (!exists("dat")) {
stop("Import your data as an object named `dat` before knitting this document.")
}
if (!all(required_variables %in% names(dat))) {
stop("The data must contain salary, salbegin, and age.")
}
reg_dat <- dat[complete.cases(dat[required_variables]), required_variables]
cat("Complete cases included in the analysis:", nrow(reg_dat))
## Complete cases included in the analysis: 473
The population structural model is:
\[ Y_i = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \varepsilon_i \]
For this analysis:
\[ \text{Salary}_i = \beta_0 + \beta_1(\text{Beginning Salary}_i) + \beta_2(\text{Age}_i) + \varepsilon_i \]
where:
Both predictors are entered into the model simultaneously. In R,
lm() includes an intercept by default.
salary_model <- lm(salary ~ salbegin + age, data = reg_dat)
model_summary <- summary(salary_model)
model_summary
##
## Call:
## lm(formula = salary ~ salbegin + age, data = reg_dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32950 -4241 -889 2396 49345
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.161e+04 1.715e+03 6.770 3.85e-11 ***
## salbegin 1.907e+00 4.549e-02 41.928 < 2e-16 ***
## age -1.976e+02 3.037e+01 -6.505 1.99e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7785 on 470 degrees of freedom
## Multiple R-squared: 0.7934, Adjusted R-squared: 0.7926
## F-statistic: 902.7 on 2 and 470 DF, p-value: < 2.2e-16
anova(salary_model)
## Analysis of Variance Table
##
## Response: salary
## Df Sum Sq Mean Sq F value Pr(>F)
## salbegin 1 1.0686e+11 1.0686e+11 1763.022 < 2.2e-16 ***
## age 1 2.5652e+09 2.5652e+09 42.321 1.986e-10 ***
## Residuals 470 2.8488e+10 6.0613e+07
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(salary_model, level = 0.95)
## 2.5 % 97.5 %
## (Intercept) 8242.595668 14984.374985
## salbegin 1.818027 1.996814
## age -257.272094 -137.905739
The estimated regression equation is:
\[ \widehat{\text{Salary}} = 11613.49 + 1.907(\text{Beginning Salary}) + -197.589(\text{Age}) \]
The intercept is 11613.49. This is the predicted current salary when beginning salary and age both equal zero. Because an age and beginning salary of zero may not be meaningful in this context, the intercept primarily anchors the regression equation.
Holding age constant, a one-unit increase in beginning salary is associated with an estimated 1.907-unit increase in predicted current salary (p < .001). If salary is measured in dollars, this means a one-dollar increase in beginning salary corresponds to an estimated $1.907 increase in current salary.
Holding beginning salary constant, each additional year of age is associated with an estimated 197.59-unit decrease in predicted current salary (p < .001).
The model had an \(R^2\) of 0.793 and an adjusted \(R^2\) of 0.793. Therefore, beginning salary and age together explained approximately 79.3% of the variance in participants’ current salaries.
The overall regression model was statistically significant, \(F(2, 470) = 902.67\), p < .001. This indicates that the model containing beginning salary and age predicted current salary significantly better than an intercept-only model. The F test establishes that at least one regression slope is not zero; the individual coefficient tests determine which predictor or predictors are statistically significant.
The following plots assess linearity, residual normality, equal residual variance, and influential observations.
par(mfrow = c(2, 2))
plot(salary_model)
par(mfrow = c(1, 1))
The table below displays cases with absolute standardized residuals of 3 or greater. An empty table indicates that no observations met this screening criterion.
diagnostics <- data.frame(
case = rownames(model.frame(salary_model)),
predicted_salary = fitted(salary_model),
residual = resid(salary_model),
standardized_residual = rstandard(salary_model),
studentized_residual = rstudent(salary_model),
leverage = hatvalues(salary_model),
cooks_distance = cooks.distance(salary_model)
)
unusual_cases <- diagnostics[
abs(diagnostics$standardized_residual) >= 3,
,
drop = FALSE
]
knitr::kable(unusual_cases, digits = 3)
| case | predicted_salary | residual | standardized_residual | studentized_residual | leverage | cooks_distance | |
|---|---|---|---|---|---|---|---|
| 127 | 127 | 45196.32 | 24803.68 | 3.192 | 3.224 | 0.004 | 0.014 |
| 128 | 128 | 63316.82 | 27308.18 | 3.526 | 3.569 | 0.010 | 0.043 |
| 146 | 146 | 33554.21 | 46445.79 | 5.975 | 6.209 | 0.003 | 0.037 |
| 151 | 151 | 44998.74 | 38751.26 | 4.987 | 5.119 | 0.004 | 0.032 |
| 221 | 221 | 59366.99 | 31883.01 | 4.111 | 4.183 | 0.008 | 0.045 |
| 274 | 274 | 69303.19 | 27696.81 | 3.581 | 3.627 | 0.013 | 0.058 |
| 299 | 299 | 54404.77 | 49345.23 | 6.357 | 6.642 | 0.006 | 0.080 |
| 317 | 317 | 87172.78 | 23452.22 | 3.057 | 3.085 | 0.029 | 0.093 |
| 330 | 330 | 91527.08 | -25527.08 | -3.337 | -3.373 | 0.034 | 0.132 |
| 375 | 375 | 99700.14 | -32950.14 | -4.337 | -4.421 | 0.048 | 0.313 |
correlation_test <- cor.test(
reg_dat$salbegin,
reg_dat$age,
method = "pearson"
)
correlation_test
##
## Pearson's product-moment correlation
##
## data: reg_dat$salbegin and reg_dat$age
## t = -0.19461, df = 471, p-value = 0.8458
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.09904778 0.08125963
## sample estimates:
## cor
## -0.008966964
The correlation between beginning salary and age was \(r = -0.009\), p = 0.846.
plot(
reg_dat$salbegin,
reg_dat$age,
xlab = "Beginning salary",
ylab = "Age",
main = "Beginning Salary and Age",
pch = 19,
col = rgb(0.2, 0.4, 0.7, 0.50)
)
abline(
lm(age ~ salbegin, data = reg_dat),
col = "firebrick",
lwd = 2
)
For each predictor, the variance inflation factor (VIF) is calculated as \(1/(1-R_j^2)\), where \(R_j^2\) comes from predicting that variable with the other predictor.
vif_values <- c(
salbegin = 1 / (1 - summary(lm(salbegin ~ age, data = reg_dat))$r.squared),
age = 1 / (1 - summary(lm(age ~ salbegin, data = reg_dat))$r.squared)
)
knitr::kable(
data.frame(Predictor = names(vif_values), VIF = unname(vif_values)),
digits = 3
)
| Predictor | VIF |
|---|---|
| salbegin | 1 |
| age | 1 |
The VIF values are 1.000 and 1.000. Values close to 1 indicate that the predictors share little variance and that multicollinearity is not a concern. VIF values greater than 5 are commonly used as a warning sign, although the chosen threshold should be justified.
The predictors are mean-centered before creating their interaction. Centering makes the lower-order coefficients easier to interpret and reduces nonessential multicollinearity between the product term and its components.
reg_dat$cent_age <- as.numeric(
scale(reg_dat$age, center = TRUE, scale = FALSE)
)
reg_dat$cent_salbegin <- as.numeric(
scale(reg_dat$salbegin, center = TRUE, scale = FALSE)
)
reg_dat$interaction <- reg_dat$cent_age * reg_dat$cent_salbegin
The first model contains the centered main effects. The second model adds the interaction term.
main_effects_model <- lm(
salary ~ cent_age + cent_salbegin,
data = reg_dat
)
interaction_model <- lm(
salary ~ cent_age + cent_salbegin + interaction,
data = reg_dat
)
summary(interaction_model)
##
## Call:
## lm(formula = salary ~ cent_age + cent_salbegin + interaction,
## data = reg_dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24136 -3862 -1038 2694 48359
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.439e+04 3.397e+02 101.216 < 2e-16 ***
## cent_age -2.912e+02 3.157e+01 -9.225 < 2e-16 ***
## cent_salbegin 2.014e+00 4.560e-02 44.168 < 2e-16 ***
## interaction -3.683e-02 5.065e-03 -7.271 1.51e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7388 on 469 degrees of freedom
## Multiple R-squared: 0.8144, Adjusted R-squared: 0.8132
## F-statistic: 685.8 on 3 and 469 DF, p-value: < 2.2e-16
interaction_change_test <- anova(main_effects_model, interaction_model)
interaction_change_test
## Analysis of Variance Table
##
## Model 1: salary ~ cent_age + cent_salbegin
## Model 2: salary ~ cent_age + cent_salbegin + interaction
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 470 2.8488e+10
## 2 469 2.5602e+10 1 2885787896 52.863 1.508e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Adding the beginning-salary-by-age interaction explained an additional 2.1% of the variance in current salary, \(\Delta R^2 = 0.021\). The change in model fit was statistically significant, \(F(1, 469) = 52.86\), p < .001. Thus, the results provide evidence that the relationship between beginning salary and current salary differs by age.
sessionInfo()
## R version 4.5.1 (2025-06-13)
## Platform: aarch64-apple-darwin20
## Running under: macOS 27.0
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## 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: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] haven_2.5.5
##
## loaded via a namespace (and not attached):
## [1] vctrs_0.7.3 cli_3.6.6 knitr_1.51 rlang_1.3.0
## [5] xfun_0.60 otel_0.2.0 forcats_1.0.1 jsonlite_2.0.0
## [9] glue_1.8.1 htmltools_0.5.9 sass_0.4.10 hms_1.1.4
## [13] rmarkdown_2.31 evaluate_1.0.5 jquerylib_0.1.4 tibble_3.3.1
## [17] tzdb_0.5.0 fastmap_1.2.0 yaml_2.3.12 lifecycle_1.0.5
## [21] compiler_4.5.1 pkgconfig_2.0.3 rstudioapi_0.19.0 digest_0.6.39
## [25] R6_2.6.1 readr_2.2.0 pillar_1.11.1 magrittr_2.0.5
## [29] bslib_0.11.0 tools_4.5.1 cachem_1.1.0