Skip to main content

7  Modeling Longitudinal Responses using Generalized Least Squares and Semiparametric Models

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

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

E

1 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 × × × × ?

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

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

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

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

6 For full efficiency, if using the working independence model

7 Or requires the user to specify one

8 For full efficiency of regression coefficient estimates

9 Unless the last observation is missing

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

11 Unless one knows how to properly do a weighted analysis

12 Or users population averages

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

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

15 If correction not applied

16 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

17 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

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.

S
Code 
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.
Code 
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
Code 
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.

Code 
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\).
Code 
v <- Variogram(a, form=~ week | uid)
plot(v)
Figure 7.5: Variogram, with assumed correlation pattern superimposed

Check constant variance and normality assumptions:

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

Code 
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
Code 
plot(anova(a))
Figure 7.7: 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.8: Estimated effects of time, baseline TWSTRS, age, and sex
Code 
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
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 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
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 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
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.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
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.10: Nomogram from GLS fit. Second axis is the baseline score.

7.8.3 Semiparametric Random Effects Models

  • Ordinal semiparametric model does away with distributional assumptions
  • Within-subject dependence can be accounted for with a Markov structure, (see later), GEE, or random effects
  • Semiparametric models’ coefficient estimates don’t change if \(Y\) is monotonically transformed
  • Start with a proportional odds (PO) model
  • Add random subject-specific intercepts to the model
  • Would have been a good fit to the correlation pattern had it been compound symmetric
  • Compare deviance of PO model to deviances of all available link functions to see which link fits best
Code 
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)
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)
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
Code 
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
  • Stay with the PO model
  • SD of random effects \(\sigma_\gamma\) is fairly large on the logit scale, indicating within-person correlation
  • orm function also fits a random effects structure with two \(\sigma\) parameters
  • \(u\) is a \(N(0,1)\) random effect
  • Within-subject time weighting function \(w(t)\)
  • Random effect contribution is \((1 - w(t))\sigma_{1}\times u + w(t)\sigma_{2}\times u\)
  • Take \(w(t)\) to be the fraction of elapsed time between the first and last follow-up visits
Code 
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)
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)
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
Code 
deviance(f)
                 intercepts                intercepts+x 
                   4078.936                    3571.416 
  intercepts+random effects intercepts+x+random effects 
                   3698.472                    3414.225
  • Deviance did not improve since \(\sigma_1\) and \(\sigma_2\) are close
  • If random effects are used, they can be the single-\(\sigma\) variety.

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
Code 
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 1 finished in 4.2 seconds.
Chain 3 finished in 4.1 seconds.
Chain 4 finished in 4.1 seconds.
Chain 2 finished in 4.5 seconds.

All 4 chains finished successfully.
Mean chain execution time: 4.2 seconds.
Total execution time: 4.6 seconds.
Code 
# 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.

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

Processing complete, no problems detected.

EBFMI: 0.763 0.825 0.735 0.724 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.001     1092     1374
2   alpha[2] 1.002      763     1473
3   alpha[3] 1.003      688     1652
4   alpha[4] 1.003      604     1412
5   alpha[5] 1.004      535     1244
6   alpha[6] 1.005      508      978
7   alpha[7] 1.005      496     1128
8   alpha[8] 1.006      487     1043
9   alpha[9] 1.006      506     1019
10 alpha[10] 1.006      491     1075
11 alpha[11] 1.007      482      975
12 alpha[12] 1.008      481     1014
13 alpha[13] 1.008      470     1039
14 alpha[14] 1.008      478     1161
15 alpha[15] 1.008      487     1083
16 alpha[16] 1.008      484      996
17 alpha[17] 1.008      475     1091
18 alpha[18] 1.008      466     1075
19 alpha[19] 1.008      459     1197
20 alpha[20] 1.008      459     1211
21 alpha[21] 1.009      451     1186
22 alpha[22] 1.009      453     1179
23 alpha[23] 1.007      454     1027
24 alpha[24] 1.007      444     1008
25 alpha[25] 1.006      453     1047
26 alpha[26] 1.006      457      975
27 alpha[27] 1.007      443      922
28 alpha[28] 1.007      437      972
29 alpha[29] 1.007      429      949
30 alpha[30] 1.006      428      949
31 alpha[31] 1.006      443      959
32 alpha[32] 1.006      442      873
33 alpha[33] 1.006      436      847
34 alpha[34] 1.006      430      914
35 alpha[35] 1.006      434     1004
36 alpha[36] 1.005      445      877
37 alpha[37] 1.005      444      839
38 alpha[38] 1.005      460      893
39 alpha[39] 1.005      469      847
40 alpha[40] 1.006      474      929
41 alpha[41] 1.006      491      975
42 alpha[42] 1.005      489     1046
43 alpha[43] 1.006      495      986
44 alpha[44] 1.006      515      986
45 alpha[45] 1.006      530     1097
46 alpha[46] 1.006      541     1170
47 alpha[47] 1.006      567     1287
48 alpha[48] 1.007      576     1214
49 alpha[49] 1.007      608     1400
50 alpha[50] 1.007      655     1326
51 alpha[51] 1.008      677     1504
52 alpha[52] 1.007      707     1493
53 alpha[53] 1.007      749     1405
54 alpha[54] 1.006      767     1728
55 alpha[55] 1.005      765     1593
56 alpha[56] 1.004      780     1470
57 alpha[57] 1.004      854     1536
58 alpha[58] 1.004      924     1884
59 alpha[59] 1.003      981     1896
60 alpha[60] 1.003     1101     2010
61 alpha[61] 1.004     1446     2397
62   beta[1] 1.002      862     1586
63   beta[2] 1.002      781     1535
64   beta[3] 1.001     1630     2600
65   beta[4] 1.000     4119     3084
66   beta[5] 1.007      681     1136
67   beta[6] 1.003      853     1674
68   beta[7] 1.007      811     1653
69   beta[8] 1.002      938     1464
70   beta[9] 1.009      891     1554
71  beta[10] 1.003      668     1723
72  beta[11] 1.002     4030     2399
73  beta[12] 1.001     3566     3060
74  beta[13] 1.000     4246     3035
75  beta[14] 1.002     3110     2354
76  beta[15] 1.002      927     1887
77  beta[16] 1.002      814     1844
78  beta[17] 1.002     1009     1797
79 sigmag[1] 1.003      747     1402
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
Obs522 LOO log L-1747.04±23.74 g3.829 [3.361, 4.374] C0.793 [0.786, 0.799]
Draws4000 LOO IC3494.07±47.49 gp0.435 [0.42, 0.449] Dxy0.585 [0.572, 0.598]
Chains4 Effective p179.57±7.91 EV0.592 [0.54, 0.638]
Time5.4s B0.149 [0.14, 0.159] v11.424 [8.357, 14.255]
p17 vp0.148 [0.134, 0.159]
Cluster onuid
Clusters108
σγ1.8894 [1.5159, 2.2534]
Mean β Median β S.E. Lower Upper P(β>0) Symmetry
y≥7   -1.7840   -1.6881  4.2337  -10.1482   6.6016  0.3402  1.01
y≥9   -2.8165   -2.7243  4.1580  -10.9672   5.4398  0.2438  0.98
y≥10   -3.9937   -3.9356  4.1069  -12.1579   4.0122  0.1630  0.97
y≥11   -4.4426   -4.3819  4.0930  -12.3525   3.8732  0.1355  0.95
y≥13   -4.6408   -4.5470  4.0862  -12.7657   3.3775  0.1227  0.96
y≥14   -4.9897   -4.9040  4.0761  -12.9207   3.1452  0.1050  0.97
y≥15   -5.2980   -5.2404  4.0777  -13.0792   3.0335  0.0950  0.97
y≥16   -5.6802   -5.6394  4.0787  -13.4537   2.6640  0.0800  0.97
y≥17   -6.4844   -6.4398  4.0698  -14.0868   1.9987  0.0528  0.97
y≥18   -6.7386   -6.6574  4.0693  -14.8344   1.1937  0.0468  0.96
y≥19   -7.0257   -6.9674  4.0676  -14.5768   1.4951  0.0402  0.98
y≥20   -7.2222   -7.1717  4.0678  -15.1630   0.9509  0.0355  0.97
y≥21   -7.4053   -7.3575  4.0680  -16.1055   0.0311  0.0322  0.97
y≥22   -7.8216   -7.7576  4.0683  -15.9216   0.1765  0.0272  0.96
y≥23   -8.0895   -8.0272  4.0707  -16.2631   -0.1560  0.0222  0.97
y≥24   -8.3737   -8.3282  4.0739  -16.6097   -0.5103  0.0198  0.97
y≥25   -8.6325   -8.6042  4.0731  -16.9457   -0.8140  0.0152  0.97
y≥26   -9.0270   -8.9584  4.0776  -17.1449   -1.0222  0.0118  0.97
y≥27   -9.3167   -9.2514  4.0794  -17.4582   -1.3196  0.0095  0.96
y≥28   -9.5614   -9.5231  4.0796  -17.7774   -1.6442  0.0088  0.96
y≥29   -9.7955   -9.7486  4.0805  -17.9787   -1.8292  0.0082  0.97
y≥30  -10.1003  -10.0752  4.0813  -17.7735   -1.6111  0.0063  0.96
y≥31  -10.3922  -10.3655  4.0849  -18.8726   -2.6706  0.0060  0.95
y≥32  -10.5083  -10.4786  4.0866  -18.7625   -2.5987  0.0053  0.95
y≥33  -10.8735  -10.8233  4.0896  -19.1521   -2.9278  0.0035  0.96
y≥34  -11.1842  -11.1411  4.0921  -19.4742   -3.2303  0.0032  0.96
y≥35  -11.4070  -11.3668  4.0934  -19.8830   -3.6724  0.0020  0.96
y≥36  -11.6527  -11.6166  4.0951  -20.0615   -3.8228  0.0018  0.96
y≥37  -11.9291  -11.8921  4.0966  -20.2361   -3.9993  0.0015  0.96
y≥38  -12.1573  -12.1244  4.0976  -20.4714   -4.2431  0.0013  0.96
y≥39  -12.4035  -12.3777  4.0987  -20.6167   -4.3564  0.0010  0.96
y≥40  -12.5885  -12.5644  4.1004  -20.7723   -4.5303  0.0010  0.95
y≥41  -12.7725  -12.7418  4.1014  -20.9607   -4.7477  0.0010  0.95
y≥42  -13.0979  -13.0776  4.1048  -21.4048   -5.1636  0.0010  0.95
y≥43  -13.3250  -13.3097  4.1062  -21.5430   -5.3319  0.0010  0.95
y≥44  -13.6650  -13.6359  4.1091  -21.8548   -5.5684  0.0008  0.94
y≥45  -13.9831  -13.9564  4.1114  -22.2197   -5.9310  0.0008  0.95
y≥46  -14.2827  -14.2367  4.1119  -22.6313   -6.3108  0.0008  0.94
y≥47  -14.7008  -14.6718  4.1136  -23.0098   -6.7562  0.0003  0.94
y≥48  -14.9911  -14.9501  4.1162  -23.2795   -6.9599  0.0003  0.94
y≥49  -15.3659  -15.3285  4.1203  -23.6278   -7.3098  0.0000  0.95
y≥50  -15.6849  -15.6543  4.1213  -23.9363   -7.6387  0.0000  0.95
y≥51  -16.2153  -16.1648  4.1252  -24.2969   -8.0110  0.0000  0.95
y≥52  -16.5795  -16.5294  4.1276  -24.7343   -8.4227  0.0000  0.95
y≥53  -17.0235  -16.9888  4.1297  -25.2565   -8.8877  0.0000  0.95
y≥54  -17.5258  -17.4819  4.1345  -26.0008   -9.6416  0.0000  0.94
y≥55  -17.9367  -17.8956  4.1358  -26.2597   -9.9551  0.0000  0.95
y≥56  -18.1892  -18.1380  4.1340  -26.3531  -10.0251  0.0000  0.94
y≥57  -18.6593  -18.6007  4.1370  -27.1446  -10.8498  0.0000  0.93
y≥58  -19.2198  -19.1340  4.1383  -27.4004  -11.0919  0.0000  0.93
y≥59  -19.5710  -19.4900  4.1389  -28.1815  -11.8497  0.0000  0.94
y≥60  -19.9007  -19.8356  4.1416  -28.3503  -11.9975  0.0000  0.94
y≥61  -20.5937  -20.4819  4.1574  -29.0448  -12.6372  0.0000  0.94
y≥62  -20.9658  -20.8511  4.1642  -29.5469  -13.1372  0.0000  0.93
y≥63  -21.3796  -21.2857  4.1787  -30.0321  -13.5796  0.0000  0.93
y≥64  -21.5205  -21.3981  4.1795  -29.8753  -13.4147  0.0000  0.94
y≥65  -22.2490  -22.1126  4.1943  -30.3133  -13.7647  0.0000  0.93
y≥66  -22.6269  -22.4614  4.1972  -30.9833  -14.5128  0.0000  0.93
y≥67  -23.0325  -22.8612  4.2043  -31.3609  -14.7547  0.0000  0.93
y≥68  -23.8077  -23.6582  4.2141  -32.1164  -15.5010  0.0000  0.94
y≥71  -24.6600  -24.5562  4.2321  -33.2293  -16.5221  0.0000  0.92
treat=5000U   0.0814   0.0846  0.7259   -1.3321   1.4545  0.5480  1.04
treat=Placebo   2.3468   2.3412  0.7362   0.9717   3.8018  0.9992  0.99
week   0.1210   0.1194  0.0781   -0.0477   0.2632  0.9415  0.96
week'   0.1922   0.1919  0.0859   0.0144   0.3449  0.9892  1.02
twstrs0   0.2281   0.2280  0.0510   0.1294   0.3281  1.0000  1.08
twstrs0'   0.1291   0.1290  0.0632   0.0062   0.2494  0.9785  0.98
age   -0.0150   -0.0146  0.0806   -0.1811   0.1360  0.4310  0.97
age'   0.1927   0.1912  0.2184   -0.2314   0.6040  0.8082  1.01
age''   -1.0716   -1.0649  0.8531   -2.7887   0.4831  0.1032  0.98
sex=M   5.1599   5.1891  6.2577   -7.4959   16.7992  0.7935  1.01
treat=5000U × week   0.0529   0.0545  0.1124   -0.1699   0.2780  0.6820  0.98
treat=Placebo × week   -0.0537   -0.0553  0.1128   -0.2797   0.1607  0.3127  1.00
treat=5000U × week'   -0.1642   -0.1649  0.1227   -0.4084   0.0693  0.0963  1.01
treat=Placebo × week'   -0.1388   -0.1375  0.1243   -0.3763   0.1039  0.1330  0.97
age × sex=M   -0.1127   -0.1134  0.1499   -0.4002   0.1828  0.2240  1.00
age' × sex=M   0.1528   0.1575  0.4203   -0.6470   0.9985  0.6418  1.00
age'' × sex=M   0.0410   0.0464  1.6383   -3.0841   3.2321  0.5122  0.98
Code 
a <- anova(bpo)
a

Relative Explained Variation for twstrs. Approximate total model Wald χ2 used in denominators of REV:252 [203.8, 326.6].

REV Lower Upper d.f.
treat (Factor+Higher Order Factors) 0.125 0.063 0.213 6
 All Interactions 0.087 0.036 0.172 4
week (Factor+Higher Order Factors) 0.584 0.457 0.689 6
 All Interactions 0.087 0.036 0.172 4
 Nonlinear (Factor+Higher Order Factors) 0.022 0.004 0.070 3
twstrs0 0.642 0.493 0.722 2
 Nonlinear 0.017 0.000 0.052 1
age (Factor+Higher Order Factors) 0.025 0.009 0.088 6
 All Interactions 0.016 0.001 0.054 3
 Nonlinear (Factor+Higher Order Factors) 0.022 0.003 0.075 4
sex (Factor+Higher Order Factors) 0.020 0.002 0.067 4
 All Interactions 0.016 0.001 0.054 3
treat × week (Factor+Higher Order Factors) 0.087 0.036 0.172 4
 Nonlinear 0.008 0.000 0.040 2
 Nonlinear Interaction : f(A,B) vs. AB 0.008 0.000 0.040 2
age × sex (Factor+Higher Order Factors) 0.016 0.001 0.054 3
 Nonlinear 0.014 0.000 0.048 2
 Nonlinear Interaction : f(A,B) vs. AB 0.014 0.000 0.048 2
TOTAL NONLINEAR 0.056 0.027 0.137 8
TOTAL INTERACTION 0.101 0.047 0.194 7
TOTAL NONLINEAR + INTERACTION 0.131 0.080 0.249 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
X
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 P(Contrast>0)
1  Week 2     2 -2.2394091 0.6019796 -3.419643 -1.0810237        0.0003
2  Week 4     4 -2.1320222 0.5334178 -3.185103 -1.0913020        0.0003
3  Week 8     8 -1.7784791 0.5865100 -2.941606 -0.6597846        0.0015
4* Week 12   12 -0.8698584 0.5299549 -1.887384  0.1802347        0.0475
5* Week 16   16  0.1775317 0.6178701 -1.060459  1.3657964        0.6192

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

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')
Running MCMC with 4 chains, at most 9 in parallel...

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

All 4 chains finished successfully.
Mean chain execution time: 2.4 seconds.
Total execution time: 2.6 seconds.
Code 
# 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: 1.003 0.856 1.041 0.993 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.002     2907     2470
2   alpha[2] 1.000     2375     2710
3   alpha[3] 1.001     1952     3111
4   alpha[4] 1.002     1747     2846
5   alpha[5] 1.002     1628     2533
6   alpha[6] 1.001     1670     2655
7   alpha[7] 1.002     1647     2616
8   alpha[8] 1.001     1809     2523
9   alpha[9] 1.001     2037     2700
10 alpha[10] 1.001     1952     3193
11 alpha[11] 1.001     1889     2734
12 alpha[12] 1.001     2012     3087
13 alpha[13] 1.001     1955     3009
14 alpha[14] 1.002     2202     2749
15 alpha[15] 1.001     2247     2879
16 alpha[16] 1.001     2381     2942
17 alpha[17] 1.001     2371     3248
18 alpha[18] 1.000     2670     2695
19 alpha[19] 1.000     2741     2718
20 alpha[20] 1.001     2965     2586
21 alpha[21] 1.001     3237     2985
22 alpha[22] 1.001     3566     2932
23 alpha[23] 1.000     3956     3151
24 alpha[24] 1.000     4053     2789
25 alpha[25] 1.000     4490     3264
26 alpha[26] 1.000     5417     3044
27 alpha[27] 1.000     5616     3224
28 alpha[28] 1.000     5356     2996
29 alpha[29] 1.000     6161     3303
30 alpha[30] 1.000     6543     3567
31 alpha[31] 1.001     6444     3369
32 alpha[32] 1.004     6639     3142
33 alpha[33] 1.002     6567     3760
34 alpha[34] 1.002     7186     3220
35 alpha[35] 1.000     7441     3479
36 alpha[36] 1.000     7829     3362
37 alpha[37] 1.000     7203     3124
38 alpha[38] 1.000     6772     3390
39 alpha[39] 1.001     6549     3132
40 alpha[40] 1.001     5729     3576
41 alpha[41] 1.001     5926     3225
42 alpha[42] 1.000     5890     2515
43 alpha[43] 1.000     5587     3164
44 alpha[44] 1.000     5136     2708
45 alpha[45] 1.001     5240     3415
46 alpha[46] 1.001     4949     3418
47 alpha[47] 1.002     4638     3228
48 alpha[48] 1.001     4581     3196
49 alpha[49] 1.000     4525     3490
50 alpha[50] 1.000     4452     3218
51 alpha[51] 1.001     4474     3221
52 alpha[52] 1.000     4420     3210
53 alpha[53] 1.000     4537     3475
54 alpha[54] 1.000     4351     3580
55 alpha[55] 1.000     4406     3318
56 alpha[56] 1.000     4365     3515
57 alpha[57] 1.001     4413     3339
58 alpha[58] 1.001     4175     3227
59 alpha[59] 1.001     4311     3079
60 alpha[60] 1.000     4379     3249
61 alpha[61] 1.000     4819     3510
62   beta[1] 1.001     9230     3214
63   beta[2] 1.001    11290     3033
64   beta[3] 1.001     5661     3306
65   beta[4] 1.000     8678     2981
66   beta[5] 1.002     2370     3172
67   beta[6] 0.999     5487     3428
68   beta[7] 1.000     6256     3140
69   beta[8] 1.001     8352     2833
70   beta[9] 1.003     9615     2846
71  beta[10] 1.001     8974     2815
72  beta[11] 0.999     9748     2808
73  beta[12] 1.001     8191     3049
74  beta[13] 1.001     8477     2894
75  beta[14] 1.001     8064     2863
76  beta[15] 1.001     8914     2872
77  beta[16] 1.003     9863     2561
78  beta[17] 1.003     9574     3013
79  beta[18] 1.002    10513     2787
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
Obs517 LOO log L-1785.74±22.31 g3.254 [3.008, 3.553] C0.828 [0.825, 0.83]
Draws4000 LOO IC3571.49±44.61 gp0.415 [0.402, 0.427] Dxy0.656 [0.651, 0.661]
Chains4 Effective p89.57±4.93 EV0.532 [0.497, 0.569]
Time3.2s B0.117 [0.113, 0.12] v8.335 [6.987, 9.747]
p18 vp0.133 [0.124, 0.143]
Mode β Mean β Median β S.E. Lower Upper P(β>0) Symmetry
treat=5000U   0.2207   0.2075   0.2173  0.5440  -0.8397   1.2575  0.6495  0.96
treat=Placebo   1.8311   1.8351   1.8334  0.5469   0.7705   2.8788  0.9995  0.99
week   0.4865   0.4870   0.4862  0.0809   0.3270   0.6467  1.0000  0.99
week'  -0.2878  -0.2878  -0.2874  0.0868  -0.4578  -0.1176  0.0003  1.05
ptwstrs   0.1998   0.2008   0.2006  0.0262   0.1505   0.2541  1.0000  1.01
ptwstrs'  -0.0621  -0.0648  -0.0651  0.0613  -0.1888   0.0527  0.1460  0.97
ptwstrs''   0.5328   0.5454   0.5469  0.2462   0.0670   1.0405  0.9880  0.99
age  -0.0295  -0.0285  -0.0288  0.0312  -0.0852   0.0365  0.1828  1.04
age'   0.1237   0.1206   0.1213  0.0887  -0.0522   0.2863  0.9072  0.96
age''  -0.5070  -0.4949  -0.4989  0.3538  -1.1742   0.1899  0.0860  1.06
sex=M  -0.4652  -0.4476  -0.4581  2.3307  -5.0263   4.1168  0.4235  1.01
treat=5000U × week  -0.0341  -0.0318  -0.0321  0.1067  -0.2447   0.1700  0.3723  1.00
treat=Placebo × week  -0.2719  -0.2733  -0.2720  0.1064  -0.4765  -0.0688  0.0058  1.00
treat=5000U × week'  -0.0340  -0.0366  -0.0342  0.1158  -0.2587   0.1931  0.3787  0.99
treat=Placebo × week'   0.1198   0.1208   0.1200  0.1159  -0.1090   0.3471  0.8570  1.01
age × sex=M   0.0113   0.0109   0.0107  0.0560  -0.0938   0.1247  0.5785  1.01
age' × sex=M  -0.0511  -0.0500  -0.0522  0.1611  -0.3540   0.2724  0.3775  0.99
age'' × sex=M   0.2619   0.2586   0.2661  0.6341  -1.0121   1.4392  0.6522  0.97
Code 
a <- anova(bmark)
a

Relative Explained Variation for twstrs. Approximate total model Wald χ2 used in denominators of REV:466.7 [391.5, 567.3].

REV Lower Upper d.f.
treat (Factor+Higher Order Factors) 0.052 0.026 0.102 6
 All Interactions 0.049 0.020 0.093 4
week (Factor+Higher Order Factors) 0.294 0.225 0.376 6
 All Interactions 0.049 0.020 0.093 4
 Nonlinear (Factor+Higher Order Factors) 0.059 0.024 0.104 3
ptwstrs 0.935 0.862 0.950 3
 Nonlinear 0.039 0.007 0.072 2
age (Factor+Higher Order Factors) 0.008 0.002 0.039 6
 All Interactions 0.001 0.000 0.018 3
 Nonlinear (Factor+Higher Order Factors) 0.005 0.001 0.028 4
sex (Factor+Higher Order Factors) 0.001 0.000 0.021 4
 All Interactions 0.001 0.000 0.018 3
treat × week (Factor+Higher Order Factors) 0.049 0.020 0.093 4
 Nonlinear 0.004 0.000 0.022 2
 Nonlinear Interaction : f(A,B) vs. AB 0.004 0.000 0.022 2
age × sex (Factor+Higher Order Factors) 0.001 0.000 0.018 3
 Nonlinear 0.001 0.000 0.014 2
 Nonlinear Interaction : f(A,B) vs. AB 0.001 0.000 0.014 2
TOTAL NONLINEAR 0.101 0.061 0.170 9
TOTAL INTERACTION 0.050 0.025 0.102 7
TOTAL NONLINEAR + INTERACTION 0.139 0.095 0.221 12
TOTAL 1.000 1.000 1.000 18
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')
Running MCMC with 4 chains, at most 9 in parallel...

Chain 1 finished in 2.7 seconds.
Chain 2 finished in 2.9 seconds.
Chain 3 finished in 3.1 seconds.
Chain 4 finished in 3.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 3.0 seconds.
Total execution time: 3.3 seconds.
Code 
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.
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.971 0.975 1.047 0.941 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.001     1815     1725
2   alpha[2] 1.002     1526     2161
3   alpha[3] 1.004     1391     2142
4   alpha[4] 1.004     1176     1927
5   alpha[5] 1.004     1048     1564
6   alpha[6] 1.003     1032     1609
7   alpha[7] 1.003      998     1442
8   alpha[8] 1.004     1006     1543
9   alpha[9] 1.003     1227     1760
10 alpha[10] 1.003     1186     1767
11 alpha[11] 1.003     1224     1840
12 alpha[12] 1.003     1213     1557
13 alpha[13] 1.003     1191     1633
14 alpha[14] 1.002     1331     1839
15 alpha[15] 1.002     1406     1937
16 alpha[16] 1.002     1469     1992
17 alpha[17] 1.002     1487     2127
18 alpha[18] 1.001     1694     2535
19 alpha[19] 1.001     1853     2406
20 alpha[20] 1.001     1836     2283
21 alpha[21] 1.000     1964     2548
22 alpha[22] 1.000     2215     2813
23 alpha[23] 1.000     2351     2733
24 alpha[24] 1.000     2410     2755
25 alpha[25] 1.000     2581     2730
26 alpha[26] 1.000     2806     2739
27 alpha[27] 1.000     2892     3067
28 alpha[28] 1.000     3057     2733
29 alpha[29] 1.001     3242     2586
30 alpha[30] 1.001     3588     2726
31 alpha[31] 1.002     3922     2843
32 alpha[32] 1.003     4090     2476
33 alpha[33] 1.003     4141     2717
34 alpha[34] 1.003     4754     3001
35 alpha[35] 1.000     4267     2722
36 alpha[36] 1.000     4180     2815
37 alpha[37] 1.001     3721     2204
38 alpha[38] 1.003     3708     2687
39 alpha[39] 1.002     3578     2645
40 alpha[40] 1.001     3314     2331
41 alpha[41] 1.001     2759     2664
42 alpha[42] 1.001     2960     2571
43 alpha[43] 1.001     2911     2930
44 alpha[44] 1.002     2663     2714
45 alpha[45] 1.001     2502     2740
46 alpha[46] 1.001     2385     2700
47 alpha[47] 1.001     2346     2017
48 alpha[48] 1.001     2446     2210
49 alpha[49] 1.001     2498     2451
50 alpha[50] 1.002     2676     2986
51 alpha[51] 1.001     2674     2813
52 alpha[52] 1.002     2518     2694
53 alpha[53] 1.001     2408     2760
54 alpha[54] 1.001     2530     2917
55 alpha[55] 1.001     2419     2781
56 alpha[56] 1.001     2404     2228
57 alpha[57] 1.001     2621     2656
58 alpha[58] 1.001     2368     2309
59 alpha[59] 1.001     2408     2441
60 alpha[60] 1.001     2476     2593
61 alpha[61] 1.001     2859     2348
62   beta[1] 1.001     5254     2785
63   beta[2] 1.003     5733     2856
64   beta[3] 1.004     3442     2735
65   beta[4] 1.002     7324     2424
66   beta[5] 1.001     1541     2641
67   beta[6] 1.000     2961     2904
68   beta[7] 1.000     4645     3273
69   beta[8] 1.001     5736     2898
70   beta[9] 1.001     6212     2992
71  beta[10] 1.000     5756     2917
72  beta[11] 1.004     6519     2615
73  beta[12] 1.001     5910     2300
74  beta[13] 1.004     5476     2811
75  beta[14] 1.001     6638     3038
76  beta[15] 1.001     6875     2819
77  beta[16] 1.001     5592     2375
78  beta[17] 1.003     6170     2650
79  beta[18] 1.002     6490     2759
80 sigmag[1] 1.006     1208     1624
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
Obs517 LOO log L-1786.83±22.3 g3.258 [3.03, 3.604] C0.828 [0.825, 0.83]
Draws4000 LOO IC3573.67±44.59 gp0.415 [0.403, 0.428] Dxy0.656 [0.65, 0.661]
Chains4 Effective p93.55±5.09 EV0.532 [0.498, 0.57]
Time4s B0.117 [0.114, 0.122] v8.354 [7.107, 10.055]
p18 vp0.133 [0.124, 0.142]
Cluster onuid
Clusters108
σγ0.1088 [0, 0.3295]
Mean β Median β S.E. Lower Upper P(β>0) Symmetry
treat=5000U   0.2131   0.2252  0.5601  -0.9279   1.2522  0.6480  0.98
treat=Placebo   1.8299   1.8429  0.5691   0.7254   2.9530  0.9995  1.05
week   0.4865   0.4860  0.0809   0.3200   0.6364  1.0000  1.09
week'  -0.2858  -0.2858  0.0875  -0.4627  -0.1193  0.0003  0.95
ptwstrs   0.1990   0.1994  0.0265   0.1470   0.2506  1.0000  1.05
ptwstrs'  -0.0626  -0.0627  0.0625  -0.1829   0.0585  0.1622  0.97
ptwstrs''   0.5387   0.5376  0.2476   0.0818   1.0415  0.9875  1.05
age  -0.0287  -0.0289  0.0330  -0.0944   0.0357  0.1872  1.05
age'   0.1227   0.1225  0.0913  -0.0540   0.3034  0.9148  0.96
age''  -0.5065  -0.5068  0.3601  -1.2313   0.1815  0.0740  1.02
sex=M  -0.4484  -0.4560  2.3833  -5.2222   4.0229  0.4252  0.99
treat=5000U × week  -0.0318  -0.0313  0.1089  -0.2528   0.1771  0.3865  1.00
treat=Placebo × week  -0.2714  -0.2720  0.1122  -0.4926  -0.0499  0.0058  1.00
treat=5000U × week'  -0.0373  -0.0368  0.1183  -0.2763   0.1830  0.3845  1.04
treat=Placebo × week'   0.1178   0.1181  0.1223  -0.1243   0.3573  0.8340  1.02
age × sex=M   0.0108   0.0111  0.0574  -0.0998   0.1232  0.5755  1.02
age' × sex=M  -0.0499  -0.0514  0.1638  -0.3680   0.2732  0.3845  0.99
age'' × sex=M   0.2583   0.2576  0.6402  -1.0470   1.4578  0.6518  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:

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 P(Contrast>0)
1  Week 2     2 -1.2884028 0.3718895 -2.0055357 -0.5779557        0.0005
2  Week 4     4 -0.7417251 0.2607277 -1.2767224 -0.2557648        0.0032
3  Week 8     8  0.2307980 0.3456263 -0.4672853  0.8873975        0.7550
4* Week 12   12  0.7199925 0.2538495  0.1871243  1.1926845        0.9980
5* Week 16   16  1.0883549 0.3890277  0.3237614  1.8310467        0.9982

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
Z
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
A
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
B
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
C
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
D
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.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. The improvement comes from canceling out a buildup of variances of \(Y\) due to using lagged raw \(Y\) instead of lag residuals as used in linear models. For this purpose, a good choice of weighting functions \(w(t)\) is an indicator function for being at the first follow-up time.

Code 
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, x=TRUE)
f

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = twstrs ~ treat * rcs(week, 3) + rcs(ptwstrs, 4) + 
    rcs(age, 4) * sex + cluster(uid) + mix_re(mix2), data = both, 
    x = TRUE)
Frequencies of Missing Values Due to Each Variable
      twstrs        treat         week      ptwstrs          age          sex 
           0            0            0            5            0            0 
cluster(uid) mix_re(mix2) 
           0            0
Model Likelihood
Ratio Test		 Discrimination
Indexes		 Rank Discrim.
Indexes
Obs517 LR χ2291.74 R20.432 ρ0.834
ESS516.7 d.f.18 R218,5170.411 Dxy0.661
Distinct Y62 P(>χ2)<0.0001 R218,516.70.411
Cluster onuid |P(Y ≥ median)-½|0.341
Clusters108
σ11.4341 [0.9403, 2.1872]
σ2-0.20520 [-0.44540, 0.03504]
Y0.542
max |∂log L/∂β|2×10-5
β S.E. Wald Z P(>|Z|)
treat=5000U   0.3063  0.6466 0.47 0.6358
treat=Placebo   1.6058  0.6526 2.46 0.0139
week   0.4910  0.0901 5.45 <0.0001
week'  -0.2870  0.0924 -3.11 0.0019
ptwstrs   0.2226  0.0288 7.74 <0.0001
ptwstrs'  -0.0728  0.0658 -1.11 0.2680
ptwstrs''   0.5705  0.2631 2.17 0.0302
age  -0.0319  0.0310 -1.03 0.3025
age'   0.1309  0.0864 1.51 0.1299
age''  -0.5317  0.3417 -1.56 0.1197
sex=M  -0.5530  2.3744 -0.23 0.8158
treat=5000U × week  -0.0485  0.1226 -0.40 0.6926
treat=Placebo × week  -0.2315  0.1250 -1.85 0.0641
treat=5000U × week'  -0.0248  0.1268 -0.20 0.8452
treat=Placebo × week'   0.0818  0.1294 0.63 0.5275
age × sex=M   0.0135  0.0572 0.24 0.8134
age' × sex=M  -0.0542  0.1621 -0.33 0.7382
age'' × sex=M   0.2631  0.6288 0.42 0.6756
Code 
deviance(f)
                 intercepts                intercepts+x 
                   4041.567                    3392.026 
  intercepts+random effects intercepts+x+random effects 
                   3674.128                    3382.387

\(\sigma_{1}\) is large; there is substantial subject-level correlation at baseline (before ptwstrs is contributing anything), consistent with a real, non-trivial random effect there.

\(\sigma_{2} ≈ -0.2\) is 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 \(\sigma_2\) is the model’s way of trimming that slight over-induction rather than needing to zero the whole random effect out just to compensate.

7.9 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?
  2. 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?
  2. In a longitudinal intervention study, what is the most typical comparison of interest? Is it best to borrow information in estimating this contrast?