---
title: "Modern Analysis of the Solomon Four-Group Design"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Modern Analysis of the Solomon Four-Group Design}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## A unified approach

The Solomon four-group design is often introduced through a sequence
of separate analyses.

A modern alternative is to represent the four-group design in a single
model and estimate the scientifically relevant questions as explicit
contrasts.

For a continuous posttest outcome, the core model in `solomonR` is

\[
Y =
\beta_0 +
\beta_T T +
\beta_P P +
\beta_{TP}(T \times P) +
\beta_X X_{\mathrm{obs}} +
\varepsilon,
\]

where:

- \(T\) indicates treatment assignment;
- \(P\) indicates whether a participant was pretested;
- \(T \times P\) represents pretest sensitization; and
- \(X_{\mathrm{obs}}\) incorporates the pretest information available
  in the pretested groups.

`solomonR` then estimates four predefined Solomon contrasts from this
model.

## Structural missingness of the pretest

One unusual feature of the Solomon design deserves special attention.

Participants in Groups 3 and 4 were deliberately **not pretested**.
Their absent pretest values are therefore structurally missing by
design.

Simply placing the raw pretest variable into an ordinary regression
formula can cause complete-case deletion of Groups 3 and 4, destroying
the four-group design.

`fit_solomon_glm()` avoids this problem internally. For the optional
pretest covariate, the model uses the observed pretest score among
pretested participants and a design-safe value among participants for
whom the pretest was never administered.

The pretest indicator remains in the model, so this coding does not
pretend that Groups 3 and 4 actually had baseline scores of zero.

## Example data

```{r example-data}
library(solomonR)

data(solomon_example)

demo_preview <- head(
  solomon_example,
  6
)

demo_preview$y_post <- round(
  demo_preview$y_post,
  2
)

demo_preview$y_pre <- round(
  demo_preview$y_pre,
  2
)

knitr::kable(
  demo_preview,
  align = "rrrr"
)
```

## Fit the model

The unified model with HC3 heteroskedasticity-consistent covariance, the
default, is a sound starting point for continuous outcomes.

```{r}
fit <- fit_solomon_glm(
  y_post, treat, pretested, y_pre,
  robust = "HC3",
  data = solomon_example
)

fit
```

The output contains model coefficients followed by the four principal
Solomon contrasts.

The contrasts can also be accessed directly:

```{r glm-effects}
effects_display <- fit$effects[
  ,
  c(
    "contrast",
    "estimate",
    "std.error",
    "statistic",
    "p.value",
    "r2"
  )
]

effects_display$estimate <- sprintf(
  "%.2f",
  effects_display$estimate
)

effects_display$std.error <- sprintf(
  "%.2f",
  effects_display$std.error
)

effects_display$statistic <- sprintf(
  "%.2f",
  effects_display$statistic
)

effects_display$p.value <- ifelse(
  effects_display$p.value < .001,
  "< .001",
  sub(
    "^0",
    "",
    sprintf(
      "%.3f",
      effects_display$p.value
    )
  )
)

effects_display$r2 <- sub(
  "^0",
  "",
  sprintf(
    "%.3f",
    effects_display$r2
  )
)

names(effects_display) <- c(
  "Contrast",
  "Estimate",
  "SE",
  "t",
  "p",
  "Wald R2"
)

knitr::kable(
  effects_display,
  align = c(
    "l",
    "r",
    "r",
    "r",
    "r",
    "r"
  )
)
```

## The four Solomon estimands

### Average treatment effect

The reported ATE is the treatment effect averaged equally over the two
pretest conditions:

\[
ATE =
\beta_T +
\frac{1}{2}\beta_{TP}.
\]

This distinction is important.

With the unpretested condition used as the reference group,
\(\beta_T\) by itself is the treatment effect among **unpretested**
participants. It is not the equal-weighted Solomon ATE.

### Pretest x Treatment

The sensitization contrast is

\[
\beta_{TP}.
\]

It asks whether the treatment effect differs according to whether the
pretest was administered.

A statistically nonsignificant interaction should not be interpreted
as proof that sensitization is absent. The equivalence test described
below addresses that stronger question.

### Treatment effect among pretested participants

The treatment effect among participants who received the pretest is

\[
\beta_T + \beta_{TP}.
\]

### Treatment effect among unpretested participants

The treatment effect among participants who did not receive the
pretest is

\[
\beta_T.
\]

Together, these contrasts describe both the treatment effect and its
possible modification by pretesting.

## Why HC3?

The `robust` argument controls covariance estimation.

```{r eval=FALSE}
fit_solomon_glm(
  y_post,
  treat,
  pretested,
  y_pre,
  robust = "none"
)

fit_solomon_glm(
  y_post,
  treat,
  pretested,
  y_pre,
  robust = "HC3"
)
```

HC3, the default, provides heteroskedasticity-consistent standard errors
(MacKinnon & White, 1985). Long and Ervin (2000) recommend HC3 when the
sample is below about 250, and Hayes and Cai (2007) recommend
heteroskedasticity-consistent standard errors as routine practice. The
Solomon design adds its own reason: adjusting for the pretest reduces
residual variance only in the pretested groups, so the unified model is
heteroskedastic whenever the pretest predicts the posttest.

For Gaussian models, tests and confidence intervals use the t distribution
with residual degrees of freedom, with or without robust covariance;
binomial and Poisson models use the normal distribution. The `df`,
`conf.low`, and `conf.high` columns of `fit$effects` record the reference
distribution and interval for each contrast, and `conf_level` sets the
confidence level.

A CR2 covariance option is also available for clustered data:

```{r eval=FALSE}
fit_solomon_glm(
  y_post,
  treat,
  pretested,
  y_pre,
  robust = "CR2",
  cluster = cluster_id
)
```

CR2 combines the bias-reduced cluster-robust covariance estimator (Bell &
McCaffrey, 2002) with Satterthwaite degrees of freedom for each
coefficient and contrast (Pustejovsky & Tipton, 2018). The degrees of
freedom are reported in the `df` column and can be small when there are
few clusters.

## Is sensitization negligible?

A nonsignificant Pretest x Treatment contrast does not show that
sensitization is absent. An equivalence test asks a different question:
can effects at least as large as the smallest effect size of interest
(SESOI) be rejected? `equivalence_solomon()` performs the two one-sided
tests (TOST) procedure (Lakens, 2017).

The bounds must be set before the data are examined. Suppose that, before
collecting data, the researchers noted that prior studies with this
measure reported a posttest standard deviation of about 10 points, and
decided that sensitization smaller than half a standard deviation
(5 points) would be negligible.

```{r glm-equivalence}
equivalence_solomon(fit, bounds = 5)
```

With 30 participants per group, the 90% confidence interval still
extends beyond the bounds, so the result is inconclusive: these data can
neither detect sensitization nor rule out sensitization of that size.
Equivalence tests need adequate power of their own (Lakens, 2017), so
equivalence bounds belong in study planning.

## Wald-based partial R-squared

For Gaussian models, `fit_solomon_glm()` reports a Wald-based partial
\(R^2\) for each one-degree-of-freedom contrast.

With conventional Gaussian OLS covariance,

\[
R^2_{\mathrm{partial}} =
\frac{t^2}{t^2 + df_{\mathrm{residual}}}.
\]

When a robust covariance estimator such as HC3 is used, `solomonR`
applies the same transformation to the robust Wald statistic. In that
case the result should be interpreted as a **descriptive Wald-based
approximation**, not an exact decomposition of model variance.

For conventional Gaussian fits (`robust = "none"`), `r2_lo` and `r2_hi`
give a confidence interval from the noncentral F distribution (Steiger,
2004). No interval is reported with robust covariance, because that pivot
no longer applies.

```{r glm-r2}
r2_display <- fit$effects[
  ,
  c(
    "contrast",
    "r2"
  )
]

r2_display$r2 <- sub(
  "^0",
  "",
  sprintf(
    "%.3f",
    r2_display$r2
  )
)

names(r2_display) <- c(
  "Contrast",
  "Wald R2"
)

knitr::kable(
  r2_display,
  align = c("l", "r")
)
```

## Diagnostics

The package provides design-relevant diagnostics through
`check_solomon_assumptions()`.

```{r}
checks <- with(
  solomon_example,
  check_solomon_assumptions(
    y_post,
    treat,
    pretested,
    y_pre
  )
)

checks
```

These include variance checks, cell-level distribution summaries, and
assessment of ANCOVA slope homogeneity.

Diagnostics should inform analysis choices rather than be treated as a
mechanical series of pass/fail gates.

## Randomization-based inference

If treatment assignment was genuinely randomized within the Solomon
pretesting conditions, treatment-label permutation provides another
inferential approach.

```{r}
perm <- perm_solomon(
  fit,
  contrast = "ATE (avg over pretest)",
  reps = 1000,
  seed = 42,
  return_dist = TRUE
)

perm
```

`perm_solomon()` shuffles treatment labels **within pretest strata**.
This preserves the number of treated and control participants within
the pretested and unpretested portions of the design.

The test uses an HC3-studentized statistic.

Its Monte Carlo p-value uses a finite-simulation correction:

\[
p =
\frac{B_{\mathrm{extreme}} + 1}
     {B_{\mathrm{valid}} + 1}.
\]

This prevents an estimated permutation p-value of exactly zero.

The null distribution can be visualized with:

```{r perm-plot, fig.width = 7, fig.height = 4.5, fig.alt = "Histogram of the permutation distribution of the test statistic, with the observed value marked."}
plot_perm(perm)
```

Randomization inference is justified by the assignment mechanism, not
merely by a small sample size. If treatment was randomized at the
cluster level, permutation must likewise occur at the cluster level. For
fits with a `cluster` variable, `perm_solomon()` permutes whole clusters:

```{r eval=FALSE}
fit_cl <- fit_solomon_glm(
  y_post, treat, pretested, y_pre,
  robust = "CR2",
  cluster = cluster_id
)

perm_solomon(fit_cl, contrast = "Treatment | unpretested")
```

Two assignment mechanisms are supported: whole clusters assigned to the
four Solomon conditions, and treatment assigned to clusters with
pretesting assigned to participants within them. The statistic is built
from covariate-adjusted cluster-level summaries, each cluster counts once,
and small designs are enumerated exactly (Gail et al., 1996; Hayes &
Moulton, 2017). Designs that assign treatment to participants within
clusters are refused. The help page, `?perm_solomon`, gives the details
and the simulation evidence.

## A likelihood-based alternative

`fit_solomon_ml()` provides a full-information maximum-likelihood
analysis inspired by van Engelenburg's treatment of the Solomon
four-group design.

```{r}
ml <- with(
  solomon_example,
  fit_solomon_ml(
    y_post,
    treat,
    pretested,
    y_pre,
    inference = "satterthwaite"
  )
)

ml
```

A useful feature of the likelihood formulation is that the structurally
absent pretests in Groups 3 and 4 are treated as part of the design
rather than as ordinary missing baseline observations.

The ML and GLM approaches target closely related Solomon effects but use
different assumptions for uncertainty estimation.

By default, `fit_solomon_ml()` uses van Engelenburg's large-sample Wald
inference. In the package's simulation validation those intervals were
too narrow with small groups, so the example above uses the small-sample
option, `inference = "satterthwaite"`, which applies Welch-Satterthwaite
degrees of freedom (Satterthwaite, 1946; Welch, 1947). The function warns
when groups are small and `inference` has not been chosen.

## Historical analysis

For teaching or reproducing the classical Tests A-I workflow, use:

```{r}
classic <- with(
  solomon_example,
  fit_solomon_classic(
    y_post,
    treat,
    pretested,
    y_pre
  )
)

classic
```

See:

```{r eval=FALSE}
vignette("classic-solomon", package = "solomonR")
```

for a detailed explanation of the historical sequence and its
limitations.

## SEM extensions

`solomonR` also includes observed- and latent-variable SEM approaches.

An observed four-group mean model can be fit with:

```{r eval=FALSE}
fit_solomon_sem(
  y_post,
  treat,
  pretested
)
```

The unrestricted four-group observed mean model is saturated, so global
fit indices are not diagnostic for that model.

For multi-item outcomes, `fit_solomon_sem_latent()` provides latent
Solomon contrasts. Latent mean comparisons require scalar measurement
invariance across the four Solomon groups.

These SEM tools extend the same basic estimands used throughout the
package:

- ATE;
- pretest sensitization;
- treatment effect among pretested participants; and
- treatment effect among unpretested participants.

## Comparing analyses

`compare_solomon_methods()` fits several analyses to the same data and
lines up their estimates of each Solomon contrast.

```{r glm-compare}
with(
  solomon_example,
  compare_solomon_methods(
    y_post,
    treat,
    pretested,
    y_pre,
    methods = c("glm", "ml", "classic")
  )
)
```

All rows target the same population contrasts. Several estimates are
identical, such as the unified GLM, maximum likelihood, and the classic
ANCOVA among pretested participants, yet their standard errors and
intervals differ because the methods make different variance
assumptions and use different reference distributions. Pretest
adjustment changes precision rather than the target (Lin, 2013).

## Choosing an analysis

For many randomized Solomon studies with continuous outcomes, a useful
workflow is:

1. check the design with `validate_solomon()` and
   `check_solomon_missing()`, and inspect the data;
2. fit the unified GLM with prespecified contrasts;
3. examine the pretest-by-treatment contrast directly;
4. inspect diagnostics;
5. use randomization inference when justified by the treatment
   assignment mechanism; and
6. use ML or SEM when the scientific question or measurement structure
   warrants those approaches.

The historical decision tree remains valuable for teaching, but a
modern analysis does not need to choose its inferential procedure based
on whether an earlier significance test crossed \(p=.05\).

## Where solomonR is heading

`ROADMAP.md` in the package repository lists the work that remains before
version 1.0.0: the CRAN pre-submission checks, the last sources of the
literature audit, and the release itself.

All works cited in solomonR are listed, with notes on how the package uses
them, on the [References](https://juhalt.github.io/solomonR/articles/references.html) page.

## References
Bell, R. M., & McCaffrey, D. F. (2002). Bias reduction in standard errors
for linear regression with multi-stage samples. *Survey Methodology, 28*(2),
169–181.

Gail, M. H., Mark, S. D., Carroll, R. J., Green, S. B., & Pee, D. (1996). On
design considerations and randomization-based inference for community
intervention trials. *Statistics in Medicine, 15*(11), 1069–1092.
<https://doi.org/10.1002/(SICI)1097-0258(19960615)15:11%3C1069::AID-SIM220%3E3.0.CO;2-Q>

Hayes, A. F., & Cai, L. (2007). Using heteroskedasticity-consistent standard
error estimators in OLS regression: An introduction and software
implementation. *Behavior Research Methods, 39*(4), 709–722.
https://doi.org/10.3758/BF03192961

Hayes, R. J., & Moulton, L. H. (2017). *Cluster randomised trials* (2nd ed.).
Chapman and Hall/CRC. https://doi.org/10.4324/9781315370286

Lakens, D. (2017). Equivalence tests: A practical primer for t tests,
correlations, and meta-analyses. *Social Psychological and Personality
Science, 8*(4), 355–362. https://doi.org/10.1177/1948550617697177

Lin, W. (2013). Agnostic notes on regression adjustments to experimental
data: Reexamining Freedman's critique. *The Annals of Applied Statistics,
7*(1), 295–318. https://doi.org/10.1214/12-AOAS583

Long, J. S., & Ervin, L. H. (2000). Using heteroscedasticity consistent
standard errors in the linear regression model. *The American Statistician,
54*(3), 217–224. https://doi.org/10.1080/00031305.2000.10474549

MacKinnon, J. G., & White, H. (1985). Some heteroskedasticity-consistent
covariance matrix estimators with improved finite sample properties.
*Journal of Econometrics, 29*(3), 305–325.
https://doi.org/10.1016/0304-4076(85)90158-7

Pustejovsky, J. E., & Tipton, E. (2018). Small-sample methods for
cluster-robust variance estimation and hypothesis testing in fixed effects
models. *Journal of Business & Economic Statistics, 36*(4), 672–683.
https://doi.org/10.1080/07350015.2016.1247004

Satterthwaite, F. E. (1946). An approximate distribution of estimates of
variance components. *Biometrics Bulletin, 2*(6), 110–114.
https://doi.org/10.2307/3002019

Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals
and tests of close fit in the analysis of variance and contrast analysis.
*Psychological Methods, 9*(2), 164–182.
https://doi.org/10.1037/1082-989X.9.2.164

van Engelenburg, G. (1999). *Statistical analysis for the Solomon four-group
design* (Research Report 99-06). University of Twente. ERIC.
https://eric.ed.gov/?id=ED435692

Welch, B. L. (1947). The generalization of "Student's" problem when several
different population variances are involved. *Biometrika, 34*(1–2), 28–35.
https://doi.org/10.1093/biomet/34.1-2.28