1 Scenario and Dataset

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

2 Questions

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

  1. Report the sample size, the help-seeking rate (overall and for each living situation), and the mean and standard deviation of 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.
  2. Plot 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?

  1. Fit a logistic regression predicting 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.
  2. In one sentence, explain what the odds ratio means.
  3. Is the effect significant? Report the \(z\) test from 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?

  1. Add stigma_beliefs (Model 2). Report the new coefficients and odds ratios.
  2. Does this improve the model? Show your work as in the lab practical. What does the sign of the stigma_beliefs coefficient mean?

Part D – Research question 3: is living situation related to help-seeking, controlling for stress and stigma?

  1. Add 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.
  2. Using Model 3, compute the predicted probability of seeking counseling for an On-Campus and an Off-Campus student, both with average stress and average stigma.

Part E – Research question 4: does the stress effect differ by living situation?

  1. Add the interaction between perceived_stress and living_situation (Model 4). Report the coefficients and test, with a deviance comparison, whether the interaction improves the model.
  2. This is where the homework differs from the lab practical. Depending on your answer to Question 13:
    • If the interaction is significant, compute and interpret the stress slope (and its odds ratio) for each living situation, as in the lab practical.
    • If it is not significant, state the correct conclusion. Should a write-up still describe the two slopes as “looking” different? Explain why or why not.
  3. Plot the predicted probability of seeking counseling against 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

  1. Write a short results paragraph for the counseling center director: what do the models support, and what do they not support? State whether these data can show that stress, or living situation, causes help-seeking.

3 Guided Walkthrough

3.1 Inspecting the data

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.

3.2 Model 1 – does perceived stress predict help-seeking?

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.

3.3 Model 2 – do stigma beliefs add to stress?

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.

3.5 Model 4 – does the stress effect differ by living situation?

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.

3.6 Reporting the results

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"
)
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.

4 Solutions

The answers below use the models fitted in the Guided Walkthrough.

Part A

  1. \(N\) = 280; overall help-seeking rate 0.425. Off-Campus: 0.319; On-Campus: 0.524. Means and SDs are in the descriptives table. \(r\)(stress, stigma) = 0.127: small, so both can go in one model, but each coefficient is “holding the other constant”.
  2. Help-seeking rises with stress in both groups, and On-Campus students are more likely to seek counseling at every stress level.

Part B

  1. Slope = 0.04 (log-odds). OR = 1.04, 95% CI [1.005, 1.078].
  2. Each extra point of perceived stress multiplies a student’s odds of seeking counseling by 1.04, about 4% higher odds per point.
  3. \(z\) = 2.22, \(p\) = 0.026. Likelihood-ratio test: null deviance 381.84 minus Model 1 deviance 376.78 = 5.06 on 1 df, \(p\) = 0.025. Both say the effect is significant.

Part C

  1. Stress OR = 1.053; stigma OR = 0.921.
  2. Model 1 deviance 376.78 minus Model 2 deviance 360.54 = \(\chi^2(1)\) = 16.24, \(p\) < .001: stigma significantly improves the model. The negative sign means that, at the same stress level, students with more stigma are less likely to seek counseling.

Part D

  1. Off-Campus is the reference (0), On-Campus is 1. The choice does not change the fit or the conclusions, only the direction in which the odds ratio is expressed. On-Campus OR = 2.456, 95% CI [1.478, 4.13]. Model 2 deviance 360.54 minus Model 3 deviance 348.39 = \(\chi^2(1)\) = 12.15, \(p\) < .001.
  2. Off-Campus: 30.8%; On-Campus: 52.3% (average stress and stigma).

Part E

  1. Interaction coefficient = -0.016 (OR = 0.984). Model 3 deviance 348.39 minus Model 4 deviance 348.21 = \(\chi^2(1)\) = 0.18, \(p\) = 0.673: not significant.
  2. The interaction is not significant, so the correct conclusion is that there is no evidence that living situation changes the stress effect. The two slopes (ORs of 1.072 and 1.055) should not be described as different: the interaction test is the test of that difference, and it was not significant. Report Model 3.
  3. No. The lines are not perfectly parallel, but that is the non-significant interaction. A write-up should describe Model 3’s result: the same stress effect in both groups, with On-Campus students more likely to seek counseling overall.

Part F

  1. See the example paragraph under Reporting the results. It should report the stress, stigma and living-situation effects, the non-significant interaction, and state that these are associations: the data cannot show that stress or living situation causes help-seeking.

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