Init

library(kirkegaard)
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── 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
## Loading required package: magrittr
## 
## 
## Attaching package: 'magrittr'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## 
## Loading required package: weights
## 
## Loading required package: Hmisc
## 
## 
## Attaching package: 'Hmisc'
## 
## 
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## 
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## 
## 
## Loading required package: assertthat
## 
## 
## Attaching package: 'assertthat'
## 
## 
## The following object is masked from 'package:tibble':
## 
##     has_name
## 
## 
## Loading required package: psych
## 
## 
## Attaching package: 'psych'
## 
## 
## The following object is masked from 'package:Hmisc':
## 
##     describe
## 
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
## 
## 
## 
## Attaching package: 'kirkegaard'
## 
## 
## The following object is masked from 'package:psych':
## 
##     rescale
## 
## 
## The following object is masked from 'package:assertthat':
## 
##     are_equal
## 
## 
## The following object is masked from 'package:purrr':
## 
##     is_logical
## 
## 
## The following object is masked from 'package:base':
## 
##     +
load_packages(
  
)

theme_set(theme_bw())

options(
    digits = 3
)

Analysis

#parameters
mz_r = .50
dz_r = .25
bw_latent_d = 1
pairs_per_group = 1e5

#matrices
mz_mat = matrix(c(1, mz_r, mz_r, 1), ncol = 2)
dz_mat = matrix(c(1, dz_r, dz_r, 1), ncol = 2)

#bias vector, random chance that gets convicted/arrested for hostile racism reasons
bias_vector = seq(0, 1, by = .025) %>% rep(3) %>% sort()

#function
run_simulations = function(bias_vector, prob_sd = 0, prob_mean = 1) {
  map_df(bias_vector, function(bias) {
  # browser()
  #simulate mz's
  mz_white_latent = MASS::mvrnorm(n = pairs_per_group, mu = c(0, 0), Sigma = mz_mat)
  #dz's
  dz_white_latent = MASS::mvrnorm(n = pairs_per_group, mu = c(0, 0), Sigma = dz_mat)
  
  #blacks
  mz_black_latent = MASS::mvrnorm(n = pairs_per_group, mu = c(bw_latent_d, bw_latent_d), Sigma = mz_mat)
  dz_black_latent = MASS::mvrnorm(n = pairs_per_group, mu = c(bw_latent_d, bw_latent_d), Sigma = dz_mat)
  
  #join vectors
  d = cbind(
    mz_white_latent,
    dz_white_latent,
    mz_black_latent,
    dz_black_latent
  ) %>% as.matrix() %>% set_colnames(c("white_mz_1", "white_mz_2", "white_dz_1", "white_dz_2", "black_mz_1", "black_mz_2", "black_dz_1", "black_dz_2"))
  
  #prob of getting caught
  d_chance_for_caught = pnorm(d, mean = prob_mean, sd = prob_sd, lower.tail = T)
  
  #convert to binary by sampling
  d_caught = matrix(rbinom(length(as.vector(d_chance_for_caught)), size = 1, prob = as.vector(d_chance_for_caught)), nrow = nrow(d))
  
  #add bias
  bias_caughts = rbinom(pairs_per_group*4, size = 1, prob = bias)
  d_caught[, 5:8] = d_caught[, 5:8] + bias_caughts
  d_caught[d_caught > 1] = 1
  
  #twin tetrachorics
  tetras = psych::tetrachoric(d_caught)
  
  #h2s from falconer
  tibble(
    bias = bias,
    white_h2 = try_else((tetras$rho[1, 2]-tetras$rho[3, 4]) * 2, else. = NA),
    black_h2 = try_else((tetras$rho[5, 6]-tetras$rho[7, 8]) * 2, else. = NA)
  )
})
}

#simulate mz twins
set.seed(1)
results = run_simulations(bias_vector)
## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted
## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted
## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted
#with randomness
results_05random = run_simulations(bias_vector, prob_sd = 0.5)
## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted
## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted
## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted

## Warning in psych::tetrachoric(d_caught): Item = had no variance and was deleted
#plot results
results %>% 
  pivot_longer(cols = -bias) %>% 
  ggplot(aes(bias, value, color = name)) +
  geom_smooth() +
  scale_color_discrete("Parameter") +
  scale_x_continuous("Random extra chance of getting caught that only applies to Blacks", labels = scales::percent) +
  scale_y_continuous("Heritability estimate from MZ-DZ design", labels = scales::percent) +
  ggtitle("Based on simulated MZ-DZ data", subtitle = "No randomness")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (`stat_smooth()`).

GG_save("figs/heritability, race, police bias.png")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (`stat_smooth()`).
results_05random %>% 
  pivot_longer(cols = -bias) %>% 
  ggplot(aes(bias, value, color = name)) +
  geom_smooth() +
  scale_color_discrete("Parameter") +
  scale_x_continuous("Random extra chance of getting caught that only applies to Blacks", labels = scales::percent) +
  scale_y_continuous("Heritability estimate from MZ-DZ design", labels = scales::percent) +
  ggtitle("Based on simulated MZ-DZ data", subtitle = "With some randomness")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (`stat_smooth()`).

GG_save("figs/heritability, race, police bias randomness.png")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (`stat_smooth()`).

Meta

#versions
write_sessioninfo()
## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Linux Mint 21.1
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_DK.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_DK.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Berlin
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kirkegaard_2023-04-25 psych_2.3.3           assertthat_0.2.1     
##  [4] weights_1.0.4         Hmisc_5.0-1           magrittr_2.0.3       
##  [7] lubridate_1.9.2       forcats_1.0.0         stringr_1.5.0        
## [10] dplyr_1.1.2           purrr_1.0.1           readr_2.1.4          
## [13] tidyr_1.3.0           tibble_3.2.1          ggplot2_3.4.2        
## [16] tidyverse_2.0.0      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.3      xfun_0.39         bslib_0.4.2       htmlwidgets_1.6.2
##  [5] lattice_0.20-45   tzdb_0.3.0        vctrs_0.6.2       tools_4.3.0      
##  [9] generics_0.1.3    parallel_4.3.0    fansi_1.0.4       highr_0.10       
## [13] cluster_2.1.4     pkgconfig_2.0.3   Matrix_1.5-1      data.table_1.14.8
## [17] checkmate_2.2.0   lifecycle_1.0.3   farver_2.1.1      compiler_4.3.0   
## [21] textshaping_0.3.6 mnormt_2.1.1      munsell_0.5.0     htmltools_0.5.5  
## [25] sass_0.4.5        yaml_2.3.7        htmlTable_2.4.1   Formula_1.2-5    
## [29] mice_3.15.0       nloptr_2.0.3      pillar_1.9.0      jquerylib_0.1.4  
## [33] MASS_7.3-58.3     gdata_2.18.0.1    cachem_1.0.7      boot_1.3-28      
## [37] rpart_4.1.19      nlme_3.1-162      gtools_3.9.4      tidyselect_1.2.0 
## [41] digest_0.6.31     stringi_1.7.12    labeling_0.4.2    splines_4.3.0    
## [45] fastmap_1.1.1     grid_4.3.0        colorspace_2.1-0  cli_3.6.1        
## [49] base64enc_0.1-3   utf8_1.2.3        broom_1.0.4       foreign_0.8-82   
## [53] withr_2.5.0       scales_1.2.1      backports_1.4.1   timechange_0.2.0 
## [57] rmarkdown_2.21    nnet_7.3-18       lme4_1.1-33       gridExtra_2.3    
## [61] ragg_1.2.5        hms_1.1.3         evaluate_0.20     knitr_1.42       
## [65] mgcv_1.8-42       rlang_1.1.0       Rcpp_1.0.10       glue_1.6.2       
## [69] minqa_1.2.5       rstudioapi_0.14   jsonlite_1.8.4    R6_2.5.1         
## [73] systemfonts_1.0.4