Rescue analysis for % change from baseline as an outcome in RCT

[feizuo]

1. If one must use percent change from baseline (% CFB) as an outcome, does the following approach rescue the analysis?

% CFB in 28-day seizure frequency = treatment + baseline seizure frequency

Understanding that the concept of change from baseline does not respect the parallel group design of an RCT, percent change from baseline is unlikely to have good distributional properties for statistical analysis, and a Poisson-type of approach (perhaps negative binomial with an offset term for drug exposure) would be better, but unfortunately the % CFB outcome is what’s been vetted by clinicians and the FDA and some of us are now stuck with using this measure as an endpoint.

2. Would the following approach be preferable over modelling it as count data using a Poisson type of approach? Why or why not?

natural log(28-day seizure rate) = treatment + baseline seizure rate

Many thanks.

[R_cubed]

I’d consider reading the paper mentioned in this thread on the same issue:

Changing the "Status Quo" analyses for Regulatory Submission - Change from Baseline

Two points: First, to anyone who tries the “FDA won’t like it argument”, show them the paper I posted above. In the paper, FDA statisticians considered 5 ways to handle baseline data - they recommeneded 2 of them, neither of which corresponds to the naive “change from baseline” analysis. While this does not decisively refute the “FDA won’t like it” argument, it shows that they may be more open minded and thoughtful than we presume. Second, I do not think technical statistical arguments shoul…

Handling Baselines in Repeated Measures Analyses with Missing Data at Random

[f2harrell]

As a side note, percent change from baseline can only be rescued by undoing the percent change, since it is not legal to do math on percent change for reasons detailed at Statistical Thinking - Why I Don’t Like Percents . So you have to take the percent change, divide it by 100, multiply by pre and add pre to get the original post value.

I would analyze seizure frequency with an ordinal semiparametric model, which allows aribtrary clumping at zero.

[feizuo]

Thank you, Prof Harrell. I’ll also check out your post on Regression Modeling Strategies - 15  Regression Models for Continuous \(Y\) and Case Study in Ordinal Regression before posting additional questions.

[feizuo]

Thank again, Prof Harrell. Please see a few follow-up questions below:

1. For the ordinal semiparametric model, would you recommend the approach below?

require(rms)

f ← orm(outcome ~ base + treatment,data)

where the outcome is the seizure rate during the follow-up period (obtained by the total number of seizure count in the follow-up period divided by the number of days on drug during the follow-up period), and base is the seizure rate during the baseline period pre-randomization. I’m using seizure rate as a continuous outcome to account for the duration of exposure for each subject.

2. How does one examine the PO assumption when the outcome is continuous?

2.1 for example, when we’re still blinded to the treatment would it make sense to include the dummy_treatment variable for the plot below? or should the dummy_treatment variable be taken out and we can examine the PO assumption for other baseline variables to start?

require(Hmisc)

(s ← summary(cumcategory(outcome) ~ base + dummy_treatment + age + sex,data))

2.2 when the outcome is a continuous variable, how are the levels specified? would you look at f$yunique to derive the qlogis(mean(y>= )) levels for the plot?

sf ← function(y){

c(‘Y>=1’=qlogis(mean(y >= the second lowest unique value)), ‘Y>=2’=qlogis(mean(y >= the third lowest unique value)),

…,

‘Y>=# of intercepts’=qlogis(mean(y >=the highest unique value)))

}

plot(s,which=1:# of intercepts,pch=1:# of intercepts,

xlab=‘logit’,vnames=‘names’,main=’ ’

3. I like your approach where the fit for the whole distribution is checked by stratifying on the predicted mean using ols and then checking the agreement between the observed and predicted exceedance probability distribution within each stratum. For this approach, would you recommend trying the different families (i.e., logistic, loglog, probit etc.) and see which one would give the best agreement? Should this method be used for determining the link for the above model f?

4. What if all of the ordinal models using different links show a significant departure of the estimated intercepts from that implied by Gaussian residuals? Is this evidence that these ordinal models would not be a good fit for the data and shouldn’t be used?

Many thanks.

[f2harrell]

I’ll leave most of the questions to further research on your part. I only have time now to comment on the initial part where you mention dividing by days on drug. This implies that seizure count is a longitudinal measurement and should be analyzed longitudinally. If treatment is not constant then interpretation is cloudy unless you have a crossover design. If you have a time-dependent treatment in a non-randomized-crossover design then you’ll need to carefully model time-dependent confounding. Longitudinal models such as a longitudinal proportional odds model (see case study in the last chapter of Regression Modeling Strategies) can do that.

Time intervals for analyzing seizure counts might be days or weeks.

[feizuo]

Thank you, Prof Harrell. I will investigate the longitudinal approach.

This is for a parallel-group RCT, seizure count is recorded on a daily basis throughout the study period, and duration of treatment is theoretically the same for all subjects but some patients may early discontinue the treatment or withdraw consent for having seizure data collected after treatment discontinuation so that not all subjects will have been exposed for the same amount of time and unfortunately the daily seizure count data is not collected after treatment discontinuation. In this case, would you still recommend analyzing the data with a longitudinal proportional odds model?

Is a negative binomial model with an offset term for treatment exposure another reasonable approach for modelling this data?

Thank you.

[f2harrell]

If you use a non-longitudinal approach you can’t have seizure counts missing for some days on-treatment. A longitudinal approach handles missing data and increases power. Don’t know what to do post treatment discontinuation. An analysis that censors at treatment discontinuation will not have that problem. Good to hear that treatment isn’t really time-dependent. If few patients withdrew consent for having seizure data collected things are in even better shape.

[feizuo]

I’m exploring one of the longitudinal ordinal analysis approaches for seizure data using the VGAM package, and started with a simple model below:

formula ← daily_count ~ p_daily_count + baseline.count + group

VGAM::vgam(formula, VGAM::cumulative(parallel = TRUE ~ group, reverse=TRUE), data=long_format)

where data is organized into a long format for the post-baseline daily seizure counts, baseline seizure count is repeated for each subject record, p_daily_count is the seizure count of the previous day. No gap effect is incorporated into the model as post-baseline seizure is recorded daily during the treatment period though there could be missed days for a small number of subjects.

Question 1: It didn’t like relaxing the PO assumption and gave the following warning, so I changed it back to parallel=TRUE for the above model and it worked okay. Do you have other advice for this problem?

Warning: fitted values close to 0 or 1Warning: 333 diagonal elements of the working weights variable ‘wz’ have been replaced by 1.819e-12Error in lm.fit(X.vlm, y = z.vlm, …) : NA/NaN/Inf in ‘x’

Question 2: How should one go about deciding to incorporate gap and/or day (i.e., day since randomization) effects?

Question 3: Both the outcome (daily_count) and the previous state (p_daily_count) should be ordered as factors, correct? I noticed that the estimated coefficients are different depending on whether these are coded as numeric, or as ordered factors.

Many thanks.

[f2harrell]

When Y has more than say 5 levels, vgam is going to choke when you allow non-PO because of so many parameters being added. A Bayesian model using the rmsb package will allow you to use the constrained partial PO model which will add only 1 parameter to the model by assuming that the departure from PO is linear in Y for example.

I’m surprised that you need baseline repeated. Most of the time we’d just let the baseline by p_daily_count for the first record. It would be interesting to see a statistical test on how much the baseline continues to add to the transition model.

Hopefully you can ignore that warning when you are assuming PO.

See COVID-19 for some examples of various gap time analyses.

For vgam Y should be an ordered factor I think. p_daily_count should be a continuous numeric variable and you may want to add a quadratic term for it.

[feizuo]

Thank you, Prof Harrell. I will investigate the examples and the Bayesian approach.

Baseline data were collected during several weeks of the baseline period, so these are also longitudinal daily seizure data.

Question 1: It sounds like for the VGAM approach you’re suggesting that we can simply incorporate the total baseline count (total=20) into the first record of p_daily_count, correct?

Subject ID	Postbaseline daily count (daily_count)	Previous count (p_daily_count)	Day (i.e., after randomization)
1	5	20	1
1	4	5	2
1	2	4	3
…	…	2	…

Question 2: I had difficulty finding an option for extracting the confidence intervals of the coefficients from VGAM.

Do we need to calculate this manually, for example: coef(VGAM_fit) +/- 1.96*sqrt(diag(vcov(VGAM_fit)))?

Thank you.

[f2harrell]

Q1: Correct
Q2: Correct unless summary(fit) gives you something

[feizuo]

For the Bayesian constrained partial PO approach I started with this:

rmsb:blrm(daily_count ~ p_daily_count + group, ppo = ~p_daily_count, cppo = function(daily_count) daily_count, data = long_format, refresh = 100)

• It appears that the outcome variable daily_count needs to be numeric. Model would not run after converting it to an ordered factor.
• Seems to take a bit of time to run (~15-20 min)

It’s mentioned that blrm also allows for random effects with PPO models to handle longitudinal data (i.e., by adding cluster(Subject_ID).

Question 1: Would you keep both p_daily_count and cluster(Subject_ID) or only one of them in the model?

• Keeping both p_daily_count and cluster(Subject_ID) resulted in the following warning, which I’m looking into:

Warning: There were 2 divergent transitions after warmup. See Runtime warnings and convergence problems to find out why this is a problem and how to eliminate them.Warning: Examine the pairs() plot to diagnose sampling problems

Question 2: There isn’t an option for incorporating random effects for the VGAM approach, correct?

Question 3: Is the difference between VGAM and VGLM that: VGAM handles smoothing, but if no smoothing is introduced, VGAM(formula, …) and VGLM(formula, …) are the same?

Many thanks.

[f2harrell]

Often what is not in proportional odds is time and not Y, but I think you have a special reason for not trusting PO for previous daily count. Yes it needs to be numeric, and the time you see is typical.

Most of the serial correlation pattern will be absorbed by the previous count so I add random effects only as a test of the Markov assumption.

Correct about VGAM.

Don’t know about the smoothing question. I advise using vglm.

[feizuo]

The Markov PO ordinal logistic model using VGLM (as below) gave a warning for detecting a Hauck-Donner effect (HDE) on the intercepts.

fit1 ← VGAM::vglm(formula, VGAM::cumulative(parallel = TRUE , reverse=TRUE), data=long_format)

Formula ← daily_count ~ p_daily_count + group

Question 1: is the warning acceptable or should we worry about the treatment estimates and p values from the above model?

Note: Using a quadratic or cubic form of p_daily_count didn’t affect the treatment group coefficients and SEs of fit1.

For fit2, I tried relaxing the PO assumption by incorporating VGAM::cumulative(parallel = TRUE ~ group, …), and VGLM was able to handle it better than VGAM but gave the following warning:

fitted values close to 0 or 1. Warning: It seems that the nonparallelism assumption has resulted in intersecting linear/additive predictors. Try propodds() or fitting a partial nonproportional odds model or choosing some other link function, etc.

Question 2: fit1 is essentially propodds() and resulted in the HDE warning on the intercepts; fit2 is a partial non-PO model; so the only remaining option is trying other links but I’m not sure if anything other than the logit link would make sense.

Question 3: is vglm preferred over vgam here because like you said: when Y has more than 5 levels, vgam will have issues when allowing for non-PO? Some examples from COVID-19 used vgam instead of vglm. I guess it didn’t matter for those examples where the Y levels were less than 5, correct?

Many thanks.

[f2harrell]

These will be the last questions I have time to answer for a while.

Ignore the warning except don’t use standard error estimates that are huge.

State the model in terms of parallel = FALSE ~ ... to better show what variables have PO relaxed. With VGAM when you don’t put constraints (can’t remember if VGAM allows contraints) the non-PO can only be with regard to discrete Y.

The VGAM author said I should be using vglm so I didn’t question that.

[feizuo]

Thank you so much for all your help, Prof Harrell!