A university counseling center wants to know what predicts whether a first-year student seeks counseling during the year: perceived stress, stigma beliefs about mental health treatment, and living situation (on or off campus) – and whether living situation changes how strongly stress predicts help-seeking.
You have anonymized survey and records data for 280 first-year
students. The file counseling_helpseeking_homework.csv has
one row per student. Use the lab practical as a model for how to reason
through each step, not as a template to copy: the variables are
different, and so is one of the conclusions.
| Variable | Type | Description |
|---|---|---|
| student_id | Integer | Anonymous student identifier, 1-280 |
| perceived_stress | Continuous | Perceived stress at the start of the spring term, 0-40 (higher = more stress) |
| stigma_beliefs | Continuous | Stigma beliefs about seeking mental health treatment, 0-30 (higher = more stigma) |
| living_situation | Categorical (2 levels) | “On-Campus” or “Off-Campus” |
| sought_counseling | Binary (0/1) | 1 = sought counseling during the first year, 0 = did not |
These questions ask for more independent judgment than the lab practical: in several places you choose and justify a step yourself.
Part A – Getting to know the data
perceived_stress and stigma_beliefs for each
living situation. Report the correlation between
perceived_stress and stigma_beliefs, and say
what it means for a model that includes both.sought_counseling against
perceived_stress, with a separate curve for each living
situation. What does the plot suggest?Part B – Research question 1: does perceived stress predict help-seeking?
sought_counseling
from perceived_stress (Model 1). Report the slope on the
log-odds scale and as an odds ratio with its 95% confidence
interval.summary() and a
likelihood-ratio test against the intercept-only model, showing the two
deviances, the degrees of freedom and the \(p\)-value.Part C – Research question 2: do stigma beliefs add to stress?
stigma_beliefs (Model 2). Report the new
coefficients and odds ratios.stigma_beliefs coefficient
mean?Part D – Research question 3: is living situation related to help-seeking, controlling for stress and stigma?
living_situation (Model 3). We use Off-Campus as
the reference category (coded 0) and On-Campus as 1. Report the odds
ratio with its confidence interval, and test whether Model 3 improves on
Model 2.Part E – Research question 4: does the stress effect differ by living situation?
perceived_stress and
living_situation (Model 4). Report the coefficients and
test, with a deviance comparison, whether the interaction improves the
model.perceived_stress, one line per living situation, with
stigma_beliefs at its mean. Given your answer to Question
14, is this plot evidence that living situation changes the stress
effect?Part F – Reporting
d <- read.csv("data/counseling_helpseeking_homework.csv", stringsAsFactors = TRUE)
d$living_situation <- relevel(d$living_situation, ref = "Off-Campus")
str(d)
## 'data.frame': 280 obs. of 5 variables:
## $ student_id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ perceived_stress : num 15.2 21.3 17.8 33.4 22.2 5.4 21.2 17.2 25.6 35.7 ...
## $ stigma_beliefs : num 14.2 15.3 14 16 23.9 21.1 9.7 17.3 13.1 11.6 ...
## $ living_situation : Factor w/ 2 levels "Off-Campus","On-Campus": 2 1 1 1 1 1 1 2 2 1 ...
## $ sought_counseling: int 0 1 0 1 1 1 1 0 1 0 ...
d %>%
summarise(n = n(), helpseeking_rate = mean(sought_counseling)) %>%
kable(digits = 3)
| n | helpseeking_rate |
|---|---|
| 280 | 0.425 |
d %>%
group_by(living_situation) %>%
summarise(n = n(),
helpseeking_rate = mean(sought_counseling),
stress_M = mean(perceived_stress), stress_SD = sd(perceived_stress),
stigma_M = mean(stigma_beliefs), stigma_SD = sd(stigma_beliefs)) %>%
kable(digits = 2)
| living_situation | n | helpseeking_rate | stress_M | stress_SD | stigma_M | stigma_SD |
|---|---|---|---|---|---|---|
| Off-Campus | 135 | 0.32 | 21.34 | 6.83 | 16.13 | 5.97 |
| On-Campus | 145 | 0.52 | 19.94 | 7.07 | 14.75 | 6.47 |
cor(d$perceived_stress, d$stigma_beliefs)
## [1] 0.1272
On-Campus students sought counseling more often (52.4% vs. 31.9%). Stress and stigma are weakly positively correlated: students with more stress tend to hold slightly more stigma. The correlation is small, so a model with both is fine, but each coefficient should be read as “holding the other constant”.
ggplot(d, aes(perceived_stress, sought_counseling, color = living_situation)) +
geom_jitter(height = 0.04, alpha = 0.6) +
geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
scale_color_manual(values = c("Off-Campus" = "grey30", "On-Campus" = ntnublue)) +
labs(x = "Perceived stress", y = "Sought counseling (0/1, jittered)",
title = "Perceived stress and help-seeking, by living situation", color = "Living situation")
Help-seeking rises with stress in both groups, and the On-Campus curve lies above the Off-Campus curve.
m1 <- glm(sought_counseling ~ perceived_stress, data = d, family = binomial(link = "logit"))
summary(m1)
##
## Call:
## glm(formula = sought_counseling ~ perceived_stress, family = binomial(link = "logit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1228 0.3913 -2.87 0.0041 **
## perceived_stress 0.0395 0.0178 2.22 0.0263 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 381.84 on 279 degrees of freedom
## Residual deviance: 376.78 on 278 degrees of freedom
## AIC: 380.8
##
## Number of Fisher Scoring iterations: 4
exp(coef(m1))
## (Intercept) perceived_stress
## 0.3254 1.0403
exp(confint(m1))
## 2.5 % 97.5 %
## (Intercept) 0.1487 0.6923
## perceived_stress 1.0051 1.0779
m0 <- glm(sought_counseling ~ 1, data = d, family = binomial(link = "logit"))
chisq_1 <- deviance(m0) - deviance(m1)
df_1 <- df.residual(m0) - df.residual(m1)
p_1 <- pchisq(chisq_1, df_1, lower.tail = FALSE)
c(null_deviance = deviance(m0), model1_deviance = deviance(m1), chisq = chisq_1, df = df_1, p = p_1)
## null_deviance model1_deviance chisq df p
## 381.83858 376.78082 5.05776 1.00000 0.02452
lrtest(m0, m1)
## Likelihood ratio test
##
## Model 1: sought_counseling ~ 1
## Model 2: sought_counseling ~ perceived_stress
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 1 -191
## 2 2 -188 1 5.06 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Each extra point of perceived stress multiplies the odds of seeking counseling by 1.04. The deviance drops from 381.84 to 376.78, a difference of 5.06 on 1 degree of freedom, \(p\) = 0.025. Stress is a significant but fairly modest predictor on its own.
m2 <- glm(sought_counseling ~ perceived_stress + stigma_beliefs, data = d, family = binomial(link = "logit"))
summary(m2)
##
## Call:
## glm(formula = sought_counseling ~ perceived_stress + stigma_beliefs,
## family = binomial(link = "logit"), data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1233 0.4711 -0.26 0.79347
## perceived_stress 0.0519 0.0188 2.76 0.00570 **
## stigma_beliefs -0.0828 0.0214 -3.87 0.00011 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 381.84 on 279 degrees of freedom
## Residual deviance: 360.54 on 277 degrees of freedom
## AIC: 366.5
##
## Number of Fisher Scoring iterations: 4
exp(coef(m2))
## (Intercept) perceived_stress stigma_beliefs
## 0.8840 1.0533 0.9206
exp(confint(m2))
## 2.5 % 97.5 %
## (Intercept) 0.3497 2.230
## perceived_stress 1.0158 1.094
## stigma_beliefs 0.8818 0.959
lrtest(m1, m2)
## Likelihood ratio test
##
## Model 1: sought_counseling ~ perceived_stress
## Model 2: sought_counseling ~ perceived_stress + stigma_beliefs
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 2 -188
## 2 3 -180 1 16.2 5.6e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Adding stigma lowers the deviance from 376.78 to 360.54: \(\chi^2(1)\) = 16.24, \(p\) < .001. The stigma coefficient is negative: each extra point of stigma multiplies the odds of seeking counseling by 0.921, holding stress constant. Stigma acts as a barrier to seeking help.
Notice that the stress odds ratio went up a little, from 1.04 to 1.053. The odds ratio for stress rose slightly, from 1.04 to 1.053, because more stressed students also tend to hold more stigma, and stigma lowers the odds of help-seeking, so in Model 1 part of the positive stress effect was cancelled out. Once stigma is held constant in Model 2, students with the same level of stigma are compared, and the direct link between stress and help-seeking shows more clearly. This is a small example of a suppression effect.
m4 <- glm(sought_counseling ~ perceived_stress * living_situation + stigma_beliefs,
data = d, family = binomial(link = "logit"))
summary(m4)
##
## Call:
## glm(formula = sought_counseling ~ perceived_stress * living_situation +
## stigma_beliefs, family = binomial(link = "logit"), data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0339 0.7163 -1.44 0.14889
## perceived_stress 0.0694 0.0295 2.35 0.01873 *
## living_situationOn-Campus 1.2435 0.8609 1.44 0.14860
## stigma_beliefs -0.0791 0.0218 -3.63 0.00028 ***
## perceived_stress:living_situationOn-Campus -0.0163 0.0386 -0.42 0.67299
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 381.84 on 279 degrees of freedom
## Residual deviance: 348.21 on 275 degrees of freedom
## AIC: 358.2
##
## Number of Fisher Scoring iterations: 4
exp(coef(m4))
## (Intercept) perceived_stress
## 0.3556 1.0719
## living_situationOn-Campus stigma_beliefs
## 3.4678 0.9239
## perceived_stress:living_situationOn-Campus
## 0.9838
exp(confint(m4))
## 2.5 % 97.5 %
## (Intercept) 0.08279 1.3979
## perceived_stress 1.01304 1.1382
## living_situationOn-Campus 0.65560 19.4450
## stigma_beliefs 0.88425 0.9634
## perceived_stress:living_situationOn-Campus 0.91147 1.0609
lrtest(m3, m4)
## Likelihood ratio test
##
## Model 1: sought_counseling ~ perceived_stress + stigma_beliefs + living_situation
## Model 2: sought_counseling ~ perceived_stress * living_situation + stigma_beliefs
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -174
## 2 5 -174 1 0.18 0.67
b <- coef(m4)
slope_off <- b[["perceived_stress"]]
slope_on <- b[["perceived_stress"]] + b[["perceived_stress:living_situationOn-Campus"]]
tibble(living_situation = c("Off-Campus", "On-Campus"),
slope_log_odds = c(slope_off, slope_on),
odds_ratio_per_point = exp(c(slope_off, slope_on))) %>%
kable(digits = 3)
| living_situation | slope_log_odds | odds_ratio_per_point |
|---|---|---|
| Off-Campus | 0.069 | 1.072 |
| On-Campus | 0.053 | 1.055 |
Look at the test before the table. Adding the interaction lowers the deviance only from 348.39 to 348.21: \(\chi^2(1)\) = 0.18, \(p\) = 0.673. This is not significant. The two slopes in the table are slightly different, but the interaction test is exactly the test of whether that difference is bigger than chance – and it is not. So the data give no evidence that living situation changes the stress effect, and Model 3 (without the interaction) remains the model to report.
newdat4 <- expand.grid(
perceived_stress = seq(min(d$perceived_stress), max(d$perceived_stress), length.out = 100),
living_situation = levels(d$living_situation),
stigma_beliefs = mean(d$stigma_beliefs)
)
newdat4$predicted <- predict(m4, newdata = newdat4, type = "response")
ggplot(newdat4, aes(perceived_stress, predicted, color = living_situation)) +
geom_line(linewidth = 1.2) +
scale_color_manual(values = c("Off-Campus" = "grey30", "On-Campus" = ntnublue)) +
labs(x = "Perceived stress", y = "Predicted probability of seeking counseling",
title = "Predicted help-seeking by stress and living situation",
subtitle = "Model 4 (interaction not significant); stigma held at the sample mean",
color = "Living situation")
The lines are not perfectly parallel, but that small difference is the non-significant interaction we just tested. The plot should not be described as showing that living situation changes the stress effect.
modelsummary(
list("Model 1" = m1, "Model 2" = m2, "Model 3" = m3, "Model 4" = m4),
exponentiate = TRUE,
statistic = "conf.int",
gof_map = c("nobs", "aic", "bic"),
title = "Odds ratios (95% CI) for Models 1-4, counseling help-seeking"
)
| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| (Intercept) | 0.325 | 0.884 | 0.435 | 0.356 |
| [0.149, 0.692] | [0.350, 2.230] | [0.153, 1.209] | [0.083, 1.398] | |
| perceived_stress | 1.040 | 1.053 | 1.062 | 1.072 |
| [1.005, 1.078] | [1.016, 1.094] | [1.023, 1.104] | [1.013, 1.138] | |
| stigma_beliefs | 0.921 | 0.924 | 0.924 | |
| [0.882, 0.959] | [0.885, 0.964] | [0.884, 0.963] | ||
| living_situationOn-Campus | 2.456 | 3.468 | ||
| [1.478, 4.130] | [0.656, 19.445] | |||
| perceived_stress × living_situationOn-Campus | 0.984 | |||
| [0.911, 1.061] | ||||
| Num.Obs. | 280 | 280 | 280 | 280 |
| AIC | 380.8 | 366.5 | 356.4 | 358.2 |
| BIC | 388.1 | 377.4 | 370.9 | 376.4 |
Higher perceived stress predicted higher odds of seeking counseling (Model 1). Stigma beliefs added a significant contribution, with more stigma predicting lower odds of seeking help (Model 2). With stress and stigma held constant, On-Campus students had higher odds of seeking counseling than Off-Campus students (Model 3). The stress-by-living-situation interaction was not significant (Model 4), so there is no evidence that the stress effect differs between On-Campus and Off-Campus students; Model 3 is the best description of these data. Because stress, stigma and living situation were measured rather than manipulated, these results describe associations only. They do not show that stress or living situation causes help-seeking: unmeasured factors, such as earlier experience with mental health services, could affect both where students live and whether they seek help.
The answers below use the models fitted in the Guided Walkthrough.
Part A
Part B
Part C
Part D
Part E
Part F
Common errors to avoid
sessionInfo()
## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS 27.0.1
##
## 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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] knitr_1.51 broom_1.0.12 caret_7.0-1 lattice_0.22-9 modelsummary_2.6.0 lmtest_0.9-40
## [7] zoo_1.8-15 lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0 dplyr_1.2.0 purrr_1.2.1
## [13] readr_2.2.0 tidyr_1.3.2 tibble_3.3.1 ggplot2_4.0.2 tidyverse_2.0.0 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] pROC_1.19.0.1 sandwich_3.1-1 rlang_1.1.7 magrittr_2.0.4 multcomp_1.4-30
## [6] otel_0.2.0 compiler_4.4.2 mgcv_1.9-4 vctrs_0.7.2 reshape2_1.4.5
## [11] pkgconfig_2.0.3 fastmap_1.2.0 backports_1.5.0 labeling_0.4.3 effectsize_1.0.2
## [16] rmarkdown_2.30 prodlim_2026.03.11 tzdb_0.5.0 xfun_0.57 cachem_1.1.0
## [21] litedown_0.9 jsonlite_2.0.0 tinytable_0.16.0 recipes_1.3.1 parallel_4.4.2
## [26] R6_2.6.1 bslib_0.10.0 tables_0.9.33 stringi_1.8.7 RColorBrewer_1.1-3
## [31] parallelly_1.46.1 rpart_4.1.24 jquerylib_0.1.4 estimability_1.5.1 Rcpp_1.1.1
## [36] iterators_1.0.14 future.apply_1.20.2 parameters_0.28.3 Matrix_1.7-5 splines_4.4.2
## [41] nnet_7.3-20 timechange_0.4.0 tidyselect_1.2.1 rstudioapi_0.18.0 yaml_2.3.12
## [46] timeDate_4052.112 codetools_0.2-20 listenv_0.10.1 plyr_1.8.9 withr_3.0.2
## [51] bayestestR_0.17.0 S7_0.2.1 coda_0.19-4.1 evaluate_1.0.5 future_1.70.0
## [56] survival_3.8-6 pillar_1.11.1 checkmate_2.3.4 foreach_1.5.2 stats4_4.4.2
## [61] insight_1.4.6 generics_0.1.4 hms_1.1.4 scales_1.4.0 xtable_1.8-8
## [66] globals_0.19.1 class_7.3-23 glue_1.8.0 emmeans_2.0.2 tools_4.4.2
## [71] data.table_1.18.2.1 ModelMetrics_1.2.2.2 gower_1.0.2 mvtnorm_1.3-6 grid_4.4.2
## [76] datawizard_1.3.0 ipred_0.9-15 nlme_3.1-168 performance_0.16.0 cli_3.6.5
## [81] lava_1.8.2 gtable_0.3.6 sass_0.4.10 digest_0.6.39 TH.data_1.1-5
## [86] farver_2.1.2 htmltools_0.5.9 lifecycle_1.0.5 hardhat_1.4.2 MASS_7.3-65