About this exercise

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.



Introduction

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):

  • knowledge about the construct to be measured (Do we know what we want to measure?)
  • complexity of the construct (What dimensions / sub-aspects does the construct have; what content needs to be covered, etc.)
  • dependency on the situation (The same tool can be “valid” in one situation (young adults) and less-valid in another (e.g. older adults with mild cognitive impairment))
  • validation of scores, not measurement instruments (the conclusions need to be valid)
  • formulation of specific hypotheses (we need to establish hypotheses and test them)
  • validation as a continuous process (we always hear “the instrument is valid”. This is kind of an exageration. Rarely are all aspects of validity tested.)



Load the data

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.

data <- rio::import("https://ndownloader.figshare.com/files/22063455")

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.

data <- janitor::clean_names(data, case="none")

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



Calculation of socers for the questionnaires

Scores of the Kessler 6 questionnaire

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.

table_k <- data %>% 
  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



Sub-scores of the DASS21

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)

table_dass <- data %>% 
  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

EQ-5D utilities

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.

df.eq5d <- data %>% 
  mutate(MO = H1,
         SC = H2, 
         UA = H3, 
         PD = H4, 
         AD = H5) %>% 
  select(MO:AD)

data$EQ5D_UtilityIndex <- eq5d::eq5d(df.eq5d, version="5L", type="VT", country="Vietnam", ignore.invalid=TRUE)



Job Content Questionnaire (JCQ)

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)



Item-total-correlation and Cronbach’s alpha

Although internal consistency is not the learning objective of this exercise , we still calculate it to reproduce parts or the Table 3:

Figure 2. Table 3 of the paper 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.

Kessler_items<-data %>% 
  select(K1:K6)

psych::alpha(Kessler_items)
## 
## 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



Testing the Hypotheses

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:

Figure 4. Table of the Hypotheses 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.

Variables<-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)")
dass_depression<-psych::describe(data$DASS_Depression)
dass_anxiety<-psych::describe(data$DASS_Anxiety)
dass_stress<-psych::describe(data$DASS_Stress)
eq5dui<-psych::describe(data$EQ5D_UtilityIndex)
eq5d_sh<-psych::describe(data$H6)
jcs_demand<-psych::describe(data$psychological_demand)
jcs_control<-psych::describe(data$decision_authority)
jcs_supervisor_support<-psych::describe(data$supervisor_support)
jcs_coworker_support<-psych::describe(data$coworker_support)
edu<-psych::describe(data$education)
personal_income<-psych::describe(data$monthly_income)

r_DASS_Anxiety<-psych::cor.ci(data[c("DASS_Anxiety", "Score_K6")], plot=FALSE, method="pearson")
r_DASS_Depression<-psych::cor.ci(data[c("DASS_Depression", "Score_K6")], plot=FALSE, method="pearson")
r_DASS_Stress<-psych::cor.ci(data[c("DASS_Stress", "Score_K6")], plot=FALSE, method="pearson")
r_EQ5D<-psych::cor.ci(data[c("EQ5D_UtilityIndex", "Score_K6")], plot=FALSE, method="pearson")
r_self_rep_health<-psych::cor.ci(data[c("H6", "Score_K6")], plot=FALSE, method="pearson")
r_job_demand<-psych::cor.ci(data[c("psychological_demand", "Score_K6")], plot=FALSE, method="pearson")
r_job_control<-psych::cor.ci(data[c("decision_authority", "Score_K6")], plot=FALSE, method="pearson")
r_job_supervisor_support<-psych::cor.ci(data[c("supervisor_support", "Score_K6")], plot=FALSE, method="pearson")
r_job_coworker_support<-psych::cor.ci(data[c("coworker_support", "Score_K6")], plot=FALSE, method="pearson")
r_education<-psych::cor.ci(data[c("education", "Score_K6")], plot=FALSE, method="pearson")
r_income<-psych::cor.ci(data[c("monthly_income", "Score_K6")], plot=FALSE, method="pearson")

N<-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 )
Average<-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 )
SD<-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 )
rs<-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])
ci_lb<-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_ub<-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])

table.df<-data.frame(Variables, N, Average, SD, rs, ci_lb, ci_ub)

table.df<-table.df %>% 
  mutate(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.

Figure 5. Table 4 of the paper

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.