Why Response Profiles?

Goal: Analyze longitudinal data with minimal structure on:

• Mean response trajectory
• Covariance among repeated measures

What Are Response Profiles?

Response profiles are sequences of group-specific mean responses over time.

Key features:

• Applies to balanced longitudinal designs
• Same measurement occasions for all subjects
• Can accommodate missing data (balanced but incomplete designs)

Computing Group-Specific Means

Suppose groups are defined by treatment or observational factors.

At each time \(t_j\), compute the group-specific mean:

\[ \hat{\mu}_g(t_j) = \frac{1}{n_g} \sum_{i \in g} Y_{ij} \]

Plotting \(\hat{\mu}_g(t_j)\) over time yields the mean response profile.

Example: TLC Trial Data

# Published TLC group means from FLW Table 5.6; full long data loaded later (Section 6)
time <- c(0, 1, 4, 6)
succimer <- c(26.5, 13.5, 15.5, 20.8)
placebo  <- c(26.3, 24.7, 24.1, 23.6)

plot(time, succimer, type="b", col="blue", pch=19, ylim=c(10,30),
     xlab="Time (weeks)", ylab="Mean Blood Lead (µg/dL)",
     main="TLC group means (FLW Table 5.6)")
lines(time, placebo, type="b", col="red", pch=17)
legend("topright", legend=c("Succimer", "Placebo"),
       col=c("blue","red"), pch=c(19,17))

TLC Trial Observations

Discussion: Handling Baseline

Baseline measurements are special:

• Typically pre-randomization in trials (independent of treatment)
• May be range-restricted or differ systematically in observational studies

Key Question: How Should Baseline Be Handled?

Implications for hypothesis testing, efficiency, and power.

Three Main Hypotheses

Given \(n\) occasions (time points) across \(G\) groups, we test:

Testing Order Principle

If \(H_{01}\) is rejected, main effects (\(H_{02}\), \(H_{03}\)) are not directly interpretable due to the presence of interaction.

Takeaway: Test parallelism first; only then consider main effects.

Randomized vs. Observational Studies

Visual: Parallel Profiles (No Interaction)

Two groups with parallel mean trajectories (constant difference over time).

set.seed(1)
t <- 1:5
mu1 <- 10 + 2*t                # Group 1
mu2 <- 12 + 2*t                # Group 2 (constant +2 offset)

plot(t, mu1, type="b", pch=19, ylim=range(c(mu1, mu2)),
     xlab="Time", ylab="Mean Response", main="(a) Parallel Mean Profiles")
lines(t, mu2, type="b", pch=17)
abline(v=t, col="gray90", lty=3)
legend("topleft", legend=c("Group 1","Group 2"), pch=c(19,17), bty="n")

Visual: Flat Profiles (No Time Effect)

Assuming parallelism holds, flat profiles test:

\[ H_{02}: \mu_1 = \mu_2 = \cdots = \mu_n \]

t <- 1:5
mu1 <- rep(15, length(t))
mu2 <- rep(20, length(t))

plot(t, mu1, type="b", pch=19, ylim=range(c(mu1, mu2)),
     xlab="Time", ylab="Mean Response", main="(b) Flat Mean Profiles (No Time Effect)")
lines(t, mu2, type="b", pch=17)
abline(v=t, col="gray90", lty=3)
legend("topleft", legend=c("Group 1","Group 2"), pch=c(19,17), bty="n")

Visual: Same Level (No Group Effect)

Assuming parallelism holds, same level tests:

\[ H_{03}: \boldsymbol{\mu}^{(1)} = \boldsymbol{\mu}^{(2)} = \cdots = \boldsymbol{\mu}^{(G)} \]

t <- 1:5
mu  <- 10 + 1.5*t

# The two group profiles coincide by construction (no group effect);
# add a tiny vertical jitter so both series are visible.
plot(t, mu, type="b", pch=19, ylim=range(mu),
     xlab="Time", ylab="Mean Response", main="(c) Same Level Across Groups (No Group Effect)")
lines(t, mu + 0.1, type="b", pch=17, col="gray40")
abline(v=t, col="gray90", lty=3)
legend("topleft", legend=c("Group 1","Group 2 (jittered)"), pch=c(19,17),
       col=c("black","gray40"), bty="n")

Formal Statement: Two Groups

For Treatment (\(T\)) vs. Control (\(C\)), define:

\[ \Delta_j = \mu_j^{(T)} - \mu_j^{(C)} \]

Test parallelism:

\[ H_{01}: \Delta_1 = \Delta_2 = \cdots = \Delta_n \quad (n-1 \text{ df}) \]

Formal Statement: More Than Two Groups

With reference group \(G\), define:

\[ \Delta_j^{(g)} = \mu_j^{(g)} - \mu_j^{(G)} \]

Test parallelism:

\[ H_{01}: \Delta_1^{(g)} = \cdots = \Delta_n^{(g)} \quad \forall g=1,\dots,G-1 \]

with \((G-1)(n-1)\) degrees of freedom.

The General Linear Model (LM) for Response Profiles

Represent mean response profiles using the general (normal-errors) linear model, the LM:

\[ E(\mathbf{Y}_i \mid \mathbf{X}_i) = \boldsymbol{\mu}_i = \mathbf{X}_i \boldsymbol{\beta} \]

Notation (Ch. 5):

Example: Two Groups, Three Occasions

Mean response profiles (column vectors):

\[ \boldsymbol{\mu}^{(1)} = \begin{pmatrix} \mu_1^{(1)} \\ \mu_2^{(1)} \\ \mu_3^{(1)} \end{pmatrix}, \quad \boldsymbol{\mu}^{(2)} = \begin{pmatrix} \mu_1^{(2)} \\ \mu_2^{(2)} \\ \mu_3^{(2)} \end{pmatrix} \]

Parameter vector:

\[ \boldsymbol{\beta} = (\beta_1, \beta_2, \ldots, \beta_6)^\top \]

Design Matrix Example

For a subject in group 1 with all three occasions (cell-means coding):

\[ \mathbf{X}_i = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix} \]

Columns 1-3 are the group-1 cell means (occasions 1, 2, 3); columns 4-6 are the group-2 cell means. So \(\mathbf{X}_i \boldsymbol{\beta}\) picks out group 1 here. The map is:

\[ \beta_1 = \mu_1^{(1)},\ \beta_2 = \mu_2^{(1)},\ \beta_3 = \mu_3^{(1)},\quad \beta_4 = \mu_1^{(2)},\ \beta_5 = \mu_2^{(2)},\ \beta_6 = \mu_3^{(2)} \]

Expressing Hypotheses via Contrasts

Hypothesis of no group \(\times\) time interaction:

\[ H_{01}: \mathbf{L}\boldsymbol{\beta} = \mathbf{0} \]

Each row of \(\mathbf{L}\) encodes one \(\Delta_j - \Delta_{j+1} = 0\) comparison, where \(\Delta_j = \mu_j^{(2)} - \mu_j^{(1)}\). With \(\boldsymbol{\beta} = (\mu_1^{(1)},\mu_2^{(1)},\mu_3^{(1)},\mu_1^{(2)},\mu_2^{(2)},\mu_3^{(2)})^\top\), one row testing \(\Delta_1 = \Delta_2\) is:

\[ \mathbf{L}_{\text{row 1}} = \begin{pmatrix} -1 & 1 & 0 & 1 & -1 & 0 \end{pmatrix} \;\Rightarrow\; \mathbf{L}_{\text{row 1}}\boldsymbol{\beta} = (\mu_2^{(2)} - \mu_2^{(1)}) - (\mu_1^{(2)} - \mu_1^{(1)}) = \Delta_2 - \Delta_1 \]

Setting all \(n-1\) such rows to \(\mathbf{0}\) forces \(\Delta_1 = \cdots = \Delta_n\) (parallelism).

Testing methods: Wald or likelihood ratio tests

Visual: Observed Means

time <- 1:3
mu_group1 <- c(10, 14, 18)
mu_group2 <- c(11, 16, 21)  # slightly different slope

plot(time, mu_group1, type="b", pch=19, col="blue", ylim=c(10,22),
     xlab="Time", ylab="Mean Response", main="Observed Means")
lines(time, mu_group2, type="b", pch=17, col="red")
legend("topleft", legend=c("Group 1","Group 2"),
       col=c("blue","red"), pch=c(19,17), bty="n")

Visual: Fitted Full Model

Full model: each group \(\times\) time cell estimated freely.

plot(time, mu_group1, type="b", pch=19, col="blue", ylim=c(10,22),
     xlab="Time", ylab="Mean Response", main="Fitted Means (Full Model)")
lines(time, mu_group2, type="b", pch=17, col="red")
legend("topleft", legend=c("Group 1","Group 2"),
       col=c("blue","red"), pch=c(19,17), bty="n")

Visual: Fitted Under the Null

Null hypothesis: group difference \(\Delta\) is constant over time.

delta <- mean(mu_group2 - mu_group1)   # average difference
mu_group2_null <- mu_group1 + delta

plot(time, mu_group1, type="b", pch=19, col="blue", ylim=c(10,22),
     xlab="Time", ylab="Mean Response", main="Fitted Means (Null: Parallel Profiles)")
lines(time, mu_group2_null, type="b", pch=17, col="red", lty=2)
legend("topleft", legend=c("Group 1","Group 2 (Null Fit)"),
       col=c("blue","red"), pch=c(19,17), lty=c(1,2), bty="n")

Handling Missing Data

Advantage of the LM formulation: Flexible design matrix \(\mathbf{X}_i\).

If a subject is missing a time point, remove the corresponding row from \(\mathbf{X}_i\).

Still valid to fit model using available observations.

Reference Group Parameterization

Alternative coding: choose one group as reference.

With \(G\) groups, define group-indicator covariates \(W_{ig}\) (the symbol \(Z\) is reserved for random effects, Ch. 8):

\[ W_{ig} = \begin{cases} 1 & \text{if subject } i \text{ belongs to group } g \\ 0 & \text{otherwise} \end{cases} \]

Reference Group Interpretation

If we set \(\mathbf{X}_i = (1, W_{i1}, \ldots, W_{i,G-1})\), then:

\[ \mu^{(g)} = \begin{cases} \beta_1 + \beta_{g+1}, & g=1,\ldots,G-1 \\ \beta_1, & g=G \quad (\text{reference group}) \end{cases} \]

R Example: Reference Coding

# Example: 3 groups, using Group 3 as reference
group <- factor(c("G1","G2","G3","G1","G2","G3"))

# Default coding in R: last group is reference
X <- model.matrix(~ group)
X

Reference Coding Interpretation

• Intercept corresponds to Group 3 mean
• Columns groupG1 and groupG2 are deviations from Group 3
• This is the default parameterization in R

Testing Main Effects

If profiles are parallel (\(H_{01}\) not rejected):

Both tested using a reduced model (main effects only, no interaction).

LM Framework: Key Takeaways

• Encodes group \(\times\) time means in \(\boldsymbol{\beta}\)
• Tests hypotheses via linear contrasts (\(\mathbf{L}\boldsymbol{\beta} = \mathbf{0}\))
• Full model: exact means per cell
• Null model: imposes parallelism (constant group difference)

Four Strategies for Baseline Handling

Strategy 1: Retain Baseline, No Equality

Fit the LM with all times, no constraint at \(t=1\):

\[ E(\mathbf{Y}_i \mid \mathbf{X}_i) = \mathbf{X}_i \boldsymbol{\beta} \]

where \(\mu_1^{(1)}, \ldots, \mu_1^{(G)}\) are free parameters.

Strategy 1: Pros and Cons

Strategy 2: Constrain Equal at Baseline

Appropriate in randomized trials (baseline pre-randomization).

Add constraint:

\[ \mu_1^{(1)} = \mu_1^{(2)} = \cdots = \mu_1^{(G)} \]

Strategy 2: Pros and Cons

Strategy 3: Change-from-Baseline

Define subject-level differences:

\[ \Delta_{ij} = Y_{ij} - Y_{i1}, \quad j = 2, \ldots, n \]

Or use a univariate summary:

\[ Y_i^* = \frac{1}{n-1}\sum_{j=2}^n (Y_{ij} - Y_{i1}) \]

Strategy 3: Pros and Cons

Strategy 4: ANCOVA

Univariate ANCOVA on average post-baseline response:

\[ Y_i^* = \alpha + \beta \cdot \text{trt}_i + \gamma \cdot Y_{i1} + \varepsilon_i \]

Or time-specific repeated-measures ANCOVA:

\[ Y_{ij} = \alpha_j + \beta_j \cdot \text{trt}_i + \gamma_j \cdot Y_{i1} + \varepsilon_{ij}, \quad j \ge 2 \]

Strategy 4: Pros and Cons

Lord’s Paradox (beyond FLW Ch. 5)

            change score says:            ANCOVA says:
   group A  ---- no change (flat) ----     conditional on baseline,
   group B  ---- no change (flat) ----     A and B differ in outcome
   (groups start at different baselines)   (regression line tilts the comparison)
   => "no group effect"                    => "group effect" (opposite conclusion)

The two methods point in different directions precisely because the groups are not comparable at baseline (no randomization).

Relative Efficiency: Change vs. ANCOVA

Under compound symmetry with variance \(\sigma^2\) and correlation \(\rho\), with one baseline and one post-baseline occasion (\(n = 2\)):

\[ \text{RE} = \frac{\text{Var}(\hat\beta_{\text{Change}})}{\text{Var}(\hat\beta_{\text{ANCOVA}})} = \frac{2}{1 + \rho} \]

This is the single-follow-up special case. The general \(n\)-occasion result (FLW eq. (5.3), Section 5.6) gives the efficiency of the change-score summary relative to ANCOVA as

\[ \frac{1}{n}\{1 + (n-1)\rho\} < 1, \]

so ANCOVA is \(\dfrac{n}{1 + (n-1)\rho}\) times as efficient (at \(n = 2\) this is \(2/(1+\rho)\)). ANCOVA’s advantage grows BOTH as \(\rho\) decreases AND as the number of occasions \(n\) increases (an \(n\)-fold gain as \(\rho \to 0\)).

Why ANCOVA wins: a change score writes \(Y_f - 1 \cdot Y_1\), forcing the baseline coefficient to exactly 1. ANCOVA fits \(Y_f - \gamma \cdot Y_1\) and estimates the optimal \(\gamma\) (here \(\gamma \approx \rho\)), the regression-to-the-mean adjustment. Pinning \(\gamma = 1\) is generally not optimal, so change scores carry extra variance.

Relative Efficiency: Implications

The table below is the single-follow-up (\(n = 2\)) special case \(\text{RE} = 2/(1+\rho)\):

RE \(\geq 1\) always, so ANCOVA is always at least as efficient as change scores. With more occasions the gap is larger still: by FLW eq. (5.3) ANCOVA is \(n/\{1+(n-1)\rho\}\) times as efficient, up to \(n\)-fold as \(\rho \to 0\).

Decision Guide: When to Use Which?

Estimand Matters

When to Use Alternatives

Rule of Thumb Summary

Mini Demo: Change vs. ANCOVA

set.seed(7)
n <- 120
group <- factor(rep(c("Control","Treatment"), each=n/2))
Y1 <- rnorm(n, 20, 4)                 # baseline
# True model: follow-up depends on baseline + treatment effect
Yf <- 5 + 0.6*Y1 + 2*(group=="Treatment") + rnorm(n, 0, 3)

# Change score and ANCOVA
change <- Yf - Y1
fit_change <- lm(change ~ group)
fit_ancova <- lm(Yf ~ group + Y1)

summary(fit_change)$coef

Takeaway: Comparing the Two Treatment-Effect SEs

se_change <- summary(fit_change)$coef["groupTreatment", "Std. Error"]
se_ancova <- summary(fit_ancova)$coef["groupTreatment", "Std. Error"]
rho_hat   <- cor(Y1, Yf)

c(SE_change = se_change, SE_ancova = se_ancova,
  ratio_SE2 = (se_change^2) / (se_ancova^2),   # empirical variance ratio
  RE_formula = 2 / (1 + rho_hat))              # 2/(1+rho)

The change-score SE is the larger of the two. This demo has a single follow-up, so the empirical variance ratio tracks the \(n = 2\) formula \(\text{RE} = 2/(1+\rho)\) (the special case of FLW eq. (5.3)): ANCOVA wins because it estimates the optimal baseline coefficient instead of forcing it to 1. The gap is largest at low \(\rho\) and vanishes as \(\rho \to 1\) (decreasing in \(\rho\)); with more occasions the gap grows toward \(n\)-fold.

Baseline Value Strategy Summary

Bottom line: Let design + scientific question dictate the strategy.

Strengths of Response Profile Analysis

• Balanced designs with common measurement times
• Arbitrary mean structure + unstructured \(\Sigma\)
• Avoids covariance misspecification by leaving \(\Sigma\) unstructured (no working-correlation assumption)
• Works with balanced-incomplete data
• Discrete covariates via the LM; contrasts/AUC

Comparison with MANOVA

Weaknesses of Response Profile Analysis

• Requires balanced timing; poor for mistimed visits
• Ignores time ordering; omnibus tests need follow-up
• Lower power for specific trends (prefer 1 df tests)
• Parameter growth: \(Gn + \frac{n(n+1)}{2}\)

Parameter Growth Example (\(G=2\))

param_counts <- function(G, n) G*n + n*(n+1)/2
data.frame(n=c(3,10), total_params=param_counts(2, c(3,10)))

Computing Goals

• Replicate PROC MIXED response-profiles analysis in R
• Model: \(Y_{ij} \sim \text{group} \times \text{time}\)
• Unstructured within-subject covariance
• Use nlme::gls() with corSymm + varIdent

Read and Prepare Data

library(dplyr); library(tidyr); library(nlme); library(ggplot2)
if (requireNamespace("emmeans", quietly = TRUE)) library(emmeans)

# Path to data file
path <- "../../data/tlc.csv"

# Read and rename columns
tlc_wide <- read.csv(path, check.names = FALSE, stringsAsFactors = FALSE)
names(tlc_wide) <- trimws(names(tlc_wide))

tlc_wide <- tlc_wide %>%
  rename(
    id  = ID,
    trt = `Treatment Group`,
    w0  = `Lead Level Week 0`,
    w1  = `Lead Level Week 1`,
    w4  = `Lead Level Week 4`,
    w6  = `Lead Level Week 6`
  ) %>%
  mutate(trt = factor(trt, levels = c("P","A"), labels = c("Placebo","Succimer")))

head(tlc_wide)

Reshape: Wide to Long

tlc_long <- tlc_wide %>%
  pivot_longer(cols = c(w0, w1, w4, w6),
               names_to = "wk", values_to = "y") %>%
  mutate(
    time = recode(wk, w0 = "0", w1 = "1", w4 = "4", w6 = "6"),
    time = factor(time, levels = c("0","1","4","6")),
    group = trt
  ) %>%
  arrange(id, time) %>%
  select(id, group, time, y)

head(tlc_long)

How gls() Builds an Unstructured Covariance

Fit Response Profiles (UN Covariance)

# Unstructured covariance = free correlations (corSymm) + per-occasion variances (varIdent)
fit_un_REML <- gls(
  y ~ group * time,
  data = tlc_long,
  correlation = corSymm(form = ~ as.numeric(time) | id),  # free correlations; sets within-id order
  weights     = varIdent(form = ~ 1 | time),               # one variance per occasion
  method = "REML",
  na.action = na.omit
)

summary(fit_un_REML)

Test Group x Time (Parallelism)

# Use ML for fixed-effect LRTs
fit_full_ML   <- update(fit_un_REML, method = "ML")
fit_no_int_ML <- update(fit_full_ML, . ~ group + time)

anova(fit_no_int_ML, fit_full_ML)  # df = (G-1)*(n-1)

Interpreting the Parallelism Test

In plain language: The likelihood ratio test compares:

• Full model: Separate mean for each group-by-time combination (allows non-parallel profiles)
• Reduced model: Main effects only (forces parallel profiles)

A significant p-value (e.g., p < 0.05) indicates that profiles are not parallel - the treatment effect differs across time points. In the TLC trial, this would mean Succimer’s effect is not constant over time (as suggested by the sharp drop then rebound).

Estimated Response Profiles

if (requireNamespace("emmeans", quietly = TRUE)) {
  emm <- emmeans::emmeans(fit_un_REML, ~ group * time)
  emm_df <- as.data.frame(emm)

  print(ggplot(emm_df, aes(x = time, y = emmean, group = group, shape = group)) +
    geom_line() + geom_point(size = 2) +
    labs(y = "Estimated mean blood lead (µg/dL)",
         title = "TLC Trial: Estimated Response Profiles (UN covariance)") +
    theme_minimal())
} else {
  # Alternative: compute fitted values from the model
  newdata <- expand.grid(group = levels(tlc_long$group),
                         time = levels(tlc_long$time))
  newdata$emmean <- predict(fit_un_REML, newdata = newdata)

  print(ggplot(newdata, aes(x = time, y = emmean, group = group, shape = group)) +
    geom_line() + geom_point(size = 2) +
    labs(y = "Estimated mean blood lead (µg/dL)",
         title = "TLC Trial: Estimated Response Profiles (UN covariance)") +
    theme_minimal())
}

Reading Off the Covariance: the \(\Sigma_i = \mathbf{S}\,\mathbf{C}\,\mathbf{S}\) Identity

gls() stores the covariance in two pieces, so we reassemble it from:

\[ \Sigma_i = \mathbf{S}\,\mathbf{C}\,\mathbf{S} \]

• \(\mathbf{S}\) = diagonal of occasion-specific standard deviations (from varIdent).
• \(\mathbf{C}\) = the correlation matrix among occasions (from corSymm).

(We use \(\mathbf{S}\)/\(\mathbf{C}\) for the SD-correlation split here so as not to collide with the reserved \(G\) = random-effects covariance and \(R_i\) = residual covariance of Ch. 8.)

So a per-occasion SD on the diagonal times the free correlations rebuilds the full unstructured \(\Sigma_i\).

Extract UN Covariance Matrix (advanced / reference)

# --- S: occasion-specific SDs (from varIdent) ---
lev <- levels(tlc_long$time)
# varIdent reports a multiplier for each NON-reference occasion; the
# reference occasion (here "0") is fixed at 1 and is NOT in this vector.
vmult <- coef(fit_un_REML$modelStruct$varStruct, unconstrained = FALSE)
sd_mult <- setNames(rep(1, length(lev)), lev)  # start all at the reference value 1
sd_mult[names(vmult)] <- vmult                 # fill in the estimated multipliers by NAME
SDdiag <- diag(fit_un_REML$sigma * sd_mult)    # S = base sigma * per-occasion multiplier

# --- C: correlation among occasions (from corSymm) ---
Rhat <- corMatrix(fit_un_REML$modelStruct$corStruct)[[1]]

# --- Sigma = S C S ---
Sigma_hat <- SDdiag %*% Rhat %*% SDdiag
dimnames(Sigma_hat) <- list(lev, lev)
round(Sigma_hat, 2)

Replicate Figure 5.1: Raw Group Means

means <- tlc_long %>%
  group_by(group, time) %>%
  summarise(mean_y = mean(y, na.rm = TRUE), .groups = "drop") %>%
  mutate(week = as.numeric(as.character(time)))

ggplot(means, aes(x = week, y = mean_y, group = group, shape = group)) +
  geom_line() + geom_point(size = 2) +
  scale_x_continuous(breaks = c(0,1,4,6)) +
  labs(x = "Time (weeks)", y = "Mean blood lead (µg/dL)",
       title = "TLC Trial: Mean Blood Lead by Group and Time") +
  theme_minimal()

Baseline-Adjusted ANCOVA (Section 5.6)

Subject-level model:

\[ Y_i^* = \frac{w1_i + w4_i + w6_i}{3} \]

\[ Y_i^* = \alpha + \beta \cdot \text{group}_i + \gamma \cdot w0_i + \varepsilon_i \]

ANCOVA Implementation

dat_subj <- tlc_wide %>%
  transmute(
    id,
    group = trt,
    w0, w1, w4, w6,
    Ystar = rowMeans(cbind(w1, w4, w6), na.rm = TRUE)
  )

fit_ancova <- lm(Ystar ~ group + w0, data = dat_subj)
summary(fit_ancova)$coef

ANCOVA Output: Interpretation

In plain language: The coefficient output shows:

• Intercept: Expected average post-baseline lead level when baseline = 0 and group = reference (Placebo)
• groupSuccimer: The treatment effect - the difference in average post-baseline lead level for Succimer vs. Placebo, adjusted for baseline
• w0: The slope relating baseline to follow-up - a positive value means higher baseline is associated with higher follow-up

A negative, significant coefficient for groupSuccimer indicates that Succimer reduces blood lead levels compared to Placebo, after accounting for baseline differences.

Same Analysis in SAS: Complete PROC MIXED Program

The R gls() fit above corresponds to this standalone SAS program (static, not run here):

/* TLC response-profiles analysis: unstructured covariance */
proc mixed data=tlc_long method=reml;
    class id group time;
    model y = group time group*time / s chisq;
    repeated time / type=un subject=id r rcorr;
run;

• class declares the categorical factors (set id, group, time as classification variables).
• model y = group time group*time is the saturated cell-means mean model.
• repeated time / type=un subject=id is the unstructured within-subject covariance (the corSymm + varIdent pair in R).
• The group*time row of the “Type 3 Tests of Fixed Effects” table is the parallelism test (\((G-1)(n-1)\) df, here \(1\times 3 = 3\) for 2 groups and 4 weeks). The chisq option (REML default, no ddfm=) reports this as a Wald chi-square (CHISQ) parallelism test, matching FLW Table 5.11. It is not a likelihood-ratio test. (No hand-written contrast is needed; the Type 3 test already carries the correct df.)
• The R workflow above instead uses an ML likelihood-ratio test (anova() on two method="ML" fits). For an exact SAS analogue of that LRT, refit the full and no-interaction models with method=ml and compare \(-2\log L\); the Wald chi-square and LR tests target the same hypothesis but are different statistics.

SAS to R Cheat Sheet

Common Mistakes to Avoid

• Interpreting main effects when interaction is significant: If parallelism (\(H_{01}\)) is rejected, main effects for group and time are misleading - the effect depends on the level of the other factor
• Forgetting to set factor levels correctly: R’s default is to use alphabetical ordering; set the reference level explicitly (e.g., Placebo before Treatment)
• Using ANCOVA in observational studies without considering Lord’s paradox (defined earlier): with unequal baselines, change-score and ANCOVA conclusions can point in opposite directions
• Ignoring the relative efficiency of baseline strategies: Change scores can be substantially less efficient than ANCOVA when baseline-outcome correlation is low
• Applying response profiles to irregular visit times: This approach assumes common measurement occasions; use parametric curves (Ch. 6) for irregular timing

Summary

• Response profiles provide a flexible framework for balanced longitudinal data
• Test parallelism (\(H_{01}\)) before interpreting main effects
• Choose baseline handling strategy based on study design and estimand
• The general linear model (LM) with unstructured covariance maximizes flexibility
• R’s nlme::gls() replicates SAS PROC MIXED functionality