A research group (Sniekers et al. 2017) yesterday published the largest GWAS of intelligence. I will calculate polygenic scores using the independent hits and compare them with polygenic scores computed from previous educational attainment GWAS and reported in Piffer (2017). You will need to load 5 external files (in this order):“Nature2017IntelligenceGWAS.csv”,“df_merged_Nature_top_positive.csv”,“Int_EA_Replicated.csv”, “Matched_18_freqs.csv”,“Replication_linked_Okbay.csv”, “FstDistances.csv” Download them from: osf.io/5yhf
Load packages
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(RMySQL)## Loading required package: DBICall to 1000 Genomes database
sql = src_mysql("population_freqs", host = "128.199.203.208", user = "anon", port = 3306)
sql_1kg = tbl(sql, "1000genomes_phase3")Load list of SNPs from file Nature2017IntelligenceGWAS.csv
Nature2017IntelligenceGWAS=read.csv(file.choose(),header=TRUE)
Nature2017_ind=Nature2017IntelligenceGWAS
snps = Nature2017_ind$rsIDFetch from SQL server
(Nature_2017_SNPs= sql_1kg %>% filter(SNP %in% snps) %>% collect())## # A tibble: 15 × 35
## CHR SNP A1 A2 ACB ASW BEB CDX CEU CHB
## <int> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 rs10191758 G A 0.1979 0.2049 0.1047 0.0000 0.3687 0.0097
## 2 7 rs10236197 C T 0.5833 0.5000 0.3372 0.4839 0.3586 0.5728
## 3 9 rs11138902 A G 0.2865 0.4016 0.4302 0.3548 0.5859 0.4612
## 4 1 rs12744310 T C 0.1042 0.1066 0.1047 0.4516 0.2172 0.3010
## 5 16 rs12928404 C T 0.6042 0.6066 0.3023 0.5806 0.3737 0.4563
## 6 2 rs13010010 T C 0.0417 0.0492 0.1163 0.4839 0.4091 0.5291
## 7 17 rs16954078 A T 0.0260 0.0492 0.1105 0.0000 0.2020 0.0049
## 8 13 rs2251499 T C 0.1562 0.1803 0.1512 0.1129 0.2727 0.1117
## 9 6 rs2490272 T C 0.1354 0.2213 0.4651 0.5591 0.6465 0.8107
## 10 22 rs36093924 T C 0.5365 0.5492 0.5174 0.2527 0.4444 0.2379
## 11 7 rs4728302 T C 0.2708 0.4508 0.3837 0.6505 0.5505 0.4757
## 12 2 rs6746731 G T 0.2448 0.3361 0.5814 0.2419 0.5909 0.3835
## 13 3 rs6779302 T G 0.4479 0.4180 0.3721 0.0107 0.3636 0.0583
## 14 3 rs7646501 G A 0.7240 0.7213 0.3721 0.1129 0.2374 0.1796
## 15 6 rs9320913 A C 0.1615 0.2295 0.2907 0.4409 0.5000 0.4369
## # ... with 25 more variables: CHS <dbl>, CLM <dbl>, ESN <dbl>, FIN <dbl>,
## # GBR <dbl>, GIH <dbl>, GWD <dbl>, IBS <dbl>, ITU <dbl>, JPT <dbl>,
## # KHV <dbl>, LWK <dbl>, MSL <dbl>, MXL <dbl>, PEL <dbl>, PJL <dbl>,
## # PUR <dbl>, STU <dbl>, TSI <dbl>, YRI <dbl>, EAS <dbl>, EUR <dbl>,
## # AFR <dbl>, AMR <dbl>, SAS <dbl>Merge file with downloaded frequencies to original input file with extra info (z scores, p value, etc)
df_merged_Nature_top=merge(Nature_2017_SNPs,Nature2017_ind,by.x="SNP",by.y="rsID")
write.csv(df_merged_Nature_top,"df_merged_Nature_top.csv")Unfortunately 3 variants were missing (rs66495454,rs113315451,rs41352752). Two indels and 1 SNP.Had to add them manually by using 1000 Genomes Ensembl browswer. Manually reversed the allele when z score was negative to get only positive alleles. Load updated file “df_merged_Nature_top_positive”
df_merged_Nature_top_positive=read.csv(file.choose(),header=TRUE)
df_pol_top=df_merged_Nature_top_positive[,1:26]Calculate polygenic scores
PS_top=apply(df_pol_top,2,mean)
PS_top## ACB ASW BEB CDX CEU CHB CHS
## 0.4753722 0.4690556 0.4312000 0.4925333 0.4747000 0.5264056 0.5071167
## CLM ESN FIN GBR GIH GWD IBS
## 0.4622167 0.4685500 0.4621556 0.4654611 0.4706389 0.4717833 0.4535333
## ITU JPT KHV LWK MSL MXL PEL
## 0.4305722 0.5315111 0.5073278 0.4567889 0.4885778 0.4422722 0.4467167
## PJL PUR STU TSI YRI
## 0.4433056 0.4551167 0.4439056 0.4579222 0.4704444Apply weights (z scores)
df_weighted_top=(df_merged_Nature_top_positive[,1:26])*(abs(df_merged_Nature_top_positive$z))
PS_top_w=apply(df_weighted_top,2,mean)
PS_top_w## ACB ASW BEB CDX CEU CHB CHS CLM
## 2.791038 2.760445 2.563153 2.922041 2.836348 3.149946 3.019268 2.760920
## ESN FIN GBR GIH GWD IBS ITU JPT
## 2.750451 2.754825 2.775733 2.799647 2.767001 2.700084 2.549037 3.160909
## KHV LWK MSL MXL PEL PJL PUR STU
## 3.014528 2.686187 2.872146 2.630699 2.655750 2.635490 2.715050 2.627435
## TSI YRI
## 2.727961 2.760959Sniekers et al. (2017) also identified 9 SNPs that reached significance both in their study and in the Okbay educational attainment GWAS. As this constitutes a quasi-replication, these SNPs are more likely to be real (or closer to) causal variants. Load file “Int_EA_Replicated.csv”
df_int_EA=read.csv(file.choose(),header=TRUE)
PS_int_EA=apply(df_int_EA[,2:27],2,mean)
PS_int_EA## ACB ASW BEB CDX CEU CHB CHS
## 0.5197333 0.5236667 0.5445667 0.6320667 0.5813111 0.6272889 0.6381778
## CLM ESN FIN GBR GIH GWD IBS
## 0.5555667 0.4988889 0.5768667 0.5689667 0.5378556 0.4970889 0.5596778
## ITU JPT KHV LWK MSL MXL PEL
## 0.5522444 0.6255556 0.6274222 0.4909889 0.5267778 0.5441889 0.5607667
## PJL PUR STU TSI YRI
## 0.5428000 0.5416667 0.5293667 0.5617000 0.5011111Compute correlation between unweighted, weighted PS, Piffer-estimated PGS, Population IQ (Piffer 2015). PGS were calculated using Okbay et al. (2016), Davies et al. (2016) and Rietveld et al. (2013) meta-analyses of educational attainment. For more details see Piffer (2017)
IQ=c(83,85,81,NA,99,105,105,83.5,71,101,100,NA,62,97,NA,105,99.4,74,64,88,85,84,83.5,79,99,71)
Piffer_2017_9= c(0.323,0.355,0.401,0.485,0.428,0.553,0.535,0.394,0.334,0.477,0.469,0.429,
0.328,0.438,0.412,0.530,0.517,0.333,0.341,0.385,0.353,0.419,0.394,0.404,0.439,
0.340)
Piffer_2017_Okbay_162SNPs=c(0.485,0.483,0.503,0.516,0.503,0.529,0.522, 0.498,0.479,0.516,0.506,
0.506,0.484,0.513,0.504,0.522,0.522,0.477,0.473, 0.495,0.474,0.505,
0.498,0.501,0.513,0.482)
df_PS=data.frame(PS_top,PS_top_w,PS_int_EA,IQ,Piffer_2017_9,Piffer_2017_Okbay_162SNPs)IQ has NAs, hence use pairwise deletion
cor(df_PS,use="pairwise.complete.obs")## PS_top PS_top_w PS_int_EA IQ
## PS_top 1.0000000 0.9936323 0.6136411 0.4088755
## PS_top_w 0.9936323 1.0000000 0.6814205 0.4966647
## PS_int_EA 0.6136411 0.6814205 1.0000000 0.8766974
## IQ 0.4088755 0.4966647 0.8766974 1.0000000
## Piffer_2017_9 0.5731021 0.6460421 0.9241116 0.8797065
## Piffer_2017_Okbay_162SNPs 0.4218000 0.4979507 0.8351598 0.8543011
## Piffer_2017_9 Piffer_2017_Okbay_162SNPs
## PS_top 0.5731021 0.4218000
## PS_top_w 0.6460421 0.4979507
## PS_int_EA 0.9241116 0.8351598
## IQ 0.8797065 0.8543011
## Piffer_2017_9 1.0000000 0.9500431
## Piffer_2017_Okbay_162SNPs 0.9500431 1.0000000Correlation between the intelligence PS and IQ,PS_Piffer2017,PS_Piffer_2017_162SNPs are r=0.496, 0.646, 0.497. Correlations between the intelligence-EA PS and IQ, PS_Piffer2017_9, PS_Piffer_2016_162SNPs are r= 0.877, 0.924, 0.835 These are lower than the correlations previously observed (Piffer, 2017)
Run Monte Carlo simulation using random SNPs First to calculate Monte Carlo p value for correlation with IQ Load random SNPs (N=13,130). “Matched_18_freqs.csv”
sql_random_freqs <- read.csv(file.choose(),header=TRUE)Keep only freqs.Transpose
t_random_freqs=t(sql_random_freqs[6:31])Compute non-overlapping polygenic scores of same size (18)
byapply <- function(t_random_freqs, by, fun, ...)
{
# Create index list
if (length(by) == 1)
{
nc <- ncol(t_random_freqs)
split.index <- rep(1:ceiling(nc / by), each = by, length.out = nc)
} else # 'by' is a vector of groups
{
nc <- length(by)
split.index <- by
}
index.list <- split(seq(from = 1, to = nc), split.index)
# Pass index list to fun using sapply() and return object
sapply(index.list, function(i)
{
do.call(fun, list(t_random_freqs[, i], ...))
})
}
# Run function
y <- byapply(t_random_freqs, 18, rowMeans)Compute corr between IQ and PGS for random sets of 18 SNPs (N=730)
cor_random18_IQ=cor(y,IQ,use="complete.obs")Calculate Monte Carlo p value (applying formula with correction) for weighted PS
r_w=length(which(cor_random18_IQ>=0.496))
n_w=length(cor_random18_IQ)
MC_p_w=(r_w+1)/(n_w+1)
MC_p_w## [1] 0.3816689The MC p value is 0.382. This means that the same results (r x IQ > 0.496) can be achieved with random SNPs roughly 38.2% of the times.
Create correlation matrix of random and real (PS 9) polygenic scores to obtain a Monte Carlo p value for the correlation between the 9 SNPs EA Polygenic score and the new Nature 2017 polygenic score.
cor_random18_PS_9=cor(y,Piffer_2017_9,use="complete.obs")Calculate p value
r_18_9=length(which(cor_random18_PS_9>=0.646))
n_18_9=length(cor_random18_PS_9)
MC_p_18_9=(r_18_9+1)/(n_18_9+1)
MC_p_18_9## [1] 0.120383Calculate PS by continent
df_weighted_top_continent=(df_merged_Nature_top_positive[,27:31])*(abs(df_merged_Nature_top_positive$z))
PS_top_cont=apply(df_weighted_top_continent,2,mean)
PS_top_cont## EAS EUR AFR AMR SAS
## 2.913380 2.736141 2.781131 2.720869 2.624233The MC p value is 0.120. This means that the same results (r x PS_9 > 0.646) can be achieved with random SNPs roughly 12% of the times. In general, the only robust finding of this study is an East Asian (EAS) advantage but the scores for the other populations are largely overlapping.
Run MC simulation for the 9 intelligence-EA pseudo-replicated SNPs.
byapply <- function(t_random_freqs, by, fun, ...)
{
# Create index list
if (length(by) == 1)
{
nc <- ncol(t_random_freqs)
split.index <- rep(1:ceiling(nc / by), each = by, length.out = nc)
} else # 'by' is a vector of groups
{
nc <- length(by)
split.index <- by
}
index.list <- split(seq(from = 1, to = nc), split.index)
# Pass index list to fun using sapply() and return object
sapply(index.list, function(i)
{
do.call(fun, list(t_random_freqs[, i], ...))
})
}
# Run function
y_9 <- byapply(t_random_freqs, 9, rowMeans)Compute corr between random sets of 9 SNPs and IQ
cor_random9_IQ=cor(y_9,IQ,use="complete.obs")Calculate p value
r_9_9=length(which(cor_random9_IQ>=0.876))
n_9_9=length(cor_random9_IQ)
MC_p_9_9=(r_9_9+1)/(n_9_9+1)
MC_p_9_9## [1] 0.00890411The MC p value for the 9 intelligence-EA quasi-replicated corr x IQ PGS is 0.009.
Factor analysis was used to extrapolate a signal of directional polygenic selection for the 18 intelligence SNPs(Piffer, 2015)
library(psych)
fa_int=fa(t(df_pol_top))
fa_int$scores## MR1
## ACB -1.27631875
## ASW -0.96102733
## BEB -0.07533494
## CDX 1.35058189
## CEU 0.84374466
## CHB 1.10943690
## CHS 1.20766955
## CLM 0.35707546
## ESN -1.65963577
## FIN 0.77070237
## GBR 0.79683553
## GIH -0.04922992
## GWD -1.35832632
## IBS 0.63069832
## ITU -0.07433635
## JPT 0.87816364
## KHV 1.26670030
## LWK -1.59864602
## MSL -1.44376294
## MXL 0.21476285
## PEL -0.05987738
## PJL 0.06571488
## PUR 0.37462562
## STU -0.39063541
## TSI 0.76435902
## YRI -1.68393985The results show clearer differences between populations, compared to the polygenic scores and look similar to the polygenic scores published by Piffer (2015, 2017)
cor(fa_int$scores,df_PS)## PS_top PS_top_w PS_int_EA IQ Piffer_2017_9
## MR1 0.3082242 0.4038636 0.8990939 NA 0.8876893
## Piffer_2017_Okbay_162SNPs
## MR1 0.8799348Compute correlation with IQ (remove NA pairwise)
cor(fa_int$scores,df_PS,use = "pairwise.complete.obs")## PS_top PS_top_w PS_int_EA IQ Piffer_2017_9
## MR1 0.3082242 0.4038636 0.8990939 0.9034576 0.8876893
## Piffer_2017_Okbay_162SNPs
## MR1 0.8799348Repeat for the intelligence-EA SNPs
fa_int_ea=fa(t(df_int_EA[,2:27]))
fa_int_ea$scores## MR1
## ACB -1.0629236351
## ASW -0.9970987866
## BEB -0.6601555423
## CDX 1.2507053503
## CEU 0.7542623611
## CHB 1.3743195824
## CHS 1.6347042174
## CLM -0.1131873487
## ESN -1.2547291601
## FIN 0.5814735810
## GBR 0.7815198753
## GIH -0.0005317457
## GWD -1.1861165277
## IBS 0.4755995807
## ITU -0.2120756729
## JPT 1.3210159053
## KHV 1.9253948560
## LWK -1.2545719258
## MSL -1.1648724693
## MXL -0.2591111303
## PEL -0.7621497853
## PJL 0.0347293802
## PUR -0.2075127835
## STU -0.4320448257
## TSI 0.6767470955
## YRI -1.2433904462cor(fa_int_ea$scores,df_PS)## PS_top PS_top_w PS_int_EA IQ Piffer_2017_9
## MR1 0.5516593 0.626515 0.9401761 NA 0.9618424
## Piffer_2017_Okbay_162SNPs
## MR1 0.9236177cor(fa_int_ea$scores,df_PS,use = "pairwise.complete.obs")## PS_top PS_top_w PS_int_EA IQ Piffer_2017_9
## MR1 0.5516593 0.626515 0.9401761 0.9061601 0.9618424
## Piffer_2017_Okbay_162SNPs
## MR1 0.9236177Indeed, as we can see above, the factor scores have very high correlations with population IQ (r=0.9) and the 162 SNPs polygenic score from Okbay et.al (2017) and the SNPs polygenic score made from the 9 pseudo-repllicated SNPs (Piffer, 2017). Thus, in this instance factor analysis outperforms simply averging allele frequencies.
Show factor loadings
fa_int$loadings##
## Loadings:
## MR1
## [1,]
## [2,] -0.382
## [3,] 0.710
## [4,] -0.778
## [5,] -0.329
## [6,] 0.929
## [7,] -0.396
## [8,] 0.238
## [9,] 0.921
## [10,] 0.568
## [11,] -0.713
## [12,] -0.422
## [13,] 0.829
## [14,] -0.977
## [15,] 0.840
## [16,] 0.726
## [17,] -0.834
## [18,] 0.310
##
## MR1
## SS loadings 7.965
## Proportion Var 0.443Calculate mean fa loadings
mean(fa_int$loadings)## [1] 0.06868852Repeat for the int-ea SNPs
fa_int_ea$loadings##
## Loadings:
## MR1
## [1,] -0.358
## [2,] -0.883
## [3,] 0.994
## [4,] 0.156
## [5,] 0.305
## [6,] 0.875
## [7,] 0.789
## [8,] 0.475
## [9,] -0.510
##
## MR1
## SS loadings 3.888
## Proportion Var 0.432mean(fa_int_ea$loadings)## [1] 0.2047414Finally, a Monte Carlo simulation with the random (N=730) factor scores. First carry out fa over groups of 18 SNPs…see first group output as test
ind=seq(1,(ncol(t_random_freqs)-18),by=18)
results_fa=lapply(ind, function(x) fa(t_random_freqs[, x + 0:17]))
results_fa[[1]]$scores## MR1
## ACB 1.11569074
## ASW 0.53001953
## BEB -0.52523418
## CDX -1.55733116
## CEU 0.47450740
## CHB -1.95284548
## CHS -1.80867302
## CLM 0.28743412
## ESN 0.67975332
## FIN 0.15299121
## GBR 0.60614566
## GIH -0.17495926
## GWD 1.03746149
## IBS 0.92244653
## ITU 0.13683089
## JPT -1.84459552
## KHV -1.81687202
## LWK 0.79571481
## MSL 1.33625037
## MXL 0.07936017
## PEL 0.31384984
## PJL -0.24407216
## PUR 0.39328959
## STU -0.23122284
## TSI 0.35191022
## YRI 0.94214976Save factor scores
results_fa_scores=lapply(results_fa, "[[" , "scores" )Convert list to df and transpose
df_fa_scores=data.frame(lapply(data.frame(t(sapply(results_fa_scores, `[`))), unlist))
df_fa_scores=t(df_fa_scores)Compute correlation between population IQ and random factor scores
cor_fa_iq=cor(df_fa_scores,IQ,use = "pairwise.complete.obs")Compute Monte Carlo p value for corr x population IQ
r_iq=length(which(cor_fa_iq>=0.9))
n_iq=length(cor_fa_iq)
MC_p_iq=(r_iq+1)/(n_iq+1)
MC_p_iq## [1] 0.004109589Repeat for intelligence-EA factor
ind=seq(1,(ncol(t_random_freqs)-9),by=9)
results_fa_int_ea=lapply(ind, function(x) fa(t_random_freqs[, x + 0:8]))
results_fa_scores_int_ea=lapply(results_fa_int_ea, "[[" , "scores" )
df_fa_scores_int_ea=data.frame(lapply(data.frame(t(sapply(results_fa_scores_int_ea, `[`))), unlist))
df_fa_scores_int_ea=t(df_fa_scores_int_ea)
cor_fa_iq_int_ea=cor(df_fa_scores_int_ea,IQ,use = "pairwise.complete.obs")Compute Monte Carlo p value for corr x population IQ
r_iq_int_ea=length(which(cor_fa_iq_int_ea>=0.9))
n_iq_int_ea=length(cor_fa_iq_int_ea)
MC_p_iq_int_ea=(r_iq_int_ea+1)/(n_iq_int_ea+1)
MC_p_iq_int_ea## [1] 0.002056203Compute correlation between PGS Okbay et al. (N=162) and random factor scores
cor_fa_Okbay=cor(df_fa_scores,Piffer_2017_Okbay_162SNPs,use = "pairwise.complete.obs")Compute Monte Carlo p value
r_Okbay=length(which(cor_fa_Okbay>=0.879))
n_Okbay=length(cor_fa_Okbay)
MC_p_Okbay=(r_Okbay+1)/(n_Okbay+1)
MC_p_Okbay## [1] 0.02191781Do the same as above for the 9 EA SNPs PS (Piffer, 2017)
cor_fa_9=cor(df_fa_scores,Piffer_2017_9,use = "pairwise.complete.obs")
r_9=length(which(cor_fa_9>=0.887))
n_9=length(cor_fa_9)
MC_p_9=(r_9+1)/(n_9+1)
MC_p_9## [1] 0.01643836The Monte Carlo simulation showed that the correlation between the factor score obtained from the intelligence GWAS and population IQ, and the previous educational attainment GWAS polygenic scores, is higher than the vast majority of the correlations produced stochastically (i.e. correlations between factor scores of random SNPs and population IQ and ed. att. GWAS PS are almost always lower). Specifically, the empirical (MC) p values are very low (0.004,0.022, 0.016 for population IQ, Okbay 162 PS and Piffer 2017 PS, respectively).
Calculate 9 SNPs EA factor score (load “Replication_linked_Okbay.csv”“)
Replication_linked_Okbay <- read.csv(file.choose(),header=TRUE)
fa_9_EA=fa(Replication_linked_Okbay[,2:10])Show factor loadings
fa_9_EA$loadings##
## Loadings:
## MR1
## rs1008078 -0.801
## rs11588857 0.831
## rs12987662 0.956
## rs148734725 -0.492
## rs11712056 0.606
## rs62263923 0.863
## rs13294439 0.910
## rs12969294 0.459
## rs9320913 0.717
##
## MR1
## SS loadings 5.153
## Proportion Var 0.573Compute mean fa loadings
mean(fa_9_EA$loadings)## [1] 0.4499548Correlation between 3 factors: 9 SNPs EA factor, Intell factor, Int-EA factor
df_fac=data.frame(fa_9_EA$scores,fa_int$scores,fa_int_ea$scores)
cor(df_fac)## MR1 MR1.1 MR1.2
## MR1 1.0000000 0.9608041 0.9459244
## MR1.1 0.9608041 1.0000000 0.9347101
## MR1.2 0.9459244 0.9347101 1.0000000Carry out Monte Carlo analysis for EA x int correlation
cor_fa_18=cor(df_fa_scores,fa_9_EA$scores,use = "pairwise.complete.obs")
r_18=length(which(cor_fa_18>=0.96))
n_18=length(cor_fa_18)
MC_p_18=(r_18+1)/(n_18+1)
MC_p_18## [1] 0.02328767The correlation between the EA and intelligence factors is 0.96 (MC p= 0.023)
Controlling for spatial autocorrelation. Load file “FstDistancesAlph_2016.csv” EA factor
FstDistances <- read.csv(file.choose(),header=TRUE)
library(caTools)
library(lm.beta)
library(car)##
## Attaching package: 'car'## The following object is masked from 'package:psych':
##
## logit## The following object is masked from 'package:dplyr':
##
## recodecombinations_EA=combs(fa_9_EA$scores, 2)
combinations_int=combs(fa_int$scores,2)
#SAC control 9 EA factor
cnames_EA=provideDimnames(combinations_EA, sep = "", base = list(LETTERS))
cnamesdf_EA=as.data.frame(cnames_EA)
distances_EA9=cnamesdf_EA$A-cnamesdf_EA$B
fa_9_dist=abs(distances_EA9)
FstDistances$fa_9_dist=fa_9_dist
fit_dist_fa_9=lm(FstDistances$IQdist ~FstDistances$Fst_Chr1+FstDistances$fa_9_dist)
summary(fit_dist_fa_9)##
## Call:
## lm(formula = FstDistances$IQdist ~ FstDistances$Fst_Chr1 + FstDistances$fa_9_dist)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.1961 -5.6657 0.0894 5.5857 17.9431
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.2076 0.8819 5.905 1.15e-08 ***
## FstDistances$Fst_Chr1 -26.1310 19.8261 -1.318 0.189
## FstDistances$fa_9_dist 10.3841 1.1716 8.864 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.683 on 250 degrees of freedom
## (72 observations deleted due to missingness)
## Multiple R-squared: 0.5027, Adjusted R-squared: 0.4987
## F-statistic: 126.3 on 2 and 250 DF, p-value: < 2.2e-16betafit_fit_dist_9=lm.beta(fit_dist_fa_9)
structure(betafit_fit_dist_9)##
## Call:
## lm(formula = FstDistances$IQdist ~ FstDistances$Fst_Chr1 + FstDistances$fa_9_dist)
##
## Standardized Coefficients::
## (Intercept) FstDistances$Fst_Chr1 FstDistances$fa_9_dist
## 0.0000000 -0.1207470 0.8120172vif(fit_dist_fa_9)## FstDistances$Fst_Chr1 FstDistances$fa_9_dist
## 4.218947 4.218947SAC control Intelligence factor
cnames_int=provideDimnames(combinations_int, sep = "", base = list(LETTERS))
cnamesdf_int=as.data.frame(cnames_int)
distances_int=cnamesdf_int$A-cnamesdf_int$B
fa_int_dist=abs(distances_int)
FstDistances$fa_int_dist=fa_int_dist
fit_dist_fa_int=lm(FstDistances$IQdist ~FstDistances$Fst_Chr1+FstDistances$fa_int_dist)
summary(fit_dist_fa_int)##
## Call:
## lm(formula = FstDistances$IQdist ~ FstDistances$Fst_Chr1 + FstDistances$fa_int_dist)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.0620 -5.0967 0.1203 4.7051 18.0239
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.5172 0.8315 5.433 1.32e-07 ***
## FstDistances$Fst_Chr1 -23.4298 16.2952 -1.438 0.152
## FstDistances$fa_int_dist 10.7256 0.9639 11.128 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.203 on 250 degrees of freedom
## (72 observations deleted due to missingness)
## Multiple R-squared: 0.5629, Adjusted R-squared: 0.5594
## F-statistic: 161 on 2 and 250 DF, p-value: < 2.2e-16betafit_fit_dist_int=lm.beta(fit_dist_fa_int)
structure(betafit_fit_dist_int)##
## Call:
## lm(formula = FstDistances$IQdist ~ FstDistances$Fst_Chr1 + FstDistances$fa_int_dist)
##
## Standardized Coefficients::
## (Intercept) FstDistances$Fst_Chr1 FstDistances$fa_int_dist
## 0.0000000 -0.1082649 0.8378747vif(fit_dist_fa_int)## FstDistances$Fst_Chr1 FstDistances$fa_int_dist
## 3.242631 3.242631The EA and Intell factor scores predict population IQ also after controlling for SAC (Beta= 0.81 and 0.84, respectively)
References:
Piffer, D. Can We Detect Polygenic Selection on Cognitive Ability Using GWAS Hits? Employing Random SNPs as a Null Model.. Preprints 2017, 2017010127 (doi: 10.20944/preprints201701.0127.v2).
Sniekers et al. Genome-wide association meta-analysis of 78,308 individuals identifies new loci and genes influencing human intelligence, Nature Genetics 2017 (doi: doi:10.1038/ng.3869)