Fitting a real-data covariance matrix: interpretation (2)

Superpopulations

covmat <- readRDS("../data/datasets/1kg_phase3_derived/1kg_phase3_derived_covmat.rds")
meta <- readRDS("../data/datasets/1kg_phase3_derived/1kg_phase3_derived_meta.rds")

plot_fl <- function(LL) {
  df <- data.frame(LL)
  colnames(df) <- paste0("Factor ", formatC(1:ncol(LL), width = 2, flag = "0"))
  df$subpop <- meta$pop
  df$superpop <- meta$super_pop
  df <- df %>% arrange(superpop, subpop)
  df$idx <- 1:nrow(df)
  gath_df <- df %>% 
    gather(K, value, -subpop, -idx, -superpop) %>%
    mutate(K = factor(K))
  med_gath_df <- gath_df %>% 
    group_by(subpop, K) %>% 
    summarise(value=median(value), idx=median(idx))
  
  p <- ggplot(gath_df, aes(x=idx, y=value, color=superpop)) + 
    geom_point() +
    facet_wrap(~K) + 
    geom_hline(yintercept = 0, linetype = "dashed") +
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank()) +
    labs(color="superpop")
  return(p)
}

plot_subpops <- function(LL) {
  df <- data.frame(LL)
  colnames(df) <- paste0("Factor ", formatC(1:ncol(LL), width = 2, flag = "0"))
  df$subpop <- meta$pop
  df$superpop <- meta$super_pop
  df <- df %>% arrange(superpop, subpop)
  df$idx <- 1:nrow(df)
  gath_df <- df %>% 
    gather(K, value, -subpop, -idx, -superpop) %>%
    mutate(K = factor(K))
  med_gath_df <- gath_df %>% 
    group_by(subpop, K) %>% 
    summarise(value=median(value), idx=median(idx))
  
  all_plots <- lapply(levels(df$superpop), function(pop) {
    p <- ggplot(filter(gath_df, superpop == pop), aes(x=idx, y=value, color=subpop)) + 
      geom_point() +
      facet_wrap(~K) + 
      geom_hline(yintercept = 0, linetype = "dashed") +
      theme(axis.title.x=element_blank(),
            axis.text.x=element_blank(),
            axis.ticks.x=element_blank()) +
      labs(color="subpop") +
      ggtitle(paste("superpop:", pop))
    return(p)  
  })
}

t_greedy <- system.time({
  ones <- matrix(1, nrow = nrow(covmat), ncol = 1)
  ls.soln <- t(solve(crossprod(ones), crossprod(ones, covmat)))
  
  covmat_diagNA <- covmat
  diag(covmat_diagNA) <- NA
  
  fl_g <- flash.init(covmat_diagNA) %>%
    flash.set.verbose(0) %>%
    flash.init.factors(EF = list(ones, ls.soln)) %>%
    flash.fix.loadings(kset = 1, mode = 1) %>%
    flash.backfit() %>%
    flash.add.greedy(Kmax = 14, 
                     prior.family = prior.point.laplace())
})

cat("Time to fit greedy factors:", round(t_greedy[3], 1), "seconds")

#> Time to fit greedy factors: 79.2 seconds

t_backfit <- system.time({
  fl_bf <- fl_g %>% flash.backfit(tol = 1)
})

cat("Time to backfit:", round(t_backfit[3] / 60, 1), "minutes")

#> Time to backfit: 5.1 minutes

LL <- fl_bf$loadings.pm[[1]] %*% diag(sqrt(fl_bf$loadings.scale))

t_drift1 <- system.time({
  dr <- init_from_covmat(covmat, LL, p = 30000, fix.EL = 1,
                         prior.family = prior.point.laplace())
  dr <- drift(dr, maxiter = 500, tol = 1e-2, verbose = FALSE)
})

#> Warning in mle_point_laplace(x_optset, s_optset, g, control): First
#> optimization attempt failed. Retrying with fewer significant digits.

#> Error in (function (f, p, ..., hessian = FALSE, typsize = rep(1, length(p)),  : 
#>   probable coding error in analytic Hessian

cat("Time to drift (point-Laplace):", round(t_drift1[3] / 60, 1), "minutes")

#> Time to drift (point-Laplace): 0.9 minutes

LL <- dr$EL %*% diag(sqrt(dr$prior_s2))
LL <- LL[, c(1, rep(2:fl_bf$n.factors, each = 2))]
LL <- t(t(LL) * c(1, rep(c(1, -1), fl_bf$n.factors - 1)))
LL <- pmax(LL, 0)

LL_scale <- apply(LL, 2, max)
LL <- t(t(LL) / LL_scale)

zero_cols <- which(LL_scale < 1e-2)
LL <- LL[, -zero_cols]
LL_scale <- LL_scale[-zero_cols]

t_drift2 <- system.time({
  dr_bm <- init_from_covmat(covmat, LL, p = 30000, prior.s2 = LL_scale^2, fix.EL = 1,
                            prior.family = prior.bimodal())
  dr_bm <- drift(dr_bm, maxiter = 500, tol = 1e-4, verbose = FALSE)
})

cat("Time to drift (bimodal):", round(t_drift2[3] / 60, 1), "minutes")

#> Time to drift (bimodal): 27.6 minutes

LL <- dr_bm$EL %*% diag(sqrt(dr_bm$prior_s2))
LL_norms <- apply(LL, 2, function(x) sum(x^2))
LL <- LL[, order(LL_norms, decreasing = TRUE)]
LL_norms <- apply(LL, 2, function(x) sum(x^2))
zero_cols <- which(LL_norms < 1e-6)
LL <- LL[, -zero_cols]

plot(plot_fl(LL))

Subpopulations

all_p <- plot_subpops(LL)
for (p in all_p) {plot(p)}

Factor descriptions

• Factor 2: The primary African factor, with varying degrees of admixture for African Americans (ASW), African Caribbeans (ACB), and Amerindians. Intriguingly, there is a small degree of admixture for many Iberians (IBR).

• Factor 3: The out-of-Africa factor. Again, there are varying degrees of admixture for African Americans, African Caribbeans, and Amerindians.

• Factor 4: The primary East Asian factor. Among Amerindians, loadings are largest among Peruvians and smallest among Puerto Ricans (PUR). Interestingly, there is a small degree of admixture for Finns (FIN) and most South Asians (most notably among Bengalis (BEB)).

• Factor 5: The primary European factor, with strong contributions from South Asian populations (especially Punjabi and Gujarati) and varying degrees of admixture for African Americans, African Caribbeans, and Amerindians. There is a single Vietnamese (KHV) individual with a nonzero loading: this can probably be ignored.

• Factor 6: The Amerindian factor, with loadings largest among Peruvians and smallest among Puerto Ricans (PUR). A few African Americans are also loaded on this factor.

• Factor 7: The primary South Asian factor. Oddly, Peruvians (and a single African Caribbean individual) are also loaded on this factor.

• Factor 8: Substructure among East Asians: large loadings for the more Southeastern Vietnamese and Chinese Dai populations, modest loadings for some Chinese Han, and mostly zero loadings for Japanese.

• Factor 9: Substructure among Europeans: largest loadings among the Finns with more modest loadings for other Northern European populations (GBR and CEU).

• Factor 10: A West African factor: large loadings for Gambians (GWD) and modest loadings for the Mende from Sierra Leone (MSL).

• Factor 11: Kenyan.

• Factor 12: Japanese.

• Factor 13: Common to some, but not all Gujarati Indians from Houston.

• Factor 14: Puerto Rican. I was surprised to see a separate factor for this population.

• Factor 15: Common to some, but not all Punjabi.

• Factor 16: Common to some, but not all Telugu Indians.

• Factor 17: Unique to several Gujarati Indians from Houston.

• Factor 18: Unique to several Sri Lankans.

• Factor 19: Unique to several African Americans.

• Factor 20: Unique to three South Asian individuals (two Telugu and one Sri Lankan).

• Factor 21: Shared among the two Nigerian populations (Esan (ESN) and Yoruba (YRI)).

• Factor 22: Shared among Japanese and most Chinese Han from Beijing.

• Factor 23: Shared among several Toscani and a few Puerto Ricans and Colombians.
  
  

  Session information

sessionInfo()

#> R version 3.5.3 (2019-03-11)
#> Platform: x86_64-apple-darwin15.6.0 (64-bit)
#> Running under: macOS Mojave 10.14.6
#> 
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] forcats_0.4.0      stringr_1.4.0      dplyr_0.8.0.1     
#>  [4] purrr_0.3.2        readr_1.3.1        tidyr_0.8.3       
#>  [7] tibble_2.1.1       ggplot2_3.2.0      tidyverse_1.2.1   
#> [10] drift.alpha_0.0.11 flashier_0.2.7    
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.4.6      lubridate_1.7.4   invgamma_1.1     
#>  [4] lattice_0.20-38   assertthat_0.2.1  rprojroot_1.3-2  
#>  [7] digest_0.6.18     truncnorm_1.0-8   R6_2.4.0         
#> [10] cellranger_1.1.0  plyr_1.8.4        backports_1.1.3  
#> [13] evaluate_0.13     httr_1.4.0        pillar_1.3.1     
#> [16] rlang_0.4.2       lazyeval_0.2.2    readxl_1.3.1     
#> [19] rstudioapi_0.10   ebnm_0.1-21       irlba_2.3.3      
#> [22] whisker_0.3-2     Matrix_1.2-15     rmarkdown_1.12   
#> [25] labeling_0.3      munsell_0.5.0     mixsqp_0.3-40    
#> [28] broom_0.5.1       compiler_3.5.3    modelr_0.1.5     
#> [31] xfun_0.6          pkgconfig_2.0.2   SQUAREM_2017.10-1
#> [34] htmltools_0.3.6   tidyselect_0.2.5  workflowr_1.2.0  
#> [37] withr_2.1.2       crayon_1.3.4      grid_3.5.3       
#> [40] nlme_3.1-137      jsonlite_1.6      gtable_0.3.0     
#> [43] git2r_0.25.2      magrittr_1.5      scales_1.0.0     
#> [46] cli_1.1.0         stringi_1.4.3     reshape2_1.4.3   
#> [49] fs_1.2.7          xml2_1.2.0        generics_0.0.2   
#> [52] tools_3.5.3       glue_1.3.1        hms_0.4.2        
#> [55] parallel_3.5.3    yaml_2.2.0        colorspace_1.4-1 
#> [58] ashr_2.2-51       rvest_0.3.4       knitr_1.22       
#> [61] haven_2.1.1