7  Modeling Longitudinal Responses using Generalized Least Squares

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

7.1 Notation

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

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

7.2.1 Common Basis Functions

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

7.2.2 Model for Mean Profile

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

7.2.3 Model Specification for Treatment Comparisons

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

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
    • 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):
    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)
  • 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]\)

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

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

    Figure 7.4: 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.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
    χ2 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.

    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.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
    W
    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.844 0.783 0.745 0.745 
    
       Parameter  Rhat ESS bulk ESS tail
    1   alpha[1] 1.003     2046     2023
    2   alpha[2] 1.003     1529     2315
    3   alpha[3] 1.005     1052     1983
    4   alpha[4] 1.006      952     1882
    5   alpha[5] 1.007      938     1893
    6   alpha[6] 1.008      805     1789
    7   alpha[7] 1.008      780     1708
    8   alpha[8] 1.009      733     1563
    9   alpha[9] 1.009      665     1293
    10 alpha[10] 1.009      636     1296
    11 alpha[11] 1.011      612     1413
    12 alpha[12] 1.010      591     1211
    13 alpha[13] 1.012      566     1111
    14 alpha[14] 1.012      537     1048
    15 alpha[15] 1.013      531     1049
    16 alpha[16] 1.011      515     1099
    17 alpha[17] 1.011      503      976
    18 alpha[18] 1.010      489     1023
    19 alpha[19] 1.010      478      975
    20 alpha[20] 1.010      474     1067
    21 alpha[21] 1.010      471      942
    22 alpha[22] 1.010      464      906
    23 alpha[23] 1.010      460      834
    24 alpha[24] 1.011      457      859
    25 alpha[25] 1.011      460      861
    26 alpha[26] 1.010      459      912
    27 alpha[27] 1.010      470      919
    28 alpha[28] 1.010      475      949
    29 alpha[29] 1.010      476      953
    30 alpha[30] 1.010      465      794
    31 alpha[31] 1.011      466      845
    32 alpha[32] 1.011      473     1012
    33 alpha[33] 1.010      475      826
    34 alpha[34] 1.009      498      956
    35 alpha[35] 1.008      500      965
    36 alpha[36] 1.007      512     1081
    37 alpha[37] 1.007      538     1281
    38 alpha[38] 1.007      548     1323
    39 alpha[39] 1.005      590     1176
    40 alpha[40] 1.005      608     1300
    41 alpha[41] 1.006      625     1251
    42 alpha[42] 1.004      693     1447
    43 alpha[43] 1.003      763     1548
    44 alpha[44] 1.003      783     1716
    45 alpha[45] 1.002      841     1812
    46 alpha[46] 1.002      931     1733
    47 alpha[47] 1.001      990     2263
    48 alpha[48] 1.001     1013     2117
    49 alpha[49] 1.001     1070     2086
    50 alpha[50] 1.001     1121     2087
    51 alpha[51] 1.001     1217     1672
    52 alpha[52] 1.001     1268     1753
    53 alpha[53] 1.001     1413     2342
    54 alpha[54] 1.000     1437     2317
    55 alpha[55] 1.000     1543     2719
    56 alpha[56] 1.000     1583     2874
    57 alpha[57] 1.001     1588     2611
    58 alpha[58] 1.001     1722     2499
    59 alpha[59] 1.001     1805     2782
    60 alpha[60] 1.001     2006     3023
    61 alpha[61] 1.000     2275     3210
    62   beta[1] 1.004      757     1387
    63   beta[2] 1.006      762     1405
    64   beta[3] 1.001     2022     2740
    65   beta[4] 1.001     3963     3075
    66   beta[5] 1.003      826     1558
    67   beta[6] 1.003      773     1543
    68   beta[7] 1.006      916     1326
    69   beta[8] 1.003      849     1238
    70   beta[9] 1.001      879     1191
    71  beta[10] 1.004      765     1320
    72  beta[11] 1.001     3836     2899
    73  beta[12] 1.001     2896     2707
    74  beta[13] 1.000     4148     3222
    75  beta[14] 1.002     3909     2921
    76  beta[15] 1.008      821     1351
    77  beta[16] 1.003      780     1463
    78  beta[17] 1.005      889     1447
    79 sigmag[1] 1.003      773     1636
    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 -1747.2±23.67 g 3.795 [3.292, 4.281] C 0.793 [0.785, 0.799]
    Draws 4000 LOO IC 3494.4±47.33 gp 0.434 [0.418, 0.447] Dxy 0.585 [0.569, 0.599]
    Chains 4 Effective p 178.76±7.92 EV 0.589 [0.539, 0.632]
    Time 9s B 0.149 [0.139, 0.16] v 11.211 [8.387, 14.176]
    p 17 vp 0.147 [0.135, 0.158]
    Cluster on uid
    Clusters 108
    σγ 1.871 [1.5275, 2.2372]
    Mean β Median β S.E. Lower Upper Pr(β>0) Symmetry
    y≥7   -1.8925   -1.9348  4.3060  -10.4593   6.4242  0.3235  1.05
    y≥9   -2.9245   -2.8859  4.1909  -11.6146   4.6908  0.2448  1.02
    y≥10   -4.0987   -4.1024  4.1327  -12.1431   3.8700  0.1608  1.04
    y≥11   -4.5443   -4.5475  4.1279  -12.2445   3.7380  0.1320  1.03
    y≥13   -4.7415   -4.7522  4.1280  -12.7626   3.1650  0.1235  1.02
    y≥14   -5.0926   -5.1122  4.1261  -13.3171   2.5673  0.1045  1.03
    y≥15   -5.3994   -5.4131  4.1229  -13.2727   2.6133  0.0912  1.04
    y≥16   -5.7832   -5.8074  4.1193  -13.6829   2.2175  0.0808  1.04
    y≥17   -6.5839   -6.5902  4.1031  -14.3847   1.2904  0.0538  1.01
    y≥18   -6.8380   -6.8208  4.0984  -14.6403   1.0348  0.0490  1.01
    y≥19   -7.1251   -7.1242  4.0945  -14.8641   0.7626  0.0418  1.02
    y≥20   -7.3212   -7.3342  4.0950  -15.0929   0.5257  0.0372  1.01
    y≥21   -7.5016   -7.5194  4.0947  -15.2810   0.2955  0.0312  1.00
    y≥22   -7.9177   -7.9096  4.0906  -15.7228   -0.0901  0.0243  1.01
    y≥23   -8.1852   -8.1709  4.0900  -16.0970   -0.3846  0.0215  1.00
    y≥24   -8.4661   -8.4382  4.0902  -16.2124   -0.5187  0.0192  1.00
    y≥25   -8.7240   -8.6958  4.0889  -16.3427   -0.6265  0.0170  1.00
    y≥26   -9.1195   -9.1033  4.0900  -16.8669   -1.1156  0.0140  1.00
    y≥27   -9.4088   -9.3790  4.0889  -17.1584   -1.3981  0.0115  1.01
    y≥28   -9.6505   -9.6273  4.0899  -17.2485   -1.4914  0.0102  1.00
    y≥29   -9.8850   -9.8647  4.0916  -17.5335   -1.8057  0.0088  1.00
    y≥30  -10.1875  -10.1761  4.0931  -17.8358   -2.1341  0.0075  1.00
    y≥31  -10.4768  -10.4712  4.0948  -18.1603   -2.4314  0.0073  1.00
    y≥32  -10.5931  -10.5991  4.0957  -18.2305   -2.5091  0.0068  1.01
    y≥33  -10.9555  -10.9465  4.0997  -18.6866   -2.8897  0.0058  1.01
    y≥34  -11.2670  -11.2598  4.1016  -18.9527   -3.1342  0.0032  1.01
    y≥35  -11.4901  -11.4901  4.1008  -19.2831   -3.4359  0.0027  1.00
    y≥36  -11.7354  -11.7072  4.1001  -19.4760   -3.6375  0.0018  1.00
    y≥37  -12.0134  -12.0031  4.1028  -19.7269   -3.9472  0.0015  0.99
    y≥38  -12.2403  -12.2209  4.1053  -20.1085   -4.3351  0.0015  0.99
    y≥39  -12.4861  -12.4689  4.1066  -20.3765   -4.5740  0.0015  1.00
    y≥40  -12.6687  -12.6621  4.1076  -20.5231   -4.7360  0.0013  1.00
    y≥41  -12.8521  -12.8450  4.1093  -20.6431   -4.8363  0.0013  1.00
    y≥42  -13.1765  -13.1852  4.1114  -20.9928   -5.1592  0.0013  1.00
    y≥43  -13.4033  -13.4016  4.1115  -21.1990   -5.3680  0.0013  1.00
    y≥44  -13.7428  -13.7425  4.1131  -21.5027   -5.6825  0.0005  0.99
    y≥45  -14.0603  -14.0603  4.1127  -21.7960   -5.8977  0.0005  1.00
    y≥46  -14.3587  -14.3558  4.1139  -22.1307   -6.2589  0.0005  1.00
    y≥47  -14.7783  -14.7739  4.1157  -22.6281   -6.7746  0.0005  0.99
    y≥48  -15.0661  -15.0521  4.1200  -22.8265   -6.9166  0.0005  0.98
    y≥49  -15.4378  -15.4149  4.1237  -23.2339   -7.2999  0.0003  0.98
    y≥50  -15.7575  -15.7567  4.1266  -23.6138   -7.6148  0.0003  0.99
    y≥51  -16.2863  -16.2865  4.1275  -24.1827   -8.2012  0.0003  1.00
    y≥52  -16.6484  -16.6180  4.1256  -24.5473   -8.5991  0.0003  0.98
    y≥53  -17.0940  -17.0594  4.1287  -25.1299   -9.1439  0.0000  0.97
    y≥54  -17.5950  -17.5497  4.1311  -25.6404   -9.6396  0.0000  0.97
    y≥55  -18.0042  -17.9506  4.1338  -26.0081  -10.0171  0.0000  0.98
    y≥56  -18.2570  -18.2166  4.1333  -26.2906  -10.2969  0.0000  0.98
    y≥57  -18.7246  -18.6750  4.1337  -26.8089  -10.8036  0.0000  0.99
    y≥58  -19.2848  -19.2290  4.1365  -27.3114  -11.2621  0.0000  1.00
    y≥59  -19.6375  -19.5835  4.1387  -27.7966  -11.7426  0.0000  0.99
    y≥60  -19.9638  -19.9158  4.1408  -27.9487  -11.9096  0.0000  1.00
    y≥61  -20.6541  -20.6027  4.1412  -28.6340  -12.7014  0.0000  0.98
    y≥62  -21.0249  -20.9845  4.1434  -29.0384  -13.0626  0.0000  0.98
    y≥63  -21.4276  -21.3862  4.1536  -29.2607  -13.1148  0.0000  0.98
    y≥64  -21.5672  -21.5031  4.1545  -29.2026  -13.1253  0.0000  0.97
    y≥65  -22.2878  -22.2548  4.1519  -29.8979  -13.6954  0.0000  0.97
    y≥66  -22.6661  -22.6505  4.1554  -30.7065  -14.3955  0.0000  0.98
    y≥67  -23.0760  -23.0377  4.1576  -31.5829  -15.2965  0.0000  0.97
    y≥68  -23.8575  -23.8083  4.1745  -32.1561  -15.9086  0.0000  0.98
    y≥71  -24.7235  -24.6678  4.1948  -32.9053  -16.6747  0.0000  0.98
    treat=5000U   0.1185   0.1304  0.7254   -1.2585   1.6167  0.5702  1.05
    treat=Placebo   2.3524   2.3516  0.7319   0.9497   3.8002  0.9998  1.00
    week   0.1218   0.1209  0.0789   -0.0335   0.2749  0.9380  1.01
    week'   0.1907   0.1902  0.0856   0.0219   0.3588  0.9862  1.00
    twstrs0   0.2273   0.2275  0.0516   0.1245   0.3250  1.0000  1.01
    twstrs0'   0.1285   0.1289  0.0632   0.0046   0.2492  0.9792  0.99
    age   -0.0130   -0.0120  0.0795   -0.1729   0.1384  0.4335  0.93
    age'   0.1921   0.1842  0.2197   -0.2576   0.6072  0.8128  1.08
    age''   -1.0722   -1.0545  0.8609   -2.6773   0.6575  0.1020  0.96
    sex=M   5.2646   5.2537  6.0906   -7.0597   16.7183  0.8082  1.01
    treat=5000U × week   0.0504   0.0500  0.1102   -0.1673   0.2585  0.6775  0.98
    treat=Placebo × week   -0.0522   -0.0528  0.1107   -0.2707   0.1658  0.3200  0.98
    treat=5000U × week'   -0.1609   -0.1603  0.1205   -0.3927   0.0767  0.0922  1.00
    treat=Placebo × week'   -0.1410   -0.1400  0.1210   -0.3905   0.0852  0.1195  1.00
    age × sex=M   -0.1157   -0.1170  0.1457   -0.4006   0.1681  0.2120  0.99
    age' × sex=M   0.1596   0.1619  0.4079   -0.6326   0.9422  0.6587  0.99
    age'' × sex=M   0.0197   -0.0055  1.5849   -3.1326   2.9931  0.4993  1.05
    Code
    a <- anova(bpo)
    a
    Relative Explained Variation for twstrs. Approximate total model Wald χ2 used in denominators of REV:262.1 [212.5, 337].
    REV Lower Upper d.f.
    treat (Factor+Higher Order Factors) 0.121 0.069 0.223 6
    All Interactions 0.083 0.030 0.166 4
    week (Factor+Higher Order Factors) 0.549 0.422 0.652 6
    All Interactions 0.083 0.030 0.166 4
    Nonlinear (Factor+Higher Order Factors) 0.021 0.000 0.070 3
    twstrs0 0.659 0.524 0.736 2
    Nonlinear 0.016 0.000 0.051 1
    age (Factor+Higher Order Factors) 0.024 0.012 0.085 6
    All Interactions 0.014 0.001 0.056 3
    Nonlinear (Factor+Higher Order Factors) 0.021 0.003 0.067 4
    sex (Factor+Higher Order Factors) 0.018 0.003 0.066 4
    All Interactions 0.014 0.001 0.056 3
    treat × week (Factor+Higher Order Factors) 0.083 0.030 0.166 4
    Nonlinear 0.008 0.000 0.045 2
    Nonlinear Interaction : f(A,B) vs. AB 0.008 0.000 0.045 2
    age × sex (Factor+Higher Order Factors) 0.014 0.001 0.056 3
    Nonlinear 0.014 0.000 0.051 2
    Nonlinear Interaction : f(A,B) vs. AB 0.014 0.000 0.051 2
    TOTAL NONLINEAR 0.053 0.022 0.127 8
    TOTAL INTERACTION 0.096 0.059 0.201 7
    TOTAL NONLINEAR + INTERACTION 0.128 0.093 0.238 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 Pr(Contrast>0)
    1  Week 2     2 -2.2478977 0.5998165 -3.439365 -1.0775322         0.0000
    2  Week 4     4 -2.1434302 0.5309202 -3.164919 -1.0687934         0.0000
    3  Week 8     8 -1.7934509 0.5785167 -2.821921 -0.6136238         0.0018
    4* Week 12   12 -0.8792943 0.5252163 -1.937805  0.1055448         0.0423
    5* Week 16   16  0.1759067 0.6120968 -1.030554  1.3625489         0.6212
    
    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:
    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')
    # 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: 0.904 0.936 1.024 1.066 
    
       Parameter  Rhat ESS bulk ESS tail
    1   alpha[1] 1.002     4210     2266
    2   alpha[2] 1.000     5448     2648
    3   alpha[3] 1.001     5708     3246
    4   alpha[4] 1.002     5310     3542
    5   alpha[5] 1.002     5412     3452
    6   alpha[6] 1.000     4983     3646
    7   alpha[7] 1.000     4537     3728
    8   alpha[8] 1.000     4396     3614
    9   alpha[9] 1.000     4659     3673
    10 alpha[10] 1.001     4294     3564
    11 alpha[11] 1.000     4062     3361
    12 alpha[12] 1.001     3907     3195
    13 alpha[13] 1.000     3974     3305
    14 alpha[14] 1.000     4083     3446
    15 alpha[15] 1.000     4109     3146
    16 alpha[16] 1.000     4009     3361
    17 alpha[17] 1.000     4082     3272
    18 alpha[18] 1.000     3929     3555
    19 alpha[19] 1.000     3912     3281
    20 alpha[20] 1.001     3940     3649
    21 alpha[21] 1.000     3899     3651
    22 alpha[22] 1.000     3811     3457
    23 alpha[23] 1.000     4187     3550
    24 alpha[24] 1.000     4191     3536
    25 alpha[25] 1.000     4334     2950
    26 alpha[26] 1.000     4628     2935
    27 alpha[27] 1.000     4786     2630
    28 alpha[28] 1.000     5216     2370
    29 alpha[29] 1.001     5496     2361
    30 alpha[30] 1.001     5947     2366
    31 alpha[31] 1.000     6171     2883
    32 alpha[32] 0.999     6520     2845
    33 alpha[33] 0.999     7052     2853
    34 alpha[34] 1.000     7180     2653
    35 alpha[35] 1.002     7991     2756
    36 alpha[36] 1.001     7469     2965
    37 alpha[37] 1.000     7865     2831
    38 alpha[38] 1.000     7703     3224
    39 alpha[39] 1.000     7350     3338
    40 alpha[40] 1.000     6548     3208
    41 alpha[41] 1.000     5862     3407
    42 alpha[42] 1.000     5915     3354
    43 alpha[43] 1.000     5789     2915
    44 alpha[44] 1.000     5532     3846
    45 alpha[45] 1.001     5671     3483
    46 alpha[46] 1.000     5546     3451
    47 alpha[47] 1.001     5877     3531
    48 alpha[48] 1.000     6120     3993
    49 alpha[49] 1.000     5717     3726
    50 alpha[50] 1.001     5591     3553
    51 alpha[51] 1.000     5507     3760
    52 alpha[52] 1.000     5689     3807
    53 alpha[53] 1.001     5765     3798
    54 alpha[54] 1.001     5980     3649
    55 alpha[55] 1.000     5454     3704
    56 alpha[56] 1.000     5654     3820
    57 alpha[57] 1.001     5843     3575
    58 alpha[58] 1.001     6060     3593
    59 alpha[59] 1.000     6152     3853
    60 alpha[60] 1.001     5753     3902
    61 alpha[61] 1.001     5715     3199
    62   beta[1] 1.001     8954     3007
    63   beta[2] 1.000    10329     2842
    64   beta[3] 1.000     6872     3509
    65   beta[4] 1.000     8459     2891
    66   beta[5] 1.000     3707     2781
    67   beta[6] 1.001     6819     3336
    68   beta[7] 1.001     8771     3188
    69   beta[8] 1.001     9635     2451
    70   beta[9] 1.001     9809     3234
    71  beta[10] 1.006     9893     2584
    72  beta[11] 1.002     8059     2929
    73  beta[12] 1.000     8871     2859
    74  beta[13] 1.001     7780     2871
    75  beta[14] 1.002    11064     2953
    76  beta[15] 1.000     8958     3145
    77  beta[16] 1.002     9525     3036
    78  beta[17] 1.003     9226     3005
    79  beta[18] 1.001     9610     2941
    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 -1787.42±22.54 g 3.258 [3.029, 3.575] C 0.828 [0.825, 0.83]
    Draws 4000 LOO IC 3574.83±45.08 gp 0.415 [0.403, 0.427] Dxy 0.655 [0.65, 0.661]
    Chains 4 Effective p 91.34±5.41 EV 0.531 [0.495, 0.566]
    Time 4.7s B 0.117 [0.113, 0.121] v 8.358 [7.076, 9.88]
    p 18 vp 0.133 [0.124, 0.141]
    Mode β Mean β Median β S.E. Lower Upper Pr(β>0) Symmetry
    treat=5000U   0.2210   0.2199   0.2249  0.5817  -0.9520   1.3083  0.6400  1.00
    treat=Placebo   1.8297   1.8301   1.8332  0.5733   0.6659   2.8994  0.9998  1.01
    week   0.4860   0.4856   0.4855  0.0849   0.3185   0.6467  1.0000  0.99
    week'  -0.2876  -0.2859  -0.2858  0.0899  -0.4581  -0.1126  0.0008  0.98
    ptwstrs   0.1996   0.2008   0.2007  0.0266   0.1473   0.2511  1.0000  0.97
    ptwstrs'  -0.0621  -0.0650  -0.0643  0.0625  -0.1866   0.0594  0.1475  0.99
    ptwstrs''   0.5323   0.5475   0.5453  0.2479   0.0566   1.0377  0.9848  1.02
    age  -0.0295  -0.0289  -0.0290  0.0316  -0.0912   0.0346  0.1735  1.04
    age'   0.1235   0.1213   0.1214  0.0881  -0.0626   0.2855  0.9115  1.01
    age''  -0.5064  -0.4971  -0.4999  0.3474  -1.1322   0.2281  0.0813  1.00
    sex=M  -0.4633  -0.4417  -0.4722  2.3794  -5.3183   3.7576  0.4280  1.02
    treat=5000U × week  -0.0341  -0.0338  -0.0332  0.1139  -0.2564   0.1807  0.3868  0.97
    treat=Placebo × week  -0.2717  -0.2714  -0.2712  0.1150  -0.4839  -0.0355  0.0098  1.01
    treat=5000U × week'  -0.0339  -0.0345  -0.0357  0.1234  -0.2701   0.2009  0.3935  1.01
    treat=Placebo × week'   0.1196   0.1186   0.1194  0.1244  -0.1358   0.3524  0.8288  1.01
    age × sex=M   0.0112   0.0107   0.0108  0.0573  -0.0894   0.1291  0.5755  1.00
    age' × sex=M  -0.0509  -0.0485  -0.0478  0.1625  -0.3634   0.2571  0.3850  0.99
    age'' × sex=M   0.2614   0.2501   0.2486  0.6318  -0.9004   1.5136  0.6492  1.02
    Code
    a <- anova(bpo)
    a
    Relative Explained Variation for twstrs. Approximate total model Wald χ2 used in denominators of REV:262.1 [215.2, 340.5].
    REV Lower Upper d.f.
    treat (Factor+Higher Order Factors) 0.121 0.064 0.214 6
    All Interactions 0.083 0.036 0.160 4
    week (Factor+Higher Order Factors) 0.549 0.424 0.656 6
    All Interactions 0.083 0.036 0.160 4
    Nonlinear (Factor+Higher Order Factors) 0.021 0.003 0.064 3
    twstrs0 0.659 0.529 0.748 2
    Nonlinear 0.016 0.000 0.048 1
    age (Factor+Higher Order Factors) 0.024 0.007 0.088 6
    All Interactions 0.014 0.001 0.057 3
    Nonlinear (Factor+Higher Order Factors) 0.021 0.006 0.076 4
    sex (Factor+Higher Order Factors) 0.018 0.005 0.072 4
    All Interactions 0.014 0.001 0.057 3
    treat × week (Factor+Higher Order Factors) 0.083 0.036 0.160 4
    Nonlinear 0.008 0.000 0.039 2
    Nonlinear Interaction : f(A,B) vs. AB 0.008 0.000 0.039 2
    age × sex (Factor+Higher Order Factors) 0.014 0.001 0.057 3
    Nonlinear 0.014 0.000 0.054 2
    Nonlinear Interaction : f(A,B) vs. AB 0.014 0.000 0.054 2
    TOTAL NONLINEAR 0.053 0.028 0.138 8
    TOTAL INTERACTION 0.096 0.054 0.190 7
    TOTAL NONLINEAR + INTERACTION 0.128 0.088 0.239 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.
    2 of 4000 (0.05%) 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: 1 0 0 1 
    
    EBFMI: 1.042 0.887 0.937 0.989 
    
       Parameter  Rhat ESS bulk ESS tail
    1   alpha[1] 1.002     3945     2235
    2   alpha[2] 1.001     4008     2617
    3   alpha[3] 1.001     3186     3143
    4   alpha[4] 1.000     2813     2535
    5   alpha[5] 1.001     2693     2591
    6   alpha[6] 1.000     2618     2471
    7   alpha[7] 1.000     2511     2560
    8   alpha[8] 1.000     2360     2421
    9   alpha[9] 1.000     2347     2700
    10 alpha[10] 1.000     2294     2807
    11 alpha[11] 1.000     2137     2154
    12 alpha[12] 1.000     2196     2518
    13 alpha[13] 1.000     2048     2404
    14 alpha[14] 1.000     2055     2446
    15 alpha[15] 1.000     2060     2622
    16 alpha[16] 1.000     2066     2785
    17 alpha[17] 1.000     2041     2487
    18 alpha[18] 1.000     1936     2547
    19 alpha[19] 1.000     1929     2202
    20 alpha[20] 1.001     1978     2372
    21 alpha[21] 1.000     1895     2345
    22 alpha[22] 1.000     1802     1946
    23 alpha[23] 1.000     1792     1769
    24 alpha[24] 1.000     1834     2317
    25 alpha[25] 1.000     1999     2355
    26 alpha[26] 1.000     2201     2248
    27 alpha[27] 1.000     2459     2700
    28 alpha[28] 1.000     2509     2849
    29 alpha[29] 1.000     2615     2714
    30 alpha[30] 1.001     2738     2644
    31 alpha[31] 1.000     3214     2496
    32 alpha[32] 1.001     3315     2250
    33 alpha[33] 1.001     3627     2668
    34 alpha[34] 1.002     3842     2425
    35 alpha[35] 1.001     4230     2976
    36 alpha[36] 1.002     4757     2939
    37 alpha[37] 1.003     4895     3213
    38 alpha[38] 1.000     4842     3148
    39 alpha[39] 1.001     4736     3356
    40 alpha[40] 1.001     4539     3264
    41 alpha[41] 1.002     4304     3365
    42 alpha[42] 1.002     4238     3251
    43 alpha[43] 1.001     3945     3324
    44 alpha[44] 1.003     3744     3227
    45 alpha[45] 1.002     3404     3017
    46 alpha[46] 1.001     3436     3195
    47 alpha[47] 1.002     3200     3122
    48 alpha[48] 1.002     3174     3111
    49 alpha[49] 1.003     3181     3028
    50 alpha[50] 1.000     3382     2772
    51 alpha[51] 1.001     3543     2832
    52 alpha[52] 1.002     3512     2833
    53 alpha[53] 1.001     3425     2904
    54 alpha[54] 1.001     3415     2931
    55 alpha[55] 1.001     3376     3152
    56 alpha[56] 1.001     3427     3079
    57 alpha[57] 1.001     3138     3054
    58 alpha[58] 1.001     3315     3150
    59 alpha[59] 1.002     2907     3059
    60 alpha[60] 1.000     3553     3026
    61 alpha[61] 1.001     4031     2978
    62   beta[1] 1.001     5933     2651
    63   beta[2] 1.000     4441     2850
    64   beta[3] 1.002     4251     2765
    65   beta[4] 1.001     5745     2538
    66   beta[5] 1.001     1745     2440
    67   beta[6] 1.000     3954     3004
    68   beta[7] 1.001     4701     2865
    69   beta[8] 1.001     5532     2605
    70   beta[9] 1.004     6115     2597
    71  beta[10] 1.001     5363     2381
    72  beta[11] 1.002     5335     2899
    73  beta[12] 1.000     6075     2504
    74  beta[13] 1.002     5991     2838
    75  beta[14] 1.001     5661     2772
    76  beta[15] 1.004     6543     2720
    77  beta[16] 1.004     5364     2611
    78  beta[17] 1.001     5527     2517
    79  beta[18] 1.001     6023     2615
    80 sigmag[1] 1.003     1264     1702
    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 -1786.26±22.28 g 3.258 [2.942, 3.537] C 0.828 [0.825, 0.83]
    Draws 4000 LOO IC 3572.52±44.55 gp 0.415 [0.401, 0.427] Dxy 0.655 [0.65, 0.66]
    Chains 4 Effective p 93.09±4.62 EV 0.531 [0.492, 0.569]
    Time 6.4s B 0.117 [0.113, 0.121] v 8.359 [6.825, 9.832]
    p 18 vp 0.133 [0.123, 0.142]
    Cluster on uid
    Clusters 108
    σγ 0.1165 [0, 0.3384]
    Mean β Median β S.E. Lower Upper Pr(β>0) Symmetry
    treat=5000U   0.2155   0.2174  0.5864  -0.9299   1.3469  0.6518  0.99
    treat=Placebo   1.8441   1.8421  0.5814   0.7542   3.0127  0.9992  1.02
    week   0.4865   0.4873  0.0801   0.3304   0.6428  1.0000  0.99
    week'  -0.2862  -0.2856  0.0836  -0.4493  -0.1219  0.0005  0.99
    ptwstrs   0.1999   0.2003  0.0269   0.1482   0.2532  1.0000  0.99
    ptwstrs'  -0.0649  -0.0637  0.0622  -0.1856   0.0537  0.1532  1.02
    ptwstrs''   0.5474   0.5462  0.2479   0.0659   1.0186  0.9875  1.04
    age  -0.0294  -0.0298  0.0316  -0.0908   0.0305  0.1767  0.98
    age'   0.1240   0.1249  0.0882  -0.0487   0.2966  0.9202  1.00
    age''  -0.5099  -0.5108  0.3477  -1.1647   0.2046  0.0722  1.01
    sex=M  -0.4235  -0.4078  2.4502  -5.1794   4.2852  0.4328  0.97
    treat=5000U × week  -0.0323  -0.0297  0.1120  -0.2718   0.1710  0.3918  0.99
    treat=Placebo × week  -0.2731  -0.2715  0.1117  -0.5055  -0.0666  0.0070  0.97
    treat=5000U × week'  -0.0360  -0.0372  0.1196  -0.2593   0.2073  0.3812  1.03
    treat=Placebo × week'   0.1197   0.1187  0.1217  -0.1148   0.3585  0.8383  1.00
    age × sex=M   0.0100   0.0104  0.0591  -0.1034   0.1242  0.5635  1.02
    age' × sex=M  -0.0477  -0.0480  0.1682  -0.3958   0.2537  0.3938  0.98
    age'' × sex=M   0.2522   0.2564  0.6546  -0.9502   1.5686  0.6450  0.99

    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.2872979 0.3788043 -2.0676854 -0.5851706         0.0003
    2  Week 4     4 -0.7444594 0.2532338 -1.2155523 -0.2251157         0.0005
    3  Week 8     8  0.2226112 0.3626920 -0.4921401  0.9308765         0.7312
    4* Week 12   12  0.7152567 0.2657876  0.1879102  1.2379469         0.9958
    5* Week 16   16  1.0892960 0.3934977  0.3237274  1.8605355         0.9970
    
    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.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?