Start-up code

### INITIAL COMMANDS TO RUN

# libs --------------------------------------------------------------------

library(pacman)
p_load(plyr, readr, nFactors, kirkegaard)

#ad hoc functions
source("scripts/functions.R")

# SNP frequency data --------------------------------------------------------------------


#pop freqs
#load from .csv files if no RDS file
if (!file.exists("data/pop_freqs.RDS")) {
  #load each file, then concatenate
  d_freqs = lapply(1:22, FUN = function(i) {
    read_csv("data/chr" + i + ".csv")
  }) %>% do.call(what = "rbind")
  
  #write to RDS file
  write_rds(x = d_freqs, path = "data/pop_freqs.RDS", compress = "gz")
} else {
  d_freqs = read_rds("data/pop_freqs.RDS")
}


# fst data ----------------------------------------------------------------

d_fst_long = read.csv("data/chr1_fst.csv", stringsAsFactors = F)

#make matrix of fst
m_fst = matrix(0, ncol = length(v_1kg_pops), nrow = length(v_1kg_pops))
rownames(m_fst) = colnames(m_fst) = v_1kg_pops

for (i in 1:nrow(d_fst_long)) {
  m_fst[d_fst_long[i, 1], d_fst_long[i, 2]] = d_fst_long[i, 3]
  m_fst[d_fst_long[i, 2], d_fst_long[i, 1]] = d_fst_long[i, 3]
}


# phenotype data ----------------------------------------------------------

d_phenotype = read.csv("data/piffer_table.csv", na.strings = "")
rownames(d_phenotype) = v_1kg_pops

Analyses

### ANALYZE THE SNPS FOUND IN OKBAY ET AL 2016

#main datafile
l_OKBAY = list()


# SNP results -------------------------------------------------------------

l_OKBAY$SNP_data = read_tsv("data/EduYears_Pooled_5000.txt") %>% as.data.frame
## Parsed with column specification:
## cols(
##   MarkerName = col_character(),
##   CHR = col_integer(),
##   POS = col_integer(),
##   A1 = col_character(),
##   A2 = col_character(),
##   EAF = col_double(),
##   Beta = col_double(),
##   SE = col_double(),
##   Pval = col_double()
## )
#rownames
rownames(l_OKBAY$SNP_data) = l_OKBAY$SNP_data$MarkerName

# LD pruned SNPs ----------------------------------------------------------

#this has both SNP and genotype
#pooled results
l_OKBAY$pooled_results = read_csv("data/OKBAY_total_pooled.csv") %>% as.data.frame()
## Parsed with column specification:
## cols(
##   SNP = col_character(),
##   Chr = col_integer(),
##   Position = col_integer(),
##   `Allele 1` = col_character(),
##   `Frequency Allele 1` = col_double(),
##   `Effect size` = col_double(),
##   `P-value` = col_double(),
##   `Heterogeneity I2` = col_double(),
##   `Heterogeneity P-value` = col_double()
## )
rownames(l_OKBAY$pooled_results) = l_OKBAY$pooled_results$SNP

#results without UKBB
l_OKBAY$discovery_results = read_csv("data/OKBAY_total_discovery.csv") %>% as.data.frame()
## Warning: Duplicated column names deduplicated: 'Effect size' => 'Effect
## size_1' [10], 'P-value' => 'P-value_1' [11], 'Effect size' => 'Effect
## size_2' [12], 'P-value' => 'P-value_2' [13]
## Parsed with column specification:
## cols(
##   SNP = col_character(),
##   Chr = col_integer(),
##   Position = col_integer(),
##   `Allele 1` = col_character(),
##   `Frequency Allele 1` = col_double(),
##   `Effect size` = col_double(),
##   `P-value` = col_double(),
##   `Heterogeneity I2` = col_double(),
##   `Heterogeneity P-value` = col_double(),
##   `Effect size_1` = col_double(),
##   `P-value_1` = col_double(),
##   `Effect size_2` = col_double(),
##   `P-value_2` = col_double()
## )
rownames(l_OKBAY$discovery_results) = l_OKBAY$discovery_results$SNP

#how many SNPs in common?
sum(l_OKBAY$pooled_results$SNP %in% d_freqs$SNP)
## [1] 154
sum(l_OKBAY$discovery_results$SNP %in% d_freqs$SNP)
## [1] 71
# subset ------------------------------------------------------------------
#subset to the ones not in LD

l_OKBAY$pop_freqs = d_freqs[d_freqs$SNP %in% l_OKBAY$pooled_results$SNP, ]


# get effect allele freqs -------------------------------------------------
l_OKBAY$effect_freqs = get_freq_effect_alleles(freq_data = l_OKBAY$pop_freqs,
                                               SNPs = l_OKBAY$pooled_results$SNP, 
                                               alleles = l_OKBAY$pooled_results$`Allele 1`#,
                                               #betas = l_OKBAY$pooled_results$`Effect size`
                                               )

#output
write.csv(l_OKBAY$effect_freqs, file = "data/effect_freqs.csv")




# factor analyze ----------------------------------------------------------
#check method variance for analysis
l_OKBAY$FA_all = fa_all_methods(l_OKBAY$effect_freqs, messages = F)
wtd.cors(l_OKBAY$FA_all$scores)
##                   regression_minres Thurstone_minres tenBerge_minres
## regression_minres         1.0000000        1.0000000       1.0000000
## Thurstone_minres          1.0000000        1.0000000       1.0000000
## tenBerge_minres           1.0000000        1.0000000       1.0000000
## Bartlett_minres           0.9981106        0.9981106       0.9981106
## regression_pa             0.9999564        0.9999564       0.9999564
## Thurstone_pa              0.9999564        0.9999564       0.9999564
## Bartlett_pa               0.9979269        0.9979269       0.9979269
## regression_minchi         1.0000000        1.0000000       1.0000000
## Thurstone_minchi          1.0000000        1.0000000       1.0000000
## tenBerge_minchi           1.0000000        1.0000000       1.0000000
## Bartlett_minchi           0.9981106        0.9981106       0.9981106
##                   Bartlett_minres regression_pa Thurstone_pa Bartlett_pa
## regression_minres       0.9981106     0.9999564    0.9999564   0.9979269
## Thurstone_minres        0.9981106     0.9999564    0.9999564   0.9979269
## tenBerge_minres         0.9981106     0.9999564    0.9999564   0.9979269
## Bartlett_minres         1.0000000     0.9982762    0.9982762   0.9999953
## regression_pa           0.9982762     1.0000000    1.0000000   0.9980969
## Thurstone_pa            0.9982762     1.0000000    1.0000000   0.9980969
## Bartlett_pa             0.9999953     0.9980969    0.9980969   1.0000000
## regression_minchi       0.9981106     0.9999564    0.9999564   0.9979269
## Thurstone_minchi        0.9981106     0.9999564    0.9999564   0.9979269
## tenBerge_minchi         0.9981106     0.9999564    0.9999564   0.9979269
## Bartlett_minchi         1.0000000     0.9982762    0.9982762   0.9999953
##                   regression_minchi Thurstone_minchi tenBerge_minchi
## regression_minres         1.0000000        1.0000000       1.0000000
## Thurstone_minres          1.0000000        1.0000000       1.0000000
## tenBerge_minres           1.0000000        1.0000000       1.0000000
## Bartlett_minres           0.9981106        0.9981106       0.9981106
## regression_pa             0.9999564        0.9999564       0.9999564
## Thurstone_pa              0.9999564        0.9999564       0.9999564
## Bartlett_pa               0.9979269        0.9979269       0.9979269
## regression_minchi         1.0000000        1.0000000       1.0000000
## Thurstone_minchi          1.0000000        1.0000000       1.0000000
## tenBerge_minchi           1.0000000        1.0000000       1.0000000
## Bartlett_minchi           0.9981106        0.9981106       0.9981106
##                   Bartlett_minchi
## regression_minres       0.9981106
## Thurstone_minres        0.9981106
## tenBerge_minres         0.9981106
## Bartlett_minres         1.0000000
## regression_pa           0.9982762
## Thurstone_pa            0.9982762
## Bartlett_pa             0.9999953
## regression_minchi       0.9981106
## Thurstone_minchi        0.9981106
## tenBerge_minchi         0.9981106
## Bartlett_minchi         1.0000000
#how many to extract?
nScree(l_OKBAY$effect_freqs)
##   noc naf nparallel nkaiser
## 1  11   2        11      11
#extract factors
l_OKBAY$fa = fa(l_OKBAY$effect_freqs, scores = "Bartlett")
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
## done
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined.
## Chi square is based upon observed residuals.
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined for the null model either.
## The Chi square is thus based upon observed correlations.
## In factor.stats, I could not find the RMSEA upper bound . Sorry about that
l_OKBAY$fa_2 = fa(l_OKBAY$effect_freqs, scores = "Bartlett", nfactors = 2)
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Loading required namespace: GPArotation
## Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
## done
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined.
## Chi square is based upon observed residuals.
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined for the null model either.
## The Chi square is thus based upon observed correlations.
## In factor.stats, I could not find the RMSEA upper bound . Sorry about that
l_OKBAY$fa_3 = fa(l_OKBAY$effect_freqs, scores = "Bartlett", nfactors = 3)
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done

## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
## done
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined.
## Chi square is based upon observed residuals.
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined for the null model either.
## The Chi square is thus based upon observed correlations.
## In factor.stats, I could not find the RMSEA upper bound . Sorry about that
l_OKBAY$fa_4 = fa(l_OKBAY$effect_freqs, scores = "Bartlett", nfactors = 4)
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in GPFoblq(L, Tmat = Tmat, normalize = normalize, eps = eps, maxit
## = maxit, : convergence not obtained in GPFoblq. 1000 iterations used.
## Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
## done
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined.
## Chi square is based upon observed residuals.
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined for the null model either.
## The Chi square is thus based upon observed correlations.
## In factor.stats, I could not find the RMSEA upper bound . Sorry about that
l_OKBAY$fa_5 = fa(l_OKBAY$effect_freqs, scores = "Bartlett", nfactors = 5)
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in GPFoblq(L, Tmat = Tmat, normalize = normalize, eps = eps, maxit
## = maxit, : convergence not obtained in GPFoblq. 1000 iterations used.
## Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
## done
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined.
## Chi square is based upon observed residuals.
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined for the null model either.
## The Chi square is thus based upon observed correlations.
## In factor.stats, I could not find the RMSEA upper bound . Sorry about that
l_OKBAY$fa_11 = fa(l_OKBAY$effect_freqs, scores = "Bartlett", nfactors = 11)
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
## done
## Warning in GPFoblq(L, Tmat = Tmat, normalize = normalize, eps = eps, maxit
## = maxit, : convergence not obtained in GPFoblq. 1000 iterations used.
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.
## Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
## done
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined.
## Chi square is based upon observed residuals.
## The determinant of the smoothed correlation was zero.
## This means the objective function is not defined for the null model either.
## The Chi square is thus based upon observed correlations.
## In factor.stats, I could not find the RMSEA upper bound . Sorry about that
#plot loadings
fa_plot_loadings(l_OKBAY$fa)

psych::describe(l_OKBAY$fa$loadings)
##     vars   n mean   sd median trimmed  mad   min  max range  skew kurtosis
## MR1    1 154 0.11 0.61   0.27    0.13 0.78 -0.96 0.96  1.91 -0.21    -1.56
##       se
## MR1 0.05
#does loadings match effect sizes?
l_OKBAY$included_SNPs = dimnames(l_OKBAY$fa$loadings)[[1]]

l_OKBAY$SNP_data[l_OKBAY$included_SNPs, "Beta"]
##   [1]  0.016 -0.016  0.018 -0.017  0.022  0.029  0.036 -0.016  0.026 -0.029
##  [11]  0.013 -0.017 -0.013  0.018 -0.014  0.014  0.020  0.018  0.029  0.016
##  [21] -0.014  0.020  0.016  0.016  0.022  0.025  0.017  0.014 -0.015 -0.015
##  [31]  0.025  0.035  0.017  0.022 -0.016  0.014 -0.014  0.016 -0.014  0.014
##  [41] -0.012  0.014  0.021 -0.015  0.023  0.023 -0.030 -0.025  0.026  0.025
##  [51] -0.024  0.021  0.038  0.022  0.018  0.012 -0.017  0.014  0.019 -0.014
##  [61]  0.015 -0.017 -0.013  0.015  0.017  0.014  0.013  0.019 -0.013 -0.017
##  [71] -0.021  0.016 -0.014 -0.013 -0.019 -0.013 -0.013 -0.018  0.013 -0.021
##  [81]  0.013  0.016 -0.013 -0.017  0.015 -0.024  0.028 -0.015  0.014  0.012
##  [91]  0.013 -0.014 -0.015 -0.014  0.016 -0.015  0.015  0.019 -0.014  0.014
## [101] -0.016  0.016  0.017 -0.016 -0.015  0.025  0.015  0.014  0.013 -0.014
## [111]  0.013 -0.014 -0.012 -0.028  0.015 -0.016 -0.017  0.016  0.014  0.014
## [121]  0.014 -0.016 -0.014  0.020  0.016 -0.019 -0.013  0.022 -0.024 -0.015
## [131] -0.016 -0.016 -0.017 -0.018 -0.015  0.017 -0.013  0.019  0.015 -0.014
## [141]  0.015  0.013 -0.013  0.018 -0.013 -0.015 -0.013  0.020 -0.029  0.016
## [151]  0.023 -0.013 -0.022 -0.014
GG_scatter(data.frame("loading" = as.vector(l_OKBAY$fa$loadings),
                      "Beta" = l_OKBAY$pooled_results[l_OKBAY$included_SNPs, "Effect size"])
           , x_var = "Beta", y_var = "loading")

d_phenotype$fs = l_OKBAY$fa$scores
d_phenotype$fs2 = l_OKBAY$fa_2$scores[, 1]
d_phenotype$fs3 = l_OKBAY$fa_3$scores[, 1]
d_phenotype$fs4 = l_OKBAY$fa_4$scores[, 1]
d_phenotype$fs5 = l_OKBAY$fa_5$scores[, 1]
d_phenotype$fs11 = l_OKBAY$fa_11$scores[, 1]

# polygenic scores --------------------------------------------------------
#score each freq according to its reported beta = weighted scores
d_phenotype$ps = calc_polygenic_score(l_OKBAY$effect_freqs, effects = l_OKBAY$pooled_results[l_OKBAY$included_SNPs, "Effect size"])


# phenotype relationships --------------------------------------------------
#plot the genomic predictions against the phenotypic values
GG_scatter(d_phenotype, x_var = "fs", y_var = "IQ")

GG_scatter(d_phenotype, x_var = "ps", y_var = "IQ")

#others
wtd.cors(cbind(d_phenotype, l_OKBAY$fa$scores, l_OKBAY$fa_2$scores, l_OKBAY$fa_3$scores, l_OKBAY$fa_4$scores, l_OKBAY$fa_5$scores))
## Warning in wtd.cors(cbind(d_phenotype, l_OKBAY$fa$scores, l_OKBAY
## $fa_2$scores, : NAs introduced by coercion
## Warning in wtd.cors(cbind(d_phenotype, l_OKBAY$fa$scores, l_OKBAY
## $fa_2$scores, : NAs introduced by coercion
##                       Population Rietveld.et.al.._2013 PS_Ed_Att_Davies
## Population                   NaN                   NaN              NaN
## Rietveld.et.al.._2013        NaN            1.00000000     0.8795920962
## PS_Ed_Att_Davies             NaN            0.87959210     1.0000000000
## PS_Ed_Att_Okbay              NaN            0.68685154     0.8497353377
## IQ                           NaN            0.91969359     0.8210296837
## fs                           NaN            0.78856747     0.5858493240
## fs2                          NaN            0.84638362     0.9115220546
## fs3                          NaN            0.87490220     0.9416790614
## fs4                          NaN            0.78643601     0.8156488855
## fs5                          NaN            0.80216597     0.9023469208
## fs11                         NaN           -0.03266627    -0.0492278077
## ps                           NaN            0.90836969     0.8926049519
## MR1                          NaN            0.78856747     0.5858493240
## MR1                          NaN            0.84638362     0.9115220546
## MR2                          NaN            0.18881588    -0.1949026832
## MR1                          NaN            0.87490220     0.9416790614
## MR2                          NaN            0.30273427    -0.1088470084
## MR3                          NaN           -0.17003340    -0.1445074609
## MR1                          NaN            0.78643601     0.8156488855
## MR2                          NaN            0.28724924    -0.1238291365
## MR3                          NaN           -0.34523478    -0.2937744911
## MR4                          NaN            0.09486742     0.1378107154
## MR1                          NaN            0.80216597     0.9023469208
## MR2                          NaN            0.06997185     0.0066447189
## MR4                          NaN            0.04598254    -0.0491776284
## MR3                          NaN           -0.30478318    -0.2423539115
## MR5                          NaN            0.05967042    -0.0001710936
##                       PS_Ed_Att_Okbay          IQ         fs         fs2
## Population                        NaN         NaN        NaN         NaN
## Rietveld.et.al.._2013      0.68685154  0.91969359 0.78856747  0.84638362
## PS_Ed_Att_Davies           0.84973534  0.82102968 0.58584932  0.91152205
## PS_Ed_Att_Okbay            1.00000000  0.65366840 0.16061663  0.64328833
## IQ                         0.65366840  1.00000000 0.75454518  0.77616263
## fs                         0.16061663  0.75454518 1.00000000  0.72108580
## fs2                        0.64328833  0.77616263 0.72108580  1.00000000
## fs3                        0.67628889  0.81969957 0.72295100  0.97146429
## fs4                        0.58983960  0.71528180 0.65483399  0.89354063
## fs5                        0.64666952  0.77325104 0.67075040  0.94027226
## fs11                      -0.10341763 -0.03686627 0.03708935  0.06322191
## ps                         0.78221849  0.83676254 0.66228301  0.79229709
## MR1                        0.16061663  0.75454518 1.00000000  0.72108580
## MR1                        0.64328833  0.77616263 0.72108580  1.00000000
## MR2                       -0.47732895  0.26428908 0.62218729 -0.07758889
## MR1                        0.67628889  0.81969957 0.72295100  0.97146429
## MR2                       -0.25446889  0.35941809 0.56002378 -0.06629354
## MR3                       -0.30321964 -0.25276637 0.10354522  0.03492828
## MR1                        0.58983960  0.71528180 0.65483399  0.89354063
## MR2                       -0.27506573  0.31181626 0.53090649 -0.07613008
## MR3                       -0.57089797 -0.35401303 0.11984993 -0.02330972
## MR4                        0.05565769  0.05585519 0.12777696  0.10545063
## MR1                        0.64666952  0.77325104 0.67075040  0.94027226
## MR2                        0.03687246  0.01361237 0.09458786  0.05699563
## MR4                       -0.17407546 -0.14462436 0.21279185  0.12814594
## MR3                       -0.50121970 -0.19630632 0.12837618 -0.04068641
## MR5                       -0.06892137  0.14341002 0.07625123 -0.08068191
##                                fs3          fs4          fs5         fs11
## Population                     NaN          NaN          NaN          NaN
## Rietveld.et.al.._2013  0.874902202  0.786436015  0.802165972 -0.032666269
## PS_Ed_Att_Davies       0.941679061  0.815648886  0.902346921 -0.049227808
## PS_Ed_Att_Okbay        0.676288891  0.589839601  0.646669523 -0.103417635
## IQ                     0.819699569  0.715281803  0.773251043 -0.036866269
## fs                     0.722950998  0.654833994  0.670750402  0.037089351
## fs2                    0.971464292  0.893540631  0.940272258  0.063221912
## fs3                    1.000000000  0.893611642  0.943955046 -0.026052238
## fs4                    0.893611642  1.000000000  0.943719307  0.021545099
## fs5                    0.943955046  0.943719307  1.000000000  0.023641276
## fs11                  -0.026052238  0.021545099  0.023641276  1.000000000
## ps                     0.831233385  0.660154818  0.716109231 -0.108740576
## MR1                    0.722950998  0.654833994  0.670750402  0.037089351
## MR1                    0.971464292  0.893540631  0.940272258  0.063221912
## MR2                   -0.078188273 -0.086656002 -0.119590537  0.047733954
## MR1                    1.000000000  0.893611642  0.943955046 -0.026052238
## MR2                   -0.054233612 -0.075914472 -0.130102925  0.042819224
## MR3                    0.052358486 -0.117120789 -0.157797793 -0.009081295
## MR1                    0.893611642  1.000000000  0.943719307  0.021545099
## MR2                   -0.086426422 -0.034301152 -0.085820984  0.082488021
## MR3                   -0.061976693  0.045146332  0.023375287  0.171904580
## MR4                    0.139773991 -0.267667156 -0.118903076 -0.104843979
## MR1                    0.943955046  0.943719307  1.000000000  0.023641276
## MR2                    0.001604441 -0.026340153 -0.001491502  0.181541772
## MR4                    0.111908633 -0.036024953 -0.124626419 -0.157540807
## MR3                   -0.111385744  0.111153985  0.112118639  0.178521796
## MR5                   -0.011840236 -0.008744385 -0.023990189 -0.165770194
##                                 ps        MR1         MR1         MR2
## Population                     NaN        NaN         NaN         NaN
## Rietveld.et.al.._2013  0.908369691 0.78856747  0.84638362  0.18881588
## PS_Ed_Att_Davies       0.892604952 0.58584932  0.91152205 -0.19490268
## PS_Ed_Att_Okbay        0.782218489 0.16061663  0.64328833 -0.47732895
## IQ                     0.836762538 0.75454518  0.77616263  0.26428908
## fs                     0.662283014 1.00000000  0.72108580  0.62218729
## fs2                    0.792297085 0.72108580  1.00000000 -0.07758889
## fs3                    0.831233385 0.72295100  0.97146429 -0.07818827
## fs4                    0.660154818 0.65483399  0.89354063 -0.08665600
## fs5                    0.716109231 0.67075040  0.94027226 -0.11959054
## fs11                  -0.108740576 0.03708935  0.06322191  0.04773395
## ps                     1.000000000 0.66228301  0.79229709  0.07688113
## MR1                    0.662283014 1.00000000  0.72108580  0.62218729
## MR1                    0.792297085 0.72108580  1.00000000 -0.07758889
## MR2                    0.076881134 0.62218729 -0.07758889  1.00000000
## MR1                    0.831233385 0.72295100  0.97146429 -0.07818827
## MR2                    0.235541566 0.56002378 -0.06629354  0.92694960
## MR3                   -0.025919694 0.10354522  0.03492828  0.10714323
## MR1                    0.660154818 0.65483399  0.89354063 -0.08665600
## MR2                    0.149334450 0.53090649 -0.07613008  0.88531467
## MR3                   -0.398580014 0.11984993 -0.02330972  0.15508404
## MR4                    0.293091339 0.12777696  0.10545063  0.07493190
## MR1                    0.716109231 0.67075040  0.94027226 -0.11959054
## MR2                    0.102745228 0.09458786  0.05699563  0.15901859
## MR4                    0.133971929 0.21279185  0.12814594  0.19662492
## MR3                   -0.445463774 0.12837618 -0.04068641  0.19325197
## MR5                   -0.009929531 0.07625123 -0.08068191  0.11917071
##                                MR1         MR2          MR3          MR1
## Population                     NaN         NaN          NaN          NaN
## Rietveld.et.al.._2013  0.874902202  0.30273427 -0.170033402  0.786436015
## PS_Ed_Att_Davies       0.941679061 -0.10884701 -0.144507461  0.815648886
## PS_Ed_Att_Okbay        0.676288891 -0.25446889 -0.303219645  0.589839601
## IQ                     0.819699569  0.35941809 -0.252766373  0.715281803
## fs                     0.722950998  0.56002378  0.103545217  0.654833994
## fs2                    0.971464292 -0.06629354  0.034928277  0.893540631
## fs3                    1.000000000 -0.05423361  0.052358486  0.893611642
## fs4                    0.893611642 -0.07591447 -0.117120789  1.000000000
## fs5                    0.943955046 -0.13010292 -0.157797793  0.943719307
## fs11                  -0.026052238  0.04281922 -0.009081295  0.021545099
## ps                     0.831233385  0.23554157 -0.025919694  0.660154818
## MR1                    0.722950998  0.56002378  0.103545217  0.654833994
## MR1                    0.971464292 -0.06629354  0.034928277  0.893540631
## MR2                   -0.078188273  0.92694960  0.107143225 -0.086656002
## MR1                    1.000000000 -0.05423361  0.052358486  0.893611642
## MR2                   -0.054233612  1.00000000  0.040795012 -0.075914472
## MR3                    0.052358486  0.04079501  1.000000000 -0.117120789
## MR1                    0.893611642 -0.07591447 -0.117120789  1.000000000
## MR2                   -0.086426422  0.91197759 -0.235635023 -0.034301152
## MR3                   -0.061976693 -0.11525079  0.339863763  0.045146332
## MR4                    0.139773991  0.07446248  0.561178817 -0.267667156
## MR1                    0.943955046 -0.13010292 -0.157797793  0.943719307
## MR2                    0.001604441  0.24847448  0.059911029 -0.026340153
## MR4                    0.111908633  0.18231481  0.783828066 -0.036024953
## MR3                   -0.111385744 -0.08519302 -0.023812979  0.111153985
## MR5                   -0.011840236  0.05345406 -0.169207851 -0.008744385
##                                MR2          MR3         MR4          MR1
## Population                     NaN          NaN         NaN          NaN
## Rietveld.et.al.._2013  0.287249242 -0.345234780  0.09486742  0.802165972
## PS_Ed_Att_Davies      -0.123829137 -0.293774491  0.13781072  0.902346921
## PS_Ed_Att_Okbay       -0.275065727 -0.570897969  0.05565769  0.646669523
## IQ                     0.311816264 -0.354013027  0.05585519  0.773251043
## fs                     0.530906494  0.119849932  0.12777696  0.670750402
## fs2                   -0.076130079 -0.023309720  0.10545063  0.940272258
## fs3                   -0.086426422 -0.061976693  0.13977399  0.943955046
## fs4                   -0.034301152  0.045146332 -0.26766716  0.943719307
## fs5                   -0.085820984  0.023375287 -0.11890308  1.000000000
## fs11                   0.082488021  0.171904580 -0.10484398  0.023641276
## ps                     0.149334450 -0.398580014  0.29309134  0.716109231
## MR1                    0.530906494  0.119849932  0.12777696  0.670750402
## MR1                   -0.076130079 -0.023309720  0.10545063  0.940272258
## MR2                    0.885314669  0.155084041  0.07493190 -0.119590537
## MR1                   -0.086426422 -0.061976693  0.13977399  0.943955046
## MR2                    0.911977594 -0.115250786  0.07446248 -0.130102925
## MR3                   -0.235635023  0.339863763  0.56117882 -0.157797793
## MR1                   -0.034301152  0.045146332 -0.26766716  0.943719307
## MR2                    1.000000000 -0.004898466 -0.12947637 -0.085820984
## MR3                   -0.004898466  1.000000000 -0.13142974  0.023375287
## MR4                   -0.129476367 -0.131429744  1.00000000 -0.118903076
## MR1                   -0.085820984  0.023375287 -0.11890308  1.000000000
## MR2                    0.141494914 -0.049762808  0.10956563 -0.001491502
## MR4                   -0.025942085  0.145390340  0.60091993 -0.124626419
## MR3                    0.041817000  0.667492155 -0.45093114  0.112118639
## MR5                    0.168418092 -0.049737361 -0.09444756 -0.023990189
##                                MR2         MR4         MR3           MR5
## Population                     NaN         NaN         NaN           NaN
## Rietveld.et.al.._2013  0.069971845  0.04598254 -0.30478318  0.0596704246
## PS_Ed_Att_Davies       0.006644719 -0.04917763 -0.24235391 -0.0001710936
## PS_Ed_Att_Okbay        0.036872460 -0.17407546 -0.50121970 -0.0689213736
## IQ                     0.013612370 -0.14462436 -0.19630632  0.1434100195
## fs                     0.094587858  0.21279185  0.12837618  0.0762512340
## fs2                    0.056995629  0.12814594 -0.04068641 -0.0806819107
## fs3                    0.001604441  0.11190863 -0.11138574 -0.0118402362
## fs4                   -0.026340153 -0.03602495  0.11115399 -0.0087443853
## fs5                   -0.001491502 -0.12462642  0.11211864 -0.0239901886
## fs11                   0.181541772 -0.15754081  0.17852180 -0.1657701938
## ps                     0.102745228  0.13397193 -0.44546377 -0.0099295312
## MR1                    0.094587858  0.21279185  0.12837618  0.0762512340
## MR1                    0.056995629  0.12814594 -0.04068641 -0.0806819107
## MR2                    0.159018591  0.19662492  0.19325197  0.1191707110
## MR1                    0.001604441  0.11190863 -0.11138574 -0.0118402362
## MR2                    0.248474480  0.18231481 -0.08519302  0.0534540602
## MR3                    0.059911029  0.78382807 -0.02381298 -0.1692078507
## MR1                   -0.026340153 -0.03602495  0.11115399 -0.0087443853
## MR2                    0.141494914 -0.02594208  0.04181700  0.1684180924
## MR3                   -0.049762808  0.14539034  0.66749215 -0.0497373607
## MR4                    0.109565626  0.60091993 -0.45093114 -0.0944475556
## MR1                   -0.001491502 -0.12462642  0.11211864 -0.0239901886
## MR2                    1.000000000  0.21320909 -0.10027722 -0.9349433730
## MR4                    0.213209085  1.00000000 -0.28034520 -0.2802294493
## MR3                   -0.100277223 -0.28034520  1.00000000  0.0587644890
## MR5                   -0.934943373 -0.28022945  0.05876449  1.0000000000
# SAC analysis ------------------------------------------------------------
#which rows has missing data?
v_keeprows = !apply(d_phenotype, MARGIN = 1, anyNA)
#the functions cannot handle missing data, so we need to subset

#is there SAC in the residuals?
#get model resids
d_phenotype$IQ_resid_fs = lm("IQ ~ fs", data = d_phenotype, na.action = na.exclude) %>% resid()
d_phenotype$IQ_resid_ps = lm("IQ ~ ps", data = d_phenotype, na.action = na.exclude) %>% resid()

#measure SAC
SAC_measures(d_phenotype[v_keeprows, ], dists = m_fst[v_keeprows, v_keeprows], vars = c("IQ", "IQ_resid_fs", "IQ_resid_ps"))
## Factors were converted to characters.
## Factors were converted to characters.
## Factors were converted to characters.
##              Morans_I    knsn_3
## IQ          0.6579384 0.9767737
## IQ_resid_fs 0.4583273 0.9553662
## IQ_resid_ps 0.2174602 0.8575929
#correct for SAC
SAC_control(d_phenotype[v_keeprows, ], dependent = "IQ", predictors = "fs", dists = m_fst[v_keeprows, v_keeprows])
## Factors were converted to characters.
##    Uncorrected KNSNR_b_k3    SLR_k3
## fs   0.7545452  0.5642874 0.6990354
SAC_control(d_phenotype[v_keeprows, ], dependent = "IQ", predictors = "fs2", dists = m_fst[v_keeprows, v_keeprows])
## Factors were converted to characters.
##     Uncorrected KNSNR_b_k3   SLR_k3
## fs2   0.7761626  0.3174107 0.380124
SAC_control(d_phenotype[v_keeprows, ], dependent = "IQ", predictors = "fs3", dists = m_fst[v_keeprows, v_keeprows])
## Factors were converted to characters.
##     Uncorrected KNSNR_b_k3    SLR_k3
## fs3   0.8196996  0.5391488 0.4565435
SAC_control(d_phenotype[v_keeprows, ], dependent = "IQ", predictors = "fs4", dists = m_fst[v_keeprows, v_keeprows])
## Factors were converted to characters.
##     Uncorrected KNSNR_b_k3     SLR_k3
## fs4   0.7152818 0.07617423 -0.3111954
SAC_control(d_phenotype[v_keeprows, ], dependent = "IQ", predictors = "fs5", dists = m_fst[v_keeprows, v_keeprows])
## Factors were converted to characters.
##     Uncorrected KNSNR_b_k3     SLR_k3
## fs5    0.773251  0.3631436 0.08125434
SAC_control(d_phenotype[v_keeprows, ], dependent = "IQ", predictors = "fs11", dists = m_fst[v_keeprows, v_keeprows])
## Factors were converted to characters.
##      Uncorrected KNSNR_b_k3   SLR_k3
## fs11 -0.03686627  0.2492557 0.507478
SAC_control(d_phenotype[v_keeprows, ], dependent = "IQ", predictors = "ps", dists = m_fst[v_keeprows, v_keeprows])
## Factors were converted to characters.
##    Uncorrected KNSNR_b_k3    SLR_k3
## ps   0.8367625  0.1902149 0.3051521
# random snps -------------------------------------------------------------
#resample random SNPs, extract their 1st factor and run analyses to get an idea of what chance results look like

set.seed(1) #reproducible
silence({
  d_simul = ldply(1:1000, .fun = function(i) {
    #pick random SNPs, same number as in study
    d_tmp = d_freqs[sample(1:nrow(d_freqs), size = length(l_OKBAY$included_SNPs)), v_1kg_pops] %>% df_t()
    
    #factor analyze them
    fa_tmp = silence(fa(d_tmp, scores = "Bartlett"))
    
    #save scores
    d_tmp_pheno = d_phenotype #copy of pheno data
    d_tmp_pheno$fs = fa_tmp$scores
    
    #do SAC control
    SAC_control(d_tmp_pheno[v_keeprows, ], dependent = "IQ", predictors = "fs", dists = m_fst[v_keeprows, v_keeprows], SLR_include_self = T)
  })
})


GG_denhist(d_simul, var = "Uncorrected")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.