6 D&H Ch6 - Statistical vs. Experimental Control: “PTSD”

Compiled: August 27, 2025

Darlington & Hayes, Chapter 6’s first example

# install.packages("remotes")
# remotes::install_github("sarbearschwartz/apaSupp")
# remotes::install_github("ddsjoberg/gtsummary")
     
library(tidyverse) 
library(flextable)
library(apaSupp)
library(car)
library(rempsyc)
library(parameters)
library(performance)
library(interactions)
library(ggResidpanel)

6.1 PURPOSE

RESEARCH QUESTION:

Does the experimental therapy reduce PTSD symptom severity more than the traditional therapy?

6.1.1 Data Description

Hypothetical study of the effect of a proposed therapy for Post-Traumatic Stress Disorder (PTSD) conducted on 14 military veterans who experienced combat. Half were randomly assigned to experience 6 weeks of proposed therapy (experimental) and the other half received 6 weeks of traditional therapy. Each was pre tested with respect to their PTSD symptoms and an identical post test assessment of their symptoms was administered upon the completion of the therapy.

pre_test symptom severity prior to therapy
post_test symptom severity after therapy
gain symptom severity gain
therapy which therapy was randomly assigned

Manually enter the data set provided on page 44 in Table 3.1

df_ptsd <- tibble::tribble(~id, ~post_test, ~pre_test, ~therapy,
                            1,  2,   1,  1,  
                            2,  4,   3,  1,
                            3,  6,   7,  1,
                            4,  6,  10,  1,
                            5,  9,  13,  1,
                            6, 10,  17,  1,
                            7, 12,  19,  1,
                            8,  6,   1,  0,
                            9,  7,   5,  0,
                           10,  9,   7,  0,
                           11,  9,   9,  0,
                           12, 12,  13,  0,
                           13, 12,  16,  0,
                           14, 15,  19,  0) %>% 
  dplyr::mutate(therapy = factor(therapy,
                                 levels = c(0, 1),
                                 labels = c("Traditional", "Experimental"))) %>% 
  dplyr::mutate(gain = post_test - pre_test)

df_ptsd %>% 
  dplyr::select("ID" = id,
                "Therapy"   = therapy,
                "Pre Test"  = pre_test,
                "Post Test" = post_test,
                "Gain"      = gain) %>% 
  flextable::flextable() %>% 
  apaSupp::theme_apa(caption = "Dataset, wide Format") %>% 
  flextable::colformat_double(digits = 0)

**Table 6.1**
Dataset, wide Format
ID	Therapy	Pre Test	Post Test	Gain
1	Experimental	1	2	1
2	Experimental	3	4	1
3	Experimental	7	6	-1
4	Experimental	10	6	-4
5	Experimental	13	9	-4
6	Experimental	17	10	-7
7	Experimental	19	12	-7
8	Traditional	1	6	5
9	Traditional	5	7	2
10	Traditional	7	9	2
11	Traditional	9	9	0
12	Traditional	13	12	-1
13	Traditional	16	12	-4
14	Traditional	19	15	-4

6.1.2 Longer Format of Repeated Measures

df_ptsd_long <- df_ptsd %>% 
  dplyr::select(-gain) %>% 
  tidyr::pivot_longer(cols = ends_with("test"),
                      names_to = c("time", ".value"),
                      names_pattern = "(.*)_(.*)") %>% 
  dplyr::mutate(time = factor(time,
                              levels = c("pre", "post"))) %>% 
  dplyr::arrange(id, time)

df_ptsd_long %>% 
  dplyr::select("ID" = id,
                "Therapy"    = therapy,
                "Time"       = time,
                "Test Score" = test) %>% 
  flextable::flextable() %>% 
  apaSupp::theme_apa(caption = "Dataset, Long Format") %>% 
  flextable::colformat_double(digits = 0)

**Table 6.2**
Dataset, Long Format
ID	Therapy	Time	Test Score
1	Experimental	pre	1
1	Experimental	post	2
2	Experimental	pre	3
2	Experimental	post	4
3	Experimental	pre	7
3	Experimental	post	6
4	Experimental	pre	10
4	Experimental	post	6
5	Experimental	pre	13
5	Experimental	post	9
6	Experimental	pre	17
6	Experimental	post	10
7	Experimental	pre	19
7	Experimental	post	12
8	Traditional	pre	1
8	Traditional	post	6
9	Traditional	pre	5
9	Traditional	post	7
10	Traditional	pre	7
10	Traditional	post	9
11	Traditional	pre	9
11	Traditional	post	9
12	Traditional	pre	13
12	Traditional	post	12
13	Traditional	pre	16
13	Traditional	post	12
14	Traditional	pre	19
14	Traditional	post	15

6.2 EXPLORATORY DATA ANALYSIS

6.2.1 Descriptive Statistics

6.2.1.1 Univariate

df_ptsd %>% 
  dplyr::select("Therapy" = therapy) %>% 
  apaSupp::tab_freq(caption = "Summary of Categorical Measures")

**Table 6.3**
Summary of Categorical Measures
		Statistic (N=14)
Therapy
Traditional		7 (50.0%)
Experimental		7 (50.0%)

df_ptsd %>% 
  dplyr::select("Post Test" = post_test,
                "Pre Test"  = pre_test,
                "Gain" = gain) %>% 
  apaSupp::tab_desc(caption = "Summary of Measures")

**Table 6.4**
Summary of Measures
	NA	M	SD	min	Q1	Mdn	Q3	max
Post Test	0	8.50	3.57	2.00	6.00	9.00	11.50	15.00
Pre Test	0	10.00	6.32	1.00	5.50	9.50	15.25	19.00
Gain	0	-1.50	3.59	-7.00	-4.00	-1.00	1.00	5.00
Note. N = 14. NA = not available or missing; Mdn = median; Q1 = 25th percentile; Q3 = 75th percentile.

6.2.1.2 Bivariate

df_ptsd %>% 
  dplyr::select(therapy,
                "Post Test" = post_test,
                "Pre Test"  = pre_test,
                "Gain" = gain) %>% 
  apaSupp::tab1(split = "therapy",
                type = list("Post Test" = "continuous",
                            "Gain"      = "continuous"),
                total_last = FALSE,
                caption = "Descriptive Summary by Therapy Group")

**Table 6.5**
Descriptive Summary by Therapy Group
	Total N = 14	Traditional n = 7 (50.0%)	Experimental n = 7 (50.0%)	p-value
Post Test	8.50 (3.57)	10.00 (3.16)	7.00 (3.51)	.119
Pre Test	10.00 (6.32)	10.00 (6.35)	10.00 (6.81)	1.000
Gain	-1.50 (3.59)	0.00 (3.32)	-3.00 (3.42)	.121
Note. Continuous variables are summarised with means (SD) and significant group differences assessed via independent t-tests. Categorical variables are summarised with counts (%) and significant group differences assessed via Chi-squared tests for independence.
* p < .05. p < .01. * p < .001.

Pearson’s correlation is only appropriate for TWO CONTINUOUS variables. The exception is when 1 variable is continuous and the other has exactly 2 levels. In this case, the binary variable needs to be converted to two numbers (numeric not factor) and the value is called the Point-Biserial Correlation ($r_{pb}$).

df_ptsd %>% 
  dplyr::mutate(therapy = as.numeric(therapy == "Experimental")) %>% 
  dplyr::select("Post Test" = post_test,
                "Pre Test"  = pre_test,
                "Gain" = gain,
                "Therapy" = therapy) %>% 
  apaSupp::tab_cor(caption = "Correlation Between Pairs of Measures",
                   general_note = "For pairs of variables with therapy, r = Point-Biserial Correlation, otherwise ") %>% 
  flextable::hline(i = 3) %>% 
  flextable::bg(i = 1:3, bg = "yellow") %>% 
  flextable::bg(i = c(4), bg = "orange")

**Table 6.6**
Correlation Between Pairs of Measures
Variable Pair		r	p
Post Test	Pre Test	.880	< .001***
Post Test	Gain	-.560	.037*
Post Test	Therapy	-.440	.119
Pre Test	Gain	-.880	< .001***
Pre Test	Therapy	.000	1.000
Gain	Therapy	-.430	.121
Note. N = 14. r = Pearson's Product-Moment correlation coefficient.For pairs of variables with therapy, r = Point-Biserial Correlation, otherwise
* p < .05. p < .01. * p < .001.

6.2.2 Visualizing Distributions

6.2.2.1 Univariate

df_ptsd %>% 
  dplyr::mutate(therapy = as.numeric(therapy == "Experimental")) %>% 
  dplyr::select(id,
                "Post Symptom Severity" = post_test,
                "Pre Symptom Severity"  = pre_test,
                "Change in Symptom Severity" = gain,
                "Therapy" = therapy) %>%  
  tidyr::pivot_longer(cols = -id) %>% 
  ggplot(aes(value)) +
  geom_histogram(binwidth = 1,
                 color = "black",
                 alpha = .25) +
  theme_bw() +
  facet_wrap(~ name,
             scale = "free_x") +
  labs(x = NULL,
       y = "Count")

Figure 6.1
Univariate Distibution of Measures

6.2.2.2 Bivariate

apaSupp::spicy_histos(df = df_ptsd,
                      var = pre_test,
                      split = therapy,
                      lab = "Symptom Severity Prior to Therapy") +
  scale_x_continuous(breaks = seq(from = -24, to = 24, by = 2))

Figure 6.2
Distribution of Symptom Severity Prior to Therapy

apaSupp::spicy_histos(df = df_ptsd,
                      var = post_test,
                      split = therapy,
                      lab = "Symptom Severity After Therapy") +
  scale_x_continuous(breaks = seq(from = -24, to = 24, by = 2))

Figure 6.3
Distribution of Symptom Severity After Therapy

apaSupp::spicy_histos(df = df_ptsd,
                      var = gain,
                      split = therapy,
                      lab = "Change in Symptom Severity") +
  scale_x_continuous(breaks = seq(from = -24, to = 24, by = 2))

Figure 6.4
Distribution of Change in Symptom Severity

df_ptsd_long %>% 
  ggplot(aes(x = time,
             y = test,
             group = id)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  labs(x = NULL,
       y = "Symptom Severity Score") +
  facet_grid(~ therapy) +
  stat_summary(aes(group = 1),
           geom = "point",
           fun = "mean",
           shape = 13,
           size = 5,
           color = "red")+
  stat_summary(aes(group = 1),
           geom = "line",
           fun = "mean",
           color = "red",
           linewidth = 1.5)

Figure 6.5
Progression in Symptom Severity, Stratified by Therapy

6.2.2.3 Multivariate

df_ptsd %>% 
  ggplot(aes(x = pre_test,
             y = post_test)) +
  geom_point(aes(color = therapy,
                 shape = therapy),
             size = 3) + 
  theme_bw() +
  geom_smooth(aes(color = therapy,
                  fill = therapy,
                  linetype = therapy),
              method = "lm",
              formula = y ~x,
              alpha = .2) +
  labs(x = "Observed Symptom Severity Prior to Therapy",
       y = "Observed Symptom Severity\nAfter Therapy",
       color    = "Therapy Type:",
       fill     = "Therapy Type:",
       shape    = "Therapy Type:",
       linetype = "Therapy Type:") +
  theme(legend.position.inside = TRUE,
        legend.position = c(0, 1),
        legend.justification = c(-.1, 1.1),
        legend.background = element_rect(color = "black"),
        legend.key.width = unit(2, "cm")) +
  scale_linetype_manual(values = c("solid", "longdash")) +
  scale_x_continuous(breaks = seq(from = 0, to = 24, by = 2)) +
  scale_y_continuous(breaks = seq(from = 0, to = 24, by = 2)) + 
  geom_abline(intercept = 0, slope = 1)

Figure 6.6
Association Between Symptom Severity Prior to and After Therapy, Stratified by Therapy Type

df_ptsd %>% 
  ggplot(aes(x = pre_test,
             y = gain)) +
  geom_point(aes(color = therapy,
                 shape = therapy),
             size = 3) + 
  theme_bw() +
  geom_smooth(aes(color = therapy,
                  fill = therapy,
                  linetype = therapy),
              method = "lm",
              formula = y ~x,
              alpha = .2) +
  labs(x = "Observed Symptom Severity Prior to Therapy",
       y = "Observed Change in Symptom Severity",
       color    = "Therapy Type:",
       fill     = "Therapy Type:",
       shape    = "Therapy Type:",
       linetype = "Therapy Type:") +
  theme(legend.position.inside = TRUE,
        legend.position = c(1, 1),
        legend.justification = c(1.1, 1.1),
        legend.background = element_rect(color = "black"),
        legend.key.width = unit(2, "cm")) +
  scale_linetype_manual(values = c("solid", "longdash")) +
  scale_x_continuous(breaks = seq(from = 0, to = 24, by = 2))+
  scale_y_continuous(breaks = seq(from = -24, to = 24, by = 2)) +
  geom_hline(yintercept = 0)

Figure 6.7
Association Between Test Gain in Score and Pre-Test, Stratified by Therapy

6.3 REGRESSION ANALYSIS

6.3.1 DV = Post Score, controlling for Pre Score

fit_ptsd_1 <- lm(post_test ~ pre_test + therapy,
                 data = df_ptsd)

apaSupp::tab_lm(fit_ptsd_1)

**Table 6.7**
Parameter Estimates for Linear Regression
	b	(SE)	p	$b^*$	$\eta^2$	$\eta^2_p$
(Intercept)	5.02	(0.39)	< .001***
pre_test	0.50	(0.03)	< .001***	0.88	.779	.963
therapy					.190	.863
Traditional	—	—
Experimental	-3.00	(0.36)	< .001***
R²	.970
Adjusted R²	.964
Note. N = 14. $b^*$ = standardize coefficient; $\eta^2$ = semi-partial correlation; $\eta^2_p$ = partial correlation; p = significance from Wald t-test for parameter estimate.
* p < .05. p < .01. * p < .001.

6.3.1.1 Equation

\[ \widehat{post} = 5.019 + 0.498(pre) -3.00(type=exp) \]

6.3.1.2 Role of Therapy

Test the significance for therapy: t(11) = -8.33, p <.001, VERY significant

After controlling for pre-test score, the experimental therapy is associated with 3 point lower post-test score compared to the traditional therapy, b = -3.00, SE = 0.36, p <.001.

Notes:

Covariate pre_test has b = -0.50, SE = 0.03, p <.001.

The entire model: SSresid = 4.998, R^2 = .970, F(2, 11) = 176.6, p < .001

6.3.2 DV = Gain Score, ignoring Pre-Score

fit_ptsd_2 <- lm(gain ~ therapy,
                 data = df_ptsd)

apaSupp::tab_lm(fit_ptsd_2)

**Table 6.8**
Parameter Estimates for Linear Regression
	b	(SE)	p	$\eta^2$	$\eta^2_p$
(Intercept)	0.00	(1.27)	1.000
therapy				.188	.188
Traditional	—	—
Experimental	-3.00	(1.80)	.121
R²	.188
Adjusted R²	.120
Note. N = 14. $b^*$ = standardize coefficient; $\eta^2$ = semi-partial correlation; $\eta^2_p$ = partial correlation; p = significance from Wald t-test for parameter estimate.
* p < .05. p < .01. * p < .001.

6.3.2.1 Equation

\[ \widehat{gain} = 0 -3.00(EXP) \]

6.3.2.2 Role of Therapy

Test the significance for therapy: t(12) = -1.67, p = .121, NOT significant

No evidence was found that the therapy type was associated with a differential change in symptom severity,b = -3.00, SE = 1.80, p =.121.

Note:

The entire model: SSresid = 136, R^2 = .188, F(1, 11) = 2.78, p =.121

6.3.2.3 Equivalent Analysis: independent groups t-test

t.test(gain ~ therapy,
       data = df_ptsd)


    Welch Two Sample t-test

data:  gain by therapy
t = 1.6672, df = 11.99, p-value = 0.1214
alternative hypothesis: true difference in means between group Traditional and group Experimental is not equal to 0
95 percent confidence interval:
 -0.9210863  6.9210863
sample estimates:
 mean in group Traditional mean in group Experimental 
                         0                         -3

6.3.3 DV = Gain Score, controlling for Pre-Score

fit_ptsd_3 <- lm(gain ~ pre_test + therapy,
                 data = df_ptsd)

apaSupp::tab_lm(fit_ptsd_3)

**Table 6.9**
Parameter Estimates for Linear Regression
	b	(SE)	p	$b^*$	$\eta^2$	$\eta^2_p$
(Intercept)	5.02	(0.39)	< .001***
pre_test	-0.50	(0.03)	< .001***	-0.88	.782	.963
therapy					.188	.863
Traditional	—	—
Experimental	-3.00	(0.36)	< .001***
R²	.970
Adjusted R²	.965
Note. N = 14. $b^*$ = standardize coefficient; $\eta^2$ = semi-partial correlation; $\eta^2_p$ = partial correlation; p = significance from Wald t-test for parameter estimate.
* p < .05. p < .01. * p < .001.

6.3.3.1 Equation

\[ \widehat{gain} = 5.019 - 0.50(pre) -3.00(EXP) \]

6.3.3.2 Role of Therapy

Test the significance for therapy type on symptom severity change: t(11) = -8.33, p <.001, VERY significant

After controlling for symptom severity prior to therapy, the experimental therapy was associated with 3 point lower change on symptom severity after therapy compared to traditional therapy, b = -3.00, SE = 0.36, p <.001.

Notes:

Covariate pre_test has b = -0.50, SE = 0.03, p <.001.

The entire model: SSresid = 4.998, R^2 = .970, F(2, 11) = 178.8, p < .001

6.3.4 Compare Models

apaSupp::tab_lm_fits(list(fit_ptsd_1, fit_ptsd_2, fit_ptsd_3))

**Table 6.10**
Comparison of Linear Model Performane Metrics
			$R^2$
Model	N	k	mult	adj	AIC	BIC	RMSE
Model 1	14	3	.970	.964	33.31	35.87	0.60
Model 2	14	2	.188	.120	77.56	79.48	3.12
Model 3	14	3	.970	.965	33.31	35.87	0.60
Note. k = number of parameters estimated in each model. Larger $R^2$ values indicated better performance. Smaller values indicated better performance for Akaike's Information Criteria (AIC), Bayesian information criteria (BIC), and Root Mean Squared Error (RMSE).

6.4 ALTERNATIVE ANALYSIS

6.4.1 2x2 mix ANOVA

fit_aov <- afex::aov_4(test ~ therapy + (time|id),
                       data = df_ptsd_long)

summary(fit_aov)


Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

              Sum Sq num Df Error SS den Df F value    Pr(>F)    
(Intercept)  2395.75      1      586     12 49.0597 1.426e-05 ***
therapy        15.75      1      586     12  0.3225    0.5806    
time           15.75      1       68     12  2.7794    0.1213    
therapy:time   15.75      1       68     12  2.7794    0.1213    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1