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(
  lavaan
)
## This is lavaan 0.6-15
## lavaan is FREE software! Please report any bugs.
## 
## Attaching package: 'lavaan'
## 
## The following object is masked from 'package:psych':
## 
##     cor2cov
theme_set(theme_bw())

options(
    digits = 3
)

Analysis

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

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

#simulate mz twins
set.seed(1)
results = 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)
  )
})
## 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)
## `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.2.3 (2023-03-15)
## 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       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lavaan_0.6-15         kirkegaard_2023-04-25 psych_2.3.3          
##  [4] assertthat_0.2.1      weights_1.0.4         Hmisc_5.0-1          
##  [7] magrittr_2.0.3        lubridate_1.9.2       forcats_1.0.0        
## [10] stringr_1.5.0         dplyr_1.1.2           purrr_1.0.1          
## [13] readr_2.1.4           tidyr_1.3.0           tibble_3.2.1         
## [16] ggplot2_3.4.2         tidyverse_2.0.0      
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.5        jsonlite_1.8.4    splines_4.2.3     gtools_3.9.4     
##  [5] bslib_0.4.2       Formula_1.2-5     highr_0.10        stats4_4.2.3     
##  [9] yaml_2.3.7        pbivnorm_0.6.0    pillar_1.9.0      backports_1.4.1  
## [13] lattice_0.20-45   quadprog_1.5-8    glue_1.6.2        digest_0.6.31    
## [17] checkmate_2.1.0   minqa_1.2.5       colorspace_2.1-0  htmltools_0.5.5  
## [21] Matrix_1.5-4      pkgconfig_2.0.3   broom_1.0.4       scales_1.2.1     
## [25] gdata_2.18.0.1    tzdb_0.3.0        lme4_1.1-32       timechange_0.2.0 
## [29] htmlTable_2.4.1   mgcv_1.8-42       farver_2.1.1      generics_0.1.3   
## [33] cachem_1.0.7      withr_2.5.0       nnet_7.3-18       cli_3.6.1        
## [37] mnormt_2.1.1      evaluate_0.20     mice_3.15.0       fansi_1.0.4      
## [41] nlme_3.1-162      MASS_7.3-58.3     foreign_0.8-82    tools_4.2.3      
## [45] data.table_1.14.8 hms_1.1.3         lifecycle_1.0.3   munsell_0.5.0    
## [49] cluster_2.1.4     compiler_4.2.3    jquerylib_0.1.4   rlang_1.1.0      
## [53] grid_4.2.3        nloptr_2.0.3      rstudioapi_0.14   htmlwidgets_1.6.2
## [57] labeling_0.4.2    base64enc_0.1-3   rmarkdown_2.21    boot_1.3-28      
## [61] gtable_0.3.3      R6_2.5.1          gridExtra_2.3     knitr_1.42       
## [65] fastmap_1.1.1     utf8_1.2.3        stringi_1.7.12    parallel_4.2.3   
## [69] Rcpp_1.0.10       vctrs_0.6.2       rpart_4.1.19      tidyselect_1.2.0 
## [73] xfun_0.39
#OSF
if (F) {
  library(osfr)
  
  #login
  osf_auth(readr::read_lines("~/.config/osf_token"))
  
  #the project we will use
  osf_proj = osf_retrieve_node("https://osf.io/XXX/")
  
  #upload all files in project
  #overwrite existing (versioning)
  osf_upload(
    osf_proj,
    path = c("data", "figures", "papers", "notebook.Rmd", "notebook.html", "sessions_info.txt"), 
    conflicts = "overwrite"
    )
}