--- title: "SRM_replication_code" max-width: 2000px format: html: toc: true toc-location: left toc-depth: 2 self-contained: true page-layout: article editor: visual code-tools: true code-overflow: wrap execute: error: false warning: false column: page-inset-right --- # Packages, Functions, and Data ## Packages ```{r} library (haven)library (broom)library (knitr)library (DT)library (psych)library (sjstats)library (pwr)library (equate)library (tidyverse)library (furrr)library (viridis)library (Cairo)library (ggrepel)library (kableExtra)``` ## Custom Functions ### Substantive Functions ```{r} # -- MOC Interpolation ----------------------------------- # See section "Interval Quantity Questions" for the intervals_df structure <- function (intervals_df, alpha){%>% mutate (moc = case_when (!= Inf ~ (lower_bound + upper_bound)/ 2 ,== Inf ~ lower_bound* (alpha/ (alpha-1 ))%>% select (moc) %>% rename_with (~ paste0 ("alpha" ,alpha))# -- OSE-RG --------------------------------------------- # Wrapper function around the equate() function from the equate package <- function (freqtab_x, freqtab_y){<- equate (freqtab_x, freqtab_y, type = "equipercentile" )$ concordance$ yx# -- Open-Ended Question Aggregation -------------------- # The function takes the vector of open-ended responses, selects values within a quantity interval, and calculates the mean quantity # Optionally, it truncates the open-ended responses first, be removing the top x values <- function (lower_bound, upper_bound, v, trim_upper = 0 ){<- v[! is.na (v)]<- sort (v) <- v[1 : floor (length (v)* (1 - trim_upper))]>= lower_bound & v <= upper_bound] %>% mean ()<- function (lower_bound, upper_bound, v, f){<- v[! is.na (v)]>= lower_bound & v <= upper_bound] %>% f ()# -- Assymetric Trimmed Mean ------------------------------------- # Only removes upper cases # in contrast to the trim argument of base::mean() <- function (x,trim_upper = 0.05 ) {<- x[! is.na (x)]<- sort (x) mean (x[1 : floor (length (x)* (1 - trim_upper))])``` ### Utility Functions ```{r} # Table styling shortcut <- function (df){%>% kable () %>% kable_styling (full_width = FALSE , position = "left" ) }# Extracts recoding information from the harmonized_equivalents_df <- function (condition, approach, rec_df = harmonized_equivalents_df){%>% filter (.data$ condition == .env$ condition) %>% pull ({{approach}})``` ## Load data ```{r} source ("SRM_replication_code_sub_population_analysis.R" )<- sub_pop_format %>% filter (sub_pop == "Full sample" )<- sub_pop_format %>% filter (sub_pop != "Full sample" )<- read_sav ("SRM_replication_data.sav" )# Dataframe only with valid answers to the quantity question: <- data %>% mutate (any_valid = coalesce (books_low_freq, books_med_freq, books_hi_freq, books_numeric)) %>% drop_na (any_valid)``` # Descriptives ## Sample size ```{r} <- data %>% nrow ()<- valid_df %>% nrow ()``` `r total_cases` cases in total. Of those, `r valid_cases` cases have valid answers to the quantity question.## Demography ### Sex ```{r} %>% group_by (sex = D01 %>% as_factor) %>% summarise (n = n ()) %>% ungroup () %>% mutate (percent = 100 * n/ sum (n),percent = round (percent, 2 )) %>% style_table ()``` ### Age ```{r} %>% mutate (age = 2020 - D02) %>% summarise (across (age, list (mean= mean, median= median, sd= sd,min= min,max= max))) %>% %>% style_table``` ### Education ```{r} %>% group_by (education %>% as_factor) %>% summarise (n = n ()) %>% ungroup () %>% mutate (percent = 100 * n/ sum (n),percent = round (percent, 2 )) %>% style_table ()``` # Preliminary analyses and instrument information ## Open ended quantity question percentiles ### Full sample result ```{r results='asis'} %>% select (quantiles_table) %>% pull () %>% 1 ]] %>% as_tibble () %>% style_table () %>% print``` ### Subsample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n " )%>% filter (sub_pop == cur_sub) %>% select (quantiles_table) %>% pull () %>% 1 ]] %>% as_tibble () %>% style_table () %>% printcat (" \n\n " )``` ## Interval quantity questions ### Full sample result ```{r} %>% select (histogram_question) %>% pull () %>% 1 ]] %>% ``` ### Subsample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n\n " )%>% filter (sub_pop == cur_sub) %>% select (histogram_question) %>% pull () %>% 1 ]] %>% cat (" \n\n\n " )``` ## Interval Quantity Question Boundaries ### Table `intervals_df` also severs as an input into the `moc_interpolation()` function. This is because the MOC interpolation requires information about the verbatim interval boundaries from the different quantity measurement instruments.```{r} <- tribble (~ condition, ~ lower_bound, ~ upper_bound,"low quantity" , 0 , 10 ,"low quantity" , 11 , 25 ,"low quantity" , 26 , 50 ,"low quantity" , 51 , Inf ,"medium quantity" , 0 , 25 ,"medium quantity" , 26 , 50 ,"medium quantity" , 51 , 100 ,"medium quantity" , 101 , Inf ,"high quantity" , 0 , 50 ,"high quantity" , 51 , 100 ,"high quantity" , 101 , 250 ,"high quantity" , 251 , Inf %>% mutate (condition = factor (condition, levels = c ("low quantity" , "medium quantity" , "high quantity" )) %>% fct_rev%>% style_table ()``` ### Plot ```{r} %>% group_by (condition) %>% mutate (row = row_number (),condition_num = condition %>% as.numeric,upper_bound_label = ifelse (upper_bound != Inf , upper_bound, NA ))%>% ggplot ()+ annotate (geom = "segment" ,x = 50 , xend = 50 , y = 0.65 , yend = 3.45 ,size = 14 , color = "#aaaaaa" , lineend = "round" )+ annotate (geom = "segment" ,x = 50 , xend = 50 , y = 0.65 , yend = 3.45 ,size = 12 , color = "#eeeeee" , lineend = "round" )+ geom_rect (aes (ymin= condition_num-0.3 , ymax= condition_num+0.3 , xmin= lower_bound, xmax = upper_bound,fill= row))+ geom_text (aes (0 , condition_num, label = condition), hjust = 1.2 , size = 4 )+ geom_text (aes (upper_bound_label, condition_num+0.4 , label = upper_bound_label), size = 4 )+ geom_text (aes (310 , condition_num+0.41 ), label = expression (phantom (0 )%->% ~~ infinity), size = 4 )+ scale_x_continuous (limits = c (- 100 , 310 ))+ scale_fill_viridis (option = "mako" , end = 0.7 , begin = 0.3 )+ theme_void ()+ theme (legend.position = "none" ,plot.background = element_rect (color = "white" , fill= "white" )``` ## Response bias demonstration ### Descriptive #### Full sample result ```{r} %>% mutate (low_50_or_less = ifelse (books_low_freq <= 3 , TRUE , FALSE ),med_50_or_less = ifelse (books_med_freq <= 2 , TRUE , FALSE ),hi_50_or_less = ifelse (books_hi_freq <= 1 , TRUE , FALSE )%>% select (contains ("_50_or_less" )) %>% pivot_longer (everything (), names_to = "condition" , values_to = "50_or_less" ) %>% drop_na (` 50_or_less ` ) %>% group_by (condition) %>% summarise (Percent_50_or_less = mean (` 50_or_less ` )* 100 ) %>% mutate (condition = factor (condition, levels = c ("low_50_or_less" , "med_50_or_less" , "hi_50_or_less" ),labels = c ("low quantity" ,"medium quantity" ,"high quantity" ),ordered = TRUE )%>% arrange (condition) %>% style_table ()``` #### Subsample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n " )%>% filter (sub_pop == cur_sub) %>% select (response_bias_table) %>% pull () %>% 1 ]] %>% as_tibble () %>% style_table () %>% cat (" \n\n " )``` ### Spearman Test #### Full sample result ```{r} <- valid_df %>% mutate (low_50_or_less = ifelse (books_low_freq <= 3 , TRUE , FALSE ),med_50_or_less = ifelse (books_low_freq <= 2 , TRUE , FALSE ),hi_50_or_less = ifelse (books_low_freq <= 1 , TRUE , FALSE )%>% select (contains ("_50_or_less" )) %>% pivot_longer (everything (), names_to = "condition" , values_to = "50_or_less" ) %>% drop_na (` 50_or_less ` ) %>% mutate (more_than_50 = as.numeric (! ` 50_or_less ` ))%>% mutate (condition = factor (condition, levels = c ("low_50_or_less" , "med_50_or_less" , "hi_50_or_less" ),ordered = TRUE ),condition_num = as.numeric (condition))cor.test (spearman_df$ condition_num, spearman_df$ more_than_50, method = "spearman" ) %>% tidy %>% mutate (across (where (is.numeric), round, 3 )) %>% t %>% style_table ()``` #### Subsample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n " )%>% filter (sub_pop == cur_sub) %>% select (response_bias_spearman) %>% pull () %>% 1 ]] %>% style_table () %>% cat (" \n\n " )``` ## Sensitivity of MOC to different Alpha values ```{r} %>% add_column (map_dfc (seq (1 ,3 ,length.out = 21 ), ~ moc_interpolation (intervals_df, .x))%>% pivot_longer (starts_with ("alpha" ), names_to = "alpha" , values_to = "moc_equivalents" ) %>% mutate (alpha = str_remove (alpha, "alpha" ) %>% as.numeric ()) %>% filter (upper_bound== Inf ,!= 1 ) %>% mutate (response_label = paste0 ("More than " , lower_bound)) %>% ggplot (aes (x= alpha, y= moc_equivalents))+ facet_grid (rows = vars (condition), scales = "free" )+ geom_point (aes (color= alpha), size = 4 )+ geom_text (aes (color= alpha, label = round (moc_equivalents, 0 )), angle = 45 , hjust = - 0.2 , size = 6 )+ geom_text (aes (x= Inf , y= Inf , label = response_label), stat = "unique" , hjust = 1.05 , vjust = 1.3 , size = 6 )+ scale_y_continuous ("MOC interp. of highest interval" , limits = c (0 , 3600 ))+ scale_x_continuous (breaks = c (1.1 , 1.5 , 2.0 , 2.5 , 3.0 ))+ theme_bw (base_size = 14 )+ theme (legend.position = "none" ,plot.background = element_rect (colour = "white" , fill= "white" )+ coord_cartesian (clip = FALSE )``` # Main Analyses ## All interval conditions harmonized towards the open ended format ```{r} # Preparing Frequency tables for OSE-RG <- freqtab (data$ books_low_freq)<- freqtab (data$ books_med_freq)<- freqtab (data$ books_hi_freq)<- freqtab (data$ books_numeric)# Extracting the vector of open responses <- data$ books_numeric# Calculate a table with equivalent values for every interval response option <- intervals_df %>% # Calculating MOC interpolated values add_column (moc_interpolation (intervals_df, alpha = 2 )%>% # Calculating OSE-RG Equivalents mutate (equip_equating = c (calculate_eq_equivalents (freqtab_low, freqtab_open),calculate_eq_equivalents (freqtab_med, freqtab_open),calculate_eq_equivalents (freqtab_hi, freqtab_open)%>% # Calculate average open-ended values within each interval mutate (# All open-ended responses open_mean_quantity = map2_dbl (lower_bound, upper_bound,~ aggregate_open_interval (.x, .y, open_vector)),# Trimmed by removing the top 5% open-ended responses open_mean_quantity_trimmed = map2_dbl (lower_bound, upper_bound,~ aggregate_open_interval (.x, .y, open_vector, trim_upper = 0.05 )),open_median_quantity = map2_dbl (lower_bound, upper_bound, ~ aggregate_open_interval_generalized (.x, .y, open_vector, median))%>% # Adding numerical interval IDs and interval labels group_by (condition) %>% mutate (interval_number = row_number (),interval_label = paste0 ("[" ,lower_bound, ", " ,upper_bound,"]" )``` ### Plot equivalent values Interval -\> Open Ended #### Full sample result ```{r fig.height= 8} %>% mutate (interval_number = (- 1 * (interval_number-1 ))+ 4 ,interval_label = paste0 (interval_number, " - " , interval_label),interval_label = factor (interval_label) %>% fct_rev) %>% pivot_longer (alpha2: open_median_quantity, values_to = "equivalents" , names_to = "approach" ) %>% mutate (condition = fct_rev (condition),approach = factor (approach,levels = c ("alpha2" ,"equip_equating" ,"open_mean_quantity" ,"open_mean_quantity_trimmed" ,"open_median_quantity" labels = c ("MOC" ,"OSE-RG" ,"Open Q. (untrimmed)" ,"Open Q. (trimmed)" ,"Open Q. (median)" %>% ggplot (aes (equivalents, interval_label, color= approach, shape = approach, linetype = approach))+ facet_wrap (vars (condition), ncol = 1 , scales = "free" )+ geom_line (aes (group = approach), size = 1 )+ geom_point (size = 4 , stroke = 2 )+ #geom_text(aes(0, label = interval_label, color = NA), hjust = 1, stat = "unique")+ scale_color_manual (values = c ("#497593" , "orange" , "#00ABA0" , "#00ABA0" , "#00ABA0" ))+ scale_linetype_manual (values = c (NA , NA , "solid" ,"dashed" , "dotted" ))+ scale_shape_manual (values = c (3 ,4 ,NA , NA , NA ))+ scale_x_continuous ("Reported books" )+ scale_y_discrete ("" )+ theme_bw (base_size = 20 )+ coord_flip ()+ #guides(color=guide_legend(ncol=2, byrow = TRUE))+ theme (panel.grid.minor = element_blank (),panel.grid.major.x = element_blank (),strip.background = element_blank (),#legend.position = "bottom", legend.title = element_blank ()``` #### Subsample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n\n " )%>% filter (sub_pop == cur_sub) %>% select (numerical_equivalents_plot) %>% pull () %>% 1 ]] %>% cat (" \n\n\n " )``` ### Values underlying the plot #### Full sample results ```{r} %>% select (interval = interval_label,MOC = alpha2,` OSE-RG ` = equip_equating,` Open mean (untrimmed) ` = open_mean_quantity,` Open mean (trimmed) ` = open_mean_quantity_trimmed,%>% mutate (across (where (is.numeric), round, 1 )) %>% ``` #### Subsample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n\n " )%>% filter (sub_pop == cur_sub) %>% select (numerical_equivalents_table) %>% pull () %>% 1 ]] %>% style_table () %>% cat (" \n\n\n " )``` ## Harmonization Performance ### Aggregated Mean bias #### Data preparation ```{r} # Here, the average quantities are calculated for untrimmed data <- data %>% select (` low quantity ` = books_low_freq,` medium quantity ` = books_med_freq,` high quantity ` = books_hi_freq,` open question ` = books_numeric) %>% # Recode response data with the equivalent values calculated earlier transmute (` open question ` = ` open question ` ,# Recode all conditions via MOC across (` low quantity ` : ` high quantity ` , ~ get_rec_vector (cur_column (), "alpha2" )[.x], .names = "{.col}_MOC" ),# Recode all conditions via OSE-RG across (` low quantity ` : ` high quantity ` , ~ get_rec_vector (cur_column (), "equip_equating" )[.x], .names = "{.col}_OSE-RG" ),%>% # calculate mean book quantity for every harmonization approach and condition summarise (across (everything (), mean, na.rm= TRUE )) %>% # Wide to long format pivot_longer (contains ("quantity" )) %>% # Calculate differences from the open-ended question expectation mutate (delta_mean = value - ` open question ` %>% # tidy up table separate (into = c ("condition" , "approach" ), sep = "_" %>% mutate (condition = factor (condition,levels = c ("low quantity" , "medium quantity" , "high quantity" )),trim = "untrimmed" )# Here, the same is done as above, but this time with trimmed data <- data %>% select (` low quantity ` = books_low_freq,` medium quantity ` = books_med_freq,` high quantity ` = books_hi_freq,` open question ` = books_numeric) %>% transmute (` open question ` = ` open question ` ,across (` low quantity ` : ` high quantity ` , ~ get_rec_vector (cur_column (), "alpha2" )[.x], .names = "{.col}_MOC" ),across (` low quantity ` : ` high quantity ` , ~ get_rec_vector (cur_column (), "equip_equating" )[.x], .names = "{.col}_OSE-RG" ),%>% summarise (across (everything (), mean_trim_upper, trim_upper = 0.05 )) %>% pivot_longer (contains ("quantity" )) %>% mutate (delta_mean = value - ` open question ` %>% separate (into = c ("condition" , "approach" ), sep = "_" %>% mutate (condition = factor (condition,levels = c ("low quantity" , "medium quantity" , "high quantity" )),trim = "trimmed 5%" )``` #### Plot ##### Full sample result ```{r} # Untrimmed and trimmed means are combined and then plotted %>% add_row (aggregated_book_estimations_trimmed_df) %>% mutate (delta_mean_percent = 100 * (delta_mean / ` open question ` ),trim = factor (trim,levels = c ("untrimmed" , "trimmed 5%" )),condition = fct_rev (condition),label_hjust = ifelse (approach == "MOC" ,- lead (delta_mean, 3 ) < 1 ,- lag (delta_mean, 3 ) < 1 )%>% ggplot (aes (delta_mean_percent, as.integer (condition), color = approach, group = approach))+ facet_grid (vars (trim))+ geom_vline (aes (xintercept= 0 ), color = "#666666" , size = 1 )+ geom_smooth (method = "lm" , se= FALSE , fullrange= TRUE , size = 2 )+ geom_smooth (method = "lm" , se= FALSE , fullrange= TRUE , color = "#ffffffaa" , size = 2 )+ geom_label (aes (label= paste0 (" " , round (delta_mean_percent), "% " ), hjust = label_hjust),size = 4 , alpha = 0.7 , label.size = 1 , show.legend = FALSE )+ geom_point (size = 3 )+ scale_x_continuous ("Average quantity difference to open question (%)" , limits = c (- 500 , 100 ))+ scale_y_continuous ("" , labels = aggregated_book_estimations_df$ condition %>% fct_rev %>% levels (), breaks = c (1 ,2 ,3 ))+ coord_cartesian (ylim = c (0.5 ,3.5 ),xlim = c (- 100 ,50 ))+ scale_color_manual (values = c ("#497593" , "orange" ))+ theme_bw (base_size = 12 )+ theme (panel.grid.minor = element_blank ()``` ##### Subsample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n\n " )%>% filter (sub_pop == cur_sub) %>% select (average_quantity_dif_plot) %>% pull () %>% 1 ]] %>% cat (" \n\n\n " )``` ## Interval to interval harmonization ```{r} # Calculating the equivalent scores <- harmonized_equivalents_df %>% select (lower_bound, upper_bound, interval_label, interval_number, alpha2) %>% ungroup () %>% mutate (equipercentile_towards_mid = c (calculate_eq_equivalents (freqtab_low, freqtab_med),calculate_eq_equivalents (freqtab_med, freqtab_med),calculate_eq_equivalents (freqtab_hi, freqtab_med)``` ### Equivalents table #### Full sample results ```{r} %>% filter (condition != "medium quantity" ) %>% transmute (condition, ` original score ` = interval_number,` OSE-RG towards medium quantity ` = equipercentile_towards_mid %>% round (2 )) %>% style_table ()``` #### Subsample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n\n " )%>% filter (sub_pop == cur_sub) %>% select (harmonization_to_medium_quantity_equivalents_table) %>% pull () %>% 1 ]] %>% style_table () %>% cat (" \n\n\n " )``` ### Equivalents plot #### Full sample results ```{r} %>% mutate (condition = factor (condition),cond_num = as.integer (condition),interval_label = ifelse (str_detect (interval_label, "Inf" ),paste0 (str_remove (interval_label, "Inf]" ), "\U221E)" ),%>% ggplot (aes (cond_num, equipercentile_towards_mid, color = interval_number))+ geom_segment (aes (x = cond_num-0.4 , xend = cond_num+0.4 ,yend = equipercentile_towards_mid), size = 2 , lineend = "round" )+ geom_label (aes (label = interval_label), label.padding = unit (0.5 , "lines" ), size = 6 , label.size = 1 )+ scale_x_continuous ("" , labels = interval_interval_equivalents$ condition %>% levels, breaks = c (1 ,2 ,3 ))+ scale_y_continuous ("OSE-RG equivalents in score format" )+ scale_color_viridis (option = "mako" , end = 0.7 , begin = 0.3 )+ theme_minimal (base_size = 12 )+ theme (legend.position = "none" ,plot.background = element_rect (color = "white" , fill = "white" )``` #### Subsample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n\n " )%>% filter (sub_pop == cur_sub) %>% select (harmonization_to_medium_quantity_plot) %>% pull () %>% 1 ]] %>% cat (" \n\n\n " )``` ## Interval to interval harmonization performance ```{r} # Calculating MOC equivalents <- intervals_df %>% add_column (moc_interpolation (intervals_df, 2 )%>% pull (alpha2)# Recoding data <- data %>% transmute (low_freq.raw = books_low_freq,med_freq.raw = books_med_freq,hi_freq.raw = books_hi_freq,# Apply OSE-RG to harmonize towards medium frequency condition low_freq.eq = calculate_eq_equivalents (freqtab_low, freqtab_med)[books_low_freq],med_freq.eq = calculate_eq_equivalents (freqtab_med, freqtab_med)[books_med_freq],hi_freq.eq = calculate_eq_equivalents (freqtab_hi, freqtab_med)[books_hi_freq],# Apply MOC to interpolate all conditions low_freq.moc = moc_low_med_hi[1 : 4 ][books_low_freq],med_freq.moc = moc_low_med_hi[5 : 8 ][books_med_freq],hi_freq.moc = moc_low_med_hi[9 : 12 ][books_hi_freq]# Caluclate <- interval_interval_data %>% summarise (across (everything (), list (mean= mean, sd = sd), na.rm= TRUE , .names = "{.col}.{.fn}" )) %>% pivot_longer (everything ()) %>% separate (name, into = c ("condition" , "approach" , "statistic" ), sep = " \\ ." )``` ### table #### Full sample results ```{r} %>% pivot_wider (names_from = statistic, values_from = value) %>% mutate (mid_mean = mean[c (2 ,2 ,2 ,5 ,5 ,5 ,8 ,8 ,8 )],mid_sd = sd[c (2 ,2 ,2 ,5 ,5 ,5 ,8 ,8 ,8 )],d = (mean- mid_mean)/ mid_sd,condition = factor (condition, levels = c ("low_freq" ,"med_freq" ,"hi_freq" ),labels = c ("low quantity" ,"medium quantity" ,"high quantity" %>% fct_rev,approach = factor (approach,levels = c ("eq" , "moc" ),labels = c ("OSE-RG" , "MOC" )) %>% fct_rev ()%>% drop_na (approach) %>% :: relocate (approach) %>% mutate (across (where (is.numeric), round, 2 )) %>% rename (` medium condition mean ` = mid_mean,` medium condition sd ` = mid_sd%>% kable () %>% kable_styling ()``` #### Sub sample results ```{r results="asis"} for (cur_sub in sub_pop_results$ sub_pop){cat ("## " ,cur_sub," \n\n\n " )%>% filter (sub_pop == cur_sub) %>% select (harmonization_to_medium_quantity_summary_table) %>% pull () %>% 1 ]] %>% style_table () %>% cat (" \n\n\n " )``` 
