Margin of error comparison.

knitr::opts_chunk$set(
  warning = TRUE, # show warnings during codebook generation
  message = TRUE, # show messages during codebook generation
  error = TRUE, # do not interrupt codebook generation in case of errors,
                # usually better for debugging
  echo = TRUE,  # show R code
  fig.width = 4,
  fig.height = 4
)
ggplot2::theme_set(ggplot2::theme_bw())

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(codebook)
library(haven)
library(labelled)

## 
## Attaching package: 'labelled'
## 
## The following object is masked from 'package:codebook':
## 
##     to_factor

Margin of error.

Comparison to pilot studies of different sizes.

Item pairs

Ns = c(10, 15, 20, 30, 40, 50, 60, 80, 100, 200)
r_distribution = seq(from = -1, to = 1, by = 0.01)
simulation_results = list()

for (e in 1:length(Ns)) {
  N = Ns[e]
  for (i in seq_along(1:1)) {
    print(paste("N =", N, "i =", i))
    se = (1 - r_distribution^2)/sqrt(N - 3)
    index = (e - 1) + i
    simulation_results[[index]] = data.frame(N = N, r = r_distribution, se = se)
  }
}

## [1] "N = 10 i = 1"
## [1] "N = 15 i = 1"
## [1] "N = 20 i = 1"
## [1] "N = 30 i = 1"
## [1] "N = 40 i = 1"
## [1] "N = 50 i = 1"
## [1] "N = 60 i = 1"
## [1] "N = 80 i = 1"
## [1] "N = 100 i = 1"
## [1] "N = 200 i = 1"

length(simulation_results)

## [1] 10

sim_res <- bind_rows(simulation_results)
nrow(sim_res)

## [1] 2010

sim_res %>%
  ggplot(aes(x = r, y = se, color = N, group = N)) +
  geom_line() +
  geom_text(aes(label = N), data = sim_res %>% group_by(N) %>% filter(r == 0), nudge_y = 0.002) +
  theme_bw()

Validation

m_lmsynth <- readRDS("ignore/m_synth_rr_r_items_lm.rds")
library(brms)

## Loading required package: Rcpp

## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').

## 
## Attaching package: 'brms'

## The following object is masked from 'package:stats':
## 
##     ar

plot_prediction_error_items <- plot(conditional_effects(m_lmsynth, dpar = "sigma"), plot = F)[[1]] + 
  theme_bw() + 
  xlab("Synthetic inter-item correlation") + 
  ylab("Prediction error (sigma) // SE") +
  geom_smooth(stat = "identity", color = "#a48500", fill = "#EDC951") +
  geom_line(aes(x = r, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res) +
  geom_text(aes(label = paste("N =", N), x = r, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res %>% group_by(N) %>% filter(r == 0), nudge_y = 0.01) +
  scale_color_gradient(guide = "none")

## Loading required package: rstan

## Loading required package: StanHeaders

## 
## rstan version 2.32.6 (Stan version 2.32.2)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)

## 
## Attaching package: 'rstan'

## The following object is masked from 'package:tidyr':
## 
##     extract

plot_prediction_error_items

ggsave("ignore/plot_prediction_error_items.png", width = 6, height = 6)

Pilot

m_lmsynth_items_pilot <- readRDS("ignore/m_synth_r_items_lm.rds")

plot_prediction_error_items_pilot <- plot(conditional_effects(m_lmsynth_items_pilot, dpar = "sigma"), plot = F)[[1]] + 
  theme_bw() + 
  xlab("Synthetic inter-item correlation") + 
  ylab("Prediction error (sigma) // SE") +
  geom_line(aes(x = r, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res) +
  geom_smooth(stat = "identity", color = "#a48500", fill = "#EDC951") +
  geom_text(aes(label = paste("N =", N), x = r, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res %>% group_by(N) %>% filter(r == 0), nudge_y = 0.01) +
  scale_color_gradient(guide = "none")

plot_prediction_error_items_pilot

ggsave("ignore/plot_prediction_error_items_pilot.png", width = 6, height = 6)

Scale pairs

Validation

m_lmsynth_scales <- readRDS("ignore/m_synth_rr_r_scales_lm8.rds")

plot_prediction_error_scales <- plot(conditional_effects(m_lmsynth_scales, dpar = "sigma"), plot = F)[[1]] + 
  theme_bw() + 
  xlab("Synthetic inter-scale correlation") + 
  ylab("Prediction error (sigma) // SE") +
  geom_smooth(stat = "identity", color = "#a48500", fill = "#EDC951") +
  geom_line(aes(x = r, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res) +
  geom_text(aes(label = paste("N =", N), x = r, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res %>% group_by(N) %>% filter(r == 0), nudge_y = 0.01) +
  scale_color_gradient(guide = "none")

plot_prediction_error_scales

ggsave("ignore/plot_prediction_error_scales.png", width = 6, height = 6)

Pilot

m_lmsynth_scales_pilot <- readRDS("ignore/m_synth_r_scales_lm8.rds")

plot_prediction_error_scales_pilot <- plot(conditional_effects(m_lmsynth_scales_pilot, dpar = "sigma"), plot = F)[[1]] + 
  theme_bw() + 
  xlab("Synthetic inter-scale correlation") + 
  ylab("Prediction error (sigma) // SE") +
  geom_line(aes(x = r, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res) +
  geom_smooth(stat = "identity", color = "#a48500", fill = "#EDC951") +
  geom_text(aes(label = paste("N =", N), x = r, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res %>% group_by(N) %>% filter(r == 0), nudge_y = 0.01) +
  scale_color_gradient(guide = "none")

plot_prediction_error_scales_pilot

ggsave("ignore/plot_prediction_error_scales_pilot.png", width = 6, height = 6)

Reliabilities

Validation

Ns = c(15, 20, 30, 40, 50, 60, 80, 100, 200)
ks = c(3, 5, 10)
a_distribution = seq(from = -1, to = 1, by = 0.01)
simulation_results = list()
rowdiff = function(x) {
  lcl = x$ALPHA - x$LCL
  ucl = x$UCL - x$ALPHA
  (lcl + ucl)/2
}

for (e in seq_along(Ns)) {
  N = Ns[e]
  for (f in seq_along(ks)) {
    k = ks[f]
    print(paste("N =", N, "k =", k))
    se = rowdiff(psychometric::alpha.CI(a_distribution, k = k, N = N, level = 0.95))/1.96
    simulation_results[[length(simulation_results) + 1]] = data.frame(N = N, k = k, a = a_distribution, se = se)
  }
}

## [1] "N = 15 k = 3"
## [1] "N = 15 k = 5"
## [1] "N = 15 k = 10"
## [1] "N = 20 k = 3"
## [1] "N = 20 k = 5"
## [1] "N = 20 k = 10"
## [1] "N = 30 k = 3"
## [1] "N = 30 k = 5"
## [1] "N = 30 k = 10"
## [1] "N = 40 k = 3"
## [1] "N = 40 k = 5"
## [1] "N = 40 k = 10"
## [1] "N = 50 k = 3"
## [1] "N = 50 k = 5"
## [1] "N = 50 k = 10"
## [1] "N = 60 k = 3"
## [1] "N = 60 k = 5"
## [1] "N = 60 k = 10"
## [1] "N = 80 k = 3"
## [1] "N = 80 k = 5"
## [1] "N = 80 k = 10"
## [1] "N = 100 k = 3"
## [1] "N = 100 k = 5"
## [1] "N = 100 k = 10"
## [1] "N = 200 k = 3"
## [1] "N = 200 k = 5"
## [1] "N = 200 k = 10"

length(simulation_results)

## [1] 27

sim_res <- bind_rows(simulation_results)
nrow(sim_res)

## [1] 5427

sim_res %>%
  ggplot(aes(x = a, y = se, color = N, linetype = factor(k), group = interaction(N, k))) +
  geom_line() +
  geom_text(aes(label = paste("N =", N, "k =", k), x = a, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res %>% group_by(N, k) %>% filter(a == -.99), nudge_y = 0.002) +
  theme_bw()

sim_res %>%
  ggplot(aes(x = a, y = se, color = N, linetype = factor(k), group = interaction(N, k))) +
  geom_line() +
  geom_text(aes(label = paste("N =", N, "k =", k), x = a, y = se, color = N, group = N, ymin = NULL, ymax = NULL), data = sim_res %>% group_by(N, k) %>% filter(a == -.99), nudge_y = 0.002) +
  theme_bw()

Validation

library(brms)
m_synth_rr_rel_lm <- readRDS("ignore/m_synth_rr_rel_lm.rds")

plot_prediction_error_rels <- plot(conditional_effects(m_synth_rr_rel_lm, dpar = "sigma"), plot = F)[[1]] + 
  theme_bw() + 
  xlab("Synthetic Cronbach's alpha") + 
  ylab("Prediction error (sigma) // SE") +
  geom_line(aes(x = a, y = se, color = N, linetype = factor(k), group = interaction(N, k), ymin = NULL, ymax = NULL), data = sim_res) +
  geom_smooth(stat = "identity", color = "#a48500", fill = "#EDC951") +
  geom_text(aes(label = paste("N =", N, "k =", k), x = a, y = se, ymin = NULL, ymax = NULL), data = sim_res %>% group_by(N, k) %>% filter(a == -.99), nudge_y = 0.002) +
  scale_color_gradient("Sample size N", guide = "legend") +
  scale_linetype_manual("Number of items k", values = c("solid", "dashed", "dotted", "twodash"),
  guide = "legend") +
  # place legend in top right corner
  theme(legend.position = c(0.75, 0.85)) +
  coord_cartesian(xlim = c(-1, 1), ylim = c(0, 0.5))

## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

plot_prediction_error_rels

ggsave("ignore/plot_prediction_error_rels.png", width = 6, height = 6)

Pilot

library(brms)
m_synth_rel_lm <- readRDS("ignore/m_synth_rel_lm.rds")

plot_prediction_error_rels_pilot <- plot(conditional_effects(m_synth_rel_lm, dpar = "sigma"), plot = F)[[1]] + 
  theme_bw() + 
  xlab("Synthetic Cronbach's alpha") + 
  ylab("Prediction error (sigma) // SE") +
  geom_line(aes(x = a, y = se, color = N, linetype = factor(k), group = interaction(N, k), ymin = NULL, ymax = NULL), data = sim_res) +
  geom_smooth(stat = "identity", color = "#a48500", fill = "#EDC951") +
  geom_text(aes(label = paste("N =", N, "k =", k), x = a, y = se, ymin = NULL, ymax = NULL), data = sim_res %>% group_by(N, k) %>% filter(a == -.99), nudge_y = 0.002) +
  scale_color_gradient("Sample size N", guide = "legend") +
  scale_linetype_manual("Number of items k", values = c("solid", "dashed", "dotted", "twodash"),
  guide = "legend") +
  # place legend in top right corner
  theme(legend.position = c(0.75, 0.85)) +
  coord_cartesian(xlim = c(-1, 1), ylim = c(0, 0.5))

plot_prediction_error_rels_pilot

ggsave("ignore/plot_prediction_error_rels_pilot.png", width = 6, height = 6)

Combined plot

library(patchwork)

(plot_prediction_error_items_pilot + ggtitle("Pilot study") +
    plot_prediction_error_rels_pilot + 
        scale_color_gradient("Sample size N", guide = "none") +
    plot_prediction_error_scales_pilot) /


(plot_prediction_error_items + ggtitle("Validation study") +
    plot_prediction_error_rels + scale_color_gradient("Sample size N", guide = "none") +
           scale_linetype_manual("Number of items k", values = c("solid", "dashed", "dotted", "twodash"),
  guide = "none") +
    plot_prediction_error_scales)

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for linetype is already present.
## Adding another scale for linetype, which will replace the existing scale.

ggsave("Figure_pred_error_vs_N.pdf", width = 8.3, height = 5, scale = 1.6, device = grDevices::cairo_pdf)
ggsave("Figure_pred_error_vs_N.png", width = 8.3, height = 5, scale = 1.6)