Precision simulations

The following precision simulations all follow the same structure.

• We randomly redraw from our holdout data to get realistic distributions of empirical estimates. SEs are estimated based on our planned N of 400.
• We repeat this random drawing process many times.
• We adjust for sampling error/the standard error of the empirical estimate to get the semi-latent accuracy.
• The target quantity is the average standard error of the semi-latent accuracy with which the synthetic estimates predicts the empirical estimates. We also report the maximal standard error across simulations.

We begin with item correlations, then reliabilities, then scale correlations.

knitr::opts_chunk$set(echo = TRUE, error = T, warning = F, message = F)

# Libraries and Settings

# Libs ---------------------------
library(tidyverse)
library(arrow)
library(glue)
library(psych)
library(lavaan)
library(ggplot2)
library(plotly)
library(gridExtra)
library(semTools)
library(semPlot)

model_name = "ItemSimilarityTraining-20240502-trial12"
#model_name = "item-similarity-20231018-122504"
pretrained_model_name = "all-mpnet-base-v2"

data_path = glue("./")
pretrained_data_path = glue("./")

set.seed(42)

number_of_items <- 246
number_of_scales <- 76
number_of_scales_with_more_than_3_items <- 56
combinations_items <- choose(number_of_items, 2)
combinations_scales <- choose(number_of_scales_with_more_than_3_items, 2)
scale_subscale_pairs <- 38
combinations_scales <- combinations_scales - scale_subscale_pairs # after eliminating scale-subscale pairs
planned_N <- 400

Precision simulation for synthetic inter-item correlations

holdout <- arrow::read_feather(file = file.path(data_path, glue("data/intermediate/{model_name}.raw.osf-bainbridge-2021-s2-0.item_correlations.feather")))

holdout_mapping_data = arrow::read_feather(
  file = file.path(data_path, glue("{model_name}.raw.osf-bainbridge-2021-s2-0.mapping2.feather"))
) %>%
  rename(scale_0 = scale0,
         scale_1 = scale1)

scales <- arrow::read_feather(file.path(data_path, glue("{model_name}.raw.osf-bainbridge-2021-s2-0.scales.feather"))
)

holdout_llm <- holdout %>%
  left_join(holdout_mapping_data %>% select(variable_1 = variable, InstrumentA = instrument, ScaleA = scale_0, SubscaleA = scale_1)) %>%
  left_join(holdout_mapping_data %>% select(variable_2 = variable, InstrumentA = instrument, ScaleA = scale_0, SubscaleA = scale_1))

sim_results <- tibble()
library(lavaan)

for(i in 1:500) {
  items <- holdout %>% select(variable_1) %>% distinct() %>% sample_n(number_of_items) %>% pull(variable_1)

  subset <- holdout %>% filter(variable_1 %in% items, variable_2 %in% items)

  N <- planned_N
  subset <- subset %>% mutate(se = (1 - empirical_r^2)/sqrt(N - 3))
  se2 <- mean(subset$se^2)

  r <- broom::tidy(cor.test(subset$empirical_r, subset$synthetic_r))

  model <- paste0('
    # Latent variables
    PearsonLatent =~ 1*empirical_r

    # Fixing error variances based on known standard errors
    empirical_r ~~ ',se2,'*empirical_r

    # Relationship between latent variables
    PearsonLatent ~~ synthetic_r
  ')

  fit <- sem(model, data = subset)

  sim_results <- bind_rows(sim_results,
    standardizedsolution(fit) %>% filter(lhs == "PearsonLatent", rhs ==  "synthetic_r")
  )
}

sim_results %>% summarise(semi_latent_r = mean(est.std), mean_se = sqrt(mean(se^2)), max_se = max(se)) %>% 
  knitr::kable()

semi_latent_r	mean_se	max_se
0.707866	0.00316	0.0035834

Precision simulation for synthetic reliabilities

source("global_functions.R")
holdout_human_data = arrow::read_feather(
  file = file.path(data_path, glue("{model_name}.raw.osf-bainbridge-2021-s2-0.human.feather"))
)
holdout_human_data <- holdout_human_data %>% haven::zap_labels()

cors_llm <- holdout_llm %>%
  select(x = variable_1, y = variable_2, r = synthetic_r) %>%
  as.data.frame() |>
  igraph::graph_from_data_frame(directed = FALSE) |>
  igraph::as_adjacency_matrix(attr = "r", sparse = FALSE)
diag(cors_llm) <- 1

cors_real <- holdout_llm %>%
  select(x = variable_1, y = variable_2, r = empirical_r) %>%
  as.data.frame() |>
  igraph::graph_from_data_frame(directed = FALSE) |>
  igraph::as_adjacency_matrix(attr = "r", sparse = FALSE)
diag(cors_real) <- 1

mapping_data <- holdout_mapping_data
items_by_scale <- bind_rows(
  scales %>% filter(scale_1 == "") %>% left_join(mapping_data %>% select(-scale_1), by = c("instrument", "scale_0")),
  scales %>% filter(scale_1 != "") %>% left_join(mapping_data, by = c("instrument", "scale_0", "scale_1"))
)

real_scales <- items_by_scale %>%
  group_by(scale) %>%
  summarise(
    items = list(variable),
    number_of_items = n_distinct(variable),
    keyed = first(keyed),
    lvn = paste(first(scale), " =~ ", paste(variable, collapse = " + "))) %>%
  group_by(scale) %>%
  mutate(reverse_items = list(find_reverse_items_by_first_item(cors_real[unlist(items), unlist(items)], keyed))) %>% 
  drop_na() %>% 
  ungroup()

random_scales <- list()
for(i in 1:1000) {
  n_items <- rpois(1, mean(real_scales$number_of_items,na.rm = T))
  n_items <- if_else(n_items < 3, 3, n_items, 3)
  random_scales[[i]] <- holdout_mapping_data %>%
    sample_n(n_items) %>%
    mutate(scale = paste0("random", i)) %>%
    group_by(scale) %>%
    summarise(
      items = list(variable),
      number_of_items = n_distinct(variable),
      lvn = paste(first(scale), " =~ ", paste(variable, collapse = " + "))) %>%
    drop_na() %>% 
    mutate(keyed = 1)
}

random_scales <- bind_rows(random_scales) %>% 
  distinct(items, .keep_all = TRUE) %>% 
  rowwise() %>% 
  mutate(
    reverse_items = list(randomly_choose_items_for_reversion(items))
    ) %>% 
  ungroup()

scales <- bind_rows(real = real_scales, random = random_scales, .id = "type")
n_distinct(scales$scale)

## [1] 1113

scales <- scales %>% filter(number_of_items >= 3)
scales <- scales %>%
  rowwise() %>%
  mutate(r_llm = list(cors_llm[items, items]),
         r_real = list(cors_real[items, items]),
         N_real = holdout %>% 
           filter(variable_1 %in% items, variable_2 %in% items) %>% 
           summarise(min_n=min(pairwise_n)) %>% pull(min_n)) %>%
  mutate(
    rel_real = list(psych::alpha(holdout_human_data[, items], keys = reverse_items, n.iter = 1000)),
    rel_llm = list(psych::alpha(r_llm, keys = reverse_items, n.obs = N_real)$feldt)) %>%
  mutate(empirical_alpha = rel_real$feldt$alpha$raw_alpha,
         synthetic_alpha = rel_llm$alpha$raw_alpha) %>%
  mutate(
    alpha_se = mean(diff(unlist(psychometric::alpha.CI(empirical_alpha, k = number_of_items, N = planned_N, level = 0.95)))/1.96)
  )

realistic_scales <- scales %>% ungroup()

sim_results <- tibble()
for(i in 1:500) {
  picked_scales <- realistic_scales %>% filter(!str_detect(scale, "random")) %>% sample_n(number_of_scales_with_more_than_3_items)
  subset <-
    bind_rows(picked_scales,
              realistic_scales %>% filter(str_detect(scale, "random")) %>% sample_n(200)
  )

  se2 <- mean(subset$alpha_se^2)

  r <- broom::tidy(cor.test(subset$empirical_alpha, subset$synthetic_alpha))
  (r$conf.high - r$conf.low)/2

  model <- paste0('
    # Latent variables
    PearsonLatent =~ 1*empirical_alpha

    # Fixing error variances based on known standard errors
    empirical_alpha ~~ ',se2,'*empirical_alpha

    # Relationship between latent variables
    PearsonLatent ~~ synthetic_alpha
  ')

  fit <- sem(model, data = subset)

  sim_results <- bind_rows(sim_results,
                           standardizedsolution(fit) %>% filter(lhs == "PearsonLatent", rhs ==  "synthetic_alpha")
  )
}
sim_results %>% summarise(semi_latent_r = mean(est.std), mean_se = sqrt(mean(se^2)), max_se = max(se)) %>% 
  knitr::kable()

semi_latent_r	mean_se	max_se
0.8798027	0.0153378	0.0191323

Precision simulation for synthetic scale correlations

manifest_scores = arrow::read_feather(file = file.path(data_path, glue("data/intermediate/{model_name}.raw.osf-bainbridge-2021-s2-0.scale_correlations.feather")))


sim_results <- tibble()
library(lavaan)

for(i in 1:500) {
  scales <- manifest_scores %>% select(scale_a) %>% distinct() %>% pull(scale_a)

  subset <- manifest_scores %>% filter(scale_a %in% scales, scale_b %in% scales) %>% 
    sample_n(combinations_scales)

  N <- planned_N
  subset <- subset %>% mutate(se = (1 - empirical_r^2)/sqrt(N - 3))
  se2 <- mean(subset$se^2)

  r <- broom::tidy(cor.test(subset$empirical_r, subset$synthetic_r))
  (r$conf.high - r$conf.low)/2

  model <- paste0('
    # Latent variables
    PearsonLatent =~ 1*empirical_r

    # Fixing error variances based on known standard errors
    empirical_r ~~ ',se2,'*empirical_r

    # Relationship between latent variables
    PearsonLatent ~~ synthetic_r
  ')

  fit <- sem(model, data = subset)

  sim_results <- bind_rows(sim_results,
                           standardizedsolution(fit) %>% filter(lhs == "PearsonLatent", rhs ==  "synthetic_r")
  )
}

sim_results %>% summarise(semi_latent_r = mean(est.std), mean_se = sqrt(mean(se^2)), max_se = max(se)) %>% 
  knitr::kable()

semi_latent_r	mean_se	max_se
0.8802321	0.006223	0.0070053