7  Modeling Longitudinal Responses using Generalized Least Squares

Some good general references on longitudinal data analysis are Davis (2002), Pinheiro & Bates (2000), Diggle et al. (2002), Venables & Ripley (2003), Hand & Crowder (1996), Verbeke & Molenberghs (2000), Lindsey (1997)

7.1 Notation

• \(N\) subjects
• Subject \(i\) (\(i=1,2,\ldots,N\)) has \(n_{i}\) responses measured at times \(t_{i1}, t_{i2}, \ldots, t_{in_{i}}\)
• Response at time \(t\) for subject \(i\): \(Y_{it}\)
• Subject \(i\) has baseline covariates \(X_{i}\)
• Generally the response measured at time \(t_{i1}=0\) is a covariate in \(X_{i}\) instead of being the first measured response \(Y_{i0}\)
• Time trend in response is modeled with \(k\) parameters so that the time “main effect” has \(k\) d.f.
• Let the basis functions modeling the time effect be \(g_{1}(t), g_{2}(t), \ldots, g_{k}(t)\)

A

7.2 Model Specification for Effects on \(E(Y)\)

7.2.1 Common Basis Functions

• \(k\) dummy variables for \(k+1\) unique times (assumes no functional form for time but may spend many d.f.)
• \(k=1\) for linear time trend, \(g_{1}(t)=t\)
• \(k\)–order polynomial in \(t\)
• \(k+1\)–knot restricted cubic spline (one linear term, \(k-1\) nonlinear terms)

B

7.2.2 Model for Mean Profile

• A model for mean time-response profile without interactions between time and any \(X\):
  \(E[Y_{it} | X_{i}] = X_{i}\beta + \gamma_{1}g_{1}(t) + \gamma_{2}g_{2}(t) + \ldots + \gamma_{k}g_{k}(t)\)
• Model with interactions between time and some \(X\)’s: add product terms for desired interaction effects
• Example: To allow the mean time trend for subjects in group 1 (reference group) to be arbitrarily different from time trend for subjects in group 2, have a dummy variable for group 2, a time “main effect” curve with \(k\) d.f. and all \(k\) products of these time components with the dummy variable for group 2
• Time should be modeled using indicator variables only when time is really discrete, e.g., when time is in weeks and subjects were followed at exactly the intended weeks. In general time should be modeled continuously (and nonlinearly if there are more than 2 followup times) using actual visit dates instead of intended dates (Donohue et al., n.d.).

C

7.2.3 Model Specification for Treatment Comparisons

• In studies comparing two or more treatments, a response is often measured at baseline (pre-randomization)
• Analyst has the option to use this measurement as \(Y_{i0}\) or as part of \(X_{i}\)

D

NoteComments from Jim Rochon

For RCTs, I draw a sharp line at the point when the intervention begins. The LHS [left hand side of the model equation] is reserved for something that is a response to treatment. Anything before this point can potentially be included as a covariate in the regression model. This includes the “baseline” value of the outcome variable. Indeed, the best predictor of the outcome at the end of the study is typically where the patient began at the beginning. It drinks up a lot of variability in the outcome; and, the effect of other covariates is typically mediated through this variable.

I treat anything after the intervention begins as an outcome. In the western scientific method, an “effect” must follow the “cause” even if by a split second.

Note that an RCT is different than a cohort study. In a cohort study, “Time 0” is not terribly meaningful. If we want to model, say, the trend over time, it would be legitimate, in my view, to include the “baseline” value on the LHS of that regression model.

Now, even if the intervention, e.g., surgery, has an immediate effect, I would include still reserve the LHS for anything that might legitimately be considered as the response to the intervention. So, if we cleared a blocked artery and then measured the MABP, then that would still be included on the LHS.

Now, it could well be that most of the therapeutic effect occurred by the time that the first repeated measure was taken, and then levels off. Then, a plot of the means would essentially be two parallel lines and the treatment effect is the distance between the lines, i.e., the difference in the intercepts.

If the linear trend from baseline to Time 1 continues beyond Time 1, then the lines will have a common intercept but the slopes will diverge. Then, the treatment effect will the difference in slopes.

One point to remember is that the estimated intercept is the value at time 0 that we predict from the set of repeated measures post randomization. In the first case above, the model will predict different intercepts even though randomization would suggest that they would start from the same place. This is because we were asleep at the switch and didn’t record the “action” from baseline to time 1. In the second case, the model will predict the same intercept values because the linear trend from baseline to time 1 was continued thereafter.

More importantly, there are considerable benefits to including it as a covariate on the RHS. The baseline value tends to be the best predictor of the outcome post-randomization, and this maneuver increases the precision of the estimated treatment effect. Additionally, any other prognostic factors correlated with the outcome variable will also be correlated with the baseline value of that outcome, and this has two important consequences. First, this greatly reduces the need to enter a large number of prognostic factors as covariates in the linear models. Their effect is already mediated through the baseline value of the outcome variable. Secondly, any imbalances across the treatment arms in important prognostic factors will induce an imbalance across the treatment arms in the baseline value of the outcome. Including the baseline value thereby reduces the need to enter these variables as covariates in the linear models.

Senn (2006) states that temporally and logically, a “baseline cannot be a response to treatment”, so baseline and response cannot be modeled in an integrated framework.

… one should focus clearly on ‘outcomes’ as being the only values that can be influenced by treatment and examine critically any schemes that assume that these are linked in some rigid and deterministic view to ‘baseline’ values. An alternative tradition sees a baseline as being merely one of a number of measurements capable of improving predictions of outcomes and models it in this way.

The final reason that baseline cannot be modeled as the response at time zero is that many studies have inclusion/exclusion criteria that include cutoffs on the baseline variable. In other words, the baseline measurement comes from a truncated distribution. In general it is not appropriate to model the baseline with the same distributional shape as the follow-up measurements. Thus the approaches recommended by Liang & Zeger (2000) and Liu et al. (2009) are problematic¹.

E

¹ In addition to this, one of the paper’s conclusions that analysis of covariance is not appropriate if the population means of the baseline variable are not identical in the treatment groups is not correct (Senn, 2006). See Kenward et al. (2010) for a rebuke of Liu et al. (2009).

7.3 Modeling Within-Subject Dependence

• Random effects and mixed effects models have become very popular
• Disadvantages:
  • Induced correlation structure for \(Y\) may be unrealistic
  • Numerically demanding
  • Require complex approximations for distributions of test statistics
• Conditional random effects vs. (subject-) marginal models:
  • Random effects are subject-conditional
  • Random effects models are needed to estimate responses for individual subjects
  • Models without random effects are marginalized with respect to subject-specific effects
  • They are natural when the interest is on group-level (i.e., covariate-specific but not patient-specific) parameters (e.g., overall treatment effect)
  • Random effects are natural when there is clustering at more than the subject level (multi-level models)
• Extended linear model (marginal; with no random effects) is a logical extension of the univariate model (e.g., few statisticians use subject random effects for univariate \(Y\))
• This was known as growth curve models and generalized least squares (Goldstein, 1989; Potthoff & Roy, 1964) and was developed long before mixed effect models became popular
• Pinheiro and Bates (Section~5.1.2) state that “in some applications, one may wish to avoid incorporating random effects in the model to account for dependence among observations, choosing to use the within-group component \(\Lambda_{i}\) to directly model variance-covariance structure of the response.”
• We will assume that \(Y_{it} | X_{i}\) has a multivariate normal distribution with mean given above and with variance-covariance matrix \(V_{i}\), an \(n_{i}\times n_{i}\) matrix that is a function of \(t_{i1}, \ldots, t_{in_{i}}\)
• We further assume that the diagonals of \(V_{i}\) are all equal
• Procedure can be generalized to allow for heteroscedasticity over time or with respect to \(X\) (e.g., males may be allowed to have a different variance than females)
• This extended linear model has the following assumptions:
  • all the assumptions of OLS at a single time point including correct modeling of predictor effects and univariate normality of responses conditional on \(X\)
  • the distribution of two responses at two different times for the same subject, conditional on \(X\), is bivariate normal with a specified correlation coefficient
  • the joint distribution of all \(n_{i}\) responses for the \(i^{th}\) subject is multivariate normal with the given correlation pattern (which implies the previous two distributional assumptions)
  • responses from any times for any two different subjects are uncorrelated

FGH

What Methods To Use for Repeated Measurements / Serial Data? ^{2 3}

	Repeated Measures ANOVA	GEE	Mixed Effects Models	GLS	Markov	LOCF	Summary Statistic4
Assumes normality	×		×	×
Assumes independence of measurements within subject	×5	×6
Assumes a correlation structure7	×	×8	×	×	×
Requires same measurement times for all subjects	×					?
Does not allow smooth modeling of time to save d.f.	×
Does not allow adjustment for baseline covariates	×
Does not easily extend to non-continuous \(Y\)	×			×
Loses information by not using intermediate measurements						×9	×
Does not allow widely varying # observations per subject	×	×10				×	×11
Does not allow for subjects to have distinct trajectories12	×	×		×	×	×
Assumes subject-specific effects are Gaussian			×
Badly biased if non-random dropouts	?	×				×
Biased in general						×
Harder to get tests & CLs			×13			×14
Requires large # subjects/clusters		×
SEs are wrong	×15					×
Assumptions are not verifiable in small samples	×	N/A	×	×		×
Does not extend to complex settings such as time-dependent covariates and dynamic 16 models	×		×	×		×	?

² Thanks to Charles Berry, Brian Cade, Peter Flom, Bert Gunter, and Leena Choi for valuable input.

³ GEE: generalized estimating equations; GLS: generalized least squares; LOCF: last observation carried forward.

⁴ E.g., compute within-subject slope, mean, or area under the curve over time. Assumes that the summary measure is an adequate summary of the time profile and assesses the relevant treatment effect.

⁵ Unless one uses the Huynh-Feldt or Greenhouse-Geisser correction

⁶ For full efficiency, if using the working independence model

⁷ Or requires the user to specify one

⁸ For full efficiency of regression coefficient estimates

⁹ Unless the last observation is missing

¹⁰ The cluster sandwich variance estimator used to estimate SEs in GEE does not perform well in this situation, and neither does the working independence model because it does not weight subjects properly.

¹¹ Unless one knows how to properly do a weighted analysis

¹² Or users population averages

¹³ Unlike GLS, does not use standard maximum likelihood methods yielding simple likelihood ratio \(\chi^2\) statistics. Requires high-dimensional integration to marginalize random effects, using complex approximations, and if using SAS, unintuitive d.f. for the various tests.

¹⁴ Because there is no correct formula for SE of effects; ordinary SEs are not penalized for imputation and are too small

¹⁵ If correction not applied

¹⁶ E.g., a model with a predictor that is a lagged value of the response variable

• Markov models use ordinary univariate software and are very flexible
• They apply the same way to binary, ordinal, nominal, and continuous Y
• They require post-fitting calculations to get probabilities, means, and quantiles that are not conditional on the previous Y value

I

Gardiner et al. (2009) compared several longitudinal data models, especially with regard to assumptions and how regression coefficients are estimated. Peters et al. (2012) have an empirical study confirming that the “use all available data” approach of likelihood–based longitudinal models makes imputation of follow-up measurements unnecessary.

J

7.4 Parameter Estimation Procedure

• Generalized least squares
• Like weighted least squares but uses a covariance matrix that is not diagonal
• Each subject can have her own shape of \(V_{i}\) due to each subject being measured at a different set of times
• Maximum likelihood
• Newton-Raphson or other trial-and-error methods used for estimating parameters
• For small number of subjects, advantages in using REML (restricted maximum likelihood) instead of ordinary MLE (Diggle et al., 2002, p. Section~5.3), (Pinheiro & Bates, 2000, p. Chapter~5), Goldstein (1989) (esp. to get more unbiased estimate of the covariance matrix)
• When imbalances are not severe, OLS fitted ignoring subject identifiers may be efficient
  • But OLS standard errors will be too small as they don’t take intra-cluster correlation into account
  • May be rectified by substituting covariance matrix estimated from Huber-White cluster sandwich estimator or from cluster bootstrap
• When imbalances are severe and intra-subject correlations are strong, OLS is not expected to be efficient because it gives equal weight to each observation
  • a subject contributing two distant observations receives \(\frac{1}{5}\) the weight of a subject having 10 tightly-spaced observations

KLM

7.5 Common Correlation Structures

• Usually restrict ourselves to isotropic correlation structures — correlation between responses within subject at two times depends only on a measure of distance between the two times, not the individual times (this will be relaxed in the last structure covered below)
• We simplify further and assume depends on \(|t_{1} - t_{2}|\)
• Can speak interchangeably of correlations of residuals within subjects or correlations between responses measured at different times on the same subject, conditional on covariates \(X\)
• Assume that the correlation coefficient for \(Y_{it_{1}}\) vs. \(Y_{it_{2}}\) conditional on baseline covariates \(X_{i}\) for subject \(i\) is \(h(|t_{1} - t_{2}|, \rho)\), where \(\rho\) is a vector (usually a scalar) set of fundamental correlation parameters
• Some commonly used structures when times are continuous and are not equally spaced (Pinheiro & Bates, 2000, Section 5.3.3) (nlme correlation function names are at the right if the structure is implemented in nlme):

NO

Table 7.1: Some longitudinal data correlation structures

Structure	nlme Function
Compound symmetry: \(h = \rho\) if \(t_{1} \neq t_{2}\), 1 if \(t_{1}=t_{2}\) 17	corCompSymm
Autoregressive-moving average lag 1: \(h = \rho^{|t_{1} - t_{2}|} = \rho^s\) where \(s = |t_{1}-t_{2}|\)	corCAR1
Exponential: \(h = \exp(-s/\rho)\)	corExp
Gaussian: \(h = \exp[-(s/\rho)^2]\)	corGaus
Linear: \(h = (1 - s/\rho)[s < \rho]\)	corLin
Rational quadratic: \(h = 1 - (s/\rho)^{2}/[1+(s/\rho)^{2}]\)	corRatio
Spherical: \(h = [1-1.5(s/\rho)+0.5(s/\rho)^{3}][s < \rho]\)	corSpher
Linear exponent AR(1): \(h = \rho^{d_{min} + \delta\frac{s - d_{min}}{d_{max} - d_{min}}}\), 1 if \(t_{1}=t_{2}\)	Simpson et al. (2010)
Exponential decline with floor and non-isotropy: \(h = k + (1-k) \exp(-\exp(a+b\frac{t_{1}+t_{2}}{2})|t_{1}-t_{2}|)\)	corFloorExp

¹⁷ Essentially what two-way ANOVA assumes

The structures 3-7 use \(\rho\) as a scaling parameter, not as something restricted to be in \([0,1]\)

The last structure is implemented in the rms package. The \(k\) parameter is restricted to be in \([0,1]\) by modeling it as \(\text{expit}(k^{*})\) where \(k^{*}\) has an unrestricted range and \(\text{expit}(x) = \frac{1}{1 + \exp(-x)}\). This structure allows for non-isotropy and for a nonzero limiting correlation as the time lag \(\rightarrow \infty\). It is a shared random effect plus serial decay mixture (a continuous-time analog of a random intercept superimposed on an Ornstein-Uhlenbeck exponential-decay process). \(k\) is the correlation contributed by a persistent, subject-level component that never fades; \(1 - k\) is the share of correlation carried by a component that decays with lag at rate \(\exp(a + b\frac{t_{1}+t_{2}}{2})\). \(b\) is the non-isotropy term, letting the rate grow or shrink with the time pair’s midpoint.

7.6 Checking Model Fit

• Constant variance assumption: usual residual plots
• Normality assumption: usual qq residual plots
• Correlation pattern: Variogram
  • Estimate correlations of all possible pairs of residuals at different time points
  • Pool all estimates at same absolute difference in time \(s\)
  • Variogram is a plot with \(y = 1 - \hat{h}(s, \rho)\) vs. \(s\) on the \(x\)-axis
  • Superimpose the theoretical variogram assumed by the model
  • This does not allow for non-isotropy

P

7.7 R Software

• Nonlinear mixed effects model package of Pinheiro & Bates
• For linear models, fitting functions are
  • lme for mixed effects models
  • gls for generalized least squares without random effects
• For this version the rms package has Gls so that many features of rms can be used:
  • anova: all partial Wald tests, test of linearity, pooled tests
  • summary: effect estimates (differences in \(\hat{Y}\)) and confidence limits, can be plotted
  • plot, ggplot, plotp: continuous effect plots
  • nomogram: nomogram
  • Function: generate R function code for fitted model
  • latex:  representation of fitted model

Q

In addition, Gls has a bootstrap option (hence you do not use rms’s bootcov for Gls fits).
To get regular gls functions named anova (for likelihood ratio tests, AIC, etc.) or summary use anova.gls or summary.gls * nlme package has many graphics and fit-checking functions * Several functions will be demonstrated in the case study

7.8 Case Study

Consider the dataset in Table~6.9 of Davis[davis-repmeas, pp. 161-163] from a multi-center, randomized controlled trial of botulinum toxin type B (BotB) in patients with cervical dystonia from nine U.S. sites.

• Randomized to placebo (\(N=36\)), 5000 units of BotB (\(N=36\)), 10,000 units of BotB (\(N=37\))
• Response variable: total score on Toronto Western Spasmodic Torticollis Rating Scale (TWSTRS), measuring severity, pain, and disability of cervical dystonia (high scores mean more impairment)
• TWSTRS measured at baseline (week 0) and weeks 2, 4, 8, 12, 16 after treatment began
• Dataset cdystonia from web site

R

7.8.1 Graphical Exploration of Data

require(rms)
require(data.table)
options(prType='html')    # for model print, summary, anova, validate
getHdata(cdystonia)
setDT(cdystonia)          # convert to data.table
cdystonia[, uid := paste(site, id)]   # unique subject ID

# Tabulate patterns of subjects' time points
g <- function(w) paste(sort(unique(w)), collapse=' ')
cdystonia[, table(tapply(week, uid, g))]

            0         0 2 4   0 2 4 12 16       0 2 4 8    0 2 4 8 12 
            1             1             3             1             1 
0 2 4 8 12 16    0 2 4 8 16   0 2 8 12 16   0 4 8 12 16      0 4 8 16 
           94             1             2             4             1 

# Plot raw data, superposing subjects
xl <- xlab('Week'); yl <- ylab('TWSTRS-total score')
ggplot(cdystonia, aes(x=week, y=twstrs, color=factor(id))) +
       geom_line() + xl + yl + facet_grid(treat ~ site) +
       guides(color=FALSE)

Figure 7.1: Time profiles for individual subjects, stratified by study site and dose

# Show quartiles
g <- function(x) {
  k <- as.list(quantile(x, (1 : 3) / 4, na.rm=TRUE))
  names(k) <- .q(Q1, Q2, Q3)
  k
}
cdys <- cdystonia[, g(twstrs), by=.(treat, week)]
ggplot(cdys, aes(x=week, y=Q2)) + xl + yl + ylim(0, 70) +
  geom_line() + facet_wrap(~ treat, nrow=2) +
  geom_ribbon(aes(ymin=Q1, ymax=Q3), alpha=0.2)

Figure 7.2: Quartiles of TWSTRS stratified by dose

# Show means with bootstrap nonparametric CLs
cdys <-  cdystonia[, as.list(smean.cl.boot(twstrs)),
                   by = list(treat, week)]
ggplot(cdys, aes(x=week, y=Mean)) + xl + yl + ylim(0, 70) +
  geom_line() + facet_wrap(~ treat, nrow=2) +
  geom_ribbon(aes(x=week, ymin=Lower, ymax=Upper), alpha=0.2)

Figure 7.3: Mean responses and nonparametric bootstrap 0.95 confidence limits for population means, stratified by dose

Model with \(Y_{i0}\) as Baseline Covariate

baseline <- cdystonia[week == 0]
baseline[, week := NULL]
setnames(baseline, 'twstrs', 'twstrs0')
followup <- cdystonia[week > 0, .(uid, week, twstrs)]
setkey(baseline, uid)
setkey(followup, uid, week)
both     <- Merge(baseline, followup, id = ~ uid)

         Vars Obs Unique IDs IDs in #1 IDs not in #1
baseline    7 109        109        NA            NA
followup    3 522        108       108             0
Merged      9 523        109       109             0

Number of unique IDs in any data frame : 109 
Number of unique IDs in all data frames: 108 

# Remove person with no follow-up record
both     <- both[! is.na(week)]
dd       <- datadist(both)
options(datadist='dd')

7.8.2 Using Generalized Least Squares

We stay with baseline adjustment and use a variety of correlation structures, with constant variance. Time is modeled as a restricted cubic spline with 3 knots, because there are only 3 unique interior values of week.

S

require(nlme)
cp <- list(corCAR1,corExp,corCompSymm,corLin,corGaus,corSpher,corFloorExp)
z  <- vector('list', length(cp))
for (k in 1:length(cp)) {
  z[[k]] <- gls(
    twstrs ~ treat * rcs(week, 3) +
      rcs(twstrs0, 3) + rcs(age, 4) * sex,
    data = both,
    correlation = cp[[k]](form = ~ week | uid)
  )
}
anova(z[[1]],z[[2]],z[[3]],z[[4]],z[[5]],z[[6]],z[[7]], test=FALSE)

       Model df      AIC      BIC    logLik
z[[1]]     1 20 3553.906 3638.357 -1756.953
z[[2]]     2 20 3553.906 3638.357 -1756.953
z[[3]]     3 20 3587.974 3672.426 -1773.987
z[[4]]     4 20 3575.079 3659.531 -1767.540
z[[5]]     5 20 3621.081 3705.532 -1790.540
z[[6]]     6 20 3570.958 3655.409 -1765.479
z[[7]]     7 22 3540.028 3632.925 -1748.014

AIC computed above is set up so that smaller values are better. From this the 3-parameter correlation structure is the best fitting. For the remainder of the analysis use corFloorExp, using Gls.

Keselman et al. (1998) did a simulation study to study the reliability of AIC for selecting the correct covariance structure in repeated measurement models. In choosing from among 11 structures, AIC selected the correct structure 47% of the time. Gurka et al. (2011) demonstrated that fixed effects in a mixed effects model can be biased, independent of sample size, when the specified covariate matrix is more restricted than the true one.

a <- Gls(twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) +
         rcs(age, 4) * sex, data=both,
         correlation=corFloorExp(form=~week | uid))
a

Generalized Least Squares Fit by REML

Gls(model = twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) + 
    rcs(age, 4) * sex, data = both, correlation = corFloorExp(form = ~week | 
    uid))

Obs522	Log-restricted-likelihood-1748.01
Clusters108	Model d.f.17
g11.481	σ8.5635
	d.f.504

	β	S.E.	t	P(>|t|)
Intercept	-2.1683	12.4446	-0.17	0.8617
treat=5000U	0.5425	2.6007	0.21	0.8348
treat=Placebo	7.2832	2.6176	2.78	0.0056
week	0.3533	0.3012	1.17	0.2414
week'	0.6672	0.2957	2.26	0.0245
twstrs0	0.8181	0.1518	5.39	<0.0001
twstrs0'	0.2180	0.1881	1.16	0.2470
age	-0.0941	0.2461	-0.38	0.7022
age'	0.6436	0.6794	0.95	0.3440
age''	-3.2761	2.6813	-1.22	0.2223
sex=M	23.6319	19.5183	1.21	0.2266
treat=5000U × week	0.0589	0.4273	0.14	0.8905
treat=Placebo × week	-0.1469	0.4297	-0.34	0.7326
treat=5000U × week'	-0.4238	0.4188	-1.01	0.3121
treat=Placebo × week'	-0.6140	0.4207	-1.46	0.1450
age × sex=M	-0.5641	0.4661	-1.21	0.2267
age' × sex=M	1.3575	1.2976	1.05	0.2960
age'' × sex=M	-3.5408	5.0401	-0.70	0.4827

Correlation Structure: corFloorExp
 Formula: ~week | uid 
 Parameter estimate(s):
          k           a           b 
 0.17571863 -1.08095556 -0.08066475 

cs <- a$modelStruct$corStruct
cs

Correlation structure of class corFloorExp representing
  corr(t1,t2) = k + (1-k)*exp(-exp(a+b*(t1+t2)/2)*|t1-t2|)
          k           a           b 
 0.17571863 -1.08095556 -0.08066475 

Plot the estimated correlation function for all observed pairs of times.

co <- coef(cs)
g <- function(t1, t2) {
  k <- plogis(co["k"])
  a <- co["a"]
  b <- co["b"]
  k + (1 - k) * exp(-exp(a + b * (t1 + t2) / 2) * abs(t1 - t2))
}
times <- sort(unique(both$week))
w <- setDT(expand.grid(t1 = times, t2 = times))
w <- w[t2 > t1]
w[, lag := t2 - t1]
w[, mid := (t1 + t2) / 2]
w[, rho := g(t1, t2)]
ggplot(w, aes(x = lag, y = rho, label=mid)) +
  geom_text(size=2.5) +
  xlab(expression(t[2] - t[1])) + ylab(expression(hat(rho)))

Figure 7.4: Fitted 3-parameter correlation \(\hat{\rho}\) as a function of time lag (\(x\)) and time midpoint (numbers)

The numbers drawn for each point are the values of \(\frac{t_{1} + t_{2}}{2}\). The degree to which points for the same time lag have different vertical values is the degree to which time midpoints figure into the correlation structure in addition to time lags, i.e., the degree of non-isotropy. We see a moderate degree of non-isotropy. We also see the possibility of greater-than-zero long-term correlation, i.e., a need for a compound symmetric correlation term and slight inadequacy of AR(1).

\(\hat{\rho} = 0.8672\) from an AR(1) fit, the estimate of the correlation between two T measurements taken one week apart on the same subject. The estimated correlation for measurements 10 weeks apart is \(0.8672^{10} = 0.24\).

v <- Variogram(a, form=~ week | uid)
plot(v)

Figure 7.5: Variogram, with assumed correlation pattern superimposed

Check constant variance and normality assumptions:

U

both$resid <- r <- resid(a); both$fitted <- fitted(a)
yl <- ylab('Residuals')
p1 <- ggplot(both, aes(x=fitted, y=resid)) + geom_point() +
      facet_grid(~ treat) + yl
p2 <- ggplot(both, aes(x=twstrs0, y=resid)) + geom_point()+yl
p3 <- ggplot(both, aes(x=week, y=resid)) + yl + ylim(-20,20) +
      stat_summary(fun.data="mean_sdl", geom='smooth')
p4 <- ggplot(both, aes(sample=resid)) + stat_qq() +
      geom_abline(intercept=mean(r), slope=sd(r)) + yl
gridExtra::grid.arrange(p1, p2, p3, p4, ncol=2)

Figure 7.6: Three residual plots to check for absence of trends in central tendency and in variability. Upper right panel shows the baseline score on the \(x\)-axis. Bottom left panel shows the mean \(\pm 2\times\) SD. Bottom right panel is the QQ plot for checking normality of residuals from the GLS fit.

Now get hypothesis tests, estimates, and graphically interpret the model.

anova(a)

Wald Statistics for twstrs

	χ2	d.f.	P
treat (Factor+Higher Order Factors)	27.50	6	0.0001
All Interactions	19.39	4	0.0007
week (Factor+Higher Order Factors)	102.07	6	<0.0001
All Interactions	19.39	4	0.0007
Nonlinear (Factor+Higher Order Factors)	5.80	3	0.1219
twstrs0	219.82	2	<0.0001
Nonlinear	1.34	1	0.2465
age (Factor+Higher Order Factors)	8.21	6	0.2231
All Interactions	4.41	3	0.2205
Nonlinear (Factor+Higher Order Factors)	6.56	4	0.1608
sex (Factor+Higher Order Factors)	5.21	4	0.2668
All Interactions	4.41	3	0.2205
treat × week (Factor+Higher Order Factors)	19.39	4	0.0007
Nonlinear	2.24	2	0.3263
Nonlinear Interaction : f(A,B) vs. AB	2.24	2	0.3263
age × sex (Factor+Higher Order Factors)	4.41	3	0.2205
Nonlinear	3.48	2	0.1753
Nonlinear Interaction : f(A,B) vs. AB	3.48	2	0.1753
TOTAL NONLINEAR	13.24	8	0.1040
TOTAL INTERACTION	23.77	7	0.0013
TOTAL NONLINEAR + INTERACTION	31.06	11	0.0011
TOTAL	332.95	17	<0.0001

plot(anova(a))

Figure 7.7: Results of anova.rms from generalized least squares fit with continuous time AR1 correlation structure

ylm <- ylim(25, 60)
p1 <- ggplot(Predict(a, week, treat, conf.int=FALSE),
             adj.subtitle=FALSE, legend.position='top') + ylm
p2 <- ggplot(Predict(a, twstrs0), adj.subtitle=FALSE) + ylm
p3 <- ggplot(Predict(a, age, sex), adj.subtitle=FALSE,
             legend.position='top') + ylm
gridExtra::grid.arrange(p1, p2, p3, ncol=2)

Figure 7.8: Estimated effects of time, baseline TWSTRS, age, and sex

summary(a)  # Shows for week 8

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
week	4	12	8	6.8290	1.061	4.750	8.909
twstrs0	39	53	14	13.7600	0.928	11.940	15.580
age	46	65	19	2.3580	2.147	-1.851	6.567
treat --- 5000U:10000U	1	2		0.5895	1.980	-3.291	4.470
treat --- Placebo:10000U	1	3		5.4940	1.986	1.602	9.386
sex --- M:F	1	2		-1.1110	1.864	-4.763	2.542

# To get results for week 8 for a different reference group
# for treatment, use e.g. summary(a, week=4, treat='Placebo')

# Compare low dose with placebo, separately at each time
k1 <- contrast(a, list(week=c(2,4,8,12,16), treat='5000U'),
                  list(week=c(2,4,8,12,16), treat='Placebo'))
options(width=80)
print(k1, digits=3)

    week twstrs0 age sex Contrast S.E.  Lower  Upper     Z P(>|z|)
1      2      46  56   F    -6.33 2.08 -10.41 -2.248 -3.04  0.0024
2      4      46  56   F    -5.92 1.77  -9.38 -2.451 -3.35  0.0008
3      8      46  56   F    -4.90 2.00  -8.82 -0.992 -2.46  0.0140
4*    12      46  56   F    -3.13 1.84  -6.74  0.484 -1.70  0.0896
5*    16      46  56   F    -1.17 2.10  -5.29  2.956 -0.55  0.5794

Redundant contrasts are denoted by *

Confidence intervals are 0.95 individual intervals

# Compare high dose with placebo
k2 <- contrast(a, list(week=c(2,4,8,12,16), treat='10000U'),
                  list(week=c(2,4,8,12,16), treat='Placebo'))
print(k2, digits=3)

    week twstrs0 age sex Contrast S.E.  Lower Upper     Z P(>|z|)
1      2      46  56   F    -6.99 2.05 -11.01 -2.97 -3.41  0.0007
2      4      46  56   F    -6.70 1.74 -10.11 -3.28 -3.84  0.0001
3      8      46  56   F    -5.49 1.99  -9.39 -1.60 -2.77  0.0057
4*    12      46  56   F    -1.84 1.83  -5.43  1.75 -1.00  0.3163
5*    16      46  56   F     2.44 2.09  -1.66  6.53  1.17  0.2435

Redundant contrasts are denoted by *

Confidence intervals are 0.95 individual intervals

k1 <- as.data.frame(k1[c('week', 'Contrast', 'Lower', 'Upper')])
p1 <- ggplot(k1, aes(x=week, y=Contrast)) + geom_point() +
      geom_line() + ylab('Low Dose - Placebo') +
      geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0)
k2 <- as.data.frame(k2[c('week', 'Contrast', 'Lower', 'Upper')])
p2 <- ggplot(k2, aes(x=week, y=Contrast)) + geom_point() +
      geom_line() + ylab('High Dose - Placebo') +
      geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0)
gridExtra::grid.arrange(p1, p2, ncol=2)

Figure 7.9: Contrasts and 0.95 confidence limits from GLS fit

Although multiple d.f. tests such as total treatment effects or treatment \(\times\) time interaction tests are comprehensive, their increased degrees of freedom can dilute power. In a treatment comparison, treatment contrasts at the last time point (single d.f. tests) are often of major interest. Such contrasts are informed by all the measurements made by all subjects (up until dropout times) when a smooth time trend is assumed.

V

n <- nomogram(a, age=c(seq(20, 80, by=10), 85))
plot(n, cex.axis=.55, cex.var=.8, lmgp=.25)  # Figure (*\ref{fig:longit-nomogram}*)

Figure 7.10: Nomogram from GLS fit. Second axis is the baseline score.

7.8.3 Semiparametric Random Effects Models

f <- orm(twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) +
         rcs(age, 4) * sex + cluster(uid), data=both, x=TRUE, y=TRUE, maxit.outer=300, trace=1)

NR iteration:1  -2LL:4078.9355  Max |gradient|:2069.607  Max |change in parameters|:10.20047
NR iteration:2  -2LL:3651.1652  Max |gradient|:1088.633  Max |change in parameters|:1.581404
NR iteration:3  -2LL:3572.9837  Max |gradient|:34.15156  Max |change in parameters|:0.7807072
NR iteration:4  -2LL:3571.4191  Max |gradient|:0.2065217  Max |change in parameters|:0.04135248
NR iteration:5  -2LL:3571.4155  Max |gradient|:0.006402436  Max |change in parameters|:0.000159527
NR iteration:6  -2LL:3571.4155  Max |gradient|:2.906289e-07  Max |change in parameters|:5.923429e-09
NR iteration:1  -2LL:4078.9355  Max |gradient|:1.190159e-12  Max |change in parameters|:9.2813e-15
nAGQ: 7  outer iter: 3  sigma: 2.262627  max|grad|: 2.508888  max|delta|: 0.6315961 
nAGQ: 7  outer iter: 6  sigma: 2.839968  max|grad|: 0.1056326  max|delta|: 0.04206848 
nAGQ: 7  outer iter: 9  sigma: 2.850974  max|grad|: 0.0110501  max|delta|: 0.002434221 
nAGQ: 7  outer iter: 12  sigma: 2.851705  max|grad|: 0.005924715  max|delta|: 0.001439909 
nAGQ: 7  outer iter: 15  sigma: 2.852235  max|grad|: 0.003328131  max|delta|: 0.0004303428 
nAGQ: 7  outer iter: 18  sigma: 2.852391  max|grad|: 0.001463528  max|delta|: 0.0002864179 
nAGQ: 7  outer iter: 21  sigma: 2.852458  max|grad|: 0.001335273  max|delta|: 0.0001509185 
nAGQ: 7  outer iter: 24  sigma: 2.852515  max|grad|: 0.0004278585  max|delta|: 6.43989e-05 
nAGQ: 7  outer iter: 3  sigma: 1.624202  max|grad|: 52.7541  max|delta|: 1.061674 
nAGQ: 7  outer iter: 6  sigma: 1.734737  max|grad|: 16.46641  max|delta|: 0.08589558 
nAGQ: 7  outer iter: 9  sigma: 1.736023  max|grad|: 0.09771613  max|delta|: 0.007871001 
nAGQ: 7  outer iter: 12  sigma: 1.735865  max|grad|: 0.02428374  max|delta|: 0.0003292693 
nAGQ: 7  outer iter: 15  sigma: 1.735855  max|grad|: 0.005470436  max|delta|: 8.469005e-05 
nAGQ: 7  outer iter: 18  sigma: 1.735857  max|grad|: 0.00014961  max|delta|: 1.014734e-05 

f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) + 
    rcs(age, 4) * sex + cluster(uid), data = both, x = TRUE, 
    y = TRUE, maxit.outer = 300, trace = 1)

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs522	LR χ2298.48	R20.436	ρ0.781
ESS521.7	d.f.17	R217,5220.417	Dxy0.600
Distinct Y62	P(>χ2)<0.0001	R217,521.70.417
Cluster onuid		|P(Y ≥ median)-½|0.350
Clusters108
σγ1.7359 [1.4321, 2.104]
Y0.542
max |∂log L/∂β|0.0001

	β	S.E.	Wald Z	P(>|Z|)
treat=5000U	0.1063	0.6978	0.15	0.8789
treat=Placebo	2.3443	0.7106	3.30	0.0010
week	0.1210	0.0795	1.52	0.1277
week'	0.1891	0.0872	2.17	0.0301
twstrs0	0.2253	0.0472	4.77	<0.0001
twstrs0'	0.1280	0.0581	2.20	0.0276
age	-0.0156	0.0741	-0.21	0.8336
age'	0.1927	0.2041	0.94	0.3452
age''	-1.0597	0.8037	-1.32	0.1873
sex=M	4.8434	5.8230	0.83	0.4055
treat=5000U × week	0.0499	0.1099	0.45	0.6502
treat=Placebo × week	-0.0554	0.1122	-0.49	0.6217
treat=5000U × week'	-0.1610	0.1203	-1.34	0.1807
treat=Placebo × week'	-0.1365	0.1226	-1.11	0.2656
age × sex=M	-0.1065	0.1392	-0.77	0.4441
age' × sex=M	0.1479	0.3872	0.38	0.7025
age'' × sex=M	0.0330	1.5021	0.02	0.9825

Olinks(f)

      link null.deviance deviance      AIC      LR    R2
1 logistic      3713.175 3414.694 3570.694 298.481 0.417
2   probit      3732.011 3443.773 3599.773 288.239 0.405
3   loglog      3748.263 3463.493 3619.493 284.770 0.401
4  cloglog      3712.715 3438.655 3594.655 274.061 0.389

Stick with the PO model.

both[, mix := (week - 2) / 14]
f <- orm(twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) +
         rcs(age, 4) * sex + cluster(uid) + mix_re(mix), data=both, maxit.outer=300, trace=1)

NR iteration:1  -2LL:4078.9355  Max |gradient|:2069.607  Max |change in parameters|:10.20047
NR iteration:2  -2LL:3651.1652  Max |gradient|:1088.633  Max |change in parameters|:1.581404
NR iteration:3  -2LL:3572.9837  Max |gradient|:34.15156  Max |change in parameters|:0.7807072
NR iteration:4  -2LL:3571.4191  Max |gradient|:0.2065217  Max |change in parameters|:0.04135248
NR iteration:5  -2LL:3571.4155  Max |gradient|:0.006402436  Max |change in parameters|:0.000159527
NR iteration:6  -2LL:3571.4155  Max |gradient|:2.906289e-07  Max |change in parameters|:5.923429e-09
NR iteration:1  -2LL:4078.9355  Max |gradient|:1.190159e-12  Max |change in parameters|:9.2813e-15
nAGQ: 7  outer iter: 3  sigma: 2.83153  sigma2: 2.120081  max|grad|: 8.09825  max|delta|: 0.8922446 
nAGQ: 7  outer iter: 6  sigma: 3.170743  sigma2: 2.374241  max|grad|: 8.639543  max|delta|: 0.04121938 
nAGQ: 7  outer iter: 9  sigma: 3.250566  sigma2: 2.393887  max|grad|: 6.872054  max|delta|: 0.04751811 
nAGQ: 7  outer iter: 12  sigma: 3.550237  sigma2: 2.443301  max|grad|: 2.376641  max|delta|: 0.04821976 
nAGQ: 7  outer iter: 15  sigma: 3.603117  sigma2: 2.449761  max|grad|: 2.356168  max|delta|: 0.04598941 
nAGQ: 7  outer iter: 18  sigma: 3.821636  sigma2: 2.483668  max|grad|: 1.105287  max|delta|: 0.03885251 
nAGQ: 7  outer iter: 21  sigma: 3.862305  sigma2: 2.492536  max|grad|: 1.230085  max|delta|: 0.03641951 
nAGQ: 7  outer iter: 24  sigma: 4.025185  sigma2: 2.533644  max|grad|: 0.6503175  max|delta|: 0.02674235 
nAGQ: 7  outer iter: 27  sigma: 4.05366  sigma2: 2.542306  max|grad|: 0.7158627  max|delta|: 0.02490802 
nAGQ: 7  outer iter: 30  sigma: 4.163885  sigma2: 2.578168  max|grad|: 0.4057794  max|delta|: 0.01692853 
nAGQ: 7  outer iter: 33  sigma: 4.368707  sigma2: 2.65622  max|grad|: 0.9428735  max|delta|: 0.009963075 
nAGQ: 7  outer iter: 36  sigma: 4.369212  sigma2: 2.658095  max|grad|: 0.004631852  max|delta|: 0.0007334673 
nAGQ: 7  outer iter: 39  sigma: 4.370057  sigma2: 2.658564  max|grad|: 0.01518074  max|delta|: 0.0005907678 
nAGQ: 7  outer iter: 42  sigma: 4.373283  sigma2: 2.659951  max|grad|: 0.01029467  max|delta|: 0.0004003906 
nAGQ: 7  outer iter: 45  sigma: 4.373784  sigma2: 2.660177  max|grad|: 0.009828257  max|delta|: 0.0003990524 
nAGQ: 7  outer iter: 48  sigma: 4.375643  sigma2: 2.660974  max|grad|: 0.00643406  max|delta|: 0.0002541781 
nAGQ: 7  outer iter: 51  sigma: 4.375933  sigma2: 2.6611  max|grad|: 0.005918061  max|delta|: 0.0002362594 
nAGQ: 7  outer iter: 54  sigma: 4.376999  sigma2: 2.661558  max|grad|: 0.003803108  max|delta|: 0.0001487913 
nAGQ: 7  outer iter: 57  sigma: 4.378704  sigma2: 2.662286  max|grad|: 0.0001066609  max|delta|: 4.404371e-05 
nAGQ: 7  outer iter: 3  sigma: 1.473977  sigma2: 1.410774  max|grad|: 80.87709  max|delta|: 1.384612 
nAGQ: 7  outer iter: 6  sigma: 1.63064  sigma2: 1.7726  max|grad|: 34.62541  max|delta|: 0.133379 
nAGQ: 7  outer iter: 9  sigma: 1.644124  sigma2: 1.841945  max|grad|: 3.75309  max|delta|: 0.05638793 
nAGQ: 7  outer iter: 12  sigma: 1.645255  sigma2: 1.847323  max|grad|: 2.001056  max|delta|: 0.006857605 
nAGQ: 7  outer iter: 15  sigma: 1.644087  sigma2: 1.850225  max|grad|: 0.2344561  max|delta|: 0.003300987 
nAGQ: 7  outer iter: 18  sigma: 1.644026  sigma2: 1.850494  max|grad|: 0.129864  max|delta|: 0.0004338185 
nAGQ: 7  outer iter: 21  sigma: 1.643918  sigma2: 1.850679  max|grad|: 0.01006982  max|delta|: 0.0001860752 
nAGQ: 7  outer iter: 24  sigma: 1.643916  sigma2: 1.850691  max|grad|: 0.005985366  max|delta|: 1.985624e-05 
nAGQ: 7  outer iter: 27  sigma: 1.64391  sigma2: 1.8507  max|grad|: 0.0002136181  max|delta|: 6.922586e-06 

f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) + 
    rcs(age, 4) * sex + cluster(uid) + mix_re(mix), data = both, 
    maxit.outer = 300, trace = 1)

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs522	LR χ2284.25	R20.420	ρ0.781
ESS521.7	d.f.17	R217,5220.401	Dxy0.600
Distinct Y62	P(>χ2)<0.0001	R217,521.70.401
Cluster onuid		|P(Y ≥ median)-½|0.351
Clusters108
σ11.6439 [1.2765, 2.1171]
σ21.8507 [1.3712, 2.3302]
Y0.542
max |∂log L/∂β|0.0002

	β	S.E.	Wald Z	P(>|Z|)
treat=5000U	0.1024	0.6816	0.15	0.8805
treat=Placebo	2.3625	0.6952	3.40	0.0007
week	0.1196	0.0796	1.50	0.1330
week'	0.1913	0.0873	2.19	0.0285
twstrs0	0.2234	0.0472	4.73	<0.0001
twstrs0'	0.1343	0.0587	2.29	0.0222
age	-0.0094	0.0751	-0.12	0.9006
age'	0.1648	0.2085	0.79	0.4292
age''	-0.9589	0.8169	-1.17	0.2405
sex=M	5.1700	5.8538	0.88	0.3771
treat=5000U × week	0.0504	0.1101	0.46	0.6472
treat=Placebo × week	-0.0571	0.1124	-0.51	0.6111
treat=5000U × week'	-0.1611	0.1203	-1.34	0.1807
treat=Placebo × week'	-0.1353	0.1227	-1.10	0.2704
age × sex=M	-0.1154	0.1400	-0.82	0.4096
age' × sex=M	0.1782	0.3896	0.46	0.6474
age'' × sex=M	-0.0785	1.5098	-0.05	0.9586

deviance(f)

                 intercepts                intercepts+x 
                   4078.936                    3571.416 
  intercepts+random effects intercepts+x+random effects 
                   3698.472                    3414.225 

Deviance did not improve, due to the estimates of \(\sigma_1\) and \(\sigma_2\) being so close.

7.8.4 Bayesian Proportional Odds Random Effects Model

• Develop a \(y\)-transformation invariant longitudinal model
• Proportional odds model with no grouping of TWSTRS scores
• Bayesian random effects model
• Random effects Gaussian with exponential prior distribution for its SD, with mean 1.0
• Compound symmetry correlation structure
• Demonstrates a large amount of patient-to-patient intercept variability

W

require(rmsb)
cmdstanr::set_cmdstan_path(cmdstan.loc)
# cmdstan.loc is defined in ~/.Rprofile
options(mc.cores=parallel::detectCores() - 1, rmsb.backend='cmdstan')
bpo <- blrm(twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) +
            rcs(age, 4) * sex + cluster(uid), data=both, file='bpo.rds')

Running MCMC with 4 chains, at most 9 in parallel...

Chain 2 finished in 4.3 seconds.
Chain 4 finished in 4.4 seconds.
Chain 1 finished in 4.5 seconds.
Chain 3 finished in 4.6 seconds.

All 4 chains finished successfully.
Mean chain execution time: 4.4 seconds.
Total execution time: 4.7 seconds.

# file= means that after the first time the model is run, it will not
# be re-run unless the data, fitting options, or underlying Stan code change
stanDx(bpo)

Iterations: 2000 on each of 4 chains, with 4000 posterior distribution samples saved

For each parameter, n_eff is a crude measure of effective sample size
and Rhat is the potential scale reduction factor on split chains
(at convergence, Rhat=1)


Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Rank-normalized split effective sample size satisfactory for all parameters.

The following parameters had rank-normalized split R-hat greater than 1.01:
  alpha[39]
Such high values indicate incomplete mixing and biased estimation.
You should consider regularizing your model with additional prior information or a more effective parameterization.

Processing complete.

EBFMI: 0.762 0.756 0.704 0.826 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.003      767      742
2   alpha[2] 1.005      703     1645
3   alpha[3] 1.004      687     1672
4   alpha[4] 1.004      610     1429
5   alpha[5] 1.005      544     1336
6   alpha[6] 1.006      500     1384
7   alpha[7] 1.006      481     1306
8   alpha[8] 1.005      478     1388
9   alpha[9] 1.005      497      991
10 alpha[10] 1.006      480     1115
11 alpha[11] 1.006      467      856
12 alpha[12] 1.006      457      914
13 alpha[13] 1.007      445      878
14 alpha[14] 1.008      436      882
15 alpha[15] 1.008      430      791
16 alpha[16] 1.008      427      805
17 alpha[17] 1.010      420      881
18 alpha[18] 1.011      409      907
19 alpha[19] 1.011      409      811
20 alpha[20] 1.012      403      861
21 alpha[21] 1.012      407      842
22 alpha[22] 1.012      404      856
23 alpha[23] 1.011      401      849
24 alpha[24] 1.012      397      823
25 alpha[25] 1.012      380      745
26 alpha[26] 1.012      356      693
27 alpha[27] 1.013      348      778
28 alpha[28] 1.013      337      860
29 alpha[29] 1.014      322      847
30 alpha[30] 1.014      320      774
31 alpha[31] 1.014      316      842
32 alpha[32] 1.015      309      712
33 alpha[33] 1.014      316      796
34 alpha[34] 1.014      318      831
35 alpha[35] 1.014      309      931
36 alpha[36] 1.015      304      970
37 alpha[37] 1.014      303      802
38 alpha[38] 1.014      300      986
39 alpha[39] 1.015      313     1041
40 alpha[40] 1.014      319      993
41 alpha[41] 1.014      327     1074
42 alpha[42] 1.014      332      940
43 alpha[43] 1.014      336      978
44 alpha[44] 1.013      350     1061
45 alpha[45] 1.013      354     1023
46 alpha[46] 1.013      369      954
47 alpha[47] 1.013      387      914
48 alpha[48] 1.012      416      951
49 alpha[49] 1.011      429     1215
50 alpha[50] 1.008      510     1330
51 alpha[51] 1.009      518     1505
52 alpha[52] 1.007      606     1656
53 alpha[53] 1.006      713     1696
54 alpha[54] 1.007      741     1377
55 alpha[55] 1.005      771     1358
56 alpha[56] 1.005      770     1425
57 alpha[57] 1.004      839     1491
58 alpha[58] 1.004      835     1305
59 alpha[59] 1.002      882     1664
60 alpha[60] 1.001     1074     1887
61 alpha[61] 1.001     1348     2036
62   beta[1] 1.004      854     1487
63   beta[2] 1.001      770     1313
64   beta[3] 1.001     1660     2323
65   beta[4] 1.002     4376     3297
66   beta[5] 1.003      744     1374
67   beta[6] 1.005      868     1304
68   beta[7] 1.004      867     1307
69   beta[8] 1.010      839     1509
70   beta[9] 1.003      872     1533
71  beta[10] 1.006      803     1522
72  beta[11] 1.001     3596     2618
73  beta[12] 1.005     4026     3315
74  beta[13] 1.001     3742     3074
75  beta[14] 1.000     4351     3246
76  beta[15] 1.008      850     1520
77  beta[16] 1.000      940     1809
78  beta[17] 1.004      963     1574
79 sigmag[1] 1.001      771     1950

print(bpo, intercepts=TRUE)

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.044 for Intercepts

blrm(formula = twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) + 
    rcs(age, 4) * sex + cluster(uid), data = both, file = "bpo.rds")

	Mixed Calibration/ Discrimination Indexes	Discrimination Indexes	Rank Discrim. Indexes
Obs522	LOO log L-1748.52±24.29	g3.844 [3.283, 4.314]	C0.793 [0.785, 0.799]
Draws4000	LOO IC3497.04±48.58	gp0.435 [0.418, 0.449]	Dxy0.585 [0.57, 0.598]
Chains4	Effective p181.9±8.8	EV0.593 [0.55, 0.644]
Time5.5s	B0.148 [0.139, 0.159]	v11.516 [8.823, 14.93]
p17		vp0.148 [0.136, 0.16]
Cluster onuid
Clusters108
σγ1.8926 [1.5782, 2.281]

	Mean β	Median β	S.E.	Lower	Upper	P(β>0)	Symmetry
y≥7	-1.7002	-1.7734	4.3562	-10.1605	6.8135	0.3368	1.04
y≥9	-2.7386	-2.7925	4.2058	-10.9503	5.5078	0.2527	0.99
y≥10	-3.9516	-4.0247	4.1612	-11.8687	4.3084	0.1720	0.99
y≥11	-4.4055	-4.4803	4.1542	-12.2860	3.8922	0.1467	0.98
y≥13	-4.6041	-4.6660	4.1496	-12.3684	3.7260	0.1353	0.99
y≥14	-4.9639	-4.9829	4.1437	-12.9558	3.2785	0.1170	0.99
y≥15	-5.2778	-5.3490	4.1344	-13.4157	2.7981	0.1003	0.99
y≥16	-5.6672	-5.7187	4.1305	-13.5616	2.6037	0.0860	0.99
y≥17	-6.4817	-6.5643	4.1327	-14.3804	1.8335	0.0605	1.00
y≥18	-6.7377	-6.7990	4.1292	-14.5789	1.6319	0.0528	1.00
y≥19	-7.0264	-7.0706	4.1307	-15.2614	0.9845	0.0462	0.99
y≥20	-7.2207	-7.2760	4.1282	-14.9876	1.2145	0.0402	1.00
y≥21	-7.4023	-7.4656	4.1270	-15.0974	1.0784	0.0362	0.99
y≥22	-7.8229	-7.8747	4.1279	-15.6106	0.5303	0.0295	0.99
y≥23	-8.0925	-8.1531	4.1298	-16.1513	0.0098	0.0243	0.99
y≥24	-8.3768	-8.4289	4.1289	-16.5752	-0.4271	0.0217	0.99
y≥25	-8.6385	-8.6934	4.1266	-16.4834	-0.3560	0.0173	0.99
y≥26	-9.0356	-9.0950	4.1285	-16.9506	-0.8271	0.0132	0.99
y≥27	-9.3280	-9.3903	4.1287	-17.2768	-1.1810	0.0118	1.00
y≥28	-9.5719	-9.6338	4.1307	-17.4612	-1.4034	0.0107	1.00
y≥29	-9.8082	-9.8490	4.1334	-17.5704	-1.4282	0.0098	1.00
y≥30	-10.1140	-10.1702	4.1356	-17.9337	-1.7710	0.0078	1.00
y≥31	-10.4073	-10.4624	4.1387	-18.1932	-2.0589	0.0065	0.99
y≥32	-10.5260	-10.5803	4.1400	-18.3682	-2.2305	0.0060	0.99
y≥33	-10.8934	-10.9444	4.1417	-18.7100	-2.5744	0.0043	0.99
y≥34	-11.2066	-11.2709	4.1435	-18.9938	-2.8572	0.0032	0.98
y≥35	-11.4298	-11.5081	4.1466	-19.2641	-3.1457	0.0025	0.99
y≥36	-11.6772	-11.7521	4.1471	-19.4919	-3.4056	0.0018	1.00
y≥37	-11.9571	-12.0224	4.1499	-19.7771	-3.6524	0.0010	0.99
y≥38	-12.1852	-12.2339	4.1503	-20.0313	-3.9027	0.0010	1.00
y≥39	-12.4337	-12.5020	4.1528	-20.4083	-4.2526	0.0008	0.99
y≥40	-12.6190	-12.6872	4.1541	-20.5246	-4.3113	0.0008	0.99
y≥41	-12.8047	-12.8779	4.1572	-20.6270	-4.4020	0.0008	0.99
y≥42	-13.1304	-13.1955	4.1604	-21.0835	-4.8443	0.0005	0.99
y≥43	-13.3602	-13.4238	4.1622	-21.4213	-5.1536	0.0000	0.99
y≥44	-13.7010	-13.7739	4.1647	-21.7736	-5.4974	0.0000	0.99
y≥45	-14.0200	-14.1077	4.1677	-22.1482	-5.8790	0.0000	0.98
y≥46	-14.3192	-14.4020	4.1700	-22.3163	-6.0467	0.0000	0.98
y≥47	-14.7406	-14.8200	4.1719	-22.8046	-6.4998	0.0000	0.99
y≥48	-15.0323	-15.0997	4.1724	-23.2039	-6.9180	0.0000	0.98
y≥49	-15.4064	-15.4866	4.1750	-23.6554	-7.3281	0.0000	0.98
y≥50	-15.7244	-15.7922	4.1772	-23.9265	-7.5581	0.0000	0.99
y≥51	-16.2551	-16.3167	4.1806	-24.3386	-7.9772	0.0000	0.98
y≥52	-16.6199	-16.6845	4.1837	-24.6298	-8.2474	0.0000	0.97
y≥53	-17.0663	-17.1230	4.1859	-25.1847	-8.7691	0.0000	0.98
y≥54	-17.5688	-17.6424	4.1901	-25.5397	-9.0972	0.0000	0.98
y≥55	-17.9802	-18.0497	4.1911	-26.0094	-9.5808	0.0000	0.98
y≥56	-18.2327	-18.3182	4.1907	-26.2549	-9.8063	0.0000	0.98
y≥57	-18.7055	-18.7569	4.1924	-26.7439	-10.2821	0.0000	0.97
y≥58	-19.2710	-19.3209	4.1941	-27.2645	-10.8867	0.0000	0.98
y≥59	-19.6241	-19.6879	4.1944	-27.6708	-11.1285	0.0000	0.98
y≥60	-19.9547	-19.9935	4.1960	-28.4045	-11.9046	0.0000	0.96
y≥61	-20.6544	-20.7051	4.2039	-28.7605	-12.2844	0.0000	0.97
y≥62	-21.0269	-21.0729	4.2062	-28.9758	-12.5107	0.0000	0.96
y≥63	-21.4328	-21.4859	4.2066	-29.6942	-13.2036	0.0000	0.96
y≥64	-21.5732	-21.6422	4.2102	-29.7582	-13.2620	0.0000	0.96
y≥65	-22.3016	-22.3586	4.2174	-30.5591	-13.9429	0.0000	0.98
y≥66	-22.6756	-22.7163	4.2189	-31.0016	-14.3725	0.0000	0.96
y≥67	-23.0863	-23.1313	4.2243	-31.3163	-14.7646	0.0000	0.97
y≥68	-23.8654	-23.9219	4.2451	-32.3933	-15.6250	0.0000	0.97
y≥71	-24.7268	-24.7299	4.2890	-33.3965	-16.6165	0.0000	0.99
treat=5000U	0.0815	0.0810	0.7121	-1.2903	1.4693	0.5498	0.97
treat=Placebo	2.3377	2.3297	0.7336	0.9406	3.8091	0.9995	1.00
week	0.1218	0.1229	0.0784	-0.0285	0.2733	0.9372	0.93
week'	0.1945	0.1938	0.0859	0.0306	0.3638	0.9890	1.07
twstrs0	0.2291	0.2296	0.0509	0.1287	0.3274	1.0000	1.01
twstrs0'	0.1304	0.1307	0.0628	0.0064	0.2521	0.9815	1.00
age	-0.0159	-0.0173	0.0805	-0.1708	0.1419	0.4215	1.03
age'	0.1928	0.1898	0.2224	-0.2297	0.6489	0.8125	0.99
age''	-1.0650	-1.0468	0.8814	-2.8618	0.6524	0.1052	0.98
sex=M	5.0797	4.9504	6.2514	-7.3130	17.2240	0.7908	1.06
treat=5000U × week	0.0541	0.0528	0.1095	-0.1762	0.2610	0.6912	1.01
treat=Placebo × week	-0.0524	-0.0522	0.1109	-0.2762	0.1562	0.3123	1.02
treat=5000U × week'	-0.1678	-0.1685	0.1203	-0.3961	0.0849	0.0867	0.98
treat=Placebo × week'	-0.1436	-0.1446	0.1214	-0.3784	0.0922	0.1200	1.00
age × sex=M	-0.1108	-0.1060	0.1489	-0.4129	0.1726	0.2285	0.94
age' × sex=M	0.1545	0.1518	0.4159	-0.6854	0.9329	0.6450	1.04
age'' × sex=M	0.0283	0.0034	1.6239	-3.1358	3.1332	0.5022	1.00

a <- anova(bpo)
a

Relative Explained Variation for twstrs. Approximate total model Wald χ² used in denominators of REV:269.1 [224.1, 346.3].

	REV	Lower	Upper	d.f.
treat (Factor+Higher Order Factors)	0.125	0.066	0.204	6
All Interactions	0.087	0.037	0.162	4
week (Factor+Higher Order Factors)	0.582	0.451	0.688	6
All Interactions	0.087	0.037	0.162	4
Nonlinear (Factor+Higher Order Factors)	0.021	0.001	0.064	3
twstrs0	0.635	0.498	0.707	2
Nonlinear	0.016	0.000	0.047	1
age (Factor+Higher Order Factors)	0.023	0.007	0.084	6
All Interactions	0.015	0.001	0.057	3
Nonlinear (Factor+Higher Order Factors)	0.020	0.002	0.069	4
sex (Factor+Higher Order Factors)	0.019	0.003	0.070	4
All Interactions	0.015	0.001	0.057	3
treat × week (Factor+Higher Order Factors)	0.087	0.037	0.162	4
Nonlinear	0.008	0.000	0.037	2
Nonlinear Interaction : f(A,B) vs. AB	0.008	0.000	0.037	2
age × sex (Factor+Higher Order Factors)	0.015	0.001	0.057	3
Nonlinear	0.013	0.000	0.049	2
Nonlinear Interaction : f(A,B) vs. AB	0.013	0.000	0.049	2
TOTAL NONLINEAR	0.053	0.025	0.130	8
TOTAL INTERACTION	0.102	0.047	0.189	7
TOTAL NONLINEAR + INTERACTION	0.133	0.088	0.245	11
TOTAL	1.000	1.000	1.000	17

plot(a)

• Show the final graphic (high dose:placebo contrast as function of time
• Intervals are 0.95 highest posterior density intervals
• \(y\)-axis: log-odds ratio

X

wks <- c(2,4,8,12,16)
k <- contrast(bpo, list(week=wks, treat='10000U'),
                   list(week=wks, treat='Placebo'),
              cnames=paste('Week', wks))
k

           week   Contrast      S.E.      Lower      Upper P(Contrast>0)
1  Week 2     2 -2.2329821 0.6045026 -3.4779977 -1.0997215        0.0000
2  Week 4     4 -2.1282800 0.5395822 -3.1836898 -1.0892759        0.0000
3  Week 8     8 -1.7752623 0.5926660 -2.9287838 -0.6109096        0.0020
4* Week 12   12 -0.8477916 0.5417945 -1.8838867  0.1892500        0.0528
5* Week 16   16  0.2232923 0.6242989 -0.9726945  1.4339875        0.6435

Redundant contrasts are denoted by *

Intervals are 0.95 highest posterior density intervals
Contrast is the posterior mean 

plot(k)

k <- as.data.frame(k[c('week', 'Contrast', 'Lower', 'Upper')])
ggplot(k, aes(x=week, y=Contrast)) + geom_point() +
  geom_line() + ylab('High Dose - Placebo') +
  geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0)

For each posterior draw compute the difference in means and get an exact (to within simulation error) 0.95 highest posterior density intervals for these differences.

M <- Mean(bpo)   # create R function that computes mean Y from X*beta
k <- contrast(bpo, list(week=wks, treat='10000U'),
                   list(week=wks, treat='Placebo'),
              fun=M, cnames=paste('Week', wks))
plot(k, which='diff') + theme(legend.position='bottom')

f <- function(x) {
  hpd <- HPDint(x, prob=0.95)   # is in rmsb
  r <- c(mean(x), median(x), hpd)
  names(r) <- c('Mean', 'Median', 'Lower', 'Upper')
  r
}
w    <- as.data.frame(t(apply(k$esta - k$estb, 2, f)))
week <- as.numeric(sub('Week ', '', rownames(w)))
ggplot(w, aes(x=week, y=Mean)) + geom_point() +
  geom_line() + ylab('High Dose - Placebo') +
  geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0) +
  scale_y_continuous(breaks=c(-8, -4, 0, 4))

7.8.5 Bayesian Markov Semiparametric Model

• First-order Markov model
• Serial correlation induced by Markov model is similar to AR(1) which we already know fits these data
• Markov model is more likely to fit the data than the random effects model, which induces a compound symmetry correlation structure
• Models state transitions
• PO model at each visit, with Y from previous visit conditioned upon just like any covariate
• Need to uncondition (marginalize) on previous Y to get the time-response profile we usually need
• Semiparametric model is especially attractive because one can easily “uncondition” a discrete Y model, and the distribution of Y for control subjects can be any shape
• Let measurement times be \(t_{1}, t_{2}, \dots, t_{m}\), and the measurement for a subject at time \(t\) be denoted \(Y(t)\)
• First-order Markov model:

Y

\[\begin{array}{ccc} \Pr(Y(t_{i}) \geq y | X, Y(t_{i-1})) &=& \mathrm{expit}(\alpha_{y} + X\beta\\ &+& g(Y(t_{i-1}), t_{i}, t_{i} - t_{i-1})) \end{array}\]

• \(g\) involves any number of regression coefficients for a main effect of \(t\), the main effect of time gap \(t_{i} - t_{i-1}\) if this is not collinear with absolute time, a main effect of the previous state, and interactions between these
• Examples of how the previous state may be modeled in \(g\):
  • linear in numeric codes for \(Y\)
  • spline function in same
  • discontinuous bi-linear relationship where there is a slope for in-hospital outcome severity, a separate slope for outpatient outcome severity, and an intercept jump at the transition from inpatient to outpatient (or vice versa)
• Markov model is quite flexible in handling time trends and serial correlation patterns
• Can allow for irregular measurement times:
  hbiostat.org/stat/irreg.html

Fit the model and run standard Stan diagnostics.

# Create a new variable to hold previous value of Y for the subject
# For week 2, previous value is the baseline value
setDT(both, key=c('uid', 'week'))
both[, ptwstrs := shift(twstrs), by=uid]
both[week == 2, ptwstrs := twstrs0]
dd <- datadist(both)
bmark <- blrm(twstrs ~  treat * rcs(week, 3) + rcs(ptwstrs, 4) +
                        rcs(age, 4) * sex,
              data=both, file='bmark.rds')

Running MCMC with 4 chains, at most 9 in parallel...

Chain 1 finished in 2.4 seconds.
Chain 4 finished in 2.2 seconds.
Chain 3 finished in 2.4 seconds.
Chain 2 finished in 2.5 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.4 seconds.
Total execution time: 2.6 seconds.

# When adding partial PO terms for week and ptwstrs, z=-1.8, 5.04
stanDx(bmark)

Iterations: 2000 on each of 4 chains, with 4000 posterior distribution samples saved

For each parameter, n_eff is a crude measure of effective sample size
and Rhat is the potential scale reduction factor on split chains
(at convergence, Rhat=1)


Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Rank-normalized split effective sample size satisfactory for all parameters.

Rank-normalized split R-hat values satisfactory for all parameters.

Processing complete, no problems detected.

EBFMI: 0.982 1.028 0.905 0.943 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.001     2681     1972
2   alpha[2] 1.001     2367     2669
3   alpha[3] 1.001     2272     2402
4   alpha[4] 1.001     2164     2571
5   alpha[5] 1.001     1995     2292
6   alpha[6] 1.002     1898     2800
7   alpha[7] 1.002     1860     2553
8   alpha[8] 1.001     1863     2462
9   alpha[9] 1.001     2020     2939
10 alpha[10] 1.001     2042     2763
11 alpha[11] 1.001     2055     2673
12 alpha[12] 1.001     2015     2768
13 alpha[13] 1.001     2039     2577
14 alpha[14] 1.002     2277     2935
15 alpha[15] 1.002     2370     3220
16 alpha[16] 1.002     2396     2833
17 alpha[17] 1.002     2544     2777
18 alpha[18] 1.001     2644     3005
19 alpha[19] 1.001     2792     2982
20 alpha[20] 1.001     2899     3070
21 alpha[21] 1.001     3058     2995
22 alpha[22] 1.001     3349     3070
23 alpha[23] 1.001     3632     3424
24 alpha[24] 1.002     3633     3142
25 alpha[25] 1.000     3981     2813
26 alpha[26] 1.002     4298     3197
27 alpha[27] 1.001     4509     3065
28 alpha[28] 1.001     4973     3355
29 alpha[29] 1.001     5333     3798
30 alpha[30] 1.001     5602     3446
31 alpha[31] 1.001     5929     3338
32 alpha[32] 1.000     6123     3412
33 alpha[33] 1.000     6054     3211
34 alpha[34] 1.000     5821     3054
35 alpha[35] 1.000     5721     3058
36 alpha[36] 1.000     5820     3094
37 alpha[37] 1.001     5601     3046
38 alpha[38] 1.002     5503     3367
39 alpha[39] 1.001     5597     3358
40 alpha[40] 1.000     5199     3152
41 alpha[41] 1.000     4769     3010
42 alpha[42] 1.001     4647     2729
43 alpha[43] 1.001     4377     3393
44 alpha[44] 1.000     4254     3128
45 alpha[45] 1.001     4179     3324
46 alpha[46] 1.000     4054     3012
47 alpha[47] 1.001     4313     3177
48 alpha[48] 1.001     4157     3065
49 alpha[49] 1.000     4066     3368
50 alpha[50] 1.001     3951     3401
51 alpha[51] 1.001     4046     3434
52 alpha[52] 1.001     4019     3472
53 alpha[53] 1.001     3908     3161
54 alpha[54] 1.001     3787     3004
55 alpha[55] 1.001     3824     3337
56 alpha[56] 1.001     3854     3233
57 alpha[57] 1.002     3570     3094
58 alpha[58] 1.001     3730     3041
59 alpha[59] 1.001     3905     3135
60 alpha[60] 1.001     4164     3272
61 alpha[61] 1.001     4531     3136
62   beta[1] 1.002     6978     2874
63   beta[2] 1.000     7070     2670
64   beta[3] 1.002     5088     3282
65   beta[4] 1.001     6620     2993
66   beta[5] 1.002     2246     2747
67   beta[6] 1.000     4494     3444
68   beta[7] 1.001     5813     3172
69   beta[8] 1.001     7584     2876
70   beta[9] 1.002     6939     3039
71  beta[10] 1.001     7208     3086
72  beta[11] 1.002     6929     3246
73  beta[12] 1.001     7752     2965
74  beta[13] 1.003     6007     3013
75  beta[14] 1.001     6862     2889
76  beta[15] 1.001     7634     2949
77  beta[16] 1.001     8525     2859
78  beta[17] 1.002     8079     2882
79  beta[18] 1.000     7843     2923

stanDxplot(bmark)

Note that posterior sampling is much more efficient without random effects.

bmark

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.044 for Intercepts

blrm(formula = twstrs ~ treat * rcs(week, 3) + rcs(ptwstrs, 4) + 
    rcs(age, 4) * sex, data = both, file = "bmark.rds")

Frequencies of Missing Values Due to Each Variable

 twstrs   treat    week ptwstrs     age     sex 
      0       0       0       5       0       0 

	Mixed Calibration/ Discrimination Indexes	Discrimination Indexes	Rank Discrim. Indexes
Obs517	LOO log L-1786.5±22.38	g3.276 [3.006, 3.597]	C0.828 [0.826, 0.831]
Draws4000	LOO IC3572.99±44.76	gp0.416 [0.404, 0.429]	Dxy0.656 [0.651, 0.661]
Chains4	Effective p90.4±4.94	EV0.533 [0.498, 0.571]
Time3.2s	B0.117 [0.114, 0.121]	v8.452 [6.96, 10.022]
p18		vp0.133 [0.124, 0.143]

	Mode β	Mean β	Median β	S.E.	Lower	Upper	P(β>0)	Symmetry
treat=5000U	0.2211	0.2085	0.2184	0.5904	-0.9714	1.3388	0.6412	0.98
treat=Placebo	1.8312	1.8284	1.8317	0.5890	0.6600	2.9640	1.0000	0.98
week	0.4864	0.4886	0.4877	0.0836	0.3305	0.6549	1.0000	1.05
week'	-0.2878	-0.2891	-0.2894	0.0892	-0.4466	-0.0971	0.0003	1.00
ptwstrs	0.1997	0.2014	0.2011	0.0274	0.1494	0.2551	1.0000	1.05
ptwstrs'	-0.0621	-0.0658	-0.0658	0.0644	-0.1908	0.0656	0.1460	1.04
ptwstrs''	0.5327	0.5523	0.5504	0.2538	0.0757	1.0785	0.9830	0.98
age	-0.0295	-0.0287	-0.0289	0.0319	-0.0901	0.0328	0.1820	1.03
age'	0.1236	0.1222	0.1237	0.0887	-0.0547	0.2945	0.9130	1.02
age''	-0.5067	-0.5036	-0.5043	0.3512	-1.2231	0.1524	0.0720	0.99
sex=M	-0.4614	-0.3673	-0.3579	2.3818	-4.7272	4.4595	0.4368	0.97
treat=5000U × week	-0.0341	-0.0317	-0.0339	0.1134	-0.2506	0.1925	0.3908	1.00
treat=Placebo × week	-0.2719	-0.2717	-0.2713	0.1136	-0.5031	-0.0551	0.0053	1.02
treat=5000U × week'	-0.0340	-0.0371	-0.0356	0.1218	-0.2671	0.2071	0.3895	0.99
treat=Placebo × week'	0.1197	0.1183	0.1185	0.1228	-0.1133	0.3609	0.8305	0.98
age × sex=M	0.0112	0.0089	0.0087	0.0572	-0.0961	0.1242	0.5618	1.03
age' × sex=M	-0.0509	-0.0457	-0.0468	0.1609	-0.3744	0.2507	0.3895	1.02
age'' × sex=M	0.2613	0.2444	0.2496	0.6265	-0.9803	1.4537	0.6575	1.01

a <- anova(bmark)
a

Relative Explained Variation for twstrs. Approximate total model Wald χ² used in denominators of REV:448 [375.5, 540.3].

	REV	Lower	Upper	d.f.
treat (Factor+Higher Order Factors)	0.053	0.027	0.104	6
All Interactions	0.050	0.017	0.093	4
week (Factor+Higher Order Factors)	0.295	0.212	0.376	6
All Interactions	0.050	0.017	0.093	4
Nonlinear (Factor+Higher Order Factors)	0.062	0.026	0.109	3
ptwstrs	0.948	0.871	0.960	3
Nonlinear	0.042	0.010	0.085	2
age (Factor+Higher Order Factors)	0.009	0.002	0.039	6
All Interactions	0.001	0.000	0.019	3
Nonlinear (Factor+Higher Order Factors)	0.005	0.001	0.032	4
sex (Factor+Higher Order Factors)	0.001	0.001	0.022	4
All Interactions	0.001	0.000	0.019	3
treat × week (Factor+Higher Order Factors)	0.050	0.017	0.093	4
Nonlinear	0.004	0.000	0.023	2
Nonlinear Interaction : f(A,B) vs. AB	0.004	0.000	0.023	2
age × sex (Factor+Higher Order Factors)	0.001	0.000	0.019	3
Nonlinear	0.001	0.000	0.015	2
Nonlinear Interaction : f(A,B) vs. AB	0.001	0.000	0.015	2
TOTAL NONLINEAR	0.106	0.071	0.175	9
TOTAL INTERACTION	0.052	0.024	0.105	7
TOTAL NONLINEAR + INTERACTION	0.144	0.104	0.232	12
TOTAL	1.000	1.000	1.000	18

plot(a)

Let’s add subject-level random effects to the model. Smallness of the standard deviation of the random effects provides support for the assumption of conditional independence that we like to make for Markov models and allows us to simplify the model by omitting random effects.

bmarkre <- blrm(twstrs ~  treat * rcs(week, 3) + rcs(ptwstrs, 4) +
                          rcs(age, 4) * sex + cluster(uid),
                data=both, file='bmarkre.rds')

Running MCMC with 4 chains, at most 9 in parallel...

Chain 4 finished in 3.3 seconds.
Chain 1 finished in 3.4 seconds.
Chain 2 finished in 3.4 seconds.
Chain 3 finished in 3.5 seconds.

All 4 chains finished successfully.
Mean chain execution time: 3.4 seconds.
Total execution time: 3.7 seconds.

stanDx(bmarkre)

Iterations: 2000 on each of 4 chains, with 4000 posterior distribution samples saved

For each parameter, n_eff is a crude measure of effective sample size
and Rhat is the potential scale reduction factor on split chains
(at convergence, Rhat=1)


Checking sampler transitions for divergences.
6 of 4000 (0.15%) transitions ended with a divergence.
These divergent transitions indicate that HMC is not fully able to explore the posterior distribution.
Try increasing adapt delta closer to 1.
If this doesn't remove all divergences, try to reparameterize the model.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Rank-normalized split effective sample size satisfactory for all parameters.

Rank-normalized split R-hat values satisfactory for all parameters.

Processing complete.
Divergent samples: 0 3 3 0 

EBFMI: 0.991 0.969 0.986 0.915 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.001     1574     1846
2   alpha[2] 1.003     1182     1848
3   alpha[3] 1.004     1226     1844
4   alpha[4] 1.004     1125     1997
5   alpha[5] 1.004      944     1855
6   alpha[6] 1.004      917     1530
7   alpha[7] 1.006      866     1536
8   alpha[8] 1.007      862     1625
9   alpha[9] 1.007     1019     1731
10 alpha[10] 1.008      979     1900
11 alpha[11] 1.008      902     1775
12 alpha[12] 1.008      838     2094
13 alpha[13] 1.009      790     2085
14 alpha[14] 1.008      838     1796
15 alpha[15] 1.009      806     2240
16 alpha[16] 1.008      922     2289
17 alpha[17] 1.007     1069     2060
18 alpha[18] 1.006     1255     1789
19 alpha[19] 1.005     1323     1601
20 alpha[20] 1.004     1393     1722
21 alpha[21] 1.003     1495     1578
22 alpha[22] 1.002     1644     2043
23 alpha[23] 1.002     1899     2147
24 alpha[24] 1.002     1889     2194
25 alpha[25] 1.002     2361     2296
26 alpha[26] 1.001     2532     2123
27 alpha[27] 1.002     2894     2750
28 alpha[28] 1.002     3144     2796
29 alpha[29] 1.002     3149     2790
30 alpha[30] 1.001     3259     3005
31 alpha[31] 1.000     3462     2579
32 alpha[32] 1.000     3408     3118
33 alpha[33] 1.001     3219     3092
34 alpha[34] 1.002     3110     2693
35 alpha[35] 1.001     3106     2814
36 alpha[36] 1.001     3022     2835
37 alpha[37] 1.001     2826     2715
38 alpha[38] 1.001     2427     2649
39 alpha[39] 1.003     2336     2412
40 alpha[40] 1.002     2263     2672
41 alpha[41] 1.003     2182     2623
42 alpha[42] 1.002     2166     2578
43 alpha[43] 1.002     2165     2644
44 alpha[44] 1.003     1956     2537
45 alpha[45] 1.002     1934     2520
46 alpha[46] 1.002     1843     2224
47 alpha[47] 1.001     1951     2473
48 alpha[48] 1.001     1966     2156
49 alpha[49] 1.001     1940     2314
50 alpha[50] 1.001     2028     2162
51 alpha[51] 1.000     1976     1935
52 alpha[52] 1.001     1961     1796
53 alpha[53] 1.002     2041     2171
54 alpha[54] 1.001     1995     1872
55 alpha[55] 1.002     2008     1958
56 alpha[56] 1.000     1946     1431
57 alpha[57] 1.002     2063     2200
58 alpha[58] 1.001     1917     1469
59 alpha[59] 1.001     2032     1889
60 alpha[60] 1.003     2132     2313
61 alpha[61] 1.001     2487     2585
62   beta[1] 1.003     4205     2579
63   beta[2] 1.000     4471     2812
64   beta[3] 1.000     2705     2299
65   beta[4] 1.001     4294     2895
66   beta[5] 1.003     1418     2267
67   beta[6] 1.002     2617     3058
68   beta[7] 1.001     2874     2181
69   beta[8] 1.000     3763     2656
70   beta[9] 1.000     4394     3010
71  beta[10] 1.000     3382     2369
72  beta[11] 1.002     3843     2357
73  beta[12] 1.001     4436     2611
74  beta[13] 1.000     3859     2830
75  beta[14] 1.001     4835     2877
76  beta[15] 1.002     5072     3362
77  beta[16] 1.003     3987     2545
78  beta[17] 1.001     3852     2414
79  beta[18] 1.001     4498     2801
80 sigmag[1] 1.002     1112     1093

bmarkre

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.044 for Intercepts

blrm(formula = twstrs ~ treat * rcs(week, 3) + rcs(ptwstrs, 4) + 
    rcs(age, 4) * sex + cluster(uid), data = both, file = "bmarkre.rds")

Frequencies of Missing Values Due to Each Variable

      twstrs        treat         week      ptwstrs          age          sex 
           0            0            0            5            0            0 
cluster(uid) 
           0 

	Mixed Calibration/ Discrimination Indexes	Discrimination Indexes	Rank Discrim. Indexes
Obs517	LOO log L-1786.09±22.2	g3.263 [2.994, 3.583]	C0.828 [0.825, 0.83]
Draws4000	LOO IC3572.17±44.4	gp0.416 [0.402, 0.429]	Dxy0.655 [0.65, 0.661]
Chains4	Effective p92.63±4.74	EV0.532 [0.493, 0.572]
Time4.4s	B0.117 [0.113, 0.121]	v8.378 [7.024, 10.043]
p18		vp0.133 [0.123, 0.143]
Cluster onuid
Clusters108
σγ0.1041 [3e-04, 0.3307]

	Mean β	Median β	S.E.	Lower	Upper	P(β>0)	Symmetry
treat=5000U	0.2009	0.1960	0.5765	-0.9476	1.3048	0.6355	1.01
treat=Placebo	1.8223	1.8180	0.5687	0.7796	2.9748	0.9988	1.02
week	0.4836	0.4834	0.0836	0.3176	0.6464	1.0000	1.03
week'	-0.2838	-0.2848	0.0879	-0.4534	-0.1079	0.0003	1.00
ptwstrs	0.2003	0.2002	0.0264	0.1526	0.2551	1.0000	1.00
ptwstrs'	-0.0658	-0.0654	0.0633	-0.2010	0.0467	0.1530	1.03
ptwstrs''	0.5492	0.5467	0.2530	0.0611	1.0336	0.9842	0.97
age	-0.0294	-0.0297	0.0323	-0.0903	0.0358	0.1767	1.00
age'	0.1234	0.1245	0.0903	-0.0543	0.2991	0.9158	1.02
age''	-0.5075	-0.5068	0.3571	-1.1996	0.1933	0.0737	0.98
sex=M	-0.4107	-0.3876	2.4665	-5.1849	4.4356	0.4412	1.02
treat=5000U × week	-0.0288	-0.0305	0.1128	-0.2377	0.2051	0.3985	0.99
treat=Placebo × week	-0.2693	-0.2671	0.1119	-0.4823	-0.0532	0.0068	0.98
treat=5000U × week'	-0.0399	-0.0390	0.1206	-0.2933	0.1829	0.3760	1.00
treat=Placebo × week'	0.1164	0.1158	0.1201	-0.0960	0.3639	0.8350	1.05
age × sex=M	0.0097	0.0090	0.0593	-0.0999	0.1321	0.5620	0.97
age' × sex=M	-0.0464	-0.0448	0.1667	-0.3667	0.2857	0.3910	1.01
age'' × sex=M	0.2468	0.2440	0.6457	-1.0782	1.4671	0.6488	1.02

The random effects SD is only 0.11 on the logit scale. Also, the standard deviations of all the regression parameter posterior distributions are virtually unchanged with the addition of random effects:

plot(sqrt(diag(vcov(bmark))), sqrt(diag(vcov(bmarkre))),
     xlab='Posterior SDs in Conditional Independence Markov Model',
     ylab='Posterior SDs in Random Effects Markov Model')
abline(a=0, b=1, col=gray(0.85))

So we will use the model omitting random effects.

Show the partial effects of all the predictors, including the effect of the previous measurement of TWSTRS. Also compute high dose:placebo treatment contrasts on these conditional estimates.

ggplot(Predict(bmark))

ggplot(Predict(bmark, week, treat))

k <- contrast(bmark, list(week=wks, treat='10000U'),
                     list(week=wks, treat='Placebo'),
              cnames=paste('Week', wks))
k

           week   Contrast      S.E.      Lower      Upper P(Contrast>0)
1  Week 2     2 -1.2850122 0.3969281 -2.0786432 -0.5439447        0.0000
2  Week 4     4 -0.7415990 0.2674335 -1.2895859 -0.2339872        0.0030
3  Week 8     8  0.2269131 0.3510888 -0.4650531  0.9082032        0.7515
4* Week 12   12  0.7221680 0.2543682  0.2167805  1.2176876        0.9980
5* Week 16   16  1.0991087 0.3952105  0.3638820  1.9214141        0.9955

Redundant contrasts are denoted by *

Intervals are 0.95 highest posterior density intervals
Contrast is the posterior mean 

plot(k)

k <- as.data.frame(k[c('week', 'Contrast', 'Lower', 'Upper')])
ggplot(k, aes(x=week, y=Contrast)) + geom_point() +
  geom_line() + ylab('High Dose - Placebo') +
  geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0)

Using posterior means for parameter values, compute the probability that at a given week twstrs will be \(\geq 40\) when at the previous visit it was 40. Also show the conditional mean twstrs when it was 40 at the previous visit.

ex <- ExProb(bmark)
ex40 <- function(lp, ...) ex(lp, y=40, ...)
ggplot(Predict(bmark, week, treat, ptwstrs=40, fun=ex40))

ggplot(Predict(bmark, week, treat, ptwstrs=40, fun=Mean(bmark)))

• Semiparametric models provide not only estimates of tendencies of Y but also estimate the whole distribution of Y
• Estimate the entire conditional distribution of Y at week 12 for high-dose patients having TWSTRS=42 at week 8
• Other covariates set to median/mode
• Use posterior mean of all the cell probabilities
• Also show pointwise 0.95 highest posterior density intervals
• To roughly approximate simultaneous confidence bands make the pointwise limits sum to 1 like the posterior means do

Z

# Get median/mode for covariates including ptwstrs (TWSTRS in previous visit)
d <- gendata(bmark)
d

   treat week ptwstrs age sex
1 10000U    8      42  56   F

d$week <- 12
p <- predict(bmark, d, type='fitted.ind')   # defaults to posterior means
yvals <- as.numeric(sub('twstrs=', '', p$y))
lo <- p$Lower / sum(p$Lower)
hi <- p$Upper / sum(p$Upper)
plot(yvals, p$Mean, type='l', xlab='TWSTRS', ylab='',
     ylim=range(c(lo, hi)))
lines(yvals, lo, col=gray(0.8))
lines(yvals, hi, col=gray(0.8))

• Repeat this showing the variation over 5 posterior draws

A

p <- predict(bmark, d, type='fitted.ind', posterior.summary='all')
cols <- adjustcolor(1 : 10, 0.7)
for(i in 1 : 5) {
  if(i == 1) plot(yvals, p[i, 1, ], type='l', col=cols[1], xlab='TWSTRS', ylab='')
  else lines(yvals, p[i, 1, ], col=cols[i])
}

• Turn to marginalized (unconditional on previous twstrs) quantities
• Capitalize on PO model being a multinomial model, just with PO restrictions
• Manipulations of conditional probabilities to get the unconditional probability that twstrs=y doesn’t need to know about PO
• Compute all cell probabilities and use the law of total probability recursively \[\Pr(Y_{t} = y | X) = \sum_{j=1}^{k} \Pr(Y_{t} = y | X, Y_{t-1} = j) \Pr(Y_{t-1} = j | X)\]
• predict.blrm method with type='fitted.ind' computes the needed conditional cell probabilities, optionally for all posterior draws at once
• Easy to get highest posterior density intervals for derived parameters such as unconditional probabilities or unconditional means
• Hmisc package soprobMarkovOrdm function (in version 4.6) computes an array of all the state occupancy probabilities for all the posterior draws

B

# Baseline twstrs to 42 in d
# For each dose, get all the posterior draws for all state occupancy
# probabilities for all visit
ylev <- sort(unique(both$twstrs))
tlev <- c('Placebo', '10000U')
R <- list()
for(trt in tlev) {   # separately by treatment
  d$treat <- trt
  u <- soprobMarkovOrdm(bmark, d, wks, ylev,
                        tvarname='week', pvarname='ptwstrs')
  R[[trt]] <- u
}
dim(R[[1]])    # posterior draws x times x distinct twstrs values

[1] 4000    5   62

# For each posterior draw, treatment, and week compute the mean TWSTRS
# Then compute posterior mean of means, and HPD interval
Rmean <- Rmeans <- list()
for(trt in tlev) {
  r <- R[[trt]]
  # Mean Y at each week and posterior draw (mean from a discrete distribution)
  m <- apply(r, 1:2, function(x) sum(ylev * x))
  Rmeans[[trt]] <- m
  # Posterior mean and median and HPD interval over draws
  u <- apply(m, 2, f)   # f defined above
  u <- rbind(week=as.numeric(colnames(u)), u)
  Rmean[[trt]] <- u
}
r <- lapply(Rmean, function(x) as.data.frame(t(x)))
for(trt in tlev) r[[trt]]$treat <- trt
r <- do.call(rbind, r)
ggplot(r, aes(x=week, y=Mean, color=treat)) + geom_line() +
  geom_ribbon(aes(ymin=Lower, ymax=Upper), alpha=0.2, linetype=0)

• Use the same posterior draws of unconditional probabilities of all values of TWSTRS to get the posterior distribution of differences in mean TWSTRS between high and low dose

C

Dif <- Rmeans$`10000U` - Rmeans$Placebo
dif <- as.data.frame(t(apply(Dif, 2, f)))
dif$week <- as.numeric(rownames(dif))
ggplot(dif, aes(x=week, y=Mean)) + geom_line() +
  geom_ribbon(aes(ymin=Lower, ymax=Upper), alpha=0.2, linetype=0) +
  ylab('High Dose - Placebo TWSTRS')

• Get posterior mean of all cell probabilities estimates at week 12
• Distribution of TWSTRS conditional high dose, median age, mode sex
• Not conditional on week 8 value

D

p <- R$`10000U`[, '12', ]   # 4000 x 62
pmean <- apply(p, 2, mean)
yvals <- as.numeric(names(pmean))
plot(yvals, pmean, type='l', xlab='TWSTRS', ylab='')

7.8.6 Frequentist Markov Semiparametric Model

Fitting a frequentist MOST model will allow us to take advantage of the orm functions 2-parameter random effects contribution, which was designed to improve fit to correlation structure for Markov-1 models.

both[, mix2 := ifelse(week == 2, 0, 1)]
f <- orm(twstrs ~  treat * rcs(week, 3) + rcs(ptwstrs, 4) +
                    rcs(age, 4) * sex + cluster(uid) + mix_re(mix2),
         data=both, trace=1, x=TRUE)

NR iteration:1  -2LL:4041.5670  Max |gradient|:2902.772  Max |change in parameters|:14.9951
NR iteration:2  -2LL:3683.9746  Max |gradient|:1970.889  Max |change in parameters|:9.546635
NR iteration:3  -2LL:3556.3319  Max |gradient|:1814.975  Max |change in parameters|:7.043506
NR iteration:4  -2LL:3553.6259  Max |gradient|:1901.59  Max |change in parameters|:4.070974
NR iteration:5  -2LL:3488.7724  Max |gradient|:1754.493  Max |change in parameters|:2.797338
NR iteration:6  -2LL:3435.9552  Max |gradient|:1216.329  Max |change in parameters|:0.9311894
NR iteration:7  -2LL:3394.2030  Max |gradient|:237.9145  Max |change in parameters|:0.5766449
NR iteration:8  -2LL:3392.0287  Max |gradient|:2.672245  Max |change in parameters|:0.03623263
NR iteration:9  -2LL:3392.0255  Max |gradient|:0.004129342  Max |change in parameters|:0.0001897399
NR iteration:10  -2LL:3392.0255  Max |gradient|:6.029551e-07  Max |change in parameters|:2.485188e-08
NR iteration:1  -2LL:4041.5670  Max |gradient|:5.933032e-13  Max |change in parameters|:2.288985e-14
nAGQ: 7  outer iter: 3  sigma: 1.257209  sigma2: 1.115926  max|grad|: 115.7443  max|delta|: 0.7554744 
nAGQ: 7  outer iter: 6  sigma: 2.020361  sigma2: 2.089614  max|grad|: 33.25373  max|delta|: 1.066224 
nAGQ: 7  outer iter: 9  sigma: 2.357603  sigma2: 2.519725  max|grad|: 11.80749  max|delta|: 0.09524813 
nAGQ: 7  outer iter: 12  sigma: 2.446416  sigma2: 2.629527  max|grad|: 7.779044  max|delta|: 0.0495836 
nAGQ: 7  outer iter: 15  sigma: 2.531955  sigma2: 2.737963  max|grad|: 4.535269  max|delta|: 0.03074439 
nAGQ: 7  outer iter: 18  sigma: 2.647757  sigma2: 2.899259  max|grad|: 0.8073045  max|delta|: 0.009915426 
nAGQ: 7  outer iter: 21  sigma: 2.652798  sigma2: 2.909892  max|grad|: 0.6091979  max|delta|: 0.008416825 
nAGQ: 7  outer iter: 24  sigma: 2.662308  sigma2: 2.938355  max|grad|: 0.1373392  max|delta|: 0.004755316 
nAGQ: 7  outer iter: 27  sigma: 2.649453  sigma2: 2.948284  max|grad|: 0.1060364  max|delta|: 0.002184838 
nAGQ: 7  outer iter: 30  sigma: 2.652174  sigma2: 2.951541  max|grad|: 0.04619817  max|delta|: 0.000363112 
nAGQ: 7  outer iter: 33  sigma: 2.654163  sigma2: 2.953963  max|grad|: 0.002952467  max|delta|: 0.0002620181 
nAGQ: 7  outer iter: 36  sigma: 2.654168  sigma2: 2.953989  max|grad|: 0.0001268465  max|delta|: 1.910208e-05 
nAGQ: 7  outer iter: 3  sigma: 1.24748  sigma2: -0.1544452  max|grad|: 13.59846  max|delta|: 0.1439604 
nAGQ: 7  outer iter: 6  sigma: 1.435999  sigma2: -0.2042124  max|grad|: 1.274882  max|delta|: 0.004703056 
nAGQ: 7  outer iter: 9  sigma: 1.434293  sigma2: -0.2052287  max|grad|: 0.01861679  max|delta|: 0.0001795966 
nAGQ: 7  outer iter: 12  sigma: 1.434105  sigma2: -0.2052022  max|grad|: 0.004508967  max|delta|: 1.150756e-05 
nAGQ: 7  outer iter: 15  sigma: 1.434104  sigma2: -0.205204  max|grad|: 3.394429e-05  max|delta|: 1.819564e-07 

deviance(f)

                 intercepts                intercepts+x 
                   4041.567                    3392.026 
  intercepts+random effects intercepts+x+random effects 
                   3674.128                    3382.387 

sigma1 large: substantial subject-level correlation at baseline (before ptwstrs is contributing anything), consistent with a real, non-trivial random effect there. sigma2 ≈ -0.2: small and negative, but not collapsed to the boundary the way the single-sigma model did. That’s exactly the pattern the additive design was meant to reveal — once the lag term is active, a plain constant-weight random effect would be adding back in some of the same correlation the lag term already explains, and a small negative sigma2 is the model’s way of trimming that slight over-induction rather than needing to zero the whole random effect out just to compensate.

## Study Questions

**Section 7.2**

1. When should one model the time-response profile using discrete time?
   
**Section 7.3**

1. What makes generalized least squares and mixed effect models
   relatively robust to non-completely-random dropouts?
1. What does the last observation carried forward method always violate?

**Section 7.4**

1. Which correlation structure do you expect to fit the data when there are rapid repetitions over a short time span?  When the follow-up time span is very long?

**Section 7.8**

1. What can go wrong if many correlation structures are tested in one dataset?
1. In a longitudinal intervention study, what is the most typical comparison of interest?  Is it best to borrow information in estimating this contrast?


::: {.cell}

:::






:::{#quarto-navigation-envelope .hidden}
[Regression Modeling Strategies]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyLXRpdGxl"}
[Regression Modeling Strategies]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1uYXZiYXItdGl0bGU="}
[<span class='chapter-number'>8</span>  <span class='chapter-title'>Case Study in Data Reduction</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1uZXh0"}
[<span class='chapter-number'>6</span>  <span class='chapter-title'>R Software</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1wcmV2"}
[Preface]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9pbmRleC5odG1sUHJlZmFjZQ=="}
[<span class='chapter-number'>1</span>  <span class='chapter-title'>Introduction</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9pbnRyby5odG1sPHNwYW4tY2xhc3M9J2NoYXB0ZXItbnVtYmVyJz4xPC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkludHJvZHVjdGlvbjwvc3Bhbj4="}
[<span class='chapter-number'>2</span>  <span class='chapter-title'>General Aspects of Fitting Regression Models</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9nZW5yZWcuaHRtbDxzcGFuLWNsYXNzPSdjaGFwdGVyLW51bWJlcic+Mjwvc3Bhbj4tLTxzcGFuLWNsYXNzPSdjaGFwdGVyLXRpdGxlJz5HZW5lcmFsLUFzcGVjdHMtb2YtRml0dGluZy1SZWdyZXNzaW9uLU1vZGVsczwvc3Bhbj4="}
[<span class='chapter-number'>3</span>  <span class='chapter-title'>Missing Data</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9taXNzaW5nLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjM8L3NwYW4+LS08c3Bhbi1jbGFzcz0nY2hhcHRlci10aXRsZSc+TWlzc2luZy1EYXRhPC9zcGFuPg=="}
[<span class='chapter-number'>4</span>  <span class='chapter-title'>Multivariable Modeling Strategies</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9tdWx0aXZhci5odG1sPHNwYW4tY2xhc3M9J2NoYXB0ZXItbnVtYmVyJz40PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPk11bHRpdmFyaWFibGUtTW9kZWxpbmctU3RyYXRlZ2llczwvc3Bhbj4="}
[<span class='chapter-number'>5</span>  <span class='chapter-title'>Describing, Resampling, Validating, and Simplifying the Model</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi92YWxpZGF0ZS5odG1sPHNwYW4tY2xhc3M9J2NoYXB0ZXItbnVtYmVyJz41PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkRlc2NyaWJpbmcsLVJlc2FtcGxpbmcsLVZhbGlkYXRpbmcsLWFuZC1TaW1wbGlmeWluZy10aGUtTW9kZWw8L3NwYW4+"}
[<span class='chapter-number'>6</span>  <span class='chapter-title'>R Software</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9zb2Z0d2FyZS5odG1sPHNwYW4tY2xhc3M9J2NoYXB0ZXItbnVtYmVyJz42PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPlItU29mdHdhcmU8L3NwYW4+"}
[<span class='chapter-number'>7</span>  <span class='chapter-title'>Modeling Longitudinal Responses using Generalized Least Squares</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9sb25nLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjc8L3NwYW4+LS08c3Bhbi1jbGFzcz0nY2hhcHRlci10aXRsZSc+TW9kZWxpbmctTG9uZ2l0dWRpbmFsLVJlc3BvbnNlcy11c2luZy1HZW5lcmFsaXplZC1MZWFzdC1TcXVhcmVzPC9zcGFuPg=="}
[<span class='chapter-number'>8</span>  <span class='chapter-title'>Case Study in Data Reduction</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9pbXByZWQuaHRtbDxzcGFuLWNsYXNzPSdjaGFwdGVyLW51bWJlcic+ODwvc3Bhbj4tLTxzcGFuLWNsYXNzPSdjaGFwdGVyLXRpdGxlJz5DYXNlLVN0dWR5LWluLURhdGEtUmVkdWN0aW9uPC9zcGFuPg=="}
[<span class='chapter-number'>9</span>  <span class='chapter-title'>Overview of Maximum Likelihood Estimation</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9tbGUuaHRtbDxzcGFuLWNsYXNzPSdjaGFwdGVyLW51bWJlcic+OTwvc3Bhbj4tLTxzcGFuLWNsYXNzPSdjaGFwdGVyLXRpdGxlJz5PdmVydmlldy1vZi1NYXhpbXVtLUxpa2VsaWhvb2QtRXN0aW1hdGlvbjwvc3Bhbj4="}
[<span class='chapter-number'>10</span>  <span class='chapter-title'>Binary Logistic Regression</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9scm0uaHRtbDxzcGFuLWNsYXNzPSdjaGFwdGVyLW51bWJlcic+MTA8L3NwYW4+LS08c3Bhbi1jbGFzcz0nY2hhcHRlci10aXRsZSc+QmluYXJ5LUxvZ2lzdGljLVJlZ3Jlc3Npb248L3NwYW4+"}
[<span class='chapter-number'>11</span>  <span class='chapter-title'>Binary Logistic Regression Case Study 1</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9scmNhc2UxLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjExPC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkJpbmFyeS1Mb2dpc3RpYy1SZWdyZXNzaW9uLUNhc2UtU3R1ZHktMTwvc3Bhbj4="}
[<span class='chapter-number'>12</span>  <span class='chapter-title'>Logistic Model Case Study: Survival of Titanic Passengers</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi90aXRhbmljLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjEyPC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkxvZ2lzdGljLU1vZGVsLUNhc2UtU3R1ZHk6LVN1cnZpdmFsLW9mLVRpdGFuaWMtUGFzc2VuZ2Vyczwvc3Bhbj4="}
[<span class='chapter-number'>13</span>  <span class='chapter-title'>Ordinal Logistic Regression</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9vcmRpbmFsLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjEzPC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPk9yZGluYWwtTG9naXN0aWMtUmVncmVzc2lvbjwvc3Bhbj4="}
[<span class='chapter-number'>14</span>  <span class='chapter-title'>Case Study in Ordinal Regression, Data Reduction, and Penalization</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9vcmRjYXNlLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjE0PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkNhc2UtU3R1ZHktaW4tT3JkaW5hbC1SZWdyZXNzaW9uLC1EYXRhLVJlZHVjdGlvbiwtYW5kLVBlbmFsaXphdGlvbjwvc3Bhbj4="}
[<span class='chapter-number'>15</span>  <span class='chapter-title'>Regression Models for Continuous $Y$ and Case Study in Ordinal Regression</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9jb255Lmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjE1PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPlJlZ3Jlc3Npb24tTW9kZWxzLWZvci1Db250aW51b3VzLSRZJC1hbmQtQ2FzZS1TdHVkeS1pbi1PcmRpbmFsLVJlZ3Jlc3Npb248L3NwYW4+"}
[<span class='chapter-number'>16</span>  <span class='chapter-title'>Transform-Both-Sides Regression</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9hcmVnLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjE2PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPlRyYW5zZm9ybS1Cb3RoLVNpZGVzLVJlZ3Jlc3Npb248L3NwYW4+"}
[<span class='chapter-number'>17</span>  <span class='chapter-title'>Introduction to Survival Analysis</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9zdXJ2Lmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjE3PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkludHJvZHVjdGlvbi10by1TdXJ2aXZhbC1BbmFseXNpczwvc3Bhbj4="}
[<span class='chapter-number'>18</span>  <span class='chapter-title'>Parametric Survival Models</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9wYXJzdXJ2Lmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjE4PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPlBhcmFtZXRyaWMtU3Vydml2YWwtTW9kZWxzPC9zcGFuPg=="}
[<span class='chapter-number'>19</span>  <span class='chapter-title'>Case Study in Parametric Survival Modeling and Model Approximation</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9wc21jYXNlLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjE5PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkNhc2UtU3R1ZHktaW4tUGFyYW1ldHJpYy1TdXJ2aXZhbC1Nb2RlbGluZy1hbmQtTW9kZWwtQXBwcm94aW1hdGlvbjwvc3Bhbj4="}
[<span class='chapter-number'>20</span>  <span class='chapter-title'>Cox Proportional Hazards Regression Model</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9jb3guaHRtbDxzcGFuLWNsYXNzPSdjaGFwdGVyLW51bWJlcic+MjA8L3NwYW4+LS08c3Bhbi1jbGFzcz0nY2hhcHRlci10aXRsZSc+Q294LVByb3BvcnRpb25hbC1IYXphcmRzLVJlZ3Jlc3Npb24tTW9kZWw8L3NwYW4+"}
[<span class='chapter-number'>21</span>  <span class='chapter-title'>Case Study in Cox Regression</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9jb3hjYXNlLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjIxPC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkNhc2UtU3R1ZHktaW4tQ294LVJlZ3Jlc3Npb248L3NwYW4+"}
[<span class='chapter-number'>22</span>  <span class='chapter-title'>Semiparametric Ordinal Longitudinal Models</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9tYXJrb3YuaHRtbDxzcGFuLWNsYXNzPSdjaGFwdGVyLW51bWJlcic+MjI8L3NwYW4+LS08c3Bhbi1jbGFzcz0nY2hhcHRlci10aXRsZSc+U2VtaXBhcmFtZXRyaWMtT3JkaW5hbC1Mb25naXR1ZGluYWwtTW9kZWxzPC9zcGFuPg=="}
[<span class='chapter-number'>23</span>  <span class='chapter-title'>Body Fat: Case Study in Linear Modeling</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9ib2R5ZmF0Lmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjIzPC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkJvZHktRmF0Oi1DYXNlLVN0dWR5LWluLUxpbmVhci1Nb2RlbGluZzwvc3Bhbj4="}
[<span class='chapter-number'>24</span>  <span class='chapter-title'>Bacteremia: Case Study in Nonlinear Data Reduction with Imputation</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9iYWN0ZXJlbWlhLmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjI0PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPkJhY3RlcmVtaWE6LUNhc2UtU3R1ZHktaW4tTm9ubGluZWFyLURhdGEtUmVkdWN0aW9uLXdpdGgtSW1wdXRhdGlvbjwvc3Bhbj4="}
[<span class='chapter-number'>25</span>  <span class='chapter-title'>Ordinal Semiparametric Regression for Survival Analysis</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9vcmRzdXJ2Lmh0bWw8c3Bhbi1jbGFzcz0nY2hhcHRlci1udW1iZXInPjI1PC9zcGFuPi0tPHNwYW4tY2xhc3M9J2NoYXB0ZXItdGl0bGUnPk9yZGluYWwtU2VtaXBhcmFtZXRyaWMtUmVncmVzc2lvbi1mb3ItU3Vydml2YWwtQW5hbHlzaXM8L3NwYW4+"}
[List of Figures]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9maWd1cmVzLmh0bWxMaXN0LW9mLUZpZ3VyZXM="}
[References]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWludC1zaWRlYmFyOi9yZWZlcmVuY2VzLmh0bWxSZWZlcmVuY2Vz"}
[<span class='chapter-number'>7</span>  <span class='chapter-title'>Modeling Longitudinal Responses using Generalized Least Squares</span>]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLWJyZWFkY3J1bWJzLTxzcGFuLWNsYXNzPSdjaGFwdGVyLW51bWJlcic+Nzwvc3Bhbj4tLTxzcGFuLWNsYXNzPSdjaGFwdGVyLXRpdGxlJz5Nb2RlbGluZy1Mb25naXR1ZGluYWwtUmVzcG9uc2VzLXVzaW5nLUdlbmVyYWxpemVkLUxlYXN0LVNxdWFyZXM8L3NwYW4+"}

:::{.hidden .quarto-markdown-envelope-contents render-id="Zm9vdGVyLWxlZnQtQ29weXJpZ2h0IDIwMjYsIEZyYW5rIEUgSGFycmVsbCBKcg=="}
Copyright 2026, Frank E Harrell Jr
:::


:::{.hidden .quarto-markdown-envelope-contents render-id="Zm9vdGVyLXJpZ2h0LWh0dHA6Ly9jcmVhdGl2ZWNvbW1vbnMub3JnL2xpY2Vuc2VzL2J5LW5jLXNhLzQuMA=="}
License
:::

:::



:::{#quarto-meta-markdown .hidden}
[[[7]{.chapter-number}  [Modeling Longitudinal Responses using Generalized Least Squares]{.chapter-title}]{#sec-long .quarto-section-identifier} – Regression Modeling Strategies]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLW1ldGF0aXRsZQ=="}
[[[7]{.chapter-number}  [Modeling Longitudinal Responses using Generalized Least Squares]{.chapter-title}]{#sec-long .quarto-section-identifier} – Regression Modeling Strategies]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLXR3aXR0ZXJjYXJkdGl0bGU="}
[[[7]{.chapter-number}  [Modeling Longitudinal Responses using Generalized Least Squares]{.chapter-title}]{#sec-long .quarto-section-identifier} – Regression Modeling Strategies]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLW9nY2FyZHRpdGxl"}
[Regression Modeling Strategies]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLW1ldGFzaXRlbmFtZQ=="}
[]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLXR3aXR0ZXJjYXJkZGVzYw=="}
[]{.hidden .quarto-markdown-envelope-contents render-id="cXVhcnRvLW9nY2FyZGRkZXNj"}
:::




<!-- -->

::: {.quarto-embedded-source-code}
```````````````````{.markdown shortcodes="false"}
```{r include=FALSE}
require(Hmisc)
options(qproject='rms', prType='html')
require(qreport)
getRs('qbookfun.r')
hookaddcap()
knitr::set_alias(w = 'fig.width', h = 'fig.height', cap = 'fig.cap', scap ='fig.scap')
knitr::read_chunk('long.R')

8 Modeling Longitudinal Responses using Generalized Least Squares

Some good general references on longitudinal data analysis are Davis (2002), Pinheiro & Bates (2000), Diggle et al. (2002), Venables & Ripley (2003), Hand & Crowder (1996), Verbeke & Molenberghs (2000), Lindsey (1997)

8.1 Notation

r mrg(sound("gls-1"))

• \(N\) subjects r ipacue()
• Subject \(i\) (\(i=1,2,\ldots,N\)) has \(n_{i}\) responses measured at times \(t_{i1}, t_{i2}, \ldots, t_{in_{i}}\)
• Response at time \(t\) for subject \(i\): \(Y_{it}\)
• Subject \(i\) has baseline covariates \(X_{i}\)
• Generally the response measured at time \(t_{i1}=0\) is a covariate in \(X_{i}\) instead of being the first measured response \(Y_{i0}\)
• Time trend in response is modeled with \(k\) parameters so that the time “main effect” has \(k\) d.f.
• Let the basis functions modeling the time effect be \(g_{1}(t), g_{2}(t), \ldots, g_{k}(t)\)

8.2 Model Specification for Effects on \(E(Y)\)

r mrg(sound("gls-2"))

8.2.1 Common Basis Functions

• \(k\) dummy variables for \(k+1\) unique times (assumes no r ipacue() functional form for time but may spend many d.f.)
• \(k=1\) for linear time trend, \(g_{1}(t)=t\)
• \(k\)–order polynomial in \(t\)
• \(k+1\)–knot restricted cubic spline (one linear term, \(k-1\) nonlinear terms)

8.2.2 Model for Mean Profile

• A model for mean time-response profile without interactions between time r ipacue() and any \(X\):
  \(E[Y_{it} | X_{i}] = X_{i}\beta + \gamma_{1}g_{1}(t) + \gamma_{2}g_{2}(t) + \ldots + \gamma_{k}g_{k}(t)\)
• Model with interactions between time and some \(X\)’s: add product terms for desired interaction effects
• Example: To allow the mean time trend for subjects in group 1 (reference group) to be arbitrarily different from time trend for subjects in group 2, have a dummy variable for group 2, a time “main effect” curve with \(k\) d.f. and all \(k\) products of these time components with the dummy variable for group 2
• Time should be modeled using indicator variables only when time is really discrete, e.g., when time is in weeks and subjects were followed at exactly the intended weeks. In general time should be modeled continuously (and nonlinearly if there are more than 2 followup times) using actual visit dates instead of intended dates (Donohue et al., n.d.).

8.2.3 Model Specification for Treatment Comparisons

r mrg(sound("gls-3"))

• In studies comparing two or more treatments, a response is often r ipacue() measured at baseline (pre-randomization)
• Analyst has the option to use this measurement as \(Y_{i0}\) or as part of \(X_{i}\)

NoteComments from Jim Rochon

For RCTs, I draw a sharp line at the point when the intervention begins. The LHS [left hand side of the model equation] is reserved for something that is a response to treatment. Anything before this point can potentially be included as a covariate in the regression model. This includes the “baseline” value of the outcome variable. Indeed, the best predictor of the outcome at the end of the study is typically where the patient began at the beginning. It drinks up a lot of variability in the outcome; and, the effect of other covariates is typically mediated through this variable.

I treat anything after the intervention begins as an outcome. In the western scientific method, an “effect” must follow the “cause” even if by a split second.

Note that an RCT is different than a cohort study. In a cohort study, “Time 0” is not terribly meaningful. If we want to model, say, the trend over time, it would be legitimate, in my view, to include the “baseline” value on the LHS of that regression model.

Now, even if the intervention, e.g., surgery, has an immediate effect, I would include still reserve the LHS for anything that might legitimately be considered as the response to the intervention. So, if we cleared a blocked artery and then measured the MABP, then that would still be included on the LHS.

Now, it could well be that most of the therapeutic effect occurred by the time that the first repeated measure was taken, and then levels off. Then, a plot of the means would essentially be two parallel lines and the treatment effect is the distance between the lines, i.e., the difference in the intercepts.

If the linear trend from baseline to Time 1 continues beyond Time 1, then the lines will have a common intercept but the slopes will diverge. Then, the treatment effect will the difference in slopes.

One point to remember is that the estimated intercept is the value at time 0 that we predict from the set of repeated measures post randomization. In the first case above, the model will predict different intercepts even though randomization would suggest that they would start from the same place. This is because we were asleep at the switch and didn’t record the “action” from baseline to time 1. In the second case, the model will predict the same intercept values because the linear trend from baseline to time 1 was continued thereafter.

More importantly, there are considerable benefits to including it as a covariate on the RHS. The baseline value tends to be the best predictor of the outcome post-randomization, and this maneuver increases the precision of the estimated treatment effect. Additionally, any other prognostic factors correlated with the outcome variable will also be correlated with the baseline value of that outcome, and this has two important consequences. First, this greatly reduces the need to enter a large number of prognostic factors as covariates in the linear models. Their effect is already mediated through the baseline value of the outcome variable. Secondly, any imbalances across the treatment arms in important prognostic factors will induce an imbalance across the treatment arms in the baseline value of the outcome. Including the baseline value thereby reduces the need to enter these variables as covariates in the linear models.

Senn (2006) states that temporally and logically, a “baseline cannot be a response to treatment”, so baseline and response cannot be modeled in an integrated framework.

r quoteit('... one should focus clearly on \'outcomes\' as being the only values that can be influenced by treatment and examine critically any schemes that assume that these are linked in some rigid and deterministic view to \'baseline\' values. An alternative tradition sees a baseline as being merely one of a number of measurements capable of improving predictions of outcomes and models it in this way.')

The final reason that baseline cannot be modeled as the response at r ipacue() time zero is that many studies have inclusion/exclusion criteria that include cutoffs on the baseline variable. In other words, the baseline measurement comes from a truncated distribution. In general it is not appropriate to model the baseline with the same distributional shape as the follow-up measurements. Thus the approaches recommended by Liang & Zeger (2000) and Liu et al. (2009) are problematic¹⁸.

¹⁸ In addition to this, one of the paper’s conclusions that analysis of covariance is not appropriate if the population means of the baseline variable are not identical in the treatment groups is not correct (Senn, 2006). See Kenward et al. (2010) for a rebuke of Liu et al. (2009).

8.3 Modeling Within-Subject Dependence

r mrg(sound("gls-4"))

• Random effects and mixed effects models have become very popular r ipacue()
• Disadvantages:
  • Induced correlation structure for \(Y\) may be unrealistic
  • Numerically demanding
  • Require complex approximations for distributions of test statistics
• Conditional random effects vs. (subject-) marginal models:
  • Random effects are subject-conditional
  • Random effects models are needed to estimate responses for individual subjects
  • Models without random effects are marginalized with respect to subject-specific effects
  • They are natural when the interest is on group-level (i.e., covariate-specific but not patient-specific) parameters (e.g., overall treatment effect)
  • Random effects are natural when there is clustering at more than the subject level (multi-level models)
• Extended linear model (marginal; with no random effects) is a logical extension of the univariate model (e.g., few statisticians use subject random effects for univariate \(Y\))
• This was known as growth curve models and generalized least squares (Goldstein, 1989; Potthoff & Roy, 1964) and was developed long before mixed effect models became popular
• Pinheiro and Bates (Section~5.1.2) state that “in some applications, one may wish to avoid incorporating random effects in the model to account for dependence among observations, choosing to use the within-group component \(\Lambda_{i}\) to directly model variance-covariance structure of the response.”
• We will assume that \(Y_{it} | X_{i}\) has a multivariate normal r ipacue() distribution with mean given above and with variance-covariance matrix \(V_{i}\), an \(n_{i}\times n_{i}\) matrix that is a function of \(t_{i1}, \ldots, t_{in_{i}}\)
• We further assume that the diagonals of \(V_{i}\) are all equal
• Procedure can be generalized to allow for heteroscedasticity over time or with respect to \(X\) (e.g., males may be allowed to have a different variance than females)
• This extended linear model has the following assumptions: r ipacue()
  • all the assumptions of OLS at a single time point including correct modeling of predictor effects and univariate normality of responses conditional on \(X\)
  • the distribution of two responses at two different times for the same subject, conditional on \(X\), is bivariate normal with a specified correlation coefficient
  • the joint distribution of all \(n_{i}\) responses for the \(i^{th}\) subject is multivariate normal with the given correlation pattern (which implies the previous two distributional assumptions)
  • responses from any times for any two different subjects are uncorrelated

What Methods To Use for Repeated Measurements / Serial Data? ^{19 20}

	Repeated Measures ANOVA	GEE	Mixed Effects Models	GLS	Markov	LOCF	Summary Statistic21
Assumes normality	×		×	×
Assumes independence of measurements within subject	×22	×23
Assumes a correlation structure24	×	×25	×	×	×
Requires same measurement times for all subjects	×					?
Does not allow smooth modeling of time to save d.f.	×
Does not allow adjustment for baseline covariates	×
Does not easily extend to non-continuous \(Y\)	×			×
Loses information by not using intermediate measurements						×26	×
Does not allow widely varying # observations per subject	×	×27				×	×28
Does not allow for subjects to have distinct trajectories29	×	×		×	×	×
Assumes subject-specific effects are Gaussian			×
Badly biased if non-random dropouts	?	×				×
Biased in general						×
Harder to get tests & CLs			×30			×31
Requires large # subjects/clusters		×
SEs are wrong	×32					×
Assumptions are not verifiable in small samples	×	N/A	×	×		×
Does not extend to complex settings such as time-dependent covariates and dynamic 33 models	×		×	×		×	?

¹⁹ Thanks to Charles Berry, Brian Cade, Peter Flom, Bert Gunter, and Leena Choi for valuable input.

²⁰ GEE: generalized estimating equations; GLS: generalized least squares; LOCF: last observation carried forward.

²¹ E.g., compute within-subject slope, mean, or area under the curve over time. Assumes that the summary measure is an adequate summary of the time profile and assesses the relevant treatment effect.

²² Unless one uses the Huynh-Feldt or Greenhouse-Geisser correction

²³ For full efficiency, if using the working independence model

²⁴ Or requires the user to specify one

²⁵ For full efficiency of regression coefficient estimates

²⁶ Unless the last observation is missing

²⁷ The cluster sandwich variance estimator used to estimate SEs in GEE does not perform well in this situation, and neither does the working independence model because it does not weight subjects properly.

²⁸ Unless one knows how to properly do a weighted analysis

²⁹ Or users population averages

³⁰ Unlike GLS, does not use standard maximum likelihood methods yielding simple likelihood ratio \(\chi^2\) statistics. Requires high-dimensional integration to marginalize random effects, using complex approximations, and if using SAS, unintuitive d.f. for the various tests.

³¹ Because there is no correct formula for SE of effects; ordinary SEs are not penalized for imputation and are too small

³² If correction not applied

³³ E.g., a model with a predictor that is a lagged value of the response variable

• Markov models use ordinary univariate software and are very flexible r ipacue()
• They apply the same way to binary, ordinal, nominal, and continuous Y
• They require post-fitting calculations to get probabilities, means, and quantiles that are not conditional on the previous Y value

Gardiner et al. (2009) compared several longitudinal data r ipacue() models, especially with regard to assumptions and how regression coefficients are estimated. Peters et al. (2012) have an empirical study confirming that the “use all available data” approach of likelihood–based longitudinal models makes imputation of follow-up measurements unnecessary.

8.4 Parameter Estimation Procedure

r mrg(sound("gls-5"))

• Generalized least squares r ipacue()
• Like weighted least squares but uses a covariance matrix that is not diagonal
• Each subject can have her own shape of \(V_{i}\) due to each subject being measured at a different set of times
• Maximum likelihood
• Newton-Raphson or other trial-and-error methods used for estimating parameters
• For small number of subjects, advantages in using REML (restricted maximum likelihood) instead of ordinary MLE (Diggle et al., 2002, p. Section~5.3), (Pinheiro & Bates, 2000, p. Chapter~5), Goldstein (1989) (esp. to get more unbiased estimate of the covariance matrix)
• When imbalances are not severe, OLS fitted ignoring subject r ipacue() identifiers may be efficient
  • But OLS standard errors will be too small as they don’t take intra-cluster correlation into account
  • May be rectified by substituting covariance matrix estimated from Huber-White cluster sandwich estimator or from cluster bootstrap
• When imbalances are severe and intra-subject correlations are r ipacue() strong, OLS is not expected to be efficient because it gives equal weight to each observation
  • a subject contributing two distant observations receives \(\frac{1}{5}\) the weight of a subject having 10 tightly-spaced observations

8.5 Common Correlation Structures

r mrg(sound("gls-6"))

• Usually restrict ourselves to isotropic correlation structures r ipacue() — correlation between responses within subject at two times depends only on a measure of distance between the two times, not the individual times (this will be relaxed in the last structure covered below)
• We simplify further and assume depends on \(|t_{1} - t_{2}|\)
• Can speak interchangeably of correlations of residuals within subjects or correlations between responses measured at different times on the same subject, conditional on covariates \(X\)
• Assume that the correlation coefficient for \(Y_{it_{1}}\) vs. \(Y_{it_{2}}\) conditional on baseline covariates \(X_{i}\) for subject \(i\) is \(h(|t_{1} - t_{2}|, \rho)\), where \(\rho\) is a vector (usually a scalar) set of fundamental correlation parameters
• Some commonly used structures when times are continuous and are r ipacue() not equally spaced (Pinheiro & Bates, 2000, Section 5.3.3) (nlme correlation function names are at the right if the structure is implemented in nlme):

Table 8.1: Some longitudinal data correlation structures

Structure	nlme Function
Compound symmetry: \(h = \rho\) if \(t_{1} \neq t_{2}\), 1 if \(t_{1}=t_{2}\) 34	corCompSymm
Autoregressive-moving average lag 1: \(h = \rho^{|t_{1} - t_{2}|} = \rho^s\) where \(s = |t_{1}-t_{2}|\)	corCAR1
Exponential: \(h = \exp(-s/\rho)\)	corExp
Gaussian: \(h = \exp[-(s/\rho)^2]\)	corGaus
Linear: \(h = (1 - s/\rho)[s < \rho]\)	corLin
Rational quadratic: \(h = 1 - (s/\rho)^{2}/[1+(s/\rho)^{2}]\)	corRatio
Spherical: \(h = [1-1.5(s/\rho)+0.5(s/\rho)^{3}][s < \rho]\)	corSpher
Linear exponent AR(1): \(h = \rho^{d_{min} + \delta\frac{s - d_{min}}{d_{max} - d_{min}}}\), 1 if \(t_{1}=t_{2}\)	Simpson et al. (2010)
Exponential decline with floor and non-isotropy: \(h = k + (1-k) \exp(-\exp(a+b\frac{t_{1}+t_{2}}{2})|t_{1}-t_{2}|)\)	corFloorExp

³⁴ Essentially what two-way ANOVA assumes

The structures 3-7 use \(\rho\) as a scaling parameter, not as something restricted to be in \([0,1]\)

The last structure is implemented in the rms package. The \(k\) parameter is restricted to be in \([0,1]\) by modeling it as \(\text{expit}(k^{*})\) where \(k^{*}\) has an unrestricted range and \(\text{expit}(x) = \frac{1}{1 + \exp(-x)}\). This structure allows for non-isotropy and for a nonzero limiting correlation as the time lag \(\rightarrow \infty\). It is a shared random effect plus serial decay mixture (a continuous-time analog of a random intercept superimposed on an Ornstein-Uhlenbeck exponential-decay process). \(k\) is the correlation contributed by a persistent, subject-level component that never fades; \(1 - k\) is the share of correlation carried by a component that decays with lag at rate \(\exp(a + b\frac{t_{1}+t_{2}}{2})\). \(b\) is the non-isotropy term, letting the rate grow or shrink with the time pair’s midpoint.

8.6 Checking Model Fit

r mrg(sound("gls-7"))

• Constant variance assumption: usual residual plots r ipacue()
• Normality assumption: usual qq residual plots
• Correlation pattern: Variogram
  • Estimate correlations of all possible pairs of residuals at different time points
  • Pool all estimates at same absolute difference in time \(s\)
  • Variogram is a plot with \(y = 1 - \hat{h}(s, \rho)\) vs. \(s\) on the \(x\)-axis
  • Superimpose the theoretical variogram assumed by the model
  • This does not allow for non-isotropy

8.7 R Software

r mrg(sound("gls-8"))

• Nonlinear mixed effects model package of Pinheiro & Bates r ipacue()
• For linear models, fitting functions are
  • lme for mixed effects models
  • gls for generalized least squares without random effects
• For this version the rms package has Gls so that many features of rms can be used:
  • anova: all partial Wald tests, test of linearity, pooled tests
  • summary: effect estimates (differences in \(\hat{Y}\)) and confidence limits, can be plotted
  • plot, ggplot, plotp: continuous effect plots
  • nomogram: nomogram
  • Function: generate R function code for fitted model
  • latex:  representation of fitted model

In addition, Gls has a bootstrap option (hence you do not use rms’s bootcov for Gls fits).
To get regular gls functions named anova (for likelihood ratio tests, AIC, etc.) or summary use anova.gls or summary.gls * nlme package has many graphics and fit-checking functions * Several functions will be demonstrated in the case study

8.8 Case Study

r mrg(sound("gls-9")) Consider the dataset in Table~6.9 of Davis[davis-repmeas, pp. 161-163] from a multi-center, randomized controlled trial of botulinum toxin type B (BotB) in patients with cervical dystonia from nine U.S. sites.

• Randomized to placebo (\(N=36\)), 5000 units of BotB (\(N=36\)), r ipacue() 10,000 units of BotB (\(N=37\))
• Response variable: total score on Toronto Western Spasmodic Torticollis Rating Scale (TWSTRS), measuring severity, pain, and disability of cervical dystonia (high scores mean more impairment)
• TWSTRS measured at baseline (week 0) and weeks 2, 4, 8, 12, 16 after treatment began
• Dataset cdystonia from web site

8.8.1 Graphical Exploration of Data

```{r spaghetti,h=5,w=7,cap=‘Time profiles for individual subjects, stratified by study site and dose’} #| label: fig-long-spaghetti require(rms) require(data.table) options(prType=‘html’) # for model print, summary, anova, validate getHdata(cdystonia) setDT(cdystonia) # convert to data.table cdystonia[, uid := paste(site, id)] # unique subject ID

9 Tabulate patterns of subjects’ time points

g <- function(w) paste(sort(unique(w)), collapse=’ ’) cdystonia[, table(tapply(week, uid, g))]

10 Plot raw data, superposing subjects

xl <- xlab(‘Week’); yl <- ylab(‘TWSTRS-total score’) ggplot(cdystonia, aes(x=week, y=twstrs, color=factor(id))) + geom_line() + xl + yl + facet_grid(treat ~ site) + guides(color=FALSE)

```{r quartiles,cap='Quartiles of `TWSTRS` stratified by dose',w=5,h=4}
#| label: fig-long-quartiles
# Show quartiles
g <- function(x) {
  k <- as.list(quantile(x, (1 : 3) / 4, na.rm=TRUE))
  names(k) <- .q(Q1, Q2, Q3)
  k
}
cdys <- cdystonia[, g(twstrs), by=.(treat, week)]
ggplot(cdys, aes(x=week, y=Q2)) + xl + yl + ylim(0, 70) +
  geom_line() + facet_wrap(~ treat, nrow=2) +
  geom_ribbon(aes(ymin=Q1, ymax=Q3), alpha=0.2)

{r bootcls,cap='Mean responses and nonparametric bootstrap 0.95 confidence limits for population means, stratified by dose',w=5,h=4} #| label: fig-long-bootcls # Show means with bootstrap nonparametric CLs cdys <- cdystonia[, as.list(smean.cl.boot(twstrs)), by = list(treat, week)] ggplot(cdys, aes(x=week, y=Mean)) + xl + yl + ylim(0, 70) + geom_line() + facet_wrap(~ treat, nrow=2) + geom_ribbon(aes(x=week, ymin=Lower, ymax=Upper), alpha=0.2)

Model with \(Y_{i0}\) as Baseline Covariate

quarto-executable-code-5450563D

baseline <- cdystonia[week == 0]
baseline[, week := NULL]
setnames(baseline, 'twstrs', 'twstrs0')
followup <- cdystonia[week > 0, .(uid, week, twstrs)]
setkey(baseline, uid)
setkey(followup, uid, week)
both     <- Merge(baseline, followup, id = ~ uid)
# Remove person with no follow-up record
both     <- both[! is.na(week)]
dd       <- datadist(both)
options(datadist='dd')

10.0.1 Using Generalized Least Squares

r mrg(sound("gls-10")) We stay with baseline adjustment and use a variety of correlation r ipacue() structures, with constant variance. Time is modeled as a restricted cubic spline with 3 knots, because there are only 3 unique interior values of week.

AIC computed above is set up so that smaller values are better. From this the 3-parameter correlation structure is the best fitting. For the remainder of the analysis use corFloorExp, using Gls.

Keselman et al. (1998) did a simulation study to study the reliability of AIC for selecting the correct covariance structure in repeated measurement models. In choosing from among 11 structures, AIC selected the correct structure 47% of the time. Gurka et al. (2011) demonstrated that fixed effects in a mixed effects model can be biased, independent of sample size, when the specified covariate matrix is more restricted than the true one.

Plot the estimated correlation function for all observed pairs of times.

{r fig-long-corFloorExp, fig.cap='Fitted 3-parameter correlation $\\hat{\\rho}$ as a function of time lag ($x$) and time midpoint (numbers)',h=3.75,w=4.25} co <- coef(cs) g <- function(t1, t2) { k <- plogis(co["k"]) a <- co["a"] b <- co["b"] k + (1 - k) * exp(-exp(a + b * (t1 + t2) / 2) * abs(t1 - t2)) } times <- sort(unique(both$week)) w <- setDT(expand.grid(t1 = times, t2 = times)) w <- w[t2 > t1] w[, lag := t2 - t1] w[, mid := (t1 + t2) / 2] w[, rho := g(t1, t2)] ggplot(w, aes(x = lag, y = rho, label=mid)) + geom_text(size=2.5) + xlab(expression(t[2] - t[1])) + ylab(expression(hat(rho)))

The numbers drawn for each point are the values of \(\frac{t_{1} + t_{2}}{2}\). The degree to which points for the same time lag have different vertical values is the degree to which time midpoints figure into the correlation structure in addition to time lags, i.e., the degree of non-isotropy. We see a moderate degree of non-isotropy. We also see the possibility of greater-than-zero long-term correlation, i.e., a need for a compound symmetric correlation term and slight inadequacy of AR(1).

\(\hat{\rho} = 0.8672\) from an AR(1) fit, the estimate of the correlation between two r ipacue() measurements taken one week apart on the same subject. The estimated correlation for measurements 10 weeks apart is \(0.8672^{10} = 0.24\).

{r fig-long-variogram,fig.cap='Variogram, with assumed correlation pattern superimposed',h=3.75,w=4.25} #| label: fig-long-variogram

Check constant variance and normality assumptions: r ipacue()

{r fig-long-resid,h=6,w=7.5,cap='Three residual plots to check for absence of trends in central tendency and in variability. Upper right panel shows the baseline score on the $x$-axis. Bottom left panel shows the mean $\\pm 2\\times$ SD. Bottom right panel is the QQ plot for checking normality of residuals from the GLS fit.'} #| label: fig-long-resid

Now get hypothesis tests, estimates, and graphically interpret the model. r mrg(sound("gls-11"))

{r fig-long-anova,cap='Results of `anova.rms` from generalized least squares fit with continuous time AR1 correlation structure',w=5,h=4} #| label: fig-long-anova

{r fig-long-pleffects,h=5.5,w=7,cap='Estimated effects of time, baseline `TWSTRS`, age, and sex'} #| label: fig-long-pleffects

{r fig-long-contrasts,h=4,w=6,cap='Contrasts and 0.95 confidence limits from GLS fit'} #| label: fig-long-contrasts

Although multiple d.f. tests such as total treatment effects or r ipacue() treatment \(\times\) time interaction tests are comprehensive, their increased degrees of freedom can dilute power. In a treatment comparison, treatment contrasts at the last time point (single d.f. tests) are often of major interest. Such contrasts are informed by all the measurements made by all subjects (up until dropout times) when a smooth time trend is assumed.

{r nomogram,h=6.5,w=7.5,cap='Nomogram from GLS fit. Second axis is the baseline score.'} #| label: fig-long-nomogram n <- nomogram(a, age=c(seq(20, 80, by=10), 85)) plot(n, cex.axis=.55, cex.var=.8, lmgp=.25) # Figure (*\ref{fig:longit-nomogram}*)

10.0.2 Semiparametric Random Effects Models

quarto-executable-code-5450563D

f <- orm(twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) +
         rcs(age, 4) * sex + cluster(uid), data=both, x=TRUE, y=TRUE, maxit.outer=300, trace=1)
f
Olinks(f)

Stick with the PO model.

quarto-executable-code-5450563D

both[, mix := (week - 2) / 14]
f <- orm(twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) +
         rcs(age, 4) * sex + cluster(uid) + mix_re(mix), data=both, maxit.outer=300, trace=1)
f
deviance(f)

Deviance did not improve, due to the estimates of \(\sigma_1\) and \(\sigma_2\) being so close.

10.0.3 Bayesian Proportional Odds Random Effects Model

• Develop a \(y\)-transformation invariant longitudinal model r ipacue()
• Proportional odds model with no grouping of TWSTRS scores
• Bayesian random effects model
• Random effects Gaussian with exponential prior distribution for its SD, with mean 1.0
• Compound symmetry correlation structure
• Demonstrates a large amount of patient-to-patient intercept variability

• Show the final graphic (high dose:placebo contrast as function of time r ipacue()
• Intervals are 0.95 highest posterior density intervals
• \(y\)-axis: log-odds ratio

For each posterior draw compute the difference in means and get an exact (to within simulation error) 0.95 highest posterior density intervals for these differences.

{r bayesfit5,w=7,h=3.75}

10.0.4 Bayesian Markov Semiparametric Model

• First-order Markov model r ipacue()
• Serial correlation induced by Markov model is similar to AR(1) which we already know fits these data
• Markov model is more likely to fit the data than the random effects model, which induces a compound symmetry correlation structure
• Models state transitions
• PO model at each visit, with Y from previous visit conditioned upon just like any covariate
• Need to uncondition (marginalize) on previous Y to get the time-response profile we usually need
• Semiparametric model is especially attractive because one can easily “uncondition” a discrete Y model, and the distribution of Y for control subjects can be any shape
• Let measurement times be \(t_{1}, t_{2}, \dots, t_{m}\), and the measurement for a subject at time \(t\) be denoted \(Y(t)\)
• First-order Markov model:

\[\begin{array}{ccc} \Pr(Y(t_{i}) \geq y | X, Y(t_{i-1})) &=& \mathrm{expit}(\alpha_{y} + X\beta\\ &+& g(Y(t_{i-1}), t_{i}, t_{i} - t_{i-1})) \end{array}\]

• \(g\) involves any number of regression coefficients for a main effect of \(t\), the main effect of time gap \(t_{i} - t_{i-1}\) if this is not collinear with absolute time, a main effect of the previous state, and interactions between these
• Examples of how the previous state may be modeled in \(g\):
  • linear in numeric codes for \(Y\)
  • spline function in same
  • discontinuous bi-linear relationship where there is a slope for in-hospital outcome severity, a separate slope for outpatient outcome severity, and an intercept jump at the transition from inpatient to outpatient (or vice versa)
• Markov model is quite flexible in handling time trends and serial correlation patterns
• Can allow for irregular measurement times:
  hbiostat.org/stat/irreg.html

Fit the model and run standard Stan diagnostics.

{r mark1,h=6,w=7.5}

Note that posterior sampling is much more efficient without random effects.

Let’s add subject-level random effects to the model. Smallness of the standard deviation of the random effects provides support for the assumption of conditional independence that we like to make for Markov models and allows us to simplify the model by omitting random effects.

The random effects SD is only 0.11 on the logit scale. Also, the standard deviations of all the regression parameter posterior distributions are virtually unchanged with the addition of random effects:

{r mark4r,w=7,h=7}

So we will use the model omitting random effects.

Show the partial effects of all the predictors, including the effect of the previous measurement of TWSTRS. Also compute high dose:placebo treatment contrasts on these conditional estimates.

Using posterior means for parameter values, compute the probability that at a given week twstrs will be \(\geq 40\) when at the previous visit it was 40. Also show the conditional mean twstrs when it was 40 at the previous visit.

• Semiparametric models provide not only estimates of tendencies of Y r ipacue() but also estimate the whole distribution of Y
• Estimate the entire conditional distribution of Y at week 12 for high-dose patients having TWSTRS=42 at week 8
• Other covariates set to median/mode
• Use posterior mean of all the cell probabilities
• Also show pointwise 0.95 highest posterior density intervals
• To roughly approximate simultaneous confidence bands make the pointwise limits sum to 1 like the posterior means do

• Repeat this showing the variation over 5 posterior draws r ipacue()

• Turn to marginalized (unconditional on previous twstrs) r ipacue() quantities
• Capitalize on PO model being a multinomial model, just with PO restrictions
• Manipulations of conditional probabilities to get the unconditional probability that twstrs=y doesn’t need to know about PO
• Compute all cell probabilities and use the law of total probability recursively \[\Pr(Y_{t} = y | X) = \sum_{j=1}^{k} \Pr(Y_{t} = y | X, Y_{t-1} = j) \Pr(Y_{t-1} = j | X)\]
• predict.blrm method with type='fitted.ind' computes the needed conditional cell probabilities, optionally for all posterior draws at once
• Easy to get highest posterior density intervals for derived parameters such as unconditional probabilities or unconditional means
• Hmisc package soprobMarkovOrdm function (in version 4.6) computes an array of all the state occupancy probabilities for all the posterior draws

• Use the same posterior draws of unconditional probabilities of all r ipacue() values of TWSTRS to get the posterior distribution of differences in mean TWSTRS between high and low dose

• Get posterior mean of all cell probabilities estimates at week 12 r ipacue()
• Distribution of TWSTRS conditional high dose, median age, mode sex
• Not conditional on week 8 value

10.0.5 Frequentist Markov Semiparametric Model

Fitting a frequentist MOST model will allow us to take advantage of the orm functions 2-parameter random effects contribution, which was designed to improve fit to correlation structure for Markov-1 models.

quarto-executable-code-5450563D

both[, mix2 := ifelse(week == 2, 0, 1)]
f <- orm(twstrs ~  treat * rcs(week, 3) + rcs(ptwstrs, 4) +
                    rcs(age, 4) * sex + cluster(uid) + mix_re(mix2),
         data=both, trace=1, x=TRUE)
deviance(f)

sigma1 large: substantial subject-level correlation at baseline (before ptwstrs is contributing anything), consistent with a real, non-trivial random effect there. sigma2 ≈ -0.2: small and negative, but not collapsed to the boundary the way the single-sigma model did. That’s exactly the pattern the additive design was meant to reveal — once the lag term is active, a plain constant-weight random effect would be adding back in some of the same correlation the lag term already explains, and a small negative sigma2 is the model’s way of trimming that slight over-induction rather than needing to zero the whole random effect out just to compensate.

## Study Questions

**Section 7.2**

1. When should one model the time-response profile using discrete time?
   
**Section 7.3**

1. What makes generalized least squares and mixed effect models
   relatively robust to non-completely-random dropouts?
1. What does the last observation carried forward method always violate?

**Section 7.4**

1. Which correlation structure do you expect to fit the data when there are rapid repetitions over a short time span?  When the follow-up time span is very long?

**Section 7.8**

1. What can go wrong if many correlation structures are tested in one dataset?
1. In a longitudinal intervention study, what is the most typical comparison of interest?  Is it best to borrow information in estimating this contrast?

```{r echo=FALSE}
saveCap('07')

``````````````````` :::

Davis, C. S. (2002). Statistical Methods for the Analysis of Repeated Measurements. Springer.

Diggle, P. J., Heagerty, P., Liang, K.-Y., & Zeger, S. L. (2002). Analysis of Longitudinal Data (second). Oxford University Press.

Donohue, M. C., Langford, O., Insel, P. S., van Dyck, C. H., Petersen, R. C., Craft, S., Sethuraman, G., Raman, R., Aisen, P. S., & Initiative, F. the A. D. N. (n.d.). Natural cubic splines for the analysis of Alzheimer’s clinical trials. Pharmaceutical Statistics, n/a(n/a). https://doi.org/10.1002/pst.2285

Gardiner, J. C., Luo, Z., & Roman, L. A. (2009). Fixed effects, random effects and GEE: What are the differences? Stat Med, 28, 221–239.

nice comparison of models; econometrics; different use of the term "fixed effects model"

Goldstein, H. (1989). Restricted unbiased iterative generalized least-squares estimation. Biometrika, 76(3), 622–623.

derivation of REML

Gurka, M. J., Edwards, L. J., & Muller, K. E. (2011). Avoiding bias in mixed model inference for fixed effects. Stat Med, 30(22), 2696–2707. https://doi.org/10.1002/sim.4293

Hand, D., & Crowder, M. (1996). Practical Longitudinal Data Analysis. Chapman & Hall.

Kenward, M. G., White, I. R., & Carpener, J. R. (2010). Should baseline be a covariate or dependent variable in analyses of change from baseline in clinical trials? (Letter to the editor). Stat Med, 29, 1455–1456.

sharp rebuke of liu09sho

Keselman, H. J., Algina, J., Kowalchuk, R. K., & Wolfinger, R. D. (1998). A comparison of two approaches for selecting covariance structures in the analysis of repeated measurements. Comm Stat - Sim Comp, 27, 591–604.

use of AIC and BIC for selecting the covariance structure in repeated measurements;serial data;longitudinal data;when chosing from 11 covariance patterns, AIC selected the correct structure 0.47 of the time; BIC was correct in 0.35

Liang, K.-Y., & Zeger, S. L. (2000). Longitudinal data analysis of continuous and discrete responses for pre-post designs. Sankhyā, 62, 134–148.

makes an error in assuming the baseline variable will have the same univariate distribution as the response except for a shift;baseline may have for example a truncated distribution based on a trial’s inclusion criteria;if correlation between baseline and response is zero, ANCOVA will be twice as efficient as simple analysis of change scores;if correlation is one they may be equally efficient

Lindsey, J. K. (1997). Models for Repeated Measurements. Clarendon Press.

Liu, G. F., Lu, K., Mogg, R., Mallick, M., & Mehrotra, D. V. (2009). Should baseline be a covariate or dependent variable in analyses of change from baseline in clinical trials? Stat Med, 28, 2509–2530.

seems to miss several important points, such as the fact that the baseline variable is often part of the inclusion/exclusion criteria and so has a truncated distribution that is different from that of the follow-up measurements;sharp rebuke in ken10sho

Peters, S. A., Bots, M. L., den Ruijter, H. M., Palmer, M. K., Grobbee, D. E., Crouse, J. R., O’Leary, D. H., Evans, G. W., Raichlen, J. S., Moons, K. G., Koffijberg, H., & METEOR study group. (2012). Multiple imputation of missing repeated outcome measurements did not add to linear mixed-effects models. J Clin Epi, 65(6), 686–695. https://doi.org/10.1016/j.jclinepi.2011.11.012

Pinheiro, J. C., & Bates, D. M. (2000). Mixed-Effects Models in S and S-PLUS. Springer.

Potthoff, R. F., & Roy, S. N. (1964). A generalized multivariate analysis of variance model useful especially for growth curve problems. Biometrika, 51, 313–326.

included an AR1 example

Senn, S. (2006). Change from baseline and analysis of covariance revisited. Stat Med, 25, 4334–4344.

shows that claims that in a 2-arm study it is not true that ANCOVA requires the population means at baseline to be identical;refutes some claims of lia00lon;problems with counterfactuals;temporal additivity ("amounts to supposing that despite the fact that groups are difference at baseline they would show the same evolution over time");causal additivity;is difficult to design trials for which simple analysis of change scores is unbiased, ANCOVA is biased, and a causal interpretation can be given;temporally and logically, a "baseline cannot be a \(<\)i\(>\)response\(<\)/i\(>\) to treatment", so baseline and response cannot be modeled in an integrated framework as Laird and Ware’s model has been used;"one should focus clearly on “outcomes” as being the only values that can be influenced by treatment and examine critically any schemes that assume that these are linked in some rigid and deterministic view to “baseline” values. An alternative tradition sees a baseline as being merely one of a number of measurements capable of improving predictions of outcomes and models it in this way.";"You cannot establish necessary conditions for an estimator to be valid by nominating a model and seeing what the model implies unless the model is universally agreed to be impeccable. On the contrary it is appropriate to start with the estimator and see what assumptions are implied by valid conclusions.";this is in distinction to lia00lon

Simpson, S. L., Edwards, L. J., Muller, K. E., Sen, P. K., & Styner, M. A. (2010). A linear exponent AR(1) family of correlation structures. Stat Med, 29, 1825–1838.

Venables, W. N., & Ripley, B. D. (2003). Modern Applied Statistics with S (Fourth). Springer-Verlag.

Verbeke, G., & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. Springer.