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)\)

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)

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.).

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}\)

Comments 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¹.

¹ 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? ² ³
	Repeated Measures ANOVA	GEE	Mixed Effects Models	GLS	Markov	LOCF	Summary Statistic⁴
Assumes normality	×		×	×
Assumes independence of measurements within subject	×⁵	×⁶
Assumes a correlation structure⁷	×	×⁸	×	×	×
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						×⁹	×
Does not allow widely varying # observations per subject	×	×¹⁰				×	×¹¹
Does not allow for subjects to have distinct trajectories¹²	×	×		×	×	×
Assumes subject-specific effects are Gaussian			×
Badly biased if non-random dropouts	?	×				×
Biased in general						×
Harder to get tests & CLs			×¹³			×¹⁴
Requires large # subjects/clusters		×
SEs are wrong	×¹⁵					×
Assumptions are not verifiable in small samples	×	N/A	×	×		×
Does not extend to complex settings such as time-dependent covariates and dynamic ¹⁶ 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

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.

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
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, sec. 5.3.3) (nlme correlation function names are at the right if the structure is implemented in nlme):

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}\) ¹⁷	`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)

¹⁷ 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]\)

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

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

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

7.8.1 Graphical Exploration of Data

Code

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

Code

# 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

Code

# 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

Code

# 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

Code

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

Code

# 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.

Code

require(nlme)
cp <- list(corCAR1,corExp,corCompSymm,corLin,corGaus,corSpher)
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]])

       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

AIC computed above is set up so that smaller values are best. From this the continuous-time AR1 and exponential structures are tied for the best. For the remainder of the analysis use corCAR1, 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.

Code

a <- Gls(twstrs ~ treat * rcs(week, 3) + rcs(twstrs0, 3) +
         rcs(age, 4) * sex, data=both,
         correlation=corCAR1(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 = corCAR1(form = ~week | 
    uid))


Obs 522	Log-restricted-likelihood -1756.95
Clusters 108	Model d.f. 17
g 11.334	σ 8.5917
	d.f. 504

	β	S.E.	t	Pr(>\|t\|)
Intercept	-0.3093	11.8804	-0.03	0.9792
treat=5000U	0.4344	2.5962	0.17	0.8672
treat=Placebo	7.1433	2.6133	2.73	0.0065
week	0.2879	0.2973	0.97	0.3334
week'	0.7313	0.3078	2.38	0.0179
twstrs0	0.8071	0.1449	5.57	<0.0001
twstrs0'	0.2129	0.1795	1.19	0.2360
age	-0.1178	0.2346	-0.50	0.6158
age'	0.6968	0.6484	1.07	0.2830
age''	-3.4018	2.5599	-1.33	0.1845
sex=M	24.2802	18.6208	1.30	0.1929
treat=5000U × week	0.0745	0.4221	0.18	0.8599
treat=Placebo × week	-0.1256	0.4243	-0.30	0.7674
treat=5000U × week'	-0.4389	0.4363	-1.01	0.3149
treat=Placebo × week'	-0.6459	0.4381	-1.47	0.1411
age × sex=M	-0.5846	0.4447	-1.31	0.1892
age' × sex=M	1.4652	1.2388	1.18	0.2375
age'' × sex=M	-4.0338	4.8123	-0.84	0.4023

Correlation Structure: Continuous AR(1)
 Formula: ~week | uid 
 Parameter estimate(s):
      Phi 
0.8666689

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

Code

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

Figure 7.4: Variogram, with assumed correlation pattern superimposed

Check constant variance and normality assumptions:

Code

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.5: 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.

Code

anova(a)

Wald Statistics for `twstrs`
	χ²	d.f.	P
treat (Factor+Higher Order Factors)	22.11	6	0.0012
All Interactions	14.94	4	0.0048
week (Factor+Higher Order Factors)	77.27	6	<0.0001
All Interactions	14.94	4	0.0048
Nonlinear (Factor+Higher Order Factors)	6.61	3	0.0852
twstrs0	233.83	2	<0.0001
Nonlinear	1.41	1	0.2354
age (Factor+Higher Order Factors)	9.68	6	0.1388
All Interactions	4.86	3	0.1826
Nonlinear (Factor+Higher Order Factors)	7.59	4	0.1077
sex (Factor+Higher Order Factors)	5.67	4	0.2252
All Interactions	4.86	3	0.1826
treat × week (Factor+Higher Order Factors)	14.94	4	0.0048
Nonlinear	2.27	2	0.3208
Nonlinear Interaction : f(A,B) vs. AB	2.27	2	0.3208
age × sex (Factor+Higher Order Factors)	4.86	3	0.1826
Nonlinear	3.76	2	0.1526
Nonlinear Interaction : f(A,B) vs. AB	3.76	2	0.1526
TOTAL NONLINEAR	15.03	8	0.0586
TOTAL INTERACTION	19.75	7	0.0061
TOTAL NONLINEAR + INTERACTION	28.54	11	0.0027
TOTAL	322.98	17	<0.0001

Code

plot(anova(a))

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

Code

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.7: Estimated effects of time, baseline `TWSTRS`, age, and sex

Code

summary(a)  # Shows for week 8

Effects Response: `twstrs`
	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
week	4	12	8	6.6910	1.1060	4.524	8.858
twstrs0	39	53	14	13.5500	0.8862	11.810	15.290
age	46	65	19	2.5030	2.0510	-1.518	6.523
treat --- 5000U:10000U	1	2		0.5917	1.9980	-3.325	4.508
treat --- Placebo:10000U	1	3		5.4930	2.0040	1.565	9.421
sex --- M:F	1	2		-1.0850	1.7790	-4.571	2.401

Code

# 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 Pr(>|z|)
1      2      46  56   F    -6.31 2.10 -10.43 -2.186 -3.00   0.0027
2      4      46  56   F    -5.91 1.82  -9.47 -2.349 -3.25   0.0011
3      8      46  56   F    -4.90 2.01  -8.85 -0.953 -2.43   0.0150
4*    12      46  56   F    -3.07 1.75  -6.49  0.361 -1.75   0.0795
5*    16      46  56   F    -1.02 2.10  -5.14  3.092 -0.49   0.6260

Redundant contrasts are denoted by *

Confidence intervals are 0.95 individual intervals

Code

# 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 Pr(>|z|)
1      2      46  56   F    -6.89 2.07 -10.96 -2.83 -3.32   0.0009
2      4      46  56   F    -6.64 1.79 -10.15 -3.13 -3.70   0.0002
3      8      46  56   F    -5.49 2.00  -9.42 -1.56 -2.74   0.0061
4*    12      46  56   F    -1.76 1.74  -5.17  1.65 -1.01   0.3109
5*    16      46  56   F     2.62 2.09  -1.47  6.71  1.25   0.2099

Redundant contrasts are denoted by *

Confidence intervals are 0.95 individual intervals

Code

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.8: 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.

Code

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.9: Nomogram from GLS fit. Second axis is the baseline score.

7.8.3 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

Code

require(rmsb)
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')
# 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 treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

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

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 0.816 0.748 0.77 0.739 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.001     1833     2060
2   alpha[2] 1.001     1438     2136
3   alpha[3] 1.003     1002     2010
4   alpha[4] 1.003      888     1792
5   alpha[5] 1.003      863     1655
6   alpha[6] 1.003      800     1394
7   alpha[7] 1.004      767     1539
8   alpha[8] 1.004      727     1333
9   alpha[9] 1.005      647     1019
10 alpha[10] 1.006      606     1188
11 alpha[11] 1.007      556     1027
12 alpha[12] 1.008      528      879
13 alpha[13] 1.008      508      836
14 alpha[14] 1.009      489      837
15 alpha[15] 1.010      468      770
16 alpha[16] 1.010      453      829
17 alpha[17] 1.010      442      786
18 alpha[18] 1.010      428      744
19 alpha[19] 1.011      407      765
20 alpha[20] 1.011      398      740
21 alpha[21] 1.011      390      799
22 alpha[22] 1.012      380      816
23 alpha[23] 1.012      382      804
24 alpha[24] 1.013      379      793
25 alpha[25] 1.012      388      862
26 alpha[26] 1.012      396      934
27 alpha[27] 1.012      397      860
28 alpha[28] 1.012      404      922
29 alpha[29] 1.013      403      933
30 alpha[30] 1.013      400      985
31 alpha[31] 1.012      409     1021
32 alpha[32] 1.010      424     1102
33 alpha[33] 1.011      429     1126
34 alpha[34] 1.012      421     1143
35 alpha[35] 1.012      446     1267
36 alpha[36] 1.010      480     1108
37 alpha[37] 1.010      519     1128
38 alpha[38] 1.011      554     1125
39 alpha[39] 1.011      587     1269
40 alpha[40] 1.010      609     1239
41 alpha[41] 1.009      675     1168
42 alpha[42] 1.010      747     1231
43 alpha[43] 1.009      878     1799
44 alpha[44] 1.008      944     1896
45 alpha[45] 1.006     1017     1992
46 alpha[46] 1.005     1037     2117
47 alpha[47] 1.004     1045     2196
48 alpha[48] 1.004     1078     2089
49 alpha[49] 1.003     1187     2467
50 alpha[50] 1.002     1299     2506
51 alpha[51] 1.001     1464     2516
52 alpha[52] 1.001     1520     2247
53 alpha[53] 1.002     1856     2807
54 alpha[54] 1.000     1944     2642
55 alpha[55] 1.001     2056     2726
56 alpha[56] 1.001     2122     2886
57 alpha[57] 1.000     2259     3302
58 alpha[58] 1.000     2186     2886
59 alpha[59] 1.000     2213     2786
60 alpha[60] 1.000     2400     2938
61 alpha[61] 1.000     2701     2575
62   beta[1] 1.001      858     1554
63   beta[2] 1.005      991     1765
64   beta[3] 1.001     2058     2701
65   beta[4] 1.000     4193     3362
66   beta[5] 1.003      763     1264
67   beta[6] 1.001      932     1506
68   beta[7] 1.005      874     1469
69   beta[8] 1.003      931     1691
70   beta[9] 1.002      945     1889
71  beta[10] 1.004      746     1185
72  beta[11] 1.000     4142     2897
73  beta[12] 1.001     3869     2825
74  beta[13] 1.001     4542     3080
75  beta[14] 1.002     4277     3004
76  beta[15] 1.002     1001     1597
77  beta[16] 1.006      750     1629
78  beta[17] 1.001     1100     1765
79 sigmag[1] 1.004      964     2141

Code

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
Obs 522	LOO log L -1746.69±23.79	g 3.852 [3.348, 4.364]	C 0.793 [0.784, 0.799]
Draws 4000	LOO IC 3493.39±47.59	g_p 0.435 [0.42, 0.45]	D_xy 0.586 [0.567, 0.598]
Chains 4	Effective p 179.06±7.94	EV 0.594 [0.546, 0.642]
Time 6.8s	B 0.149 [0.139, 0.16]	v 11.552 [8.791, 14.901]
p 17		vp 0.148 [0.135, 0.16]
Cluster on `uid`
Clusters 108
σ_γ 1.886 [1.5042, 2.2337]

	Mean β	Median β	S.E.	Lower	Upper	Pr(β>0)	Symmetry
y≥7	-1.7647	-1.6959	4.2708	-9.9855	6.6421	0.3368	1.00
y≥9	-2.8034	-2.7432	4.1455	-11.1474	5.0856	0.2445	0.98
y≥10	-3.9889	-3.8958	4.0876	-12.0870	3.8710	0.1628	0.97
y≥11	-4.4318	-4.3584	4.0793	-12.8734	3.0406	0.1360	0.97
y≥13	-4.6367	-4.5643	4.0754	-13.0770	2.8095	0.1248	0.97
y≥14	-4.9895	-4.9260	4.0750	-13.0793	2.8983	0.1102	0.96
y≥15	-5.2964	-5.2053	4.0771	-13.6845	2.3406	0.0930	0.98
y≥16	-5.6813	-5.5724	4.0820	-13.9316	2.0855	0.0785	0.96
y≥17	-6.4831	-6.4038	4.0717	-14.6176	1.3199	0.0515	0.96
y≥18	-6.7350	-6.6636	4.0696	-14.8838	0.9887	0.0470	0.96
y≥19	-7.0248	-6.9500	4.0690	-14.8893	1.0117	0.0388	0.96
y≥20	-7.2204	-7.1596	4.0674	-15.1670	0.7139	0.0348	0.97
y≥21	-7.4008	-7.3141	4.0659	-15.2731	0.5973	0.0320	0.96
y≥22	-7.8184	-7.7331	4.0688	-15.9628	-0.0909	0.0250	0.96
y≥23	-8.0849	-7.9820	4.0685	-15.9963	-0.1457	0.0208	0.96
y≥24	-8.3697	-8.2505	4.0705	-16.4004	-0.5330	0.0187	0.96
y≥25	-8.6289	-8.5086	4.0731	-16.6059	-0.7542	0.0147	0.96
y≥26	-9.0232	-8.9109	4.0766	-16.9433	-1.0599	0.0110	0.96
y≥27	-9.3152	-9.2130	4.0792	-17.2624	-1.3258	0.0098	0.96
y≥28	-9.5575	-9.4585	4.0804	-17.5376	-1.5986	0.0085	0.96
y≥29	-9.7927	-9.6829	4.0839	-18.2100	-2.2951	0.0075	0.96
y≥30	-10.0949	-9.9969	4.0857	-18.2703	-2.2738	0.0068	0.97
y≥31	-10.3847	-10.2777	4.0896	-18.5375	-2.5767	0.0053	0.96
y≥32	-10.5011	-10.3934	4.0906	-18.6534	-2.6876	0.0048	0.96
y≥33	-10.8656	-10.7823	4.0942	-19.1080	-3.1488	0.0040	0.96
y≥34	-11.1760	-11.0798	4.0980	-19.0144	-2.9934	0.0035	0.95
y≥35	-11.3990	-11.3176	4.0992	-19.8350	-3.8229	0.0035	0.96
y≥36	-11.6455	-11.5728	4.1008	-20.0436	-4.1122	0.0030	0.96
y≥37	-11.9232	-11.8413	4.1027	-20.2815	-4.2485	0.0018	0.96
y≥38	-12.1514	-12.0605	4.1033	-20.5559	-4.5531	0.0018	0.96
y≥39	-12.3980	-12.3241	4.1043	-20.7292	-4.7922	0.0018	0.96
y≥40	-12.5818	-12.5198	4.1054	-20.8705	-4.9402	0.0013	0.96
y≥41	-12.7666	-12.7100	4.1074	-21.0671	-5.1490	0.0013	0.95
y≥42	-13.0941	-13.0431	4.1092	-21.3128	-5.4084	0.0010	0.96
y≥43	-13.3227	-13.2659	4.1094	-21.5872	-5.6627	0.0010	0.95
y≥44	-13.6631	-13.6120	4.1132	-21.9258	-5.9704	0.0008	0.96
y≥45	-13.9824	-13.9105	4.1160	-22.2163	-6.2164	0.0005	0.95
y≥46	-14.2817	-14.2075	4.1157	-22.6467	-6.6222	0.0005	0.94
y≥47	-14.7041	-14.6384	4.1196	-23.0632	-7.0158	0.0003	0.95
y≥48	-14.9935	-14.9302	4.1228	-23.3128	-7.2711	0.0003	0.95
y≥49	-15.3681	-15.3035	4.1228	-23.7191	-7.6533	0.0003	0.95
y≥50	-15.6892	-15.6296	4.1255	-24.0342	-7.9849	0.0000	0.95
y≥51	-16.2168	-16.1576	4.1338	-24.4611	-8.3577	0.0000	0.94
y≥52	-16.5811	-16.5205	4.1383	-24.8637	-8.7343	0.0000	0.95
y≥53	-17.0314	-16.9719	4.1382	-25.3071	-9.1679	0.0000	0.94
y≥54	-17.5314	-17.4605	4.1373	-25.8438	-9.7846	0.0000	0.95
y≥55	-17.9428	-17.8766	4.1388	-25.7649	-9.7153	0.0000	0.95
y≥56	-18.1954	-18.1190	4.1406	-26.3139	-10.2436	0.0000	0.96
y≥57	-18.6678	-18.5867	4.1393	-26.9418	-10.9034	0.0000	0.95
y≥58	-19.2278	-19.1415	4.1426	-27.0661	-11.0320	0.0000	0.95
y≥59	-19.5804	-19.5015	4.1438	-27.3972	-11.4257	0.0000	0.93
y≥60	-19.9076	-19.8430	4.1451	-27.7742	-11.8254	0.0000	0.94
y≥61	-20.6049	-20.5241	4.1445	-28.4710	-12.4467	0.0000	0.94
y≥62	-20.9791	-20.9024	4.1494	-29.2339	-13.1191	0.0000	0.94
y≥63	-21.3863	-21.3057	4.1510	-29.8729	-13.7611	0.0000	0.94
y≥64	-21.5280	-21.4520	4.1503	-29.9242	-13.8288	0.0000	0.95
y≥65	-22.2525	-22.1612	4.1586	-30.5168	-14.5523	0.0000	0.93
y≥66	-22.6317	-22.5752	4.1615	-30.9429	-14.9052	0.0000	0.94
y≥67	-23.0497	-23.0046	4.1541	-31.2547	-15.1964	0.0000	0.94
y≥68	-23.8148	-23.7682	4.1832	-31.8628	-15.8025	0.0000	0.94
y≥71	-24.6983	-24.6118	4.2138	-33.4522	-17.0710	0.0000	0.96
treat=5000U	0.1269	0.1155	0.7040	-1.2744	1.4606	0.5712	1.01
treat=Placebo	2.3441	2.3496	0.7304	0.8907	3.7464	0.9990	0.98
week	0.1222	0.1226	0.0790	-0.0272	0.2795	0.9432	0.99
week'	0.1920	0.1911	0.0871	0.0265	0.3707	0.9830	0.95
twstrs0	0.2286	0.2285	0.0527	0.1331	0.3400	0.9998	1.02
twstrs0'	0.1293	0.1287	0.0648	0.0041	0.2544	0.9795	1.03
age	-0.0159	-0.0174	0.0798	-0.1738	0.1374	0.4168	1.04
age'	0.1967	0.1999	0.2186	-0.2417	0.6222	0.8202	0.97
age''	-1.0838	-1.0802	0.8542	-2.7619	0.5937	0.1015	1.06
sex=M	5.0446	5.0870	6.2370	-7.2400	16.6571	0.7832	0.98
treat=5000U × week	0.0512	0.0512	0.1081	-0.1729	0.2480	0.6862	0.99
treat=Placebo × week	-0.0523	-0.0507	0.1126	-0.2818	0.1552	0.3262	1.00
treat=5000U × week'	-0.1640	-0.1643	0.1192	-0.3998	0.0663	0.0845	1.01
treat=Placebo × week'	-0.1415	-0.1437	0.1228	-0.3868	0.1010	0.1280	0.99
age × sex=M	-0.1120	-0.1136	0.1483	-0.3789	0.1913	0.2288	1.03
age' × sex=M	0.1652	0.1716	0.4120	-0.6473	0.9593	0.6675	0.96
age'' × sex=M	-0.0286	-0.0627	1.5988	-3.1547	3.1811	0.4812	1.05

Code

a <- anova(bpo)
a

Relative Explained Variation for `twstrs`. Approximate total model Wald χ² used in denominators of REV:247.7 [199.3, 321.8].
	REV	Lower	Upper	d.f.
treat (Factor+Higher Order Factors)	0.137	0.076	0.231	6
All Interactions	0.096	0.038	0.181	4
week (Factor+Higher Order Factors)	0.588	0.463	0.683	6
All Interactions	0.096	0.038	0.181	4
Nonlinear (Factor+Higher Order Factors)	0.022	0.001	0.068	3
twstrs0	0.666	0.522	0.739	2
Nonlinear	0.016	0.000	0.047	1
age (Factor+Higher Order Factors)	0.027	0.012	0.095	6
All Interactions	0.016	0.000	0.059	3
Nonlinear (Factor+Higher Order Factors)	0.023	0.005	0.076	4
sex (Factor+Higher Order Factors)	0.019	0.002	0.070	4
All Interactions	0.016	0.000	0.059	3
treat × week (Factor+Higher Order Factors)	0.096	0.038	0.181	4
Nonlinear	0.009	0.000	0.040	2
Nonlinear Interaction : f(A,B) vs. AB	0.009	0.000	0.040	2
age × sex (Factor+Higher Order Factors)	0.016	0.000	0.059	3
Nonlinear	0.014	0.000	0.052	2
Nonlinear Interaction : f(A,B) vs. AB	0.014	0.000	0.052	2
TOTAL NONLINEAR	0.058	0.030	0.145	8
TOTAL INTERACTION	0.110	0.058	0.212	7
TOTAL NONLINEAR + INTERACTION	0.143	0.100	0.260	11
TOTAL	1.000	1.000	1.000	17

Code

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

Code

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 Pr(Contrast>0)
1  Week 2     2 -2.2395042 0.5947754 -3.4241087 -1.1069360         0.0000
2  Week 4     4 -2.1349275 0.5247592 -3.1827973 -1.1505865         0.0000
3  Week 8     8 -1.7842982 0.5792123 -2.9541916 -0.6459085         0.0008
4* Week 12   12 -0.8677657 0.5277858 -1.8638679  0.1898032         0.0510
5* Week 16   16  0.1902426 0.6115432 -0.9667285  1.4257672         0.6218

Redundant contrasts are denoted by *

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

Code

plot(k)

Code

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.

Code

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')

Code

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.4 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:

\[\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.

Code

# 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')
# 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 treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

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

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 1.026 0.966 0.953 1.002 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.001     4924     2630
2   alpha[2] 1.000     5723     2966
3   alpha[3] 1.000     5235     3081
4   alpha[4] 1.000     5421     2955
5   alpha[5] 1.000     5202     3199
6   alpha[6] 1.000     5125     3364
7   alpha[7] 1.001     5219     3449
8   alpha[8] 1.001     4966     3312
9   alpha[9] 1.001     4604     3169
10 alpha[10] 1.001     4391     3351
11 alpha[11] 1.001     4500     3389
12 alpha[12] 1.001     4389     3502
13 alpha[13] 1.000     4500     3725
14 alpha[14] 1.001     4153     3473
15 alpha[15] 1.001     4201     3480
16 alpha[16] 1.000     4122     3253
17 alpha[17] 1.000     4056     3513
18 alpha[18] 1.000     4087     3479
19 alpha[19] 1.001     4166     3517
20 alpha[20] 1.001     4164     3297
21 alpha[21] 1.000     4062     3360
22 alpha[22] 1.001     4242     3582
23 alpha[23] 1.000     4418     3200
24 alpha[24] 1.000     4464     2937
25 alpha[25] 1.001     4895     3404
26 alpha[26] 1.002     5273     3420
27 alpha[27] 1.002     5430     2891
28 alpha[28] 1.001     6314     2722
29 alpha[29] 1.002     6377     3173
30 alpha[30] 1.002     6832     3435
31 alpha[31] 1.000     7291     3022
32 alpha[32] 1.002     7379     3081
33 alpha[33] 1.001     7633     3004
34 alpha[34] 1.002     8702     2917
35 alpha[35] 1.003     8743     3246
36 alpha[36] 1.003     8590     3033
37 alpha[37] 1.002     8372     3261
38 alpha[38] 1.000     7809     3338
39 alpha[39] 1.002     7081     2946
40 alpha[40] 1.002     6466     3100
41 alpha[41] 1.001     6343     3239
42 alpha[42] 1.001     6182     3478
43 alpha[43] 1.000     5851     3302
44 alpha[44] 1.000     5779     3540
45 alpha[45] 1.000     5319     3666
46 alpha[46] 1.001     5001     3593
47 alpha[47] 1.000     4841     3825
48 alpha[48] 0.999     4721     3886
49 alpha[49] 1.000     4809     3728
50 alpha[50] 1.001     5010     3442
51 alpha[51] 1.000     4944     3677
52 alpha[52] 1.001     5216     3474
53 alpha[53] 1.001     5312     3610
54 alpha[54] 1.000     5526     3762
55 alpha[55] 1.000     5303     3486
56 alpha[56] 1.000     5299     3816
57 alpha[57] 1.000     5508     3467
58 alpha[58] 0.999     5808     3593
59 alpha[59] 0.999     5895     3332
60 alpha[60] 1.002     5698     3028
61 alpha[61] 1.002     5788     3618
62   beta[1] 1.001    10368     2769
63   beta[2] 1.001     9238     2290
64   beta[3] 1.001     6724     3024
65   beta[4] 0.999     9205     2695
66   beta[5] 1.000     3601     3330
67   beta[6] 1.000     6792     3560
68   beta[7] 1.001     7427     3102
69   beta[8] 1.000     9182     2600
70   beta[9] 1.001     9970     2986
71  beta[10] 1.003    11084     2966
72  beta[11] 1.002    10443     2604
73  beta[12] 1.000    11349     2615
74  beta[13] 1.002     8442     2660
75  beta[14] 1.000     9386     2703
76  beta[15] 1.001    11514     2825
77  beta[16] 1.000    10240     3115
78  beta[17] 1.002    10828     3104
79  beta[18] 1.002    11247     2599

Code

stanDxplot(bmark)

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

Code

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
Obs 517	LOO log L -1786.11±22.27	g 3.267 [2.995, 3.56]	C 0.828 [0.825, 0.83]
Draws 4000	LOO IC 3572.22±44.54	g_p 0.416 [0.402, 0.43]	D_xy 0.656 [0.65, 0.661]
Chains 4	Effective p 89.99±4.7	EV 0.533 [0.491, 0.573]
Time 3.6s	B 0.117 [0.113, 0.12]	v 8.402 [6.973, 9.847]
p 18		vp 0.133 [0.123, 0.144]

	Mode β	Mean β	Median β	S.E.	Lower	Upper	Pr(β>0)	Symmetry
treat=5000U	0.2210	0.2221	0.2194	0.5808	-0.8371	1.4448	0.6435	1.01
treat=Placebo	1.8315	1.8437	1.8480	0.5861	0.6460	2.9468	0.9998	1.01
week	0.4865	0.4891	0.4882	0.0862	0.3094	0.6498	1.0000	1.01
week'	-0.2878	-0.2899	-0.2897	0.0914	-0.4751	-0.1177	0.0008	0.99
ptwstrs	0.1998	0.2014	0.2014	0.0268	0.1481	0.2520	1.0000	1.03
ptwstrs'	-0.0621	-0.0652	-0.0655	0.0650	-0.1969	0.0587	0.1578	0.98
ptwstrs''	0.5329	0.5449	0.5440	0.2598	0.0345	1.0414	0.9825	1.01
age	-0.0295	-0.0286	-0.0287	0.0317	-0.0903	0.0318	0.1880	1.00
age'	0.1237	0.1208	0.1209	0.0873	-0.0556	0.2800	0.9128	1.02
age''	-0.5070	-0.4968	-0.4948	0.3416	-1.1775	0.1460	0.0717	1.00
sex=M	-0.4644	-0.4281	-0.4606	2.4046	-5.0914	4.1596	0.4260	1.00
treat=5000U × week	-0.0341	-0.0342	-0.0340	0.1126	-0.2570	0.1762	0.3845	0.97
treat=Placebo × week	-0.2719	-0.2744	-0.2754	0.1170	-0.4856	-0.0293	0.0105	0.94
treat=5000U × week'	-0.0341	-0.0341	-0.0342	0.1201	-0.2612	0.2017	0.3872	1.04
treat=Placebo × week'	0.1197	0.1218	0.1200	0.1268	-0.1224	0.3655	0.8310	1.03
age × sex=M	0.0112	0.0103	0.0112	0.0577	-0.1025	0.1225	0.5720	0.99
age' × sex=M	-0.0511	-0.0476	-0.0490	0.1608	-0.3618	0.2725	0.3735	1.01
age'' × sex=M	0.2618	0.2488	0.2520	0.6210	-0.9583	1.4411	0.6640	1.00

Code

a <- anova(bpo)
a

Relative Explained Variation for `twstrs`. Approximate total model Wald χ² used in denominators of REV:247.7 [204.9, 333.7].
	REV	Lower	Upper	d.f.
treat (Factor+Higher Order Factors)	0.137	0.080	0.229	6
All Interactions	0.096	0.046	0.168	4
week (Factor+Higher Order Factors)	0.588	0.461	0.675	6
All Interactions	0.096	0.046	0.168	4
Nonlinear (Factor+Higher Order Factors)	0.022	0.001	0.072	3
twstrs0	0.666	0.534	0.742	2
Nonlinear	0.016	0.000	0.050	1
age (Factor+Higher Order Factors)	0.027	0.010	0.094	6
All Interactions	0.016	0.000	0.061	3
Nonlinear (Factor+Higher Order Factors)	0.023	0.006	0.077	4
sex (Factor+Higher Order Factors)	0.019	0.003	0.075	4
All Interactions	0.016	0.000	0.061	3
treat × week (Factor+Higher Order Factors)	0.096	0.046	0.168	4
Nonlinear	0.009	0.000	0.042	2
Nonlinear Interaction : f(A,B) vs. AB	0.009	0.000	0.042	2
age × sex (Factor+Higher Order Factors)	0.016	0.000	0.061	3
Nonlinear	0.014	0.000	0.053	2
Nonlinear Interaction : f(A,B) vs. AB	0.014	0.000	0.053	2
TOTAL NONLINEAR	0.058	0.031	0.151	8
TOTAL INTERACTION	0.110	0.062	0.205	7
TOTAL NONLINEAR + INTERACTION	0.143	0.095	0.256	11
TOTAL	1.000	1.000	1.000	17

Code

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.

Code

bmarkre <- blrm(twstrs ~  treat * rcs(week, 3) + rcs(ptwstrs, 4) +
                          rcs(age, 4) * sex + cluster(uid),
                data=both, file='bmarkre.rds')
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 treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
1 of 4000 (0.03%) 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.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete.
Divergent samples: 0 1 0 0 

EBFMI: 0.925 0.9 0.899 1.032 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.002     3244     2049
2   alpha[2] 1.001     3567     2602
3   alpha[3] 1.000     3194     2506
4   alpha[4] 1.000     2803     2688
5   alpha[5] 1.000     2749     2696
6   alpha[6] 1.000     2609     2901
7   alpha[7] 1.000     2552     2683
8   alpha[8] 1.000     2455     2459
9   alpha[9] 1.001     2214     2305
10 alpha[10] 1.001     2169     2405
11 alpha[11] 1.001     2181     2463
12 alpha[12] 1.001     2110     2283
13 alpha[13] 1.001     2087     2344
14 alpha[14] 1.001     2093     2441
15 alpha[15] 1.001     2092     1927
16 alpha[16] 1.001     2115     2375
17 alpha[17] 1.002     2092     2322
18 alpha[18] 1.003     2071     2365
19 alpha[19] 1.003     2016     2635
20 alpha[20] 1.003     2025     2769
21 alpha[21] 1.002     2072     2620
22 alpha[22] 1.002     2034     2706
23 alpha[23] 1.002     2056     2831
24 alpha[24] 1.002     2080     2703
25 alpha[25] 1.001     2355     2765
26 alpha[26] 1.001     2531     2579
27 alpha[27] 1.002     2585     2912
28 alpha[28] 1.001     2761     2552
29 alpha[29] 1.001     2840     2630
30 alpha[30] 1.001     3094     2653
31 alpha[31] 1.000     3116     2384
32 alpha[32] 1.000     3256     2416
33 alpha[33] 1.000     3290     2873
34 alpha[34] 1.000     3639     2618
35 alpha[35] 1.001     4027     3058
36 alpha[36] 1.000     4360     2935
37 alpha[37] 1.001     4670     3190
38 alpha[38] 1.001     4883     3300
39 alpha[39] 1.002     5106     3039
40 alpha[40] 1.001     5200     3060
41 alpha[41] 1.001     4671     3524
42 alpha[42] 1.000     4595     3088
43 alpha[43] 1.000     4454     3075
44 alpha[44] 1.001     3919     2960
45 alpha[45] 1.000     3731     3059
46 alpha[46] 1.000     3447     3335
47 alpha[47] 1.000     3591     3022
48 alpha[48] 1.000     3585     3168
49 alpha[49] 1.002     3737     2989
50 alpha[50] 1.000     3596     2898
51 alpha[51] 1.001     3657     3138
52 alpha[52] 1.001     3639     3249
53 alpha[53] 1.000     3493     2963
54 alpha[54] 1.000     3062     3090
55 alpha[55] 1.001     2927     3290
56 alpha[56] 1.001     3137     3106
57 alpha[57] 1.001     3307     3092
58 alpha[58] 1.001     3162     3232
59 alpha[59] 1.001     3252     3262
60 alpha[60] 1.001     3595     3107
61 alpha[61] 1.001     4092     3336
62   beta[1] 1.001     5219     2725
63   beta[2] 1.002     5380     2739
64   beta[3] 1.000     5055     2736
65   beta[4] 1.001     5871     2746
66   beta[5] 1.002     2034     3080
67   beta[6] 1.001     4348     3043
68   beta[7] 1.003     4437     3030
69   beta[8] 1.001     5679     2571
70   beta[9] 1.002     5978     3107
71  beta[10] 1.000     5561     2937
72  beta[11] 1.001     4398     2969
73  beta[12] 1.001     6881     3130
74  beta[13] 1.001     6225     2707
75  beta[14] 1.000     5381     2674
76  beta[15] 1.001     6027     3053
77  beta[16] 1.001     5592     2872
78  beta[17] 1.001     5026     2985
79  beta[18] 1.001     6155     3098
80 sigmag[1] 1.003     1273     1622

Code

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
Obs 517	LOO log L -1787.2±22.47	g 3.249 [2.968, 3.521]	C 0.828 [0.825, 0.83]
Draws 4000	LOO IC 3574.41±44.95	g_p 0.415 [0.403, 0.428]	D_xy 0.656 [0.65, 0.66]
Chains 4	Effective p 94.26±5	EV 0.531 [0.491, 0.567]
Time 5.3s	B 0.117 [0.113, 0.121]	v 8.309 [6.795, 9.7]
p 18		vp 0.133 [0.123, 0.142]
Cluster on `uid`
Clusters 108
σ_γ 0.1155 [1e-04, 0.3447]

	Mean β	Median β	S.E.	Lower	Upper	Pr(β>0)	Symmetry
treat=5000U	0.2185	0.2254	0.5671	-0.8612	1.3535	0.6558	1.02
treat=Placebo	1.8307	1.8308	0.5757	0.7178	2.9904	0.9990	1.01
week	0.4856	0.4856	0.0831	0.3293	0.6514	1.0000	0.99
week'	-0.2844	-0.2853	0.0872	-0.4510	-0.1142	0.0005	1.02
ptwstrs	0.2004	0.2003	0.0273	0.1473	0.2546	1.0000	0.98
ptwstrs'	-0.0651	-0.0663	0.0638	-0.1920	0.0549	0.1560	1.02
ptwstrs''	0.5448	0.5502	0.2531	0.0180	1.0140	0.9855	0.98
age	-0.0291	-0.0293	0.0319	-0.0872	0.0354	0.1830	1.00
age'	0.1239	0.1244	0.0888	-0.0442	0.3022	0.9228	1.01
age''	-0.5112	-0.5101	0.3514	-1.2099	0.1686	0.0698	0.99
sex=M	-0.4264	-0.4685	2.4366	-5.2520	4.2666	0.4232	1.00
treat=5000U × week	-0.0326	-0.0318	0.1079	-0.2334	0.1817	0.3855	1.00
treat=Placebo × week	-0.2705	-0.2712	0.1118	-0.4773	-0.0403	0.0073	0.98
treat=5000U × week'	-0.0367	-0.0379	0.1157	-0.2693	0.1862	0.3748	1.01
treat=Placebo × week'	0.1170	0.1160	0.1206	-0.1094	0.3625	0.8352	1.00
age × sex=M	0.0102	0.0119	0.0587	-0.1055	0.1237	0.5715	1.01
age' × sex=M	-0.0490	-0.0506	0.1663	-0.3579	0.2901	0.3777	1.01
age'' × sex=M	0.2585	0.2616	0.6442	-0.9888	1.5106	0.6600	0.98

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:

Code

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.

Code

ggplot(Predict(bmark))

Code

ggplot(Predict(bmark, week, treat))

Code

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

           week   Contrast      S.E.      Lower      Upper Pr(Contrast>0)
1  Week 2     2 -1.2949780 0.3894362 -2.0466496 -0.5127775         0.0000
2  Week 4     4 -0.7462096 0.2633752 -1.2528826 -0.2128366         0.0022
3  Week 8     8  0.2294987 0.3714486 -0.4870849  0.9864448         0.7312
4* Week 12   12  0.7178925 0.2645893  0.1918738  1.2360777         0.9975
5* Week 16   16  1.0844575 0.3890726  0.3315238  1.8575108         0.9972

Redundant contrasts are denoted by *

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

Code

plot(k)

Code

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.

Code

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

Code

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

Code

# 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

Code

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

Code

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

Code

# 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

Code

# 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

Code

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

Code

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.9 Study Questions

Section 7.2

When should one model the time-response profile using discrete time?

Section 7.3

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

Section 7.4

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

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

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.

7.1 Notation

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

7.2.1 Common Basis Functions

7.2.2 Model for Mean Profile

7.2.3 Model Specification for Treatment Comparisons

7.3 Modeling Within-Subject Dependence

7.4 Parameter Estimation Procedure

7.5 Common Correlation Structures

7.6 Checking Model Fit

7.7 R Software

7.8 Case Study

7.8.1 Graphical Exploration of Data

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

7.8.2 Using Generalized Least Squares

7.8.3 Bayesian Proportional Odds Random Effects Model

7.8.4 Bayesian Markov Semiparametric Model

7.9 Study Questions

7.7 `R` Software