2
$\begingroup$

I am conducting a study of a number of patients with a disease, and using an ordinal scale assessment of functional status at 3 different time points.

In this example we can say that I have 100 patients, and each one is measured at the three time points (no need for imputation).

Please see the graphical representation of the change in the frequencies of each of these 6 ordinal values across the three time points (top group has no patients with ordinal score 3,5,6):

enter image description here

I want to compare the overall difference between time point 1 and time point 2, time point 1 and time point 3, and time point 2 and timepoint 3 to see if the shifts observed are significant.

Additionally, if there is a way to determine changes between the specific ordinal levels that may be helpful as well.

I am planning on doing the analysis in R, but would like to make sure that I am using the appropriate test. Would some form of ordinal regression be appropriate here? Or is some version of a Wilcoxin rank sum better for this study paradigm?

Thank you for any assistance.

*UPDATE WITH DATASET

I have still been going through working to analyze this the best way.

Code for my dataset is below for reference:

mrsstat <-tibble(
  Patient = c(1:21),
  pMRS = c(0, 0,    0,  0,  0,  0,  4,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  1,  0,  2,  0),
  dMRS = c(4, 2,    2,  4,  3,  6,  4,  1,  5,  0,  3,  0,  6,  1,  5,  3,  3,  1,  3,  4,  3),
  fMRS = c(1, 2,    2,  3,  2,  6,  3,  1,  2,  0,  3,  0,  6,  1,  5,  2,  0,  1,  3,  4,  0))

mrsstat<- mrsstat %>% 
  pivot_longer(cols = c("pMRS", "dMRS", "fMRS"), names_to = "timepoint") %>% 
  mutate(timepoint = factor(timepoint, 
                            levels= c("fMRS", 
                                      "dMRS",
                                      "pMRS")),
         value=  factor(value, levels = c("0", "1","2","3","4","5","6")),
         Patient= factor(Patient, levels = c(1:21)))

> mrsstat
# A tibble: 63 × 3
   Patient timepoint value
   <fct>   <fct>     <fct>
 1 1       pMRS      0    
 2 1       dMRS      4    
 3 1       fMRS      1    
 4 2       pMRS      0    
 5 2       dMRS      2    
 6 2       fMRS      2    
 7 3       pMRS      0    
 8 3       dMRS      2    
 9 3       fMRS      2    
10 4       pMRS      0    
# … with 53 more rows

I attempted to do an ordinal regression analysis with methods using the ordinal package et al. in R as described here:

Two-way Ordinal Regression with CLM

I am not sure that this is feasible with my dataset as each patient has only 1 recording of the functional outcome score (variable labled "value") at each timepoint, so the when attempting to get to an ANOVA analysis of the data as seen in the link above for two-way ordinal regression, it wont work as there is obviously no variance with one measure per patient per timepoint.

I still feel like there should be some test to assess for changes in significance between the "value" proportions between the 3 timepoints and for multiple comparisons between the ordinal "value" groups.

Improve this question
$\endgroup$
7
  • $\begingroup$ Consider a random effects proportional odds model with time represented with 2 indicator variables. You can form time differences using these indicators to assess what you want. $\endgroup$
    – Frank Harrell
    Jan 24, 2022 at 17:27
  • $\begingroup$ Thank you for this suggestion- Im still learning about some of these more in depth statistical analyses- based on the dataset above do you feel this would still be an option given the single observation for each patient at one of the three time points? $\endgroup$
    – XFrost
    Jan 28, 2022 at 14:42
  • $\begingroup$ You originally stated that each patient is measured at all 3 time points so I'm not following this. $\endgroup$
    – Frank Harrell
    Jan 28, 2022 at 14:47
  • $\begingroup$ Sorry Frank, yes each patient is measured at each time point. I meant that each time point only has one measurement per patient. Using the link above for the two-way ordinal regression, there isnt a way I was finding to analyze the data using an ANOVA as it was set up a as follows: ``` model <- clm(value ~ Patient + timepoint + Patient:timepoint, data = mrsstat, threshold = "symmetric") anova(model, type = "II") ``` . Thank you for your assistance and insight $\endgroup$
    – XFrost
    Jan 28, 2022 at 14:48
  • $\begingroup$ clm does not allow for random effects. Try the clmm function from the ordinal package. Something like this: clmm(value~timepoint + (1|Patient), data = mrsstat). $\endgroup$
    – COOLSerdash
    Jan 28, 2022 at 15:00

2 Answers 2

Reset to default
3
$\begingroup$

Meta comment 1: This thread is essentially about using the ordinal package in R, not about statistics per se.

Meta comment 2: This post is in response to the discussion under the original post as well as the response posted by @XFrost .

@XFrost , it sounds like you are trying to use the summary() function in ways it's not really designed for. Instead, you can use the anova() function and emmeans package to get the information I think you're looking for.

BTW, since you cited my webpage earlier, I'll note that the relevant example for the clmm model is here: rcompanion.org/handbook/G_08.html (As always, with the caveat that I wrote it.)

To continue the analysis in the response by @XFrost :

A p-value for the effect of timepoint can be determined with the anova() function:

    model.null = clmm(value ~ 1 + (1|Patient), data = mrsstat)

    anova(model, model.null)

       ### Likelihood ratio tests of cumulative link models:
       ###
       ###            no.par    AIC   logLik LR.stat df Pr(>Chisq)    
       ### model.null      7 233.74 -109.870                          
       ### model           9 200.27  -91.135  37.472  2  7.296e-09 ***

And the pairwise comparisons of timepoint can be determined with the emmeans package:

    library(emmeans)
    
    marginal = emmeans(model, ~ timepoint)
    
    pairs(marginal, adjust="tukey")
    
    ### contrast    estimate    SE  df z.ratio p.value
    ###      pMRS - dMRS    -4.81 1.029 Inf  -4.672  <.0001
    ###      pMRS - fMRS    -3.75 0.931 Inf  -4.031  0.0002
    ###      dMRS - fMRS     1.06 0.584 Inf   1.808  0.1668
    ###     
    ###     P value adjustment: tukey method for comparing a family 
    ###      of 3 estimates
Improve this answer
$\endgroup$
6
  • $\begingroup$ @ SalMagniafico -thank you so much for your comments and assistance, the use of the ordinal and emmeans package as you showed above was incredibly helpful. I also greatly appreciated your excellent website. I know that this is certainly an R heavy question, but I also would appreciate any guidance on the appropriate statistical analysis of this question as indicated in the question title, I.e., should analysis of this ordinal data across 3 timepoints be done with an ordinal regression as weve done here or with some other method? Also are multiple comparisons possible btw the Score groups? $\endgroup$
    – XFrost
    Jan 31, 2022 at 14:10
  • 1
    $\begingroup$ Thanks for the kind words. ... Yes, with an ordinal dependent variable with seven levels, ordinal regression is probably the best approach. Using the mixed effects approach, with Patient as a random effect makes sense. ... I don't understand the last question: Also are multiple comparisons possible btw the Score groups?. Does the emmeans output not address this question? $\endgroup$
    – Sal Mangiafico
    Jan 31, 2022 at 17:10
  • 1
    $\begingroup$ As a side note, the design the OP presents is an unreplicated complete block design. In this example, Friedman's test with an appropriate post-hoc test works just as well as the ordinal regression approach. I believe Friedman's test works for an ordinal dependent variable. ... If the dependent variable were assumed to be interval in nature, rather than merely ordinal, other approaches would be viable as well. $\endgroup$
    – Sal Mangiafico
    Jan 31, 2022 at 17:23
  • 1
    $\begingroup$ @XFrost , If I understand your question, I think you are switching from treating Value as an ordinal variable to a categorical variable. You could do an analysis this way as well. But I think a plot like what you presented in your original question, or using bar plots or plots commonly used for Likert-type items might tell the story. (See next comment for examples). If you want to look at if people changed from one category of response to another, you might use something like McNemar's test between two or multiple response categories. $\endgroup$
    – Sal Mangiafico
    Feb 1, 2022 at 10:50
  • 1
    $\begingroup$ Bar plot of ordinal responses example: rcompanion.org/handbook/images/image052.png. I don't know if this plot has a specific name, but it is commonly used for Likert-type items: reganmian.net/blog/images/2013-10-02-likert-graphs-in-r-embedding-metadata-for-easier-plotting-whole-04.png. $\endgroup$
    – Sal Mangiafico
    Feb 1, 2022 at 10:52
0
$\begingroup$

So Ive continued to look at this further based on @FrankHarrell and @COOLSerdash's comments above.

I set the data analysis up attempting to use a Cumulative Link Mixed Model approach to analyze.

The Code below in R is what Ive got so far using data set above:

mrsstat <-tibble(
  Patient = c(1:21),
  pMRS = c(0, 0,    0,  0,  0,  0,  4,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  1,  0,  2,  0),
  dMRS = c(4, 2,    2,  4,  3,  6,  4,  1,  5,  0,  3,  0,  6,  1,  5,  3,  3,  1,  3,  4,  3),
  fMRS = c(1, 2,    2,  3,  2,  6,  3,  1,  2,  0,  3,  0,  6,  1,  5,  2,  0,  1,  3,  4,  0))

mrsstat<- mrsstat %>% 
  pivot_longer(cols = c("pMRS", "dMRS", "fMRS"), names_to = "timepoint") %>% 
  mutate(timepoint = factor(timepoint, 
                            levels= c("pMRS", 
                                      "dMRS",
                                      "fMRS")),
         value=  factor(value, levels = c("0", "1","2","3","4","5","6")),
         Patient= factor(Patient, levels = c(1:21)))


mrs_numeric <- mrsstat %>% 
  mutate(value= as.integer(value))

Summarize(value ~ Patient, 
          data= mrs_numeric, digits=3)

Sum= Summarize(value ~ Patient + timepoint, 
          data= mrs_numeric, digits=3)


model <- clmm(value~timepoint + (1|Patient), data = mrsstat)
summary(model)

which yields:

Cumulative Link Mixed Model fitted with the Laplace approximation

formula: value ~ timepoint + (1 | Patient)
data:    mrsstat

 link  threshold nobs logLik AIC    niter     max.grad cond.H 
 logit flexible  63   -91.13 200.27 476(1772) 5.03e-06 1.9e+02

Random effects:
 Groups  Name        Variance Std.Dev.
 Patient (Intercept) 2.374    1.541   
Number of groups:  Patient 21 

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
timepointdMRS    4.809      1.029   4.672 2.98e-06 ***
timepointfMRS    3.753      0.931   4.031 5.56e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
    Estimate Std. Error z value
0|1   1.8433     0.7538   2.445
1|2   3.1043     0.8823   3.519
2|3   4.1262     0.9839   4.194
3|4   5.5329     1.1465   4.826
4|5   6.6601     1.2936   5.149
5|6   7.4873     1.4107   5.307

If I change the reference (first timepoint) to be "fMRS" instead of "pMRS", then the same approach gives the following:

Cumulative Link Mixed Model fitted with the Laplace approximation

formula: value ~ timepoint + (1 | Patient)
data:    mrsstat

 link  threshold nobs logLik AIC    niter     max.grad cond.H 
 logit flexible  63   -91.13 200.27 448(1671) 5.04e-06 4.9e+01

Random effects:
 Groups  Name        Variance Std.Dev.
 Patient (Intercept) 2.374    1.541   
Number of groups:  Patient 21 

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
timepointdMRS   1.0564     0.5842   1.808   0.0706 .  
timepointpMRS  -3.7528     0.9310  -4.031 5.56e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
    Estimate Std. Error z value
0|1  -1.9095     0.6643  -2.874
1|2  -0.6485     0.5744  -1.129
2|3   0.3734     0.5692   0.656
3|4   1.7801     0.6638   2.682
4|5   2.9073     0.7934   3.664
5|6   3.7345     0.9149   4.082

which from what I can surmise indicates that "pMRS" is significantly different than both "dMRS" (p<2.98e-06) and "fMRS" (5.56e-05), while "dMRS" and "fMRS" are not significantly different (p= 0.0706).

Still unclear on if multiple comparisons are possible between the timepoints for each score, or if this approach appropriately addresses the question (namely, that the summary is comparing the timepoints to single reference without possibly comparing all group?)

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ When you are posting code examples, please include some indication of the packages that need to be loaded. As far as I can tell, this code uses the following packages: tibble, dplyr, tidyr, ordinal, FSA. $\endgroup$
    – Sal Mangiafico
    Jan 30, 2022 at 15:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.