This is an exercise on construct validity and hypotheses testing. In the course of this exercise, you will see a lot of new R-code. We suspect that portions of this code are not fully traceable at this time. Nevertheless, we show the code so that you can use it in due time (e.g. for the master thesis). So the main goal of this exercise is not that you can do all the calculations yourself in R, but that you can better understand the procedure to evaluate construct validity. You can reproduce the calculations by simply copy and paste the R-codes provided. Make sure, that the following packages are installed and loaded:
library(rio)
library(tidyverse)
library(psych)
library(eq5d)
library(kableExtra)
For this exercise, we follow the paper from Kawakami and colleagues, entiteled Internal consistency reliability, construct validity, and item response characteristics of the Kessler 6 scale among hospital nurses in Vietnam. The Kessler Psychological Distress Scale (K6+) is a 6-item self-report measure of psychological distress intended to be used as a quick tool to assess risk for serious mental illness in the general population. For more information, see here.
The COSMIN Panel defines validity as “the degree to which an instrument truly measures the construct(s) it puports to measure” Mokkink 2010. However, we need to be aware, that not the instrument itself is valid, but that its use and the conclusion made based on its results can be more or less valid. A blood pressure cuff might give perfectly valid measures if used on a normaly sized person, but might be biased if used on a person with a larger arm circumference.
As you have or will read in the COSMIN book Measurement in Medicine, we need to consider the following points (p. 150):
Use our now well known package rio
to load the data
directly from the internet. The link is “https://ndownloader.figshare.com/files/22063455”. Store
the data in a data frame called data
.
You might want to have a look at the data-web-page.
<- rio::import("https://ndownloader.figshare.com/files/22063455") data
As it is often the case, we first need to do a bit of data cleaning.
When you look at the data, you will notice that there are some variables
(columns) with problematic names (income and education have spaces and
special characters in their name). We can replace all spaces and special
characters by underlines using the clean_names
function
from the janitor
package.
<- janitor::clean_names(data, case="none") data
As the names for these variables are still not practical, we change
them. We also create factors for the variables “income” and “education”.
Here, we use the tidyverse
package to do this. For more
information on the code, see the
cheat sheet.
<- data %>%
data mutate(monthly_income.factor = factor(Monthly_income_5_million_220_USD_1_5_10_million_220_432_USD_2_10_million_432_USD_3, levels=c(1,2,3), labels=c("≤ 5 million (220 USD)", "5- 10 million (220-432 USD)", "≥ 10 million (432 USD)"))) %>%
mutate(monthly_income = Monthly_income_5_million_220_USD_1_5_10_million_220_432_USD_2_10_million_432_USD_3) %>%
mutate(education.factor = factor(Education_Vocation_school_1_Colleges_2_University_undergraduate_3_Postgraduate_4, levels=c(1,2,3,4), labels=c("Vocation School", "Colleges", "University undergraduate", "Postgraduate"))) %>%
mutate(education = Education_Vocation_school_1_Colleges_2_University_undergraduate_3_Postgraduate_4)
Lastly, we exclude patients who had missing values (as they did in the paper).
<- data %>%
data filter(!is.na(K1) & !is.na(K2) & !is.na(K3) & !is.na(K4) & !is.na(K5) & !is.na(K6))
The next step is to calculate the total score of the Kessler 6
questionnaire. This is simply the sum of the variables k1
to k6
.
<- data %>%
data mutate(Score_K6 = select(., K1:K6) %>%
rowSums(na.rm = TRUE))
We can just check what we did by printing only the variables we used to calculate the score and the score. This makes it easier to check for errors.
<- data %>%
table_k select(contains("k"))
slice(data.frame(table_k), 1:10) %>%
kable() %>%
kable_styling(full_width = TRUE)
K1 | K2 | K3 | K4 | K5 | K6 | Score_K6 |
---|---|---|---|---|---|---|
2 | 1 | 1 | 2 | 1 | 1 | 8 |
1 | 0 | 0 | 1 | 1 | 0 | 3 |
2 | 1 | 1 | 1 | 1 | 1 | 7 |
1 | 0 | 0 | 1 | 1 | 0 | 3 |
1 | 0 | 0 | 0 | 0 | 0 | 1 |
2 | 1 | 0 | 1 | 0 | 0 | 4 |
1 | 0 | 0 | 0 | 0 | 0 | 1 |
1 | 0 | 0 | 1 | 0 | 0 | 2 |
2 | 2 | 2 | 1 | 1 | 1 | 9 |
2 | 1 | 2 | 2 | 1 | 0 | 8 |
You can find more information on the Depression Anxiety and Stress Scale 21 here.
Again, the score of a sub-scale is the sum of all variables belonging
to this sub-scale. For example, the score of the sub-scale anxiety is
the sum of the variables DASSA1
to
DASSA7
.
As the columns in the data frame are so nicely named, we can use the
c_across()
function to calculate the sums for each
sub-scale.
<- data %>%
data rowwise() %>%
mutate(DASS_Anxiety=sum(c_across(contains("DASSA")))) %>%
mutate(DASS_Depression=sum(c_across(contains("DASSD")))) %>%
mutate(DASS_Stress=sum(c_across(contains("DASSS"))))
Let’s again check what we did (because normally the first few attempts might contain errors in our code that might produce wrong results. Therefore, it is always necessary to check)
<- data %>%
table_dass select(DASS_Stress, contains("DASSS"))
slice(data.frame(table_dass), 1:10) %>%
kable() %>%
kable_styling(full_width = TRUE)
DASS_Stress | DASSS1 | DASSS2 | DASSS3 | DASSS4 | DASSS5 | DASSS6 | DASSS7 |
---|---|---|---|---|---|---|---|
6 | 1 | 0 | 0 | 0 | 1 | 1 | 3 |
4 | 1 | 1 | 0 | 0 | 0 | 1 | 1 |
8 | 1 | 1 | 1 | 1 | 1 | 1 | 2 |
7 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
2 | 1 | 1 | 0 | 0 | 0 | 0 | 0 |
1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 |
7 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
7 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
4 | 1 | 1 | 1 | 0 | 0 | 0 | 1 |
8 | 1 | 1 | 2 | 1 | 1 | 1 | 1 |
The EQ-5D-5L estimates health-related quality of life (HR-QOL) See here with five items covering
mobility, self care, usual activities, pain/discomfort and
anxiety/depression. Each item has five response options from no problems
(1) to extreme problems (5). The item H6 is the self-reported health
status. The authors used the standard value set for Vietnam to calculate
participants’ HR-QOL scores. The calculation of this score (also called
utiliy index) is a bit more complicated, therefore we use an additional
package for this called eq5d.
We also need to create a specific data frame for the function
eq5d
.
<- data %>%
df.eq5d mutate(MO = H1,
SC = H2,
UA = H3,
PD = H4,
AD = H5) %>%
select(MO:AD)
$EQ5D_UtilityIndex <- eq5d::eq5d(df.eq5d, version="5L", type="VT", country="Vietnam", ignore.invalid=TRUE) data
The authors used the JCQ to assess the psychosocial work environment with four subscales: five-item psychological demand scale, nine-item decision latitude scale, four-item supervisor support scale, four-item co-worker support scale. Items are likert-type items with four response options from 1 strongly disagree to 4 strongly agree.
Information on the sub-scale scoring was derived from here, however the items are differently named.
<- data %>%
data mutate(psychological_demand=(JCQ19+JCQ20)*3+(15-(JCQ22+JCQ23+JCQ26))*2,
skill_discretion=(JCQ3+JCQ5+JCQ7+JCQ9+JCQ11+(5-JCQ4))*2,
decision_authority=(JCQ6+JCQ10+(5-JCQ8))*4,
decision_latitute=skill_discretion+decision_authority,
supervisor_support=JCQ48+JCQ49+JCQ51+JCQ52,
coworker_support=JCQ53+JCQ54+JCQ56+JCQ58)
Although internal consistency is not the learning objective of this exercise , we still calculate it to reproduce parts or the Table 3:
We will calculate the item-total
correlations and the Cronbach’s alpha for the Kessler 6 items. For the
item-total correlation, you need to drop the item for which you
calculate the item-total-correlation. If you do not, the value is biased
because the item is on both sides of the correlation, once as item, and
once in the score. Therefore you need to take the r.drop value (see here for more
information on this issue). We do this by using the
alpha
function from the psych
package.
<-data %>%
Kessler_itemsselect(K1:K6)
::alpha(Kessler_items) psych
##
## Reliability analysis
## Call: psych::alpha(x = Kessler_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.86 0.87 0.85 0.52 6.5 0.0068 0.76 0.58 0.53
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.85 0.86 0.88
## Duhachek 0.85 0.86 0.88
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## K1 0.86 0.87 0.84 0.56 6.5 0.0069 0.0023 0.57
## K2 0.83 0.83 0.82 0.50 5.1 0.0087 0.0104 0.50
## K3 0.84 0.84 0.82 0.52 5.3 0.0083 0.0107 0.53
## K4 0.83 0.84 0.82 0.51 5.1 0.0086 0.0087 0.53
## K5 0.83 0.83 0.81 0.50 5.0 0.0086 0.0072 0.52
## K6 0.85 0.85 0.83 0.53 5.7 0.0077 0.0045 0.55
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## K1 948 0.70 0.68 0.58 0.54 1.66 0.83
## K2 948 0.82 0.81 0.77 0.72 0.61 0.78
## K3 948 0.79 0.78 0.72 0.68 0.68 0.75
## K4 948 0.81 0.81 0.76 0.70 0.77 0.81
## K5 948 0.81 0.82 0.79 0.72 0.53 0.71
## K6 948 0.72 0.74 0.68 0.62 0.31 0.63
##
## Non missing response frequency for each item
## 0 1 2 3 4 miss
## K1 0.07 0.34 0.44 0.14 0.01 0
## K2 0.55 0.31 0.12 0.02 0.00 0
## K3 0.48 0.37 0.14 0.01 0.00 0
## K4 0.44 0.38 0.15 0.03 0.00 0
## K5 0.59 0.30 0.11 0.01 0.00 0
## K6 0.76 0.18 0.05 0.01 0.00 0
# item-total correlation:
#item_analyis$item.stats[5]
# Cronbach's alpha total score
#item_analyis$total
# Cronbach's alpha per item
# item_analyis$alpha.drop
What hypotheses to test? Here is what they wrote in the method section about their eleven hypotheses:
Figure 3, Method section on Hypotheses testing
In order not to loose the overview, we write these hypotheses in a table:
We will now calculate all the
values needed for the Table 4, with the exception of the p-value. The
p-value is not really helpful in this context, as he tests the
hypothesis whether the correlation is zero. We want to see the 95%
confidence intervals and if they include the postulated values in the
table. So let’s calculate the correlation coefficients and their
corresponding 95% CI’s.
Note: The following code is primarily long enough to make a nice table. Calculating the correlations alone is much easier.
<-c("DASS depression score (0–21)", "DASS anxiety score (0–21)", "DASS stress score (0–21)", "HR-QOL score (0–1)", "Self-rated health (0–100)", "JCQ job demands score (12–48)", "JCQ job control score (24–96)", "JCQ supervisor support score (4–16)", "JCQ coworker support score (4–16)", "Education (1–4)", "Personal income (1–3)")
Variables<-psych::describe(data$DASS_Depression)
dass_depression<-psych::describe(data$DASS_Anxiety)
dass_anxiety<-psych::describe(data$DASS_Stress)
dass_stress<-psych::describe(data$EQ5D_UtilityIndex)
eq5dui<-psych::describe(data$H6)
eq5d_sh<-psych::describe(data$psychological_demand)
jcs_demand<-psych::describe(data$decision_authority)
jcs_control<-psych::describe(data$supervisor_support)
jcs_supervisor_support<-psych::describe(data$coworker_support)
jcs_coworker_support<-psych::describe(data$education)
edu<-psych::describe(data$monthly_income)
personal_income
<-psych::cor.ci(data[c("DASS_Anxiety", "Score_K6")], plot=FALSE, method="pearson")
r_DASS_Anxiety<-psych::cor.ci(data[c("DASS_Depression", "Score_K6")], plot=FALSE, method="pearson")
r_DASS_Depression<-psych::cor.ci(data[c("DASS_Stress", "Score_K6")], plot=FALSE, method="pearson")
r_DASS_Stress<-psych::cor.ci(data[c("EQ5D_UtilityIndex", "Score_K6")], plot=FALSE, method="pearson")
r_EQ5D<-psych::cor.ci(data[c("H6", "Score_K6")], plot=FALSE, method="pearson")
r_self_rep_health<-psych::cor.ci(data[c("psychological_demand", "Score_K6")], plot=FALSE, method="pearson")
r_job_demand<-psych::cor.ci(data[c("decision_authority", "Score_K6")], plot=FALSE, method="pearson")
r_job_control<-psych::cor.ci(data[c("supervisor_support", "Score_K6")], plot=FALSE, method="pearson")
r_job_supervisor_support<-psych::cor.ci(data[c("coworker_support", "Score_K6")], plot=FALSE, method="pearson")
r_job_coworker_support<-psych::cor.ci(data[c("education", "Score_K6")], plot=FALSE, method="pearson")
r_education<-psych::cor.ci(data[c("monthly_income", "Score_K6")], plot=FALSE, method="pearson")
r_income
<-c(dass_depression$n,dass_anxiety$n,dass_stress$n, eq5dui$n,eq5d_sh$n, jcs_demand$n, jcs_control$n, jcs_supervisor_support$n,jcs_coworker_support$n, edu$n, personal_income$n )
N<-c(dass_depression$mean,dass_anxiety$mean,dass_stress$mean, eq5dui$mean,eq5d_sh$mean, jcs_demand$mean, jcs_control$mean, jcs_supervisor_support$mean,jcs_coworker_support$mean, edu$mean, personal_income$mean )
Average<-c(dass_depression$sd,dass_anxiety$sd,dass_stress$sd, eq5dui$sd,eq5d_sh$sd, jcs_demand$sd, jcs_control$sd, jcs_supervisor_support$sd,jcs_coworker_support$sd, edu$sd, personal_income$sd )
SD<-c(r_DASS_Depression$rho[1,2],r_DASS_Anxiety$rho[1,2], r_DASS_Stress$rho[1,2],r_EQ5D$rho[1,2],r_self_rep_health$rho[1,2],r_job_demand$rho[1,2], r_job_control$rho[1,2],r_job_supervisor_support$rho[1,2],r_job_coworker_support$rho[1,2],r_education$rho[1,2] ,r_income$rho[1,2])
rs<-c(r_DASS_Depression$ci[1,1],r_DASS_Anxiety$ci[1,1], r_DASS_Stress$ci[1,1],r_EQ5D$ci[1,1],r_self_rep_health$ci[1,1],r_job_demand$ci[1,1], r_job_control$ci[1,1],r_job_supervisor_support$ci[1,1],r_job_coworker_support$ci[1,1],r_education$ci[1,1] ,r_income$ci[1,1])
ci_lb
<-c(r_DASS_Depression$ci[1,3],r_DASS_Anxiety$ci[1,3], r_DASS_Stress$ci[1,3],r_EQ5D$ci[1,3],r_self_rep_health$ci[1,3],r_job_demand$ci[1,3], r_job_control$ci[1,3],r_job_supervisor_support$ci[1,3],r_job_coworker_support$ci[1,3],r_education$ci[1,3] ,r_income$ci[1,3])
ci_ub
<-data.frame(Variables, N, Average, SD, rs, ci_lb, ci_ub)
table.df
<-table.df %>%
table.dfmutate(across(where(is.numeric), round, 2)) # if you need help for this line: https://stackoverflow.com/questions/43314328/how-to-use-dplyrmutate-all-for-rounding-selected-columns/44686406
%>%
table.df kable() %>%
kable_styling(full_width = TRUE)
Variables | N | Average | SD | rs | ci_lb | ci_ub |
---|---|---|---|---|---|---|
DASS depression score (0–21) | 932 | 3.04 | 2.95 | 0.51 | 0.46 | 0.56 |
DASS anxiety score (0–21) | 935 | 3.85 | 3.11 | 0.54 | 0.48 | 0.59 |
DASS stress score (0–21) | 935 | 5.57 | 3.56 | 0.54 | 0.49 | 0.60 |
HR-QOL score (0–1) | 948 | 0.94 | 0.08 | -0.40 | -0.46 | -0.34 |
Self-rated health (0–100) | 946 | 85.67 | 11.60 | -0.36 | -0.41 | -0.30 |
JCQ job demands score (12–48) | 933 | 31.49 | 4.43 | 0.28 | 0.22 | 0.34 |
JCQ job control score (24–96) | 940 | 33.51 | 4.11 | -0.14 | -0.22 | -0.07 |
JCQ supervisor support score (4–16) | 943 | 11.99 | 1.87 | -0.20 | -0.26 | -0.15 |
JCQ coworker support score (4–16) | 946 | 12.20 | 1.49 | -0.16 | -0.22 | -0.10 |
Education (1–4) | 944 | 1.94 | 0.95 | -0.08 | -0.14 | -0.01 |
Personal income (1–3) | 938 | 2.21 | 0.60 | -0.06 | -0.11 | 0.00 |
Let’s compare the Table from the paper and our Table.
Again, The p-value is not helpful as he tests the hypothesis whether the correlation is zero.
Let’s look in the paper what they wrote about these results:
Figure 5, Results on eleven hypotheses
I am not sure why the see all their hypotheses accepted. And we would have chosen different correlation coefficients for the education and the income, as these two variables were coded as ordinal variables. So either Spearman or polychoric correlations.