1 Scenario and Dataset

An outpatient clinic offers a course of cognitive behavioral therapy (CBT) for depression and anxiety, delivered either in person or by teletherapy. Some clients complete the course; others drop out. The clinic wants to know what predicts completion: i) symptom severity at intake, ii) social support, iii) therapy format, and iv) whether format changes how much severity matters.

You have anonymized records for 300 clients. The file therapy_completion_lab.csv has one row per client.

Variable	Type	Description
client_id	Integer	Anonymous client identifier, 1-300
baseline_severity	Continuous	Symptom severity at intake, 0-40 (higher = more severe)
support_contacts	Continuous	Weekly contacts with a supportive friend, family member or peer
therapy_format	Categorical (2 levels)	“In-Person” or “Teletherapy”
completed	Binary (0/1)	1 = completed the course, 0 = dropped out

2 Questions

Work through the tasks in order. Each model builds on the one before it.

Part A – Getting to know the data

Report the sample size, the completion rate (overall and for each therapy format), and the mean and standard deviation of baseline_severity and support_contacts for each therapy format.
Plot completed against baseline_severity, with a separate curve for each therapy format. What does the plot suggest? Also report the correlation between baseline_severity and support_contacts, and say what it means for models that include both.

Part B – Research question 1: does baseline severity predict completion?

Fit a logistic regression predicting completed from baseline_severity (Model 1). Report the intercept and slope on the log-odds scale, and the odds ratio for severity with its 95% confidence interval.
In one sentence a clinician would understand, explain what the odds ratio means.
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: does social support add to severity?

Add support_contacts (Model 2). Report the new coefficients and odds ratios.
Does adding support_contacts improve the model? Show the two deviances, the degrees of freedom, the \(\chi^2\) statistic and the \(p\)-value.
A classmate says: “Once you control for social support, baseline severity no longer matters.” Use Model 2 to explain why this is wrong.

Part D – Research question 3: is therapy format related to completion, controlling for severity and support?

Add therapy_format (Model 3). Which category is coded 0 and which is coded 1? What does the intercept now represent?
Report the odds ratio for therapy_format with its confidence interval, and test whether Model 3 improves on Model 2, as in Question 7.
Using Model 3, compute the predicted probability of completion for an In-Person client and a Teletherapy client, both with average severity and average support.

Part E – Checking the model

Which client has the largest Cook’s distance in Model 3? Is this client unduly influential? Use a threshold to decide.
Using a threshold of 0.5, report Model 3’s accuracy, sensitivity and specificity. Compare the accuracy with always guessing the most common outcome.

Part F – Research question 4: does the effect of severity depend on therapy format?

Add the interaction between baseline_severity and therapy_format (Model 4). Does it improve the model over Model 3? Show the deviance comparison as in Question 7.
From Model 4, compute the slope of baseline_severity separately for In-Person and Teletherapy clients, and convert each to an odds ratio. Explain what the two odds ratios mean for a clinic director.
Plot the predicted probability of completion against baseline_severity, one line per therapy format, with support_contacts at its mean. Would Model 3 (no interaction) have told the director everything worth knowing?

Part G – Reporting

Write a short results paragraph for the clinic director: what can and cannot be concluded about severity, support and therapy format? Remember that clients were not randomly assigned to a therapy format.

3 Guided Walkthrough

3.1 Understanding the problem

The outcome, completed, is yes/no (1/0), so we use logistic regression, exactly as in the lecture. The model is a straight line on the log-odds scale, and every result can be carried to two other scales:

odds \(=\exp(\text{log-odds})\)
probability \(=\text{odds}/(1+\text{odds})\)
odds ratio \(=\exp(B)\): the number the odds are multiplied by when a predictor increases by one unit

We answer the clinic’s questions one at a time, adding one thing to the model at each step and testing whether it helps with a likelihood-ratio test.

3.2 Inspecting the data

d <- read.csv("data/therapy_completion_lab.csv", stringsAsFactors = TRUE)
d$therapy_format <- relevel(d$therapy_format, ref = "In-Person")
str(d)

## 'data.frame':    300 obs. of  5 variables:
##  $ client_id        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ baseline_severity: num  21.9 25.6 32.2 36.7 23.4 12.6 20.9 20.6 25.7 19.9 ...
##  $ support_contacts : num  4.5 4.2 2.1 0.6 4.5 6.8 6.9 4.4 4.4 3.7 ...
##  $ therapy_format   : Factor w/ 2 levels "In-Person","Teletherapy": 1 2 1 2 1 1 2 2 1 2 ...
##  $ completed        : int  0 0 1 0 1 1 0 1 1 1 ...

d %>%
  summarise(n = n(), completion_rate = mean(completed)) %>%
  kable(digits = 3)

n	completion_rate
300	0.647

d %>%
  group_by(therapy_format) %>%
  summarise(n = n(),
            completion_rate = mean(completed),
            severity_M = mean(baseline_severity), severity_SD = sd(baseline_severity),
            support_M = mean(support_contacts), support_SD = sd(support_contacts)) %>%
  kable(digits = 2)

therapy_format	n	completion_rate	severity_M	severity_SD	support_M	support_SD
In-Person	160	0.84	20.08	8.05	3.59	1.70
Teletherapy	140	0.42	20.43	8.31	3.48	1.64

The two formats are very similar in severity and support, but very different in completion: 84.4% of In-Person clients completed, against 42.1% of Teletherapy clients. Clients were not randomly assigned to a format, so this difference alone does not show that the format itself is responsible.

ggplot(d, aes(baseline_severity, completed, color = therapy_format)) +
  geom_jitter(height = 0.04, alpha = 0.6) +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  scale_color_manual(values = c("In-Person" = "grey30", "Teletherapy" = ntnublue)) +
  labs(x = "Baseline severity", y = "Completed (0/1, jittered)",
       title = "Baseline severity and completion, by therapy format", color = NULL)

Completion falls as severity rises, and the fall looks steeper for Teletherapy. Model 4 will test that formally.

cor(d$baseline_severity, d$support_contacts)

## [1] -0.1199

The correlation is small, so severity and support carry mostly separate information. Adding one to a model should not change much about the other.

3.3 Model 1 – does baseline severity predict completion?

m1 <- glm(completed ~ baseline_severity, data = d, family = binomial(link = "logit"))
summary(m1)

## 
## Call:
## glm(formula = completed ~ baseline_severity, family = binomial(link = "logit"), 
##     data = d)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         2.5290     0.4027    6.28  3.4e-10 ***
## baseline_severity  -0.0912     0.0175   -5.22  1.8e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 389.69  on 299  degrees of freedom
## Residual deviance: 357.87  on 298  degrees of freedom
## AIC: 361.9
## 
## Number of Fisher Scoring iterations: 3

exp(coef(m1))

##       (Intercept) baseline_severity 
##           12.5408            0.9128

exp(confint(m1))

##                   2.5 %  97.5 %
## (Intercept)       5.875 28.6071
## baseline_severity 0.881  0.9436

The slope is -0.091 on the log-odds scale. Its odds ratio is 0.913: each extra point of severity multiplies the odds of completing by 0.913, a 8.7% drop in the odds per point.

To see the model on all three scales, take a client with severity 10 and one with severity 30:

tibble(baseline_severity = c(10, 30)) %>%
  mutate(log_odds    = predict(m1, newdata = ., type = "link"),
         odds        = exp(log_odds),
         probability = odds / (1 + odds)) %>%
  kable(digits = 3)

baseline_severity	log_odds	odds	probability
10	1.617	5.037	0.834
30	-0.207	0.813	0.448

Likelihood-ratio test. We compare the null model (no predictors) with Model 1. The test statistic is the drop in deviance, compared with a \(\chi^2\) distribution whose degrees of freedom equal the number of added parameters.

m0 <- glm(completed ~ 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 
##       3.897e+02       3.579e+02       3.182e+01       1.000e+00       1.693e-08

lrtest(m0, m1)

## Likelihood ratio test
## 
## Model 1: completed ~ 1
## Model 2: completed ~ baseline_severity
##   #Df LogLik Df Chisq Pr(>Chisq)    
## 1   1   -195                        
## 2   2   -179  1  31.8    1.7e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The deviance drops from 389.69 to 357.87, a difference of 31.82 on 1 degree of freedom, \(p\) < .001. Severity significantly predicts completion, in agreement with the \(z\) test in summary(m1).

3.4 Model 2 – does social support add to severity?

m2 <- glm(completed ~ baseline_severity + support_contacts, data = d, family = binomial(link = "logit"))
summary(m2)

## 
## Call:
## glm(formula = completed ~ baseline_severity + support_contacts, 
##     family = binomial(link = "logit"), data = d)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.5285     0.4889    3.13  0.00177 ** 
## baseline_severity  -0.0884     0.0177   -4.98  6.3e-07 ***
## support_contacts    0.2753     0.0822    3.35  0.00081 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 389.69  on 299  degrees of freedom
## Residual deviance: 345.98  on 297  degrees of freedom
## AIC: 352
## 
## Number of Fisher Scoring iterations: 4

exp(coef(m2))

##       (Intercept) baseline_severity  support_contacts 
##            4.6112            0.9154            1.3169

exp(confint(m2))

##                   2.5 %  97.5 %
## (Intercept)       1.803 12.3313
## baseline_severity 0.883  0.9468
## support_contacts  1.124  1.5532

lrtest(m1, m2)

## Likelihood ratio test
## 
## Model 1: completed ~ baseline_severity
## Model 2: completed ~ baseline_severity + support_contacts
##   #Df LogLik Df Chisq Pr(>Chisq)    
## 1   2   -179                        
## 2   3   -173  1  11.9    0.00056 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Adding support lowers the deviance from 357.87 to 345.98: \(\chi^2(1)\) = 11.89, \(p\) < .001. Each extra weekly support contact multiplies the odds of completing by 1.317.

The severity odds ratio barely moved, from 0.913 to 0.915, and it is still highly significant. Severity does not “stop mattering” once support is added; the two predictors are nearly unrelated, so each keeps its own effect.

3.5 Model 3 – is therapy format related to completion?

therapy_format is coded 0 = In-Person (the reference category) and 1 = Teletherapy. The intercept is now the log-odds of completion for an In-Person client with severity 0 and 0 support contacts – outside the typical data, like most raw intercepts.

m3 <- glm(completed ~ baseline_severity + support_contacts + therapy_format,
          data = d, family = binomial(link = "logit"))
summary(m3)

## 
## Call:
## glm(formula = completed ~ baseline_severity + support_contacts + 
##     therapy_format, family = binomial(link = "logit"), data = d)
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 3.1122     0.6142    5.07  4.0e-07 ***
## baseline_severity          -0.1138     0.0213   -5.34  9.2e-08 ***
## support_contacts            0.3505     0.0961    3.65  0.00027 ***
## therapy_formatTeletherapy  -2.4073     0.3280   -7.34  2.1e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 389.69  on 299  degrees of freedom
## Residual deviance: 277.42  on 296  degrees of freedom
## AIC: 285.4
## 
## Number of Fisher Scoring iterations: 5

exp(coef(m3))

##               (Intercept)         baseline_severity          support_contacts therapy_formatTeletherapy 
##                  22.47083                   0.89244                   1.41979                   0.09006

exp(confint(m3))

##                             2.5 %  97.5 %
## (Intercept)               7.02550 78.6404
## baseline_severity         0.85431  0.9290
## support_contacts          1.18158  1.7247
## therapy_formatTeletherapy 0.04603  0.1673

lrtest(m2, m3)

## Likelihood ratio test
## 
## Model 1: completed ~ baseline_severity + support_contacts
## Model 2: completed ~ baseline_severity + support_contacts + therapy_format
##   #Df LogLik Df Chisq Pr(>Chisq)    
## 1   3   -173                        
## 2   4   -139  1  68.6     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Holding severity and support constant, a Teletherapy client’s odds of completing are 0.09 times an In-Person client’s (95% CI 0.046 to 0.167) – put the other way round, an In-Person client’s odds are about 11 times higher.

An odds ratio below 1, like 0.09, is hard to picture. To compare the other way round, In-Person against Teletherapy, take the reciprocal: 1 / 0.09 = 11.1. So an In-Person client’s odds of completing are about 11 times a Teletherapy client’s. That is the same finding, just expressed from the other group’s side.

Adding format lowers the deviance from 345.98 to 277.42: \(\chi^2(1)\) = 68.56, \(p\) < .001.

newdat3 <- tibble(
  baseline_severity = mean(d$baseline_severity),
  support_contacts  = mean(d$support_contacts),
  therapy_format    = factor(levels(d$therapy_format), levels = levels(d$therapy_format))
) %>%
  mutate(log_odds    = predict(m3, newdata = ., type = "link"),
         odds        = exp(log_odds),
         probability = odds / (1 + odds))
newdat3 %>% kable(digits = 3)

baseline_severity	support_contacts	therapy_format	log_odds	odds	probability
20.24	3.542	In-Person	2.050	7.768	0.886
20.24	3.542	Teletherapy	-0.357	0.700	0.412

For an average client, the predicted probability of completing is 88.6% In-Person and 41.2% by Teletherapy.

3.6 Checking the model before trusting it

Influential clients. Cook’s distance measures how much the model would change if one client were removed. A common rule is that values above 1 deserve a closer look.

m3_aug <- augment(m3, data = d)
m3_aug %>%
  arrange(desc(.cooksd)) %>%
  select(client_id, baseline_severity, support_contacts, therapy_format, completed, .cooksd) %>%
  head(5) %>%
  kable(digits = 4)

client_id	baseline_severity	support_contacts	therapy_format	completed	.cooksd
298	6.2	1.3	In-Person	0	0.0499
42	8.1	4.2	In-Person	0	0.0473
253	36.5	3.0	Teletherapy	1	0.0377
75	23.1	0.0	Teletherapy	1	0.0333
130	12.2	1.5	In-Person	0	0.0282

The largest Cook’s distance is 0.0499 (client 298), far below 1. No single client drives the results.

Classification.

predicted_class <- factor(as.integer(predict(m3, type = "response") >= 0.5), levels = c(0, 1))
cm3 <- confusionMatrix(predicted_class, factor(d$completed, levels = c(0, 1)), positive = "1")
cm3

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0  65  26
##          1  41 168
##                                         
##                Accuracy : 0.777         
##                  95% CI : (0.725, 0.823)
##     No Information Rate : 0.647         
##     P-Value [Acc > NIR] : 7.38e-07      
##                                         
##                   Kappa : 0.495         
##                                         
##  Mcnemar's Test P-Value : 0.0872        
##                                         
##             Sensitivity : 0.866         
##             Specificity : 0.613         
##          Pos Pred Value : 0.804         
##          Neg Pred Value : 0.714         
##              Prevalence : 0.647         
##          Detection Rate : 0.560         
##    Detection Prevalence : 0.697         
##       Balanced Accuracy : 0.740         
##                                         
##        'Positive' Class : 1             
##

Model 3 classifies 77.7% of clients correctly, compared with 64.7% from always guessing “completed” (the more common outcome). Sensitivity (completers correctly identified) is 86.6%, and specificity (dropouts correctly identified) is 61.3%. These are measured on the same clients the model was fitted to, so they are somewhat optimistic.

3.7 Model 4 – does the effect of severity depend on therapy format?

m4 <- glm(completed ~ baseline_severity * therapy_format + support_contacts,
          data = d, family = binomial(link = "logit"))
summary(m4)

## 
## Call:
## glm(formula = completed ~ baseline_severity * therapy_format + 
##     support_contacts, family = binomial(link = "logit"), data = d)
## 
## Coefficients:
##                                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                   1.4844     0.7728    1.92  0.05476 .  
## baseline_severity                            -0.0473     0.0301   -1.57  0.11595    
## therapy_formatTeletherapy                     0.2222     0.9657    0.23  0.81799    
## support_contacts                              0.3750     0.0993    3.78  0.00016 ***
## baseline_severity:therapy_formatTeletherapy  -0.1226     0.0444   -2.76  0.00579 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 389.69  on 299  degrees of freedom
## Residual deviance: 269.55  on 295  degrees of freedom
## AIC: 279.6
## 
## Number of Fisher Scoring iterations: 5

exp(coef(m4))

##                                 (Intercept)                           baseline_severity 
##                                      4.4124                                      0.9538 
##                   therapy_formatTeletherapy                            support_contacts 
##                                      1.2489                                      1.4550 
## baseline_severity:therapy_formatTeletherapy 
##                                      0.8847

exp(confint(m4))

##                                              2.5 %  97.5 %
## (Intercept)                                 1.0359 21.8594
## baseline_severity                           0.8970  1.0101
## therapy_formatTeletherapy                   0.1825  8.2886
## support_contacts                            1.2043  1.7801
## baseline_severity:therapy_formatTeletherapy 0.8091  0.9639

lrtest(m3, m4)

## Likelihood ratio test
## 
## Model 1: completed ~ baseline_severity + support_contacts + therapy_format
## Model 2: completed ~ baseline_severity * therapy_format + support_contacts
##   #Df LogLik Df Chisq Pr(>Chisq)   
## 1   4   -139                       
## 2   5   -135  1  7.87      0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Adding the interaction lowers the deviance from 277.42 to 269.55: \(\chi^2(1)\) = 7.87, \(p\) = 0.005. The interaction significantly improves the model.

With an interaction, the other coefficients change meaning. baseline_severity is now the severity slope for In-Person clients (the reference group), and therapy_formatTeletherapy is the format difference at severity 0. The severity slope for Teletherapy clients is the In-Person slope plus the interaction:

b <- coef(m4)
slope_inperson    <- b[["baseline_severity"]]
slope_teletherapy <- b[["baseline_severity"]] + b[["baseline_severity:therapy_formatTeletherapy"]]
tibble(therapy_format = c("In-Person", "Teletherapy"),
       slope_log_odds = c(slope_inperson, slope_teletherapy),
       odds_ratio_per_point = exp(c(slope_inperson, slope_teletherapy))) %>%
  kable(digits = 3)

therapy_format	slope_log_odds	odds_ratio_per_point
In-Person	-0.047	0.954
Teletherapy	-0.170	0.844

For In-Person clients, each extra point of severity multiplies the odds of completing by 0.954 – a small decline, and not significant on its own (\(p\) = 0.116). For Teletherapy clients it multiplies the odds by 0.844 – a much steeper decline. Severity matters far more under teletherapy.

newdat4 <- expand.grid(
  baseline_severity = seq(min(d$baseline_severity), max(d$baseline_severity), length.out = 100),
  therapy_format    = levels(d$therapy_format),
  support_contacts  = mean(d$support_contacts)
)
newdat4$predicted <- predict(m4, newdata = newdat4, type = "response")

ggplot(newdat4, aes(baseline_severity, predicted, color = therapy_format)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = c("In-Person" = "grey30", "Teletherapy" = ntnublue)) +
  labs(x = "Baseline severity", y = "Predicted probability of completion",
       title = "Predicted completion by severity and therapy format",
       subtitle = "Support contacts held at the sample mean", color = NULL)

The two lines start close together at low severity and move apart as severity rises. Model 3 would have said “Teletherapy clients complete less often”; Model 4 adds that the gap is small for mild cases and large for severe ones.

3.8 Interpreting the findings

Severity predicts completion on its own (Model 1). Social support adds a separate contribution (Model 2). Therapy format is strongly related to completion even with severity and support held constant (Model 3). And the format gap depends on severity: under teletherapy, more severe clients are much less likely to finish (Model 4).

3.9 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, therapy completion"
)

Odds ratios (95% CI) for Models 1-4, therapy completion
	Model 1	Model 2	Model 3	Model 4
(Intercept)	12.541	4.611	22.471	4.412
	[5.875, 28.607]	[1.803, 12.331]	[7.025, 78.640]	[1.036, 21.859]
baseline_severity	0.913	0.915	0.892	0.954
	[0.881, 0.944]	[0.883, 0.947]	[0.854, 0.929]	[0.897, 1.010]
support_contacts		1.317	1.420	1.455
		[1.124, 1.553]	[1.182, 1.725]	[1.204, 1.780]
therapy_formatTeletherapy			0.090	1.249
			[0.046, 0.167]	[0.182, 8.289]
baseline_severity × therapy_formatTeletherapy				0.885
				[0.809, 0.964]
Num.Obs.	300	300	300	300
AIC	361.9	352.0	285.4	279.6
BIC	369.3	363.1	300.2	298.1

Higher baseline severity predicted lower odds of completing CBT (Model 1), and more weekly support contacts predicted higher odds of completion (Model 2). With both held constant, Teletherapy clients had much lower odds of completion than In-Person clients (Model 3). A significant severity-by-format interaction (Model 4) showed that this gap depended on severity: the two formats performed similarly for clients with mild symptoms, but under Teletherapy the odds of completion fell much more steeply as severity increased. Because clients were not randomly assigned to a format, these results show associations, not proof that the format caused the difference in completion.

4 Solutions

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

Part A

\(N\) = 300; overall completion rate 0.647. In-Person: completion 0.844; Teletherapy: completion 0.421. The means and SDs of severity and support are in the descriptives table: very similar in the two formats.
Completion falls as severity rises, and more steeply for Teletherapy. \(r\) = -0.12: small, so the two predictors carry mostly separate information and each coefficient in Models 2-4 is not distorted by the other.

Part B

Intercept = 2.529, slope = -0.091 (log-odds). OR = 0.913, 95% CI [0.881, 0.944].
Each extra point of severity at intake multiplies a client’s odds of completing the course by 0.913, about 8.7% lower odds per point.
\(z\) = -5.22, \(p\) < .001. Likelihood-ratio test: null deviance 389.69 minus Model 1 deviance 357.87 = 31.82 on 1 df, \(p\) < .001. Both say severity is significant.

Part C

Severity OR = 0.915; support OR = 1.317.
Model 1 deviance 357.87 (\(df\) = 298) minus Model 2 deviance 345.98 (\(df\) = 297) = \(\chi^2(1)\) = 11.89, \(p\) < .001. Support significantly improves the model.
The severity OR is almost unchanged (0.913 in Model 1, 0.915 in Model 2) and still highly significant. Severity still matters; support is an additional, largely separate predictor.

Part D

In-Person = 0 (reference), Teletherapy = 1. The intercept is the log-odds of completion for an In-Person client with severity 0 and 0 support contacts.
Teletherapy OR = 0.09, 95% CI [0.046, 0.167]. Model 2 deviance 345.98 minus Model 3 deviance 277.42 = \(\chi^2(1)\) = 68.56, \(p\) < .001.
In-Person: 88.6%; Teletherapy: 41.2% (average severity and support).

Part E

Client 298, Cook’s distance = 0.0499, far below the threshold of 1: not unduly influential.
Accuracy = 77.7%, sensitivity = 86.6%, specificity = 61.3%. Always guessing “completed” would give 64.7%, so the model does better.

Part F

Model 3 deviance 277.42 minus Model 4 deviance 269.55 = \(\chi^2(1)\) = 7.87, \(p\) = 0.005. The interaction significantly improves the model.
In-Person: slope -0.047, OR = 0.954 per point (small, and not significant on its own). Teletherapy: slope -0.17, OR = 0.844 per point. For the director: severity makes little difference to completion in person, but under teletherapy each extra point of severity lowers the odds of completing by about 16%.
The lines start close together and move apart as severity rises. Model 3 would have missed this: it assumes the same format gap at every severity level, when the gap is in fact small for mild cases and large for severe ones.

Part G

See the example paragraph under Reporting the results. It should state all four findings and that, because format was not randomly assigned, the results do not show that format caused the difference (for example, clients choosing teletherapy may differ in ways that were not measured).

Common errors to avoid

Treating the raw difference in completion rates between formats as evidence of a format effect, before controlling for severity and support.
Reporting only the \(p\)-value of a likelihood-ratio test, without the two deviances and the degrees of freedom.
After adding the interaction, still reading the therapy_format coefficient as “the” format effect; it is the format difference at severity 0.
Reading a log-odds coefficient as a change in probability. Exponentiate it to get an odds ratio, or compute predicted probabilities.
Concluding that therapy format causes differences in completion.

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/Zurich
## 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           e1071_1.7-17         compiler_4.4.2       mgcv_1.9-4           vctrs_0.7.2         
## [11] reshape2_1.4.5       pkgconfig_2.0.3      fastmap_1.2.0        backports_1.5.0      labeling_0.4.3      
## [16] effectsize_1.0.2     rmarkdown_2.30       prodlim_2026.03.11   tzdb_0.5.0           xfun_0.57           
## [21] cachem_1.1.0         litedown_0.9         jsonlite_2.0.0       tinytable_0.16.0     recipes_1.3.1       
## [26] parallel_4.4.2       R6_2.6.1             bslib_0.10.0         tables_0.9.33        stringi_1.8.7       
## [31] RColorBrewer_1.1-3   parallelly_1.46.1    rpart_4.1.24         jquerylib_0.1.4      estimability_1.5.1  
## [36] Rcpp_1.1.1           iterators_1.0.14     future.apply_1.20.2  parameters_0.28.3    Matrix_1.7-5        
## [41] splines_4.4.2        nnet_7.3-20          timechange_0.4.0     tidyselect_1.2.1     rstudioapi_0.18.0   
## [46] yaml_2.3.12          timeDate_4052.112    codetools_0.2-20     listenv_0.10.1       plyr_1.8.9          
## [51] withr_3.0.2          bayestestR_0.17.0    S7_0.2.1             coda_0.19-4.1        evaluate_1.0.5      
## [56] future_1.70.0        survival_3.8-6       proxy_0.4-29         pillar_1.11.1        checkmate_2.3.4     
## [61] foreach_1.5.2        stats4_4.4.2         insight_1.4.6        generics_0.1.4       hms_1.1.4           
## [66] scales_1.4.0         xtable_1.8-8         globals_0.19.1       class_7.3-23         glue_1.8.0          
## [71] emmeans_2.0.2        tools_4.4.2          data.table_1.18.2.1  ModelMetrics_1.2.2.2 gower_1.0.2         
## [76] mvtnorm_1.3-6        grid_4.4.2           datawizard_1.3.0     ipred_0.9-15         nlme_3.1-168        
## [81] performance_0.16.0   cli_3.6.5            lava_1.8.2           gtable_0.3.6         sass_0.4.10         
## [86] digest_0.6.39        TH.data_1.1-5        farver_2.1.2         htmltools_0.5.9      lifecycle_1.0.5     
## [91] hardhat_1.4.2        MASS_7.3-65

Lab Practical, Topic 2: Logistic Regression

Scenario and Dataset, Questions, Guided Walkthrough, and Solutions

Frederick Anyan – PSY3100

Teaching Assistant: Hanna Pay Eriksen

30 September 2026