Received: 9 July 2016
Revised: 23 January 2017
Accepted: 24 January 2017
DOI 10.1111/jep.12729
ORIGINAL ARTICLE
Evaluating therapeutic effect on WOMAC subscales in osteoarthritis RCTs: When model choice matters Giulia Lorenzoni MA1
|
Danila Azzolina MA1
|
Nicola Soriani PhD2
|
Dario Gregori PhD3
1
PhD student, Unit of Biostatistics, Epidemiology and Public Health, Department of Cardiac, Thoracic and Vascular Sciences, University of Padova, Padova, Italy 2 Post-Doc Fellow, Unit of Biostatistics, Epidemiology and Public Health, Department of Cardiac, Thoracic and Vascular Sciences, University of Padova, Padova, Italy 3 Full Professor, Unit of Biostatistics, Epidemiology and Public Health, Department of Cardiac, Thoracic and Vascular Sciences, University of Padova, Padova, Italy Correspondence Prof Dario Gregori, Unit of Biostatistics, Epidemiology and Public Health, Department of Cardiac, Thoracic and Vascular Sciences, University of Padova, Via Loredan, 18, 35121 Padova, Italy. Email: [email protected]
Abstract Rationale, Aims, and Objectives:
The study aimed at developing a method for modelling
the Western Ontario and McMaster Universities index (WOMAC), accounting for correlation between its subscales and for heterogeneity of treatment effect (HTE), using data from 2 twin trials on knee osteoarthritis.
Method:
Two randomized, double‐blind, placebo‐controlled clinical trials (twin trials). Studies
aimed at investigating the effectiveness of a pharmacological treatment on clinical outcomes of knee osteoarthritis, measured using WOMAC index. To take into account that the WOMAC subscales are correlated and skewed, we proposed and compared multivariate gamma and Gaussian approaches with latent variable capturing correlation between outcomes. Besides the latent term, the interaction between the latent term and treatment, accounting for HTE, was further estimated.
Results:
Modelling the subscales by using a gamma approach accounting for skewness of data,
we found out different results compared with Gaussian models. The main difference regarded the latent variable interacting with treatment (accounting for unobserved heterogeneity), which is not significant for the Gaussian approach (P value = .102) and significant in the gamma model (P value < .002). Thus, indicating that unobserved covariates affect treatment's performance. Additionally, plotting the observed and the estimated values of WOMAC index of the Gaussian and gamma models, we showed that, compared with the Gaussian, the gamma one best fits the data, especially among poor responders.
Conclusion:
Multivariate gamma approach accounting for correlation between outcomes and
for HTE has been demonstrated to be more suitable to model WOMAC subscales and to provide more information on effect of therapy. KEY W ORDS
gamma approach, Gaussian approach, latent variable, WOMAC
1
|
I N T RO D U CT I O N
rheumatic diseases, and comorbidities are significantly associated with higher risk of developing osteoarthritis.
The most prevalent joint disorder in the United States is osteoarthritis,
Regarding symptoms, those that occur more often are pain, joint
with more than 27 million people suffering of this degenerative
stiffness, and impaired physical function, which are assessed using
joint disease.1
the Western Ontario and McMaster Universities index (WOMAC),6 2–5
Factors associated with osteoarthritis onset are heterogeneous
an instrument consisting on 3 subscales accounting for the symptoms
and are commonly categorized as modifiable and nonmodifiable.
mentioned above. Higher scores on the WOMAC indicate worse pain,
Nonmodifiable risk factors include genetic profile, female gender,
stiffness, and functional limitations.
aging, and ethnicity. Among modifiable risk factors, studies have found
Different treatment options are available for osteoarthritis
that obesity, joint injuries (sport, occupational or traumatic), job type
management,7 and several studies have demonstrated that they
(causing excessive and continuous mechanical stresses), other
possibly have different impact on osteoarthritis' symptoms (therefore,
J Eval Clin Pract. 2017;1–8.
wileyonlinelibrary.com/journal/jep
© 2017 John Wiley & Sons, Ltd.
1
2
LORENZONI
on WOMAC subscales), because of different mechanisms of action and
concomitant
patients' heterogeneous characteristics. Regarding pharmacological
osteoarthritis (years) were recorded.
medications,
and
time‐from‐diagnosis
of
ET AL.
knee
treatments, Bannuru et al8 performed a meta‐analysis demonstrating
Signs and symptoms of knee osteoarthritis were also deepening
that the interventions significantly differ from placebo on pain but
analyzed. At baseline, physical examination of the knees was
not on stiffness; Wang et al8 have shown that nonpharmacological
performed, including knees' degree of flextion, presence of joint
treatments better perform on physical function but not on stiffness.
inflammation (temperature increasing, palpable effusion), presence of
Given the heterogeneous characteristics of subjects suffering of osteoarthritis affecting the performance of treatments on WOMAC
crepitus on passive and active motion, and knees' alignment (valgus vs varo).
subscales, authors claim the need of developing personalized strategies for the management of such diseases,8 identifying the best treatment option according to patients' features. However, studies
2.2
|
Outcomes' assessment
investigating the effectiveness of treatments for knee osteoarthritis generally evaluate separately interventions' effect on each symptom, thus assuming independence of such outcomes. This is of particular interest when treatment has heterogeneous effect on 1 or more symptoms according to some characteristics of the patient (this phenomenon, according to which different patients respond differently to the same treatment, is known as heterogeneity of treatment effect). In this case, the independent modelling approach usually adopted in literature foresees a set of tests for interaction of treatment with 1 covariate and only 1 outcome at a time. This approach does not take into account outcomes' correlation and not allow for the identification of subgroups of patients, which, based on their characteristics,
Symptoms' modification during the observation period was assessed by administering (at baseline and at follow‐up visits) the WOMAC index.6 It is a tridimensional, disease‐specific, health status measure, assessing symptoms in the areas of pain (5 questions), stiffness (2 questions), and physical function (17 questions) in patients with osteoarthritis of the hip and/or knee. The visual analogue scale version of the index was used in both studies (100 mm visual analogue scale for each question, and the total score is represented by the sum of the 24 items scores). A higher WOMAC score represents worse symptoms, with 2400 mm being the worst possible total score. For the purpose of the study, we considered the absolute levels of
show a heterogeneous behavior after treatment on more than 1
WOMAC subscales (pain, stiffness, and physical function), which
outcome simultaneously. The aim of the study is to develop a method for modelling
represent often a long‐term benchmark for therapeutic approach.
WOMAC subscales accounting for correlation between the considered outcomes and for heterogeneity of treatment effect (HTE), using data
2.3
from 2 twin trials on knee osteoarthritis.
|
Statistical analysis
Analyses were performed on a matched‐pair cohort resulting from the
2
pooled twin trials. The matching procedure was performed on the basis
METHODS AND MATERIALS
|
of a propensity score built on the variables collected at baseline using
Two randomized, double‐blind, placebo‐controlled clinical trials (twin trials). Studies aimed at investigating the effectiveness of a pharmacological treatment on clinical outcomes of knee osteoarthritis. Inclusion criteria were age over 50 years and diagnosis of primary knee osteoar-
k‐nearest algorithm. It has been demonstrated that matched‐pair analysis can be effectively used to reduce residual variability and improve stability of estimates10,11 even in an RCT setting. Analyses were performed using R system12 and Hglm13 package.
thritis according to the clinical and radiological criteria of the American College of Rheumatology.9 The
pharmacological
treatment
and
the
placebo
were
2.4
|
WOMAC subscales' density function
manufactured and packaged so that they were indistinguishable.
Several studies have demonstrated that WOMAC subscales are
Patients
random
skewed distributed (asymmetrical distributed with a tail on the right
numbers. The treatment code of each patient was kept in a single
corresponding to poor responders).14,15 Assuming that the subscales
opaque sealed envelope to be opened only in case of emergency.
are normally distributed, when they are not (using improper modelling
were
randomized
using
computer‐generated
approach), may affect the power and increases the chance of Type I
2.1
|
error,16 by failing in the identification of subjects with extreme
Variables collected at baseline
values at the tail of the distribution (with consequently difficulties in
In both Randomized Clinical Trial (RCTs) patients' medical history was
characterizing poor responders in the case of WOMAC subscales).
assessed at baseline. Besides demographic characteristics (age and
After examining WOMAC subscales for normality using kernel
gender), patients' habits were also assessed, referring to: smoking habit
density estimation based on Epanechnikov kernel,17 we tested for
(classified as smokers vs nonsmokers), alcohol consumption (consumers
nonnormality using the Henze‐Zirkler's Multivariate Normality Test18
vs nonconsumers), caffeine consumption (consumers vs nonconsumers),
on the WOMAC distribution over the entire follow‐up.
and diet type (normal vs special, eg, hypocaloric diet).
Eventually, we propose a gamma‐based model for WOMAC
Referring to clinical status, patients' BMI, blood pressure, heart rate,
concomitant
infectious
diseases
(presence
vs
absence),
subscale, comparing it with the most widely used (in the context of osteoarthritis research) Gaussian approach.
LORENZONI
3
ET AL.
2.5 | Modelling WOMAC subscales' taking into account their potential correlation Since WOMAC subscales measure related quantities (pain, stiffness,
the latent variable. Such latent variable is assumed to have a normal distribution with mean 0 and variance σ2.
to be correlated to each other. However, generally, studies on
2.6 | Modelling WOMAC subscales taking into account HTE
osteoarthritis ignore the correlation between the outcomes and fit a
Given the highly heterogeneous characteristics of patients suffering
separate model to each response variable. Despite that this approach
from knee osteoarthritis affecting treatment performance on different
seems to be simple to apply, it completely not accounts for the
WOMAC subscales, we estimate a further model including an interac-
correlation between outcomes, determining a reduction of power in
tion term between the latent variable and treatment (representing the
detecting
misleading
unobserved covariates influencing treatment performance). This is to
estimates.19 To overcome this limitation, after testing WOMAC
account for unobserved heterogeneity affecting the likelihood of
subscales for correlation, we develop a multivariate model, which is
responding to treatment and identifying subsets of patients who are
an extension of Sammel latent variable model.20 Conditional on this
poor or high responders, sharing the same features.
and physical function) in the same individuals, probably outcomes tend
treatment
effectiveness
and
leading
to
latent variable, the outcomes are assumed to be independent and are
Summarizing, after testing the subscales for normality and
modelled as functions of fixed covariates and a subject‐specific latent
correlation, both gamma and Gaussian models (to compare each other)
variable. This latent variable establishes the link between the
were estimated through the following steps:
regression equations capturing correlations among the outcomes. Assuming that the latent variable completely specifies the correlation
• Gamma and Gaussian models were estimated ignoring the
among the outcomes, this permits examination of the outcomes as
correlation between the outcomes and fit a separate model to
independent of each other by accounting for the correlation through
each response variable (independent).
TABLE 1
Distribution of the 2 twin trials combined
Age Sex: M BMI Diatolic pressure Heart rate Systolic pressure
P value
Placebo
Treatment
Combined
(N = 207)
(N = 207)
(N = 414)
60.00/65.00/69.00
58.71/64.00/68.50
59.00/64.00/69.00
23% (47)
23% (48)
23% (95)
.907
25/27/28.625
25/27/28.075
25/27/28.340
.849
80.00/85.00/90.00
75.00/80.00/90.00
76.25/80.00/90.00
.611
64/70/80
64/72/80
64/72/80
.819 .277
.187
130/140/150
125/140/150
130/140/150
Smoking: Y
19% (39)
14% (30)
17% (69)
.235
Alcohol: Y
47% (98)
48% (100)
48% (198)
.844
Caffeine: Y
86% (179)
86% (177)
86% (356)
.777
Diet
17% (35)
17% (36)
17% (71)
.896
Infectious diseases: Y
26% (53)
16% (34)
21% (87)
.022
Concomitant medications: Y
62% (129)
65% (134)
64% (263)
.61
Presence of knee osteoarthritis (years)
4/8/12
4/7/10
4/ 7/11
.596
Increased temperature (left knee): Y
21% (43)
21% (44)
21% (87)
.904
Effusion (left knee): Y
26% (54)
22% (46)
24% (100)
.358
Effusion righ knee: Y
29% (59)
33% (69)
31% (128)
.288
Alignment (left knee): normal
66% (136)
74% (154)
70% (290)
.01
20% (42)
20% (42)
20% (84)
valgus varus
14% (29)
5% (11)
10% (40)
70% (145)
72% (150)
71% (295)
valgus
16% (34)
19% (40)
18% (74)
varus
14% (28)
8% (17)
11% (45)
Crepitus active (left knee): Y
77% (160)
78% (161)
78% (321)
.906
Crepitu active (right knee): Y
75% (155)
83% (172)
79% (327)
.04
Crepitus passive (left knee): Y
72% (150)
70% (145)
71% (295)
.587
Alignment (right knee): normal
.196
Crepitu passive (right knee): Y
68% (141)
75% (156)
72% (297)
.102
Degree of flexion (left knee)
70/90/135
68/95/140
70/90/140
.71
Degree of flexion (right knee)
70.00/95.00/135.00
70.00/90.00/140.00
70.00/90.00/138.75
.845
Data are I quartile/median/III quartile for continuous variables and percentage (absolute numbers) for categorical variables. P value refers to the difference between placebo and treatment.
4
LORENZONI
ET AL.
• Multigamma and Gaussian models were estimated including a
stiffness and physical function subscales, (higher scores represent
latent term (latent only) accounting for correlation between
worse symptoms, a reduction of WOMAC scores corresponds to an
outcomes.
improvement of symptoms) over the time period considered
• Multigamma and Gaussian models were estimated including,
(Figure 1). The analysis of such subscales revealed that they were
besides the latent term, the interaction between the latent term
skewed (Figure 2) and highly correlated between each other, with
and treatment (latent interaction) accounting for HTE.
correlations coefficients ranging from 0.71 up to 0.87. Table 2 shows the effect of treatment on the considered including
outcomes in the “independent”, “latent only,” and “latent interaction”
covariates resulting significant by inspecting AIC at a significance level
gamma and Gaussian models. Modelling the subscales using a gamma
of 0.10 and by clinical judgement. Variables considered as candidates
approach accounting for skewness of data, we found out different
for the predictor part of the model are listed in Table 1.
results compared with Gaussian models. For independent models, in
21,22
Multivariate models were fitted using H‐likelihood,
the Gaussian approach, the treatment resulted not significant on all WOMAC subscales while, using the gamma approach, it resulted
3
|
RESULTS
to be significant for pain (P value = .021) and physical function (P value = .024). The main difference regarded the latent
Subjects included in the analysis from the pooled twin trails were 414
variable interacting with treatment (accounting for unobserved
(Table 1). Scores of WOMAC subscales improved, especially for
heterogeneity), which is not significant for the Gaussian approach
FIGURE 1
Trends of Western Ontario and McMaster Universities index (WOMAC) subscales over time. Line is a loss fit for estimating trend
over time
FIGURE 2
Density estimation of Western Ontario and McMaster Universities index (WOMAC) subscales distribution. Henze‐Zirkler's multivariate normality test (from left to right): 69.51 (P value < .001), 68.86 (P value < .001), and 71.36 (P value < .001)
LORENZONI
5
ET AL.
TABLE 2 Gaussian and gamma models, by fitting (1) three independent equations, (2) a multivariate model with 1 latent variable, and (3) a latent and a latent‐by‐treatment interaction
Independent
Latent only
Latent interaction
Gaussian model
Some further insights can be given by examining the estimated latent variable, as emerging from the gamma model with latent treatment interaction. First, as shown in Figure 6, the distribution of the latent variable is shifted to the right for treated vs placebo patients, as perhaps expected by the goodness‐of‐fit analysis discussed above. Second, by examining correlation between latent and the other
Coeff.
P value Coeff.
P value Coeff.
P value
variables, it clearly emerges (Figure 7) that the 2 groups are characterized
Pain
−0.110 .104
−0.082
.181
−0.115 .0007
by different variables. Particularly, lifestyle variables (such as caffeine
Physical function
−0.082 .157
−0.101
.099
−0.128 .002
consumption and smoking habit) differently correlate with latent variable
Stiffness
−0.095 .157
−0.055
.363
−0.101 .015
in the 2 groups of patients (smoking resulted to be directly correlated
.102
with latent variable in treated subjects and inversely correlated in the
—
Latent Bayesian information criterion (BIC)
placebo group, conversely caffeine is inversely correlated with latent variable in the treated group and directly correlated in the placebo group).
Gamma model Coeff.
P value Coeff.
P value Coeff.
Pain
−0.157 .021
−0.173
.038
−0.173 .006
Physical function
−0.139 .024
−0.157
.061
−0.17
Stiffness
−0.111 .095
−0.127
.128
−0.123 .057
Latent BIC
—
In summary,
P value
.007
Revised: 23 January 2017
Accepted: 24 January 2017
DOI 10.1111/jep.12729
ORIGINAL ARTICLE
Evaluating therapeutic effect on WOMAC subscales in osteoarthritis RCTs: When model choice matters Giulia Lorenzoni MA1
|
Danila Azzolina MA1
|
Nicola Soriani PhD2
|
Dario Gregori PhD3
1
PhD student, Unit of Biostatistics, Epidemiology and Public Health, Department of Cardiac, Thoracic and Vascular Sciences, University of Padova, Padova, Italy 2 Post-Doc Fellow, Unit of Biostatistics, Epidemiology and Public Health, Department of Cardiac, Thoracic and Vascular Sciences, University of Padova, Padova, Italy 3 Full Professor, Unit of Biostatistics, Epidemiology and Public Health, Department of Cardiac, Thoracic and Vascular Sciences, University of Padova, Padova, Italy Correspondence Prof Dario Gregori, Unit of Biostatistics, Epidemiology and Public Health, Department of Cardiac, Thoracic and Vascular Sciences, University of Padova, Via Loredan, 18, 35121 Padova, Italy. Email: [email protected]
Abstract Rationale, Aims, and Objectives:
The study aimed at developing a method for modelling
the Western Ontario and McMaster Universities index (WOMAC), accounting for correlation between its subscales and for heterogeneity of treatment effect (HTE), using data from 2 twin trials on knee osteoarthritis.
Method:
Two randomized, double‐blind, placebo‐controlled clinical trials (twin trials). Studies
aimed at investigating the effectiveness of a pharmacological treatment on clinical outcomes of knee osteoarthritis, measured using WOMAC index. To take into account that the WOMAC subscales are correlated and skewed, we proposed and compared multivariate gamma and Gaussian approaches with latent variable capturing correlation between outcomes. Besides the latent term, the interaction between the latent term and treatment, accounting for HTE, was further estimated.
Results:
Modelling the subscales by using a gamma approach accounting for skewness of data,
we found out different results compared with Gaussian models. The main difference regarded the latent variable interacting with treatment (accounting for unobserved heterogeneity), which is not significant for the Gaussian approach (P value = .102) and significant in the gamma model (P value < .002). Thus, indicating that unobserved covariates affect treatment's performance. Additionally, plotting the observed and the estimated values of WOMAC index of the Gaussian and gamma models, we showed that, compared with the Gaussian, the gamma one best fits the data, especially among poor responders.
Conclusion:
Multivariate gamma approach accounting for correlation between outcomes and
for HTE has been demonstrated to be more suitable to model WOMAC subscales and to provide more information on effect of therapy. KEY W ORDS
gamma approach, Gaussian approach, latent variable, WOMAC
1
|
I N T RO D U CT I O N
rheumatic diseases, and comorbidities are significantly associated with higher risk of developing osteoarthritis.
The most prevalent joint disorder in the United States is osteoarthritis,
Regarding symptoms, those that occur more often are pain, joint
with more than 27 million people suffering of this degenerative
stiffness, and impaired physical function, which are assessed using
joint disease.1
the Western Ontario and McMaster Universities index (WOMAC),6 2–5
Factors associated with osteoarthritis onset are heterogeneous
an instrument consisting on 3 subscales accounting for the symptoms
and are commonly categorized as modifiable and nonmodifiable.
mentioned above. Higher scores on the WOMAC indicate worse pain,
Nonmodifiable risk factors include genetic profile, female gender,
stiffness, and functional limitations.
aging, and ethnicity. Among modifiable risk factors, studies have found
Different treatment options are available for osteoarthritis
that obesity, joint injuries (sport, occupational or traumatic), job type
management,7 and several studies have demonstrated that they
(causing excessive and continuous mechanical stresses), other
possibly have different impact on osteoarthritis' symptoms (therefore,
J Eval Clin Pract. 2017;1–8.
wileyonlinelibrary.com/journal/jep
© 2017 John Wiley & Sons, Ltd.
1
2
LORENZONI
on WOMAC subscales), because of different mechanisms of action and
concomitant
patients' heterogeneous characteristics. Regarding pharmacological
osteoarthritis (years) were recorded.
medications,
and
time‐from‐diagnosis
of
ET AL.
knee
treatments, Bannuru et al8 performed a meta‐analysis demonstrating
Signs and symptoms of knee osteoarthritis were also deepening
that the interventions significantly differ from placebo on pain but
analyzed. At baseline, physical examination of the knees was
not on stiffness; Wang et al8 have shown that nonpharmacological
performed, including knees' degree of flextion, presence of joint
treatments better perform on physical function but not on stiffness.
inflammation (temperature increasing, palpable effusion), presence of
Given the heterogeneous characteristics of subjects suffering of osteoarthritis affecting the performance of treatments on WOMAC
crepitus on passive and active motion, and knees' alignment (valgus vs varo).
subscales, authors claim the need of developing personalized strategies for the management of such diseases,8 identifying the best treatment option according to patients' features. However, studies
2.2
|
Outcomes' assessment
investigating the effectiveness of treatments for knee osteoarthritis generally evaluate separately interventions' effect on each symptom, thus assuming independence of such outcomes. This is of particular interest when treatment has heterogeneous effect on 1 or more symptoms according to some characteristics of the patient (this phenomenon, according to which different patients respond differently to the same treatment, is known as heterogeneity of treatment effect). In this case, the independent modelling approach usually adopted in literature foresees a set of tests for interaction of treatment with 1 covariate and only 1 outcome at a time. This approach does not take into account outcomes' correlation and not allow for the identification of subgroups of patients, which, based on their characteristics,
Symptoms' modification during the observation period was assessed by administering (at baseline and at follow‐up visits) the WOMAC index.6 It is a tridimensional, disease‐specific, health status measure, assessing symptoms in the areas of pain (5 questions), stiffness (2 questions), and physical function (17 questions) in patients with osteoarthritis of the hip and/or knee. The visual analogue scale version of the index was used in both studies (100 mm visual analogue scale for each question, and the total score is represented by the sum of the 24 items scores). A higher WOMAC score represents worse symptoms, with 2400 mm being the worst possible total score. For the purpose of the study, we considered the absolute levels of
show a heterogeneous behavior after treatment on more than 1
WOMAC subscales (pain, stiffness, and physical function), which
outcome simultaneously. The aim of the study is to develop a method for modelling
represent often a long‐term benchmark for therapeutic approach.
WOMAC subscales accounting for correlation between the considered outcomes and for heterogeneity of treatment effect (HTE), using data
2.3
from 2 twin trials on knee osteoarthritis.
|
Statistical analysis
Analyses were performed on a matched‐pair cohort resulting from the
2
pooled twin trials. The matching procedure was performed on the basis
METHODS AND MATERIALS
|
of a propensity score built on the variables collected at baseline using
Two randomized, double‐blind, placebo‐controlled clinical trials (twin trials). Studies aimed at investigating the effectiveness of a pharmacological treatment on clinical outcomes of knee osteoarthritis. Inclusion criteria were age over 50 years and diagnosis of primary knee osteoar-
k‐nearest algorithm. It has been demonstrated that matched‐pair analysis can be effectively used to reduce residual variability and improve stability of estimates10,11 even in an RCT setting. Analyses were performed using R system12 and Hglm13 package.
thritis according to the clinical and radiological criteria of the American College of Rheumatology.9 The
pharmacological
treatment
and
the
placebo
were
2.4
|
WOMAC subscales' density function
manufactured and packaged so that they were indistinguishable.
Several studies have demonstrated that WOMAC subscales are
Patients
random
skewed distributed (asymmetrical distributed with a tail on the right
numbers. The treatment code of each patient was kept in a single
corresponding to poor responders).14,15 Assuming that the subscales
opaque sealed envelope to be opened only in case of emergency.
are normally distributed, when they are not (using improper modelling
were
randomized
using
computer‐generated
approach), may affect the power and increases the chance of Type I
2.1
|
error,16 by failing in the identification of subjects with extreme
Variables collected at baseline
values at the tail of the distribution (with consequently difficulties in
In both Randomized Clinical Trial (RCTs) patients' medical history was
characterizing poor responders in the case of WOMAC subscales).
assessed at baseline. Besides demographic characteristics (age and
After examining WOMAC subscales for normality using kernel
gender), patients' habits were also assessed, referring to: smoking habit
density estimation based on Epanechnikov kernel,17 we tested for
(classified as smokers vs nonsmokers), alcohol consumption (consumers
nonnormality using the Henze‐Zirkler's Multivariate Normality Test18
vs nonconsumers), caffeine consumption (consumers vs nonconsumers),
on the WOMAC distribution over the entire follow‐up.
and diet type (normal vs special, eg, hypocaloric diet).
Eventually, we propose a gamma‐based model for WOMAC
Referring to clinical status, patients' BMI, blood pressure, heart rate,
concomitant
infectious
diseases
(presence
vs
absence),
subscale, comparing it with the most widely used (in the context of osteoarthritis research) Gaussian approach.
LORENZONI
3
ET AL.
2.5 | Modelling WOMAC subscales' taking into account their potential correlation Since WOMAC subscales measure related quantities (pain, stiffness,
the latent variable. Such latent variable is assumed to have a normal distribution with mean 0 and variance σ2.
to be correlated to each other. However, generally, studies on
2.6 | Modelling WOMAC subscales taking into account HTE
osteoarthritis ignore the correlation between the outcomes and fit a
Given the highly heterogeneous characteristics of patients suffering
separate model to each response variable. Despite that this approach
from knee osteoarthritis affecting treatment performance on different
seems to be simple to apply, it completely not accounts for the
WOMAC subscales, we estimate a further model including an interac-
correlation between outcomes, determining a reduction of power in
tion term between the latent variable and treatment (representing the
detecting
misleading
unobserved covariates influencing treatment performance). This is to
estimates.19 To overcome this limitation, after testing WOMAC
account for unobserved heterogeneity affecting the likelihood of
subscales for correlation, we develop a multivariate model, which is
responding to treatment and identifying subsets of patients who are
an extension of Sammel latent variable model.20 Conditional on this
poor or high responders, sharing the same features.
and physical function) in the same individuals, probably outcomes tend
treatment
effectiveness
and
leading
to
latent variable, the outcomes are assumed to be independent and are
Summarizing, after testing the subscales for normality and
modelled as functions of fixed covariates and a subject‐specific latent
correlation, both gamma and Gaussian models (to compare each other)
variable. This latent variable establishes the link between the
were estimated through the following steps:
regression equations capturing correlations among the outcomes. Assuming that the latent variable completely specifies the correlation
• Gamma and Gaussian models were estimated ignoring the
among the outcomes, this permits examination of the outcomes as
correlation between the outcomes and fit a separate model to
independent of each other by accounting for the correlation through
each response variable (independent).
TABLE 1
Distribution of the 2 twin trials combined
Age Sex: M BMI Diatolic pressure Heart rate Systolic pressure
P value
Placebo
Treatment
Combined
(N = 207)
(N = 207)
(N = 414)
60.00/65.00/69.00
58.71/64.00/68.50
59.00/64.00/69.00
23% (47)
23% (48)
23% (95)
.907
25/27/28.625
25/27/28.075
25/27/28.340
.849
80.00/85.00/90.00
75.00/80.00/90.00
76.25/80.00/90.00
.611
64/70/80
64/72/80
64/72/80
.819 .277
.187
130/140/150
125/140/150
130/140/150
Smoking: Y
19% (39)
14% (30)
17% (69)
.235
Alcohol: Y
47% (98)
48% (100)
48% (198)
.844
Caffeine: Y
86% (179)
86% (177)
86% (356)
.777
Diet
17% (35)
17% (36)
17% (71)
.896
Infectious diseases: Y
26% (53)
16% (34)
21% (87)
.022
Concomitant medications: Y
62% (129)
65% (134)
64% (263)
.61
Presence of knee osteoarthritis (years)
4/8/12
4/7/10
4/ 7/11
.596
Increased temperature (left knee): Y
21% (43)
21% (44)
21% (87)
.904
Effusion (left knee): Y
26% (54)
22% (46)
24% (100)
.358
Effusion righ knee: Y
29% (59)
33% (69)
31% (128)
.288
Alignment (left knee): normal
66% (136)
74% (154)
70% (290)
.01
20% (42)
20% (42)
20% (84)
valgus varus
14% (29)
5% (11)
10% (40)
70% (145)
72% (150)
71% (295)
valgus
16% (34)
19% (40)
18% (74)
varus
14% (28)
8% (17)
11% (45)
Crepitus active (left knee): Y
77% (160)
78% (161)
78% (321)
.906
Crepitu active (right knee): Y
75% (155)
83% (172)
79% (327)
.04
Crepitus passive (left knee): Y
72% (150)
70% (145)
71% (295)
.587
Alignment (right knee): normal
.196
Crepitu passive (right knee): Y
68% (141)
75% (156)
72% (297)
.102
Degree of flexion (left knee)
70/90/135
68/95/140
70/90/140
.71
Degree of flexion (right knee)
70.00/95.00/135.00
70.00/90.00/140.00
70.00/90.00/138.75
.845
Data are I quartile/median/III quartile for continuous variables and percentage (absolute numbers) for categorical variables. P value refers to the difference between placebo and treatment.
4
LORENZONI
ET AL.
• Multigamma and Gaussian models were estimated including a
stiffness and physical function subscales, (higher scores represent
latent term (latent only) accounting for correlation between
worse symptoms, a reduction of WOMAC scores corresponds to an
outcomes.
improvement of symptoms) over the time period considered
• Multigamma and Gaussian models were estimated including,
(Figure 1). The analysis of such subscales revealed that they were
besides the latent term, the interaction between the latent term
skewed (Figure 2) and highly correlated between each other, with
and treatment (latent interaction) accounting for HTE.
correlations coefficients ranging from 0.71 up to 0.87. Table 2 shows the effect of treatment on the considered including
outcomes in the “independent”, “latent only,” and “latent interaction”
covariates resulting significant by inspecting AIC at a significance level
gamma and Gaussian models. Modelling the subscales using a gamma
of 0.10 and by clinical judgement. Variables considered as candidates
approach accounting for skewness of data, we found out different
for the predictor part of the model are listed in Table 1.
results compared with Gaussian models. For independent models, in
21,22
Multivariate models were fitted using H‐likelihood,
the Gaussian approach, the treatment resulted not significant on all WOMAC subscales while, using the gamma approach, it resulted
3
|
RESULTS
to be significant for pain (P value = .021) and physical function (P value = .024). The main difference regarded the latent
Subjects included in the analysis from the pooled twin trails were 414
variable interacting with treatment (accounting for unobserved
(Table 1). Scores of WOMAC subscales improved, especially for
heterogeneity), which is not significant for the Gaussian approach
FIGURE 1
Trends of Western Ontario and McMaster Universities index (WOMAC) subscales over time. Line is a loss fit for estimating trend
over time
FIGURE 2
Density estimation of Western Ontario and McMaster Universities index (WOMAC) subscales distribution. Henze‐Zirkler's multivariate normality test (from left to right): 69.51 (P value < .001), 68.86 (P value < .001), and 71.36 (P value < .001)
LORENZONI
5
ET AL.
TABLE 2 Gaussian and gamma models, by fitting (1) three independent equations, (2) a multivariate model with 1 latent variable, and (3) a latent and a latent‐by‐treatment interaction
Independent
Latent only
Latent interaction
Gaussian model
Some further insights can be given by examining the estimated latent variable, as emerging from the gamma model with latent treatment interaction. First, as shown in Figure 6, the distribution of the latent variable is shifted to the right for treated vs placebo patients, as perhaps expected by the goodness‐of‐fit analysis discussed above. Second, by examining correlation between latent and the other
Coeff.
P value Coeff.
P value Coeff.
P value
variables, it clearly emerges (Figure 7) that the 2 groups are characterized
Pain
−0.110 .104
−0.082
.181
−0.115 .0007
by different variables. Particularly, lifestyle variables (such as caffeine
Physical function
−0.082 .157
−0.101
.099
−0.128 .002
consumption and smoking habit) differently correlate with latent variable
Stiffness
−0.095 .157
−0.055
.363
−0.101 .015
in the 2 groups of patients (smoking resulted to be directly correlated
.102
with latent variable in treated subjects and inversely correlated in the
—
Latent Bayesian information criterion (BIC)
placebo group, conversely caffeine is inversely correlated with latent variable in the treated group and directly correlated in the placebo group).
Gamma model Coeff.
P value Coeff.
P value Coeff.
Pain
−0.157 .021
−0.173
.038
−0.173 .006
Physical function
−0.139 .024
−0.157
.061
−0.17
Stiffness
−0.111 .095
−0.127
.128
−0.123 .057
Latent BIC
—
In summary,
P value
.007
