# How To Simulate Traits

    library(smer)
    library(genio)

For a simple illustration, we generate a synthetic allele count matrix with minor allele frequency `maf` greater than 5%. The function `simulate_traits()` requires the input genotype data to be stored in the PLINK format (`.bim`, `.bed`, and `.fam` files). Here we use the package `genio` to write the allele count matrix to a PLINK file.

    plink_file <- tempfile()
    n_samples <- 3000
    n_snps <- 5000
    maf <- 0.05 + 0.45 * runif(n_snps)
    random_minor_allele_counts   <- (runif(n_samples * n_snps) < maf) +
                                    (runif(n_samples * n_snps) < maf)
    allele_count_matrix <- matrix(random_minor_allele_counts,
      nrow = n_samples,
      ncol = n_snps,
      byrow = TRUE,
    )

    # Create .fam and .bim data frames
    fam <- data.frame(
      fam = sprintf("F%d", 1:n_samples),       # Family ID
      id = sprintf("I%d", 1:n_samples),       # Individual ID
      pat = 0,                        # Paternal ID
      mat = 0,                        # Maternal ID
      sex = sample(1:2, n_samples, replace = TRUE),
      pheno = -9
    )

    bim <- data.frame(
      chr = 1,               # Chromosome number
      id = sprintf("rs%d", 1:n_snps),   # SNP name
      posg = 0,                         # Genetic distance (cM)
      pos = 1:n_snps,          # Base-pair position
      ref = sample(c("A", "C", "G", "T"), n_snps, replace = T), # Minor allele
      alt = sample(c("A", "C", "G", "T"), n_snps, replace = T)  # Major allele
    )

    # Write to .bed, .fam, and .bim files
    write_plink(
      file = plink_file,
      X = t(allele_count_matrix),
      fam = fam,
      bim = bim
    )

To simulate traits, we need to specify the genetic architecture. E.g., the function requires us to specify

- the narrow sense heritability \(h^2\)
- the indices SNPs contributing to the narrow sense heritability
- the heritability due to epistasis
- the indices of SNPs contributing to the epistatic trait variance
- the path to the output file.

The simulated trait is written to file in the PLINK phenotype format.

    pheno_file <- tempfile()
    additive_heritability <- 0.3
    gxg_heritability <- 0.25
    n_snps_gxg_group <- 5
    n_snps_additive <- 100
    additive_snps <- sort(sample(1:n_snps, n_snps_additive, replace = F))
    gxg_group_1 <- sort(sample(additive_snps, n_snps_gxg_group, replace = F))
    gxg_group_2 <- sort(sample(setdiff(additive_snps, gxg_group_1),
                               n_snps_gxg_group, replace = F))

    simulate_traits(
      plink_file,
      pheno_file,
      additive_heritability,
      gxg_heritability,
      additive_snps,
      gxg_group_1,
      gxg_group_2
    )

The variance components of the traits are constructed to match the input. See Crawford et al. (2017) and mvMAPIT Documentation: Simulate Traits for a full description of the simulation scheme.

## SessionInfo

    sessionInfo()
    #> R version 4.5.1 (2025-06-13)
    #> Platform: aarch64-apple-darwin24.5.0
    #> Running under: macOS Sequoia 15.6.1
    #>
    #> Matrix products: default
    #> BLAS:   /opt/homebrew/Cellar/openblas/0.3.30/lib/libopenblasp-r0.3.30.dylib
    #> LAPACK: /opt/homebrew/Cellar/r/4.5.1/lib/R/lib/libRlapack.dylib;  LAPACK version 3.12.1
    #>
    #> locale:
    #> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
    #>
    #> time zone: Europe/Berlin
    #> tzcode source: internal
    #>
    #> attached base packages:
    #> [1] stats4    stats     graphics  grDevices utils     datasets  methods
    #> [8] base
    #>
    #> other attached packages:
    #>  [1] genio_1.1.2          knitr_1.50           tidyr_1.3.1
    #>  [4] GenomicRanges_1.60.0 GenomeInfoDb_1.44.0  IRanges_2.42.0
    #>  [7] S4Vectors_0.46.0     BiocGenerics_0.54.0  generics_0.1.4
    #> [10] ggplot2_3.5.2        dplyr_1.1.4          smer_0.0.2
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] sass_0.4.10             digest_0.6.37           magrittr_2.0.3
    #>  [4] evaluate_1.0.4          grid_4.5.1              RColorBrewer_1.1-3
    #>  [7] iterators_1.0.14        CompQuadForm_1.4.4      mvtnorm_1.3-3
    #> [10] fastmap_1.2.0           foreach_1.5.2           jsonlite_2.0.0
    #> [13] backports_1.5.0         httr_1.4.7              purrr_1.1.0
    #> [16] UCSC.utils_1.4.0        scales_1.4.0            codetools_0.2-20
    #> [19] jquerylib_0.1.4         cli_3.6.5               rlang_1.1.6
    #> [22] XVector_0.48.0          withr_3.0.2             cachem_1.1.0
    #> [25] yaml_2.3.10             FMStable_0.1-4          tools_4.5.1
    #> [28] parallel_4.5.1          checkmate_2.3.2         GenomeInfoDbData_1.2.14
    #> [31] vctrs_0.6.5             R6_2.6.1                lifecycle_1.0.4
    #> [34] mvMAPIT_2.0.3           pkgconfig_2.0.3         pillar_1.11.0
    #> [37] bslib_0.9.0             gtable_0.3.6            glue_1.8.0
    #> [40] Rcpp_1.1.0              harmonicmeanp_3.0.1     xfun_0.52
    #> [43] tibble_3.3.0            tidyselect_1.2.1        farver_2.1.2
    #> [46] htmltools_0.5.8.1       rmarkdown_2.29          labeling_0.4.3
    #> [49] compiler_4.5.1