Heritability Analysis

As an application of the variance component model, this note demonstrates the workflow for heritability analysis in genetics, using a sample data set cg10k with 6,670 individuals and 630,860 SNPs. Person IDs and phenotype names are masked for privacy. cg10k.bed, cg10k.bim, and cg10k.fam is a set of Plink files in binary format. cg10k_traits.txt contains 13 phenotypes of the 6,670 individuals.


In [31]:
;ls cg10k.bed cg10k.bim cg10k.fam cg10k_traits.txt


cg10k.bed
cg10k.bim
cg10k.fam
cg10k_traits.txt

Machine information:


In [32]:
versioninfo()


Julia Version 0.6.0
Commit 903644385b (2017-06-19 13:05 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin13.4.0)
  CPU: Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

Read in binary SNP data

We will use the SnpArrays.jl package to read in binary SNP data and compute the empirical kinship matrix. Issue

Pkg.clone("https://github.com/OpenMendel/SnpArrays.jl.git")

within Julia to install the SnpArrays package.


In [33]:
using SnpArrays

In [34]:
# read in genotype data from Plink binary file (~50 secs on my laptop)
@time cg10k = SnpArray("cg10k")


 22.902730 seconds (51.62 k allocations: 1005.845 MiB, 0.11% gc time)
Out[34]:
6670×630860 SnpArrays.SnpArray{2}:
 (false, true)   (false, true)   …  (true, true)    (true, true) 
 (true, true)    (true, true)       (false, true)   (true, false)
 (true, true)    (true, true)       (true, true)    (true, true) 
 (true, true)    (true, true)       (false, true)   (true, true) 
 (true, true)    (true, true)       (true, true)    (false, true)
 (false, true)   (false, true)   …  (true, true)    (true, true) 
 (false, false)  (false, false)     (true, true)    (true, true) 
 (true, true)    (true, true)       (true, true)    (false, true)
 (true, true)    (true, true)       (true, true)    (true, true) 
 (true, true)    (true, true)       (false, true)   (true, true) 
 (true, true)    (true, true)    …  (true, true)    (true, true) 
 (false, true)   (false, true)      (true, true)    (false, true)
 (true, true)    (true, true)       (true, true)    (false, true)
 ⋮                               ⋱                               
 (false, true)   (false, true)      (false, true)   (false, true)
 (false, true)   (false, true)      (false, true)   (true, true) 
 (true, true)    (true, true)    …  (false, true)   (true, true) 
 (false, true)   (false, true)      (true, true)    (false, true)
 (true, true)    (true, true)       (false, true)   (true, true) 
 (true, true)    (true, true)       (false, false)  (false, true)
 (true, true)    (true, true)       (true, true)    (false, true)
 (true, true)    (true, true)    …  (true, true)    (true, true) 
 (true, true)    (true, true)       (false, true)   (true, true) 
 (true, true)    (true, true)       (true, true)    (false, true)
 (false, true)   (false, true)      (true, true)    (true, true) 
 (true, true)    (true, true)       (true, true)    (true, true) 

Summary statistics of SNP data


In [35]:
people, snps = size(cg10k)


Out[35]:
(6670, 630860)

In [36]:
# summary statistics (~50 secs on my laptop)
@time maf, _, missings_by_snp, = summarize(cg10k);


 24

In [37]:
# 5 number summary and average MAF (minor allele frequencies)
quantile(maf, [0.0 .25 .5 .75 1.0]), mean(maf)


Out[37]:
([0.00841726 0.124063 … 0.364253 0.5], 0.24536516625042462)

In [38]:
# Pkg.add("Plots")
# Pkg.add("PyPlot")
using Plots
pyplot()

histogram(maf, xlab = "Minor Allele Frequency (MAF)", label = "MAF")


Out[38]:

In [39]:
# proportion of missing genotypes
sum(missings_by_snp) / length(cg10k)


Out[39]:
0.0013128198764010824

In [40]:
# proportion of rare SNPs with maf < 0.05
countnz(maf .< 0.05) / length(maf)


Out[40]:
0.07228069619249913

Empirical kinship matrix

We estimate empirical kinship based on all SNPs by the genetic relation matrix (GRM). Missing genotypes are imputed on the fly by drawing according to the minor allele frequencies.


In [41]:
# GRM using SNPs with maf > 0.01 (default) (~10 mins on my laptop)
srand(123)
@time Φgrm = grm(cg10k; method = :GRM)


396.943890 seconds (8.43 G allocations: 127.378 GiB, 4.38% gc time)
Out[41]:
6670×6670 Array{Float64,2}:
  0.503024      0.00335505   -0.000120075  …  -5.45185e-5   -0.00278072 
  0.00335505    0.498958     -0.00195952       0.000868471   0.0034285  
 -0.000120075  -0.00195952    0.493828         0.000174648  -0.000381467
  0.000923828  -0.00329169   -0.00194166      -0.00223595   -0.00123508 
 -8.39649e-5   -0.00353358    0.0018709        0.00222858   -0.00171176 
  0.00204208    0.000572952   0.00254025   …   0.000861385   2.99785e-5 
  0.000569323   0.0024786    -0.00185743       0.00117649   -0.00118027 
 -0.000642144   0.00317992   -0.00099777       0.00354182   -0.000260645
 -0.00102913   -0.00123475   -0.00061138       0.00173885    0.00177727 
 -0.00139442    0.00208423    0.000124525     -0.00145156   -0.001011   
 -0.00204555    0.00011055   -0.000419398  …  -0.000198235  -0.00110353 
  0.000947587   0.00167346    0.00184451      -0.000690143  -0.00304087 
  0.000322759  -0.000899805   0.00303981       0.000739331  -0.00118835 
  ⋮                                        ⋱                            
  0.00298012    0.00130003    0.000998861      4.18454e-6    0.00303991 
 -0.00207748    0.00274717   -0.00191741      -0.00107073    0.00368267 
  0.000545569  -0.00244439   -0.00299578   …  -0.000669885   0.00221027 
 -0.00423186   -0.00208514   -0.00108833      -0.000622127  -0.000567483
 -0.00325644   -0.000781353   0.0030423        0.000501423  -0.00010267 
  0.00041055   -0.00200772    0.00274867      -0.00624933   -0.00521365 
  0.00210519    0.000879889  -0.00107817      -0.000797878  -0.000557352
 -0.00230058   -0.000119132   0.000116817  …   0.000867087  -0.00233512 
 -0.0020119     0.00230772   -0.00128837       0.00194798   -0.00048733 
 -0.000944942  -0.000928073  -0.000175096      0.00126911   -0.00303766 
 -5.45185e-5    0.000868471   0.000174648      0.500829      0.000469478
 -0.00278072    0.0034285    -0.000381467      0.000469478   0.500627   

Phenotypes

Read in the phenotype data and compute descriptive statistics.


In [42]:
# Pkg.add("DataFrames")
using DataFrames

cg10k_trait = readtable(
    "cg10k_traits.txt"; 
    separator = ' ',
    names = [:FID; :IID; :Trait1; :Trait2; :Trait3; :Trait4; :Trait5; :Trait6; 
             :Trait7; :Trait8; :Trait9; :Trait10; :Trait11; :Trait12; :Trait13],  
    eltypes = [String; String; Float64; Float64; Float64; Float64; Float64; 
               Float64; Float64; Float64; Float64; Float64; Float64; Float64; Float64]
    )
# do not display FID and IID for privacy
cg10k_trait[:, 3:end]


Out[42]:
Trait1Trait2Trait3Trait4Trait5Trait6Trait7Trait8Trait9Trait10Trait11Trait12Trait13
1-1.81573145026234-0.946150461472831.11363077580442-2.098671211191590.7444166141117480.001391718840801310.934732480409667-1.226773154181031.1160784277875-0.44362803350290.824465656443384-1.02852542216546-0.394049201727681
2-1.244400943787290.1096599925471790.467119394241789-1.621313040975891.05667583556830.9789469794191811.000146339460470.324874271402281.162321752196962.69227069487053.082636724610471.090649547860130.0256616415357438
31.455669145023051.538669329232431.094029593765550.586655272226893-0.32796454430367-0.30337709778827-0.0334354881314741-0.464463064285437-0.3319396273436-0.486839089635991-1.10648681564373-1.42015780427231-0.687463456644413
4-0.7688092766985480.5134908855142490.244263028382142-1.317402544756911.193937743268451.173441277342881.087374266752320.5360225837322610.8027592407620680.2341594117498150.394174866891074-0.7673658924760290.0635385761884935
5-0.264415132547719-0.348240421825694-0.02390650834136060.004739158022449481.256191917121931.20388836676311.298007390426270.3101136602473110.6261598610593520.8992891298312240.549967833508120.5406878095420480.179675416046033
6-1.37617270917293-1.471919677445640.291179894254146-0.803110740704731-0.264239977442213-0.260573027836772-0.165372266287781-0.2192572941183621.04702422290318-0.09858155346164820.9473934380684480.5940148120314380.245407436348479
70.1009416296374-0.191615722103455-0.5674213215966770.378571487240382-0.246656179817904-0.6088107500538580.189081058215596-1.27077787326519-0.4524761991439650.7025628772977240.3326362189571790.00269165036261810.317117176705358
8-0.3198182763674641.357744806572830.818689545938528-1.155655316443520.634483681022590.2914619086346790.933323714954726-0.7410832896824920.647477683507572-0.9708776270779660.2208611654113040.852512250237764-0.225904624283945
9-0.2883341733420320.5660825380907520.254958336116175-0.6525783028697140.6689215592773470.9783091991705580.1228629660419381.47909263782140.06721324241734490.07959039175278270.1675324552432320.2469155794421390.539932616458363
10-1.15759732583991-0.781198583545165-0.595807759833517-1.005549802604020.7898288859333210.5710584133790440.951304176233755-0.2959629829848160.990420024797070.5613093669889830.733100030623233-1.73467772245684-1.35278484330654
110.7405691504590311.408738467554150.7346899994400880.0208322841295094-0.337440968561619-0.458304040611395-0.142582512772326-0.580392297464107-0.684684998101516-0.00785381461893456-0.712244337518008-0.313345561230878-0.345419463162219
12-0.6758924864549950.2798926138296820.267915996308248-1.041036653929850.9107417156458880.8660276185131711.074144317020050.03817510035383020.766355377018601-0.340118016143495-0.8090139585050590.548521663785885-0.0201828675962336
13-0.795410435603455-0.6999899397627380.3991295030063-0.5104762619007361.515522454168441.287430329394671.537723932509030.1339891601177021.020257368860370.499018733899186-0.36948273277931-1.10153460436318-0.598132438886619
14-0.193483122930324-0.286021160323518-0.6914942252629950.01315816787006991.523374706867821.40106380722621.531146204518960.3330664834780751.043724803810990.163206783570466-0.422883765001728-0.383527976713573-0.489221907788158
150.1512462033797182.091851089936142.03800472474384-1.124747171435311.665570243907131.625356751095761.587510704836550.6358521860437760.8425777846059790.450761870778952-1.39479033623028-0.5609841075677680.289349776549287
16-0.4646087408127120.361276947723031.2327673928287-0.8260337310863831.434752247099831.744518238188460.2110968874846382.648164251405481.025114331460960.119757316031840.0596832073448267-0.631231612661616-0.207878671782927
17-0.732977488012215-0.5262234258897790.61657871336593-0.554479743325930.9474848590251040.9368332141381730.9725168063355240.2902510138652271.012853597257230.516207422283291-0.03006891719881940.87873225245830.450254629309513
18-0.1673264596221190.1753271654872370.287467725892572-0.4026525320842460.5511815094180560.5222047432909750.4368376600946530.2995649338455790.583109520896067-0.704415820005353-0.730810367994577-1.95140580379896-0.933504665700164
191.411594857874181.787224079010170.843976395853640.481278083772991-0.0887673728508268-0.499577574268580.304195897924847-1.23884208383369-0.153475724036624-0.8704861027883290.0955473331150403-0.983708050882817-0.3563445644514
20-1.42997091652825-0.4901470450342130.272730237607695-1.610299929541530.9907878171977480.7116875326081841.1885836012715-0.3712291880756381.24703459239952-0.03891623322715160.8834957490728722.589880263210173.33539552370368
21-0.1472472881767650.123284304156520.617549051912237-0.187130771782620.2564381075866940.177949837350830.412611806463263-0.2448091245597370.09476248061364920.723017223849532-0.6839483546334360.0873751276309269-0.262209652750371
22-0.187112676773894-0.270777264595619-1.015568185516060.06028505686002330.2724197577579780.869133161879197-0.6575194614142342.32388522018189-0.9999360115250341.446718441783060.971157886040772-0.358747904241515-0.439657942096136
23-1.82434047163768-0.9334804460680671.29474003766977-1.945452211510360.335846511896540.3592016543028440.513652924365886-0.0731976966969581.571390428120051.533293713267281.820768218595282.227403018678291.50063347195857
24-2.29344084351335-2.491618423444180.40383988742336-2.364880747529481.41052548319561.422441171477921.170241662721720.844766501768551.790268754324950.648181858970515-0.0857231057403538-1.027895352926170.491288088952859
25-0.4341359328883050.7408819890346520.699576357578518-1.024055431877750.7595292239837130.9566561108952880.6332995686565890.7707339322685160.8249885117145261.842874376347691.91045942063443-0.5023172078693660.132670133448219
26-2.1920969546557-2.494656642722710.354854763893431-1.931558486357140.9419794002899380.9789171014141060.8948600972897360.4632394028318731.125371333171631.705284461919550.7177927144791230.6458880491082610.783968250169388
27-1.46602269088422-1.249216771018970.307977693653039-1.550973646609890.6189084944747980.6625081716620420.4759571739060780.4847186745977070.4015648920282490.55987973254026-0.376938143754217-0.9339826292282180.390013151672955
28-1.83317744236881-1.532687878287012.55674262685865-1.518277457838350.7894096017464550.9087477997285880.6499719229414790.6683736499316671.200583035199030.2779632560756371.25049531982753.313704450716382.22035828885342
29-0.7845466282431780.2765825795439313.01104958800057-1.119788432067580.9208238584227070.7502176898861511.26153730009639-0.4033638829224170.400667296857811-0.217597941303479-0.724669537565068-0.391945338467193-0.650023936358253
300.4644559163451351.3326356122229-1.23059563374303-0.3579759589374141.182497469771041.54315938069757-0.603390411540623.383088459584220.823740765148641-0.129951318508883-0.657979878422938-0.499534924074273-0.414476569095651
&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip

In [43]:
describe(cg10k_trait[:, 3:end])


Trait1
Summary Stats:
Mean:           0.002211
Minimum:        -3.204128
1st Quartile:   -0.645771
Median:         0.125010
3rd Quartile:   0.723315
Maximum:        3.479398
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait2
Summary Stats:
Mean:           0.001353
Minimum:        -3.511659
1st Quartile:   -0.642621
Median:         0.033517
3rd Quartile:   0.657467
Maximum:        4.913423
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait3
Summary Stats:
Mean:           -0.001296
Minimum:        -3.938436
1st Quartile:   -0.640907
Median:         -0.000782
3rd Quartile:   0.637108
Maximum:        7.916299
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait4
Summary Stats:
Mean:           0.002309
Minimum:        -3.608403
1st Quartile:   -0.546086
Median:         0.228165
3rd Quartile:   0.715291
Maximum:        3.127688
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait5
Summary Stats:
Mean:           -0.001790
Minimum:        -4.148749
1st Quartile:   -0.690765
Median:         0.031034
3rd Quartile:   0.734916
Maximum:        2.717184
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait6
Summary Stats:
Mean:           -0.001196
Minimum:        -3.824792
1st Quartile:   -0.662796
Median:         0.036242
3rd Quartile:   0.741176
Maximum:        2.589728
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait7
Summary Stats:
Mean:           -0.001989
Minimum:        -4.272455
1st Quartile:   -0.638923
Median:         0.069801
3rd Quartile:   0.710423
Maximum:        2.653779
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait8
Summary Stats:
Mean:           0.000614
Minimum:        -5.625488
1st Quartile:   -0.601575
Median:         -0.038630
3rd Quartile:   0.527342
Maximum:        5.805702
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait9
Summary Stats:
Mean:           -0.001810
Minimum:        -5.381968
1st Quartile:   -0.601429
Median:         0.106571
3rd Quartile:   0.698567
Maximum:        2.571936
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait10
Summary Stats:
Mean:           -0.000437
Minimum:        -3.548506
1st Quartile:   -0.633641
Median:         -0.096651
3rd Quartile:   0.498610
Maximum:        6.537820
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait11
Summary Stats:
Mean:           -0.000616
Minimum:        -3.264910
1st Quartile:   -0.673685
Median:         -0.068044
3rd Quartile:   0.655486
Maximum:        4.262410
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait12
Summary Stats:
Mean:           -0.000589
Minimum:        -8.851909
1st Quartile:   -0.539686
Median:         -0.141099
3rd Quartile:   0.350779
Maximum:        13.211402
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000

Trait13
Summary Stats:
Mean:           -0.000151
Minimum:        -5.592104
1st Quartile:   -0.492289
Median:         -0.141022
3rd Quartile:   0.324804
Maximum:        24.174436
Length:         6670
Type:           Float64
Number Missing: 0
% Missing:      0.000000


In [44]:
Y = convert(Matrix{Float64}, cg10k_trait[:, 3:15])
histogram(Y, layout = 13)


Out[44]:

Pre-processing data for heritability analysis

To prepare variance component model fitting, we form an instance of VarianceComponentVariate. The two variance components are $(2\Phi, I)$.


In [45]:
using VarianceComponentModels

# form data as VarianceComponentVariate
cg10kdata = VarianceComponentVariate(Y, (2Φgrm, eye(size(Y, 1))))
fieldnames(cg10kdata)


Out[45]:
3-element Array{Symbol,1}:
 :Y
 :X
 :V

In [46]:
cg10kdata


Out[46]:
VarianceComponentModels.VarianceComponentVariate{Float64,2,Array{Float64,2},Array{Float64,2},Array{Float64,2}}([-1.81573 -0.94615 … -1.02853 -0.394049; -1.2444 0.10966 … 1.09065 0.0256616; … ; 0.886626 0.487408 … -0.636874 -0.439825; -1.24394 0.213697 … 0.299931 0.392809], Array{Float64}(6670,0), ([1.00605 0.00671009 … -0.000109037 -0.00556144; 0.00671009 0.997916 … 0.00173694 0.00685701; … ; -0.000109037 0.00173694 … 1.00166 0.000938955; -0.00556144 0.00685701 … 0.000938955 1.00125], [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0]))

Before fitting the variance component model, we pre-compute the eigen-decomposition of $2\Phi_{\text{GRM}}$, the rotated responses, and the constant part in log-likelihood, and store them as a TwoVarCompVariateRotate instance, which is re-used in various variane component estimation procedures.


In [47]:
# pre-compute eigen-decomposition (~50 secs on my laptop)
@time cg10kdata_rotated = TwoVarCompVariateRotate(cg10kdata)
fieldnames(cg10kdata_rotated)


 48.837361 seconds (39 allocations: 1021.427 MiB, 0.57% gc time)
Out[47]:
5-element Array{Symbol,1}:
 :Yrot    
 :Xrot    
 :eigval  
 :eigvec  
 :logdetV2

Save intermediate results

We don't want to re-compute SnpArray and empirical kinship matrices again and again for heritibility analysis.


In [48]:
# # Pkg.add("JLD")
# using JLD
# @save "cg10k.jld"
# whos()

To load workspace


In [49]:
using SnpArrays, JLD, DataFrames, VarianceComponentModels, Plots
pyplot()
@load "cg10k.jld"
whos()


                          Base               Module
                       BinDeps  41348 KB     Module
                         Blosc  41202 KB     Module
                    ColorTypes  41457 KB     Module
                        Colors  41480 KB     Module
                        Compat  41196 KB     Module
                         Conda  41205 KB     Module
                          Core               Module
                    DataArrays  41456 KB     Module
                    DataFrames  41684 KB     Module
                DataStructures  41356 KB     Module
                        FileIO  41310 KB     Module
             FixedPointNumbers  41695 KB     Module
                          GZip  41181 KB     Module
                          HDF5  41403 KB     Module
                        IJulia 4185781 KB     Module
                         Ipopt  41172 KB     Module
                           JLD  41376 KB     Module
                          JSON  41245 KB     Module
                  LaTeXStrings   4058 bytes  Module
                 LegacyStrings  41212 KB     Module
                    LinearMaps     22 KB     Module
                    MacroTools  41606 KB     Module
                          Main               Module
                  MathProgBase  41353 KB     Module
                       MbedTLS  41269 KB     Module
                      Measures  41175 KB     Module
                       NaNMath  41200 KB     Module
                    PlotThemes  41167 KB     Module
                     PlotUtils  41332 KB     Module
                         Plots  42960 KB     Module
                        PyCall  41711 KB     Module
                        PyPlot  41771 KB     Module
                   RecipesBase  41283 KB     Module
                      Reexport  41160 KB     Module
                      Requires  41172 KB     Module
                           SHA     62 KB     Module
                       Showoff  41163 KB     Module
                     SnpArrays  41218 KB     Module
             SortingAlgorithms  41178 KB     Module
              SpecialFunctions  41252 KB     Module
                  StaticArrays  41744 KB     Module
                     StatsBase  41810 KB     Module
                     URIParser  41171 KB     Module
       VarianceComponentModels  41278 KB     Module
                             Y    677 KB     6670×13 Array{Float64,2}
                           ZMQ  41223 KB     Module
                             _     77 KB     630860-element BitArray{1}
                         cg10k 1027303 KB     6670×630860 SnpArrays.SnpArray{2}
                   cg10k_trait    978 KB     6670×15 DataFrames.DataFrame
                     cg10kdata 695816 KB     VarianceComponentModels.VarianceCo…
             cg10kdata_rotated 348299 KB     VarianceComponentModels.TwoVarComp…
                             h     24 bytes  3-element Array{Float64,1}
                           hST    104 bytes  13-element Array{Float64,1}
                        hST_se    104 bytes  13-element Array{Float64,1}
                           hse     24 bytes  3-element Array{Float64,1}
                           maf   4928 KB     630860-element Array{Float64,1}
               missings_by_snp   4928 KB     630860-element Array{Int64,1}
                        people      8 bytes  Int64
                          snps      8 bytes  Int64
                  trait57_data 347778 KB     VarianceComponentModels.TwoVarComp…
                 trait57_model    232 bytes  VarianceComponentModels.VarianceCo…
                traitall_model   2792 bytes  VarianceComponentModels.VarianceCo…
                      traitidx     16 bytes  3-element UnitRange{Int64}
                            Σa   3848 bytes  13×13 Array{Array{Float64,2},2}
                          Σcov   2592 bytes  18×18 Array{Float64,2}
                            Σe   3848 bytes  13×13 Array{Array{Float64,2},2}
                          Φgrm 347569 KB     6670×6670 Array{Float64,2}
                           σ2a    104 bytes  13-element Array{Float64,1}
                           σ2e    104 bytes  13-element Array{Float64,1}

Heritability of single traits

We use Fisher scoring algorithm to fit variance component model for each single trait.


In [50]:
# heritability from single trait analysis
hST = zeros(13)
# standard errors of estimated heritability
hST_se = zeros(13)
# additive genetic effects
σ2a = zeros(13)
# enviromental effects
σ2e = zeros(13)

tic()
for trait in 1:13
    println(names(cg10k_trait)[trait + 2])
    # form data set for trait j
    traitj_data = TwoVarCompVariateRotate(cg10kdata_rotated.Yrot[:, trait], cg10kdata_rotated.Xrot, 
        cg10kdata_rotated.eigval, cg10kdata_rotated.eigvec, cg10kdata_rotated.logdetV2)
    # initialize model parameters
    traitj_model = VarianceComponentModel(traitj_data)
    # estimate variance components
    _, _, _, Σcov, _, _ = mle_fs!(traitj_model, traitj_data; solver=:Ipopt, verbose=false)
    σ2a[trait] = traitj_model.Σ[1][1]
    σ2e[trait] = traitj_model.Σ[2][1]
    @show σ2a[trait], σ2e[trait]
    h, hse = heritability(traitj_model.Σ, Σcov)
    hST[trait] = h[1]
    hST_se[trait] = hse[1]
end
toc()


Trait1
(σ2a[trait], σ2e[trait]) = (0.25978160614793233, 0.7369535197912689)
Trait2
(σ2a[trait], σ2e[trait]) = (0.18647130348299173, 0.8129591079735827)
Trait3
(σ2a[trait], σ2e[trait]) = (0.3188368159422607, 0.6798809726936244)
Trait4
(σ2a[trait], σ2e[trait]) = (0.2651357653143703, 0.7308007669086968)
Trait5
(σ2a[trait], σ2e[trait]) = (0.28083388108246, 0.7172036435586534)
Trait6
(σ2a[trait], σ2e[trait]) = (0.2824159905728832, 0.7170988773569172)
Trait7
(σ2a[trait], σ2e[trait]) = (0.2155274336968625, 0.7815346282986375)
Trait8
(σ2a[trait], σ2e[trait]) = (0.194687807263945, 0.8049690651320599)
Trait9
(σ2a[trait], σ2e[trait]) = (0.24706855916591713, 0.7512942998567308)
Trait10
(σ2a[trait], σ2e[trait]) = (0.098712236297271, 0.9011756660217387)
Trait11
(σ2a[trait], σ2e[trait]) = (0.1664264642608195, 0.8322427413046204)
Trait12
(σ2a[trait], σ2e[trait]) = (0.0834296761650666, 0.9153609794266608)
Trait13
(σ2a[trait], σ2e[trait]) = (0.05893968504298988, 0.940270012443928)
elapsed time: 0.160999612 seconds
Out[50]:
0.160999612

In [51]:
# heritability and standard errors
[hST'; hST_se']


Out[51]:
2×13 Array{Float64,2}:
 0.260633   0.186578   0.319246   …  0.166648  0.0835307  0.0589863
 0.0799732  0.0869002  0.0741007     0.08862   0.0944407  0.0953238

Pairwise traits

Joint analysis of multiple traits is subject to intensive research recently. Following code snippet does joint analysis of all pairs of traits, a total of 78 bivariate variane component models.


In [52]:
# additive genetic effects (2x2 psd matrices) from bavariate trait analysis;
Σa = Array{Matrix{Float64}}(13, 13)
# environmental effects (2x2 psd matrices) from bavariate trait analysis;
Σe = Array{Matrix{Float64}}(13, 13)

tic()
for i in 1:13
    for j in (i+1):13
        println(names(cg10k_trait)[i + 2], names(cg10k_trait)[j + 2])
        # form data set for (trait1, trait2)
        traitij_data = TwoVarCompVariateRotate(cg10kdata_rotated.Yrot[:, [i;j]], cg10kdata_rotated.Xrot, 
            cg10kdata_rotated.eigval, cg10kdata_rotated.eigvec, cg10kdata_rotated.logdetV2)
        # initialize model parameters
        traitij_model = VarianceComponentModel(traitij_data)
        # estimate variance components
        mle_fs!(traitij_model, traitij_data; solver=:Ipopt, verbose=false)
        Σa[i, j] = traitij_model.Σ[1]
        Σe[i, j] = traitij_model.Σ[2]
        @show Σa[i, j], Σe[i, j]
    end
end
toc()


Trait1Trait2
(Σa[i, j], Σe[i, j]) = ([0.258822 0.174358; 0.174358 0.185108], [0.737892 0.585751; 0.585751 0.814301])
Trait1Trait3
(Σa[i, j], Σe[i, j]) = ([0.260236 -0.0144726; -0.0144726 0.319245], [0.736512 -0.11979; -0.11979 0.679488])
Trait1Trait4
(Σa[i, j], Σe[i, j]) = ([0.259615 0.222203; 0.222203 0.265149], [0.737116 0.599854; 0.599854 0.730788])
Trait1Trait5
(Σa[i, j], Σe[i, j]) = ([0.259574 -0.146827; -0.146827 0.28153], [0.737153 -0.254777; -0.254777 0.71653])
Trait1Trait6
(Σa[i, j], Σe[i, j]) = ([0.259476 -0.129115; -0.129115 0.282688], [0.73725 -0.23161; -0.23161 0.716837])
Trait1Trait7
(Σa[i, j], Σe[i, j]) = ([0.259115 -0.140455; -0.140455 0.215297], [0.737606 -0.197616; -0.197616 0.781774])
Trait1Trait8
(Σa[i, j], Σe[i, j]) = ([0.259778 -0.0327756; -0.0327756 0.194698], [0.736957 -0.127026; -0.127026 0.804959])
Trait1Trait9
(Σa[i, j], Σe[i, j]) = ([0.261858 -0.204589; -0.204589 0.246027], [0.734961 -0.307734; -0.307734 0.75232])
Trait1Trait10
(Σa[i, j], Σe[i, j]) = ([0.259649 -0.0994858; -0.0994858 0.0956585], [0.737083 -0.303942; -0.303942 0.904218])
Trait1Trait11
(Σa[i, j], Σe[i, j]) = ([0.25947 -0.138603; -0.138603 0.164709], [0.737257 -0.359557; -0.359557 0.83395])
Trait1Trait12
(Σa[i, j], Σe[i, j]) = ([0.261779 -0.145414; -0.145414 0.0807748], [0.735076 -0.041823; -0.041823 0.9181])
Trait1Trait13
(Σa[i, j], Σe[i, j]) = ([0.261125 -0.108774; -0.108774 0.0538214], [0.735674 -0.114123; -0.114123 0.945416])
Trait2Trait3
(Σa[i, j], Σe[i, j]) = ([0.186541 0.144056; 0.144056 0.320627], [0.812888 0.0995944; 0.0995944 0.678167])
Trait2Trait4
(Σa[i, j], Σe[i, j]) = ([0.186131 0.0746032; 0.0746032 0.265122], [0.813293 0.221109; 0.221109 0.730814])
Trait2Trait5
(Σa[i, j], Σe[i, j]) = ([0.186442 -0.0118093; -0.0118093 0.280842], [0.812987 -0.0365191; -0.0365191 0.717195])
Trait2Trait6
(Σa[i, j], Σe[i, j]) = ([0.18649 -0.00366533; -0.00366533 0.282471], [0.812941 -0.0206271; -0.0206271 0.717046])
Trait2Trait7
(Σa[i, j], Σe[i, j]) = ([0.186104 -0.030665; -0.030665 0.215304], [0.81332 -0.000667009; -0.000667009 0.781755])
Trait2Trait8
(Σa[i, j], Σe[i, j]) = ([0.187023 0.0331783; 0.0331783 0.195259], [0.812421 -0.0326343; -0.0326343 0.804415])
Trait2Trait9
(Σa[i, j], Σe[i, j]) = ([0.185032 -0.085334; -0.085334 0.245909], [0.814386 -0.0809638; -0.0809638 0.752433])
Trait2Trait10
(Σa[i, j], Σe[i, j]) = ([0.186587 -0.123303; -0.123303 0.0987387], [0.812872 -0.273083; -0.273083 0.901229])
Trait2Trait11
(Σa[i, j], Σe[i, j]) = ([0.185484 -0.117256; -0.117256 0.167776], [0.81393 -0.296772; -0.296772 0.830934])
Trait2Trait12
(Σa[i, j], Σe[i, j]) = ([0.185907 -0.0909104; -0.0909104 0.0827171], [0.813555 0.0457924; 0.0457924 0.916135])
Trait2Trait13
(Σa[i, j], Σe[i, j]) = ([0.185979 -0.0720811; -0.0720811 0.0568238], [0.8135 0.0751703; 0.0751703 0.942424])
Trait3Trait4
(Σa[i, j], Σe[i, j]) = ([0.3188 -0.154562; -0.154562 0.264323], [0.679917 -0.303223; -0.303223 0.731591])
Trait3Trait5
(Σa[i, j], Σe[i, j]) = ([0.319216 0.183527; 0.183527 0.282063], [0.679514 0.33724; 0.33724 0.716008])
Trait3Trait6
(Σa[i, j], Σe[i, j]) = ([0.319776 0.165672; 0.165672 0.284448], [0.678972 0.298667; 0.298667 0.715124])
Trait3Trait7
(Σa[i, j], Σe[i, j]) = ([0.318838 0.166283; 0.166283 0.215261], [0.67988 0.347706; 0.347706 0.781796])
Trait3Trait8
(Σa[i, j], Σe[i, j]) = ([0.320718 0.0566397; 0.0566397 0.197764], [0.678063 0.0451569; 0.0451569 0.801956])
Trait3Trait9
(Σa[i, j], Σe[i, j]) = ([0.319001 0.137699; 0.137699 0.246142], [0.679722 0.266704; 0.266704 0.752197])
Trait3Trait10
(Σa[i, j], Σe[i, j]) = ([0.31908 -0.076513; -0.076513 0.0996001], [0.679646 -0.142905; -0.142905 0.900298])
Trait3Trait11
(Σa[i, j], Σe[i, j]) = ([0.318094 -0.0177494; -0.0177494 0.16629], [0.6806 -0.1144; -0.1144 0.832376])
Trait3Trait12
(Σa[i, j], Σe[i, j]) = ([0.321164 0.0843842; 0.0843842 0.0874609], [0.677639 0.0341558; 0.0341558 0.911368])
Trait3Trait13
(Σa[i, j], Σe[i, j]) = ([0.323273 0.109443; 0.109443 0.0634295], [0.675635 -0.0060525; -0.0060525 0.935819])
Trait4Trait5
(Σa[i, j], Σe[i, j]) = ([0.26525 -0.215125; -0.215125 0.282572], [0.73068 -0.377406; -0.377406 0.715518])
Trait4Trait6
(Σa[i, j], Σe[i, j]) = ([0.265715 -0.199714; -0.199714 0.283942], [0.730231 -0.347732; -0.347732 0.715619])
Trait4Trait7
(Σa[i, j], Σe[i, j]) = ([0.26407 -0.18238; -0.18238 0.214324], [0.731843 -0.32655; -0.32655 0.782733])
Trait4Trait8
(Σa[i, j], Σe[i, j]) = ([0.266229 -0.0965381; -0.0965381 0.196655], [0.729739 -0.151461; -0.151461 0.803044])
Trait4Trait9
(Σa[i, j], Σe[i, j]) = ([0.269627 -0.226931; -0.226931 0.247265], [0.726443 -0.416085; -0.416085 0.751086])
Trait4Trait10
(Σa[i, j], Σe[i, j]) = ([0.265098 -0.0352926; -0.0352926 0.0981462], [0.730847 -0.226248; -0.226248 0.901736])
Trait4Trait11
(Σa[i, j], Σe[i, j]) = ([0.265178 -0.0970634; -0.0970634 0.164885], [0.73076 -0.272291; -0.272291 0.833762])
Trait4Trait12
(Σa[i, j], Σe[i, j]) = ([0.267732 -0.140985; -0.140985 0.081029], [0.728323 -0.0834791; -0.0834791 0.917815])
Trait4Trait13
(Σa[i, j], Σe[i, j]) = ([0.265695 -0.0970238; -0.0970238 0.0564809], [0.730259 -0.226115; -0.226115 0.942736])
Trait5Trait6
(Σa[i, j], Σe[i, j]) = ([0.281198 0.280259; 0.280259 0.281764], [0.716855 0.661013; 0.661013 0.717735])
Trait5Trait7
(Σa[i, j], Σe[i, j]) = ([0.280442 0.231918; 0.231918 0.211837], [0.717597 0.674491; 0.674491 0.785172])
Trait5Trait8
(Σa[i, j], Σe[i, j]) = ([0.280958 0.163168; 0.163168 0.193315], [0.717089 0.221817; 0.221817 0.806314])
Trait5Trait9
(Σa[i, j], Σe[i, j]) = ([0.283544 0.243884; 0.243884 0.240564], [0.714585 0.509072; 0.509072 0.757631])
Trait5Trait10
(Σa[i, j], Σe[i, j]) = ([0.281378 -0.0454427; -0.0454427 0.100081], [0.716678 -0.0579778; -0.0579778 0.899822])
Trait5Trait11
(Σa[i, j], Σe[i, j]) = ([0.280066 0.0195669; 0.0195669 0.165607], [0.71795 -0.0345589; -0.0345589 0.833047])
Trait5Trait12
(Σa[i, j], Σe[i, j]) = ([0.28101 0.0592641; 0.0592641 0.0831831], [0.717036 0.0552788; 0.0552788 0.915608])
Trait5Trait13
(Σa[i, j], Σe[i, j]) = ([0.281854 0.0680641; 0.0680641 0.0591899], [0.716223 0.0551992; 0.0551992 0.940027])
Trait6Trait7
(Σa[i, j], Σe[i, j]) = ([0.282435 0.220236; 0.220236 0.213997], [0.71708 0.581507; 0.581507 0.783041])
Trait6Trait8
(Σa[i, j], Σe[i, j]) = ([0.282435 0.18375; 0.18375 0.192999], [0.717081 0.436932; 0.436932 0.80663])
Trait6Trait9
(Σa[i, j], Σe[i, j]) = ([0.284516 0.233768; 0.233768 0.242478], [0.715071 0.477502; 0.477502 0.755765])
Trait6Trait10
(Σa[i, j], Σe[i, j]) = ([0.283087 -0.0427658; -0.0427658 0.100634], [0.716449 -0.0599491; -0.0599491 0.899275])
Trait6Trait11
(Σa[i, j], Σe[i, j]) = ([0.281046 0.0272144; 0.0272144 0.165044], [0.71843 -0.0516242; -0.0516242 0.833601])
Trait6Trait12
(Σa[i, j], Σe[i, j]) = ([0.28256 0.0548537; 0.0548537 0.083133], [0.716961 0.0502064; 0.0502064 0.915658])
Trait6Trait13
(Σa[i, j], Σe[i, j]) = ([0.283231 0.0585667; 0.0585667 0.0592752], [0.716314 0.055827; 0.055827 0.939942])
Trait7Trait8
(Σa[i, j], Σe[i, j]) = ([0.213998 0.0875641; 0.0875641 0.192993], [0.78304 -0.055939; -0.055939 0.806635])
Trait7Trait9
(Σa[i, j], Σe[i, j]) = ([0.219039 0.216925; 0.216925 0.243338], [0.778156 0.463024; 0.463024 0.754935])
Trait7Trait10
(Σa[i, j], Σe[i, j]) = ([0.216296 -0.0412106; -0.0412106 0.100663], [0.780785 -0.0868086; -0.0868086 0.899246])
Trait7Trait11
(Σa[i, j], Σe[i, j]) = ([0.2142 0.0204227; 0.0204227 0.165077], [0.782833 -0.0478727; -0.0478727 0.833569])
Trait7Trait12
(Σa[i, j], Σe[i, j]) = ([0.215054 0.0738562; 0.0738562 0.0814228], [0.782012 0.0366272; 0.0366272 0.917365])
Trait7Trait13
(Σa[i, j], Σe[i, j]) = ([0.216093 0.0728515; 0.0728515 0.0570272], [0.781006 0.0409945; 0.0409945 0.942189])
Trait8Trait9
(Σa[i, j], Σe[i, j]) = ([0.195154 0.111756; 0.111756 0.246453], [0.804528 0.185842; 0.185842 0.751896])
Trait8Trait10
(Σa[i, j], Σe[i, j]) = ([0.195015 -0.015506; -0.015506 0.0990776], [0.804651 0.0118538; 0.0118538 0.900815])
Trait8Trait11
(Σa[i, j], Σe[i, j]) = ([0.194421 0.0215044; 0.0215044 0.166226], [0.805231 -0.026247; -0.026247 0.83244])
Trait8Trait12
(Σa[i, j], Σe[i, j]) = ([0.194491 -0.00425152; -0.00425152 0.0832711], [0.805162 0.0349872; 0.0349872 0.915518])
Trait8Trait13
(Σa[i, j], Σe[i, j]) = ([0.19448 0.00235501; 0.00235501 0.0589351], [0.805173 0.0396048; 0.0396048 0.940275])
Trait9Trait10
(Σa[i, j], Σe[i, j]) = ([0.246455 -0.00257997; -0.00257997 0.0984563], [0.751895 0.0743439; 0.0743439 0.901429])
Trait9Trait11
(Σa[i, j], Σe[i, j]) = ([0.247001 0.0303415; 0.0303415 0.166421], [0.75136 0.153765; 0.153765 0.832248])
Trait9Trait12
(Σa[i, j], Σe[i, j]) = ([0.249421 0.0829968; 0.0829968 0.0890874], [0.749007 0.109331; 0.109331 0.909778])
Trait9Trait13
(Σa[i, j], Σe[i, j]) = ([0.24861 0.0916799; 0.0916799 0.0602352], [0.749811 0.100027; 0.100027 0.939032])
Trait10Trait11
(Σa[i, j], Σe[i, j]) = ([0.0914658 0.100613; 0.100613 0.166501], [0.908376 0.473847; 0.473847 0.83217])
Trait10Trait12
(Σa[i, j], Σe[i, j]) = ([0.0951392 0.0588424; 0.0588424 0.0796744], [0.904735 0.0828862; 0.0828862 0.919115])
Trait10Trait13
(Σa[i, j], Σe[i, j]) = ([0.0995192 -0.0257171; -0.0257171 0.0598595], [0.900397 0.163778; 0.163778 0.939368])
Trait11Trait12
(Σa[i, j], Σe[i, j]) = ([0.165386 0.0579914; 0.0579914 0.0796005], [0.83327 0.144637; 0.144637 0.919166])
Trait11Trait13
(Σa[i, j], Σe[i, j]) = ([0.166417 -0.000985185; -0.000985185 0.0595681], [0.832265 0.200012; 0.200012 0.939646])
Trait12Trait13
(Σa[i, j], Σe[i, j]) = ([0.085082 0.0696185; 0.0696185 0.0569655], [0.913729 0.572041; 0.572041 0.942247])
elapsed time: 3.587337102 seconds
Out[52]:
3.587337102

3-trait analysis

Researchers want to jointly analyze traits 5-7. Our strategy is to try both Fisher scoring and MM algorithm with different starting point, and choose the best local optimum. We first form the data set and run Fisher scoring, which yields a final objective value -1.4700991+04.


In [53]:
traitidx = 5:7
# form data set
trait57_data = TwoVarCompVariateRotate(cg10kdata_rotated.Yrot[:, traitidx], cg10kdata_rotated.Xrot, 
    cg10kdata_rotated.eigval, cg10kdata_rotated.eigvec, cg10kdata_rotated.logdetV2)
# initialize model parameters
trait57_model = VarianceComponentModel(trait57_data)
# estimate variance components
@time mle_fs!(trait57_model, trait57_data; solver=:Ipopt, verbose=true)
trait57_model


This is Ipopt version 3.12.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       78

Total number of variables............................:       12
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  3.0247565e+04 0.00e+00 1.00e+02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0 
   5  1.6835078e+04 0.00e+00 4.08e+02 -11.0 3.64e-01    -  1.00e+00 1.00e+00f  1 MaxS
  10  1.4742941e+04 0.00e+00 1.10e+02 -11.0 2.35e-01    -  1.00e+00 1.00e+00f  1 MaxS
  15  1.4701394e+04 0.00e+00 1.16e+01 -11.0 7.78e-02  -4.5 1.00e+00 1.00e+00f  1 MaxS
  20  1.4701019e+04 0.00e+00 5.75e-01 -11.0 1.51e-04  -6.9 1.00e+00 1.00e+00f  1 MaxS
  25  1.4701018e+04 0.00e+00 2.40e-02 -11.0 6.38e-06  -9.2 1.00e+00 1.00e+00f  1 MaxS
  30  1.4701018e+04 0.00e+00 9.98e-04 -11.0 2.66e-07 -11.6 1.00e+00 1.00e+00f  1 MaxS
  35  1.4701018e+04 0.00e+00 4.15e-05 -11.0 1.10e-08 -14.0 1.00e+00 1.00e+00h  1 MaxS
  40  1.4701018e+04 0.00e+00 1.72e-06 -11.0 4.59e-10 -16.4 1.00e+00 1.00e+00f  1 MaxSA
  45  1.4701018e+04 0.00e+00 7.17e-08 -11.0 1.91e-11 -18.8 1.00e+00 1.00e+00h  1 MaxSA
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls

Number of Iterations....: 49

                                   (scaled)                 (unscaled)
Objective...............:   4.4720359684330265e+02    1.4701017692082147e+04
Dual infeasibility......:   5.6081357364421780e-09    1.8435742302386474e-07
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   5.6081357364421780e-09    1.8435742302386474e-07


Number of objective function evaluations             = 50
Number of objective gradient evaluations             = 50
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 49
Total CPU secs in IPOPT (w/o function evaluations)   =      0.014
Total CPU secs in NLP function evaluations           =      0.076

EXIT: Optimal Solution Found.
  0.097955 seconds (55.15 k allocations: 5.632 MiB)
Out[53]:
VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}(Array{Float64}(0,3), ([0.280777 0.279441 0.232208; 0.279441 0.28422 0.219831; 0.232208 0.219831 0.212832], [0.717266 0.66183 0.674206; 0.66183 0.715287 0.581891; 0.674206 0.581891 0.784183]), Array{Float64}(0,0), Char[], Float64[], -Inf, Inf)

We then run the MM algorithm, starting from the Fisher scoring answer. MM finds an improved solution with objective value 8.955397e+03.


In [54]:
# trait59_model contains the fitted model by Fisher scoring now
@time mle_mm!(trait57_model, trait57_data; verbose=true)
trait57_model


     MM Algorithm
  Iter      Objective  
--------  -------------
       0  -1.470102e+04
       1  -1.470102e+04

  0.003006 seconds (21.01 k allocations: 1.551 MiB)
Out[54]:
VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}(Array{Float64}(0,3), ([0.280777 0.279441 0.232208; 0.279441 0.28422 0.219831; 0.232208 0.219831 0.212832], [0.717266 0.66183 0.674206; 0.66183 0.715287 0.581891; 0.674206 0.581891 0.784183]), Array{Float64}(0,0), Char[], Float64[], -Inf, Inf)

Do another run of MM algorithm from default starting point. It leads to a slightly better local optimum -1.470104e+04, slighly worse than the Fisher scoring result. Follow up anlaysis should use the Fisher scoring result.


In [55]:
# default starting point
trait57_model = VarianceComponentModel(trait57_data)
@time _, _, _, Σcov, = mle_mm!(trait57_model, trait57_data; verbose=true)
trait57_model


     MM Algorithm
  Iter      Objective  
--------  -------------
       0  -3.024757e+04
       1  -2.040300e+04
       2  -1.656070e+04
       3  -1.528529e+04
       4  -1.490986e+04
       5  -1.480638e+04
       6  -1.477811e+04
       7  -1.476968e+04
       8  -1.476639e+04
       9  -1.476444e+04
      10  -1.476286e+04
      20  -1.475000e+04
      30  -1.474011e+04
      40  -1.473248e+04
      50  -1.472658e+04
      60  -1.472200e+04
      70  -1.471840e+04
      80  -1.471555e+04
      90  -1.471328e+04
     100  -1.471145e+04
     110  -1.470997e+04
     120  -1.470875e+04
     130  -1.470775e+04
     140  -1.470691e+04
     150  -1.470621e+04
     160  -1.470562e+04
     170  -1.470511e+04
     180  -1.470469e+04
     190  -1.470432e+04
     200  -1.470400e+04
     210  -1.470372e+04
     220  -1.470348e+04
     230  -1.470326e+04
     240  -1.470308e+04
     250  -1.470291e+04
     260  -1.470276e+04
     270  -1.470263e+04
     280  -1.470251e+04
     290  -1.470241e+04
     300  -1.470231e+04
     310  -1.470223e+04
     320  -1.470215e+04
     330  -1.470208e+04
     340  -1.470201e+04
     350  -1.470195e+04
     360  -1.470190e+04
     370  -1.470185e+04
     380  -1.470180e+04
     390  -1.470176e+04
     400  -1.470172e+04
     410  -1.470168e+04
     420  -1.470165e+04
     430  -1.470162e+04
     440  -1.470159e+04
     450  -1.470156e+04
     460  -1.470153e+04
     470  -1.470151e+04
     480  -1.470149e+04
     490  -1.470147e+04
     500  -1.470145e+04
     510  -1.470143e+04
     520  -1.470141e+04
     530  -1.470139e+04
     540  -1.470138e+04
     550  -1.470136e+04
     560  -1.470135e+04
     570  -1.470133e+04
     580  -1.470132e+04
     590  -1.470131e+04
     600  -1.470130e+04
     610  -1.470129e+04
     620  -1.470128e+04
     630  -1.470127e+04
     640  -1.470126e+04
     650  -1.470125e+04
     660  -1.470124e+04
     670  -1.470123e+04
     680  -1.470122e+04
     690  -1.470122e+04
     700  -1.470121e+04
     710  -1.470120e+04
     720  -1.470120e+04
     730  -1.470119e+04
     740  -1.470118e+04
     750  -1.470118e+04
     760  -1.470117e+04
     770  -1.470117e+04
     780  -1.470116e+04
     790  -1.470116e+04
     800  -1.470115e+04
     810  -1.470115e+04
     820  -1.470114e+04
     830  -1.470114e+04
     840  -1.470114e+04
     850  -1.470113e+04
     860  -1.470113e+04
     870  -1.470112e+04
     880  -1.470112e+04
     890  -1.470112e+04
     900  -1.470111e+04
     910  -1.470111e+04
     920  -1.470111e+04
     930  -1.470111e+04
     940  -1.470110e+04
     950  -1.470110e+04
     960  -1.470110e+04
     970  -1.470109e+04
     980  -1.470109e+04
     990  -1.470109e+04
    1000  -1.470109e+04
    1010  -1.470109e+04
    1020  -1.470108e+04
    1030  -1.470108e+04
    1040  -1.470108e+04
    1050  -1.470108e+04
    1060  -1.470108e+04
    1070  -1.470107e+04
    1080  -1.470107e+04
    1090  -1.470107e+04
    1100  -1.470107e+04
    1110  -1.470107e+04
    1120  -1.470107e+04

  0.794377 seconds (168.12 k allocations: 15.640 MiB, 0.80% gc time)
Out[55]:
VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}(Array{Float64}(0,3), ([0.2808 0.279454 0.232256; 0.279454 0.284312 0.219977; 0.232256 0.219977 0.213052], [0.717243 0.661816 0.674158; 0.661816 0.715193 0.581746; 0.674158 0.581746 0.783965]), Array{Float64}(0,0), Char[], Float64[], -Inf, Inf)

Heritability from 3-variate estimate and their standard errors.


In [56]:
h, hse = heritability(trait57_model.Σ, Σcov)
[h'; hse']


Out[56]:
2×3 Array{Float64,2}:
 0.281351   0.284453  0.213689
 0.0778252  0.077378  0.084084

13-trait joint analysis

In some situations, such as studying the genetic covariance, we need to jointly analyze 13 traits. We first try the Fisher scoring algorithm.


In [57]:
# initialize model parameters
traitall_model = VarianceComponentModel(cg10kdata_rotated)
# estimate variance components using Fisher scoring algorithm
@time mle_fs!(traitall_model, cg10kdata_rotated; solver=:Ipopt, verbose=true)


This is Ipopt version 3.12.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:    16653

Total number of variables............................:      182
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.3113371e+05 0.00e+00 1.00e+02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0 
   5  8.2233766e+04 0.00e+00 6.03e+02 -11.0 2.32e+00    -  1.00e+00 1.00e+00f  1 MaxS
  10  1.1960260e+05 0.00e+00 8.76e+02 -11.0 6.20e+01  -5.4 1.00e+00 1.00e+00h  1 MaxS
  15  2.4416551e+05 0.00e+00 2.50e+02 -11.0 8.69e+02  -7.8 1.00e+00 1.00e+00f  1 MaxS
DomainError:
log will only return a complex result if called with a complex argument. Try log(complex(x)).

Stacktrace:
 [1] nan_dom_err at ./math.jl:300 [inlined]
 [2] log at ./math.jl:419 [inlined]
 [3] logdet(::Array{Float64,2}) at ./linalg/generic.jl:1244
 [4] VarianceComponentModels.TwoVarCompModelRotate(::VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}) at /Users/huazhou/.julia/v0.6/VarianceComponentModels/src/VarianceComponentModels.jl:127
 [5] eval_f(::VarianceComponentModels.TwoVarCompOptProb{VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}},VarianceComponentModels.TwoVarCompVariateRotate{Float64,Array{Float64,2},Array{Float64,2}},Array{Float64,2},Array{Float64,1},VarianceComponentModels.VarianceComponentAuxData{Array{Float64,2},Array{Float64,1}}}, ::Array{Float64,1}) at /Users/huazhou/.julia/v0.6/VarianceComponentModels/src/two_variance_component.jl:683
 [6] (::Ipopt.#eval_f_cb#4{VarianceComponentModels.TwoVarCompOptProb{VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}},VarianceComponentModels.TwoVarCompVariateRotate{Float64,Array{Float64,2},Array{Float64,2}},Array{Float64,2},Array{Float64,1},VarianceComponentModels.VarianceComponentAuxData{Array{Float64,2},Array{Float64,1}}}})(::Array{Float64,1}) at /Users/huazhou/.julia/v0.6/Ipopt/src/IpoptSolverInterface.jl:53
 [7] eval_f_wrapper(::Int32, ::Ptr{Float64}, ::Int32, ::Ptr{Float64}, ::Ptr{Void}) at /Users/huazhou/.julia/v0.6/Ipopt/src/Ipopt.jl:89
 [8] solveProblem(::Ipopt.IpoptProblem) at /Users/huazhou/.julia/v0.6/Ipopt/src/Ipopt.jl:304
 [9] optimize!(::Ipopt.IpoptMathProgModel) at /Users/huazhou/.julia/v0.6/Ipopt/src/IpoptSolverInterface.jl:120
 [10] #mle_fs!#29(::Int64, ::Symbol, ::Symbol, ::Bool, ::Function, ::VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}, ::VarianceComponentModels.TwoVarCompVariateRotate{Float64,Array{Float64,2},Array{Float64,2}}) at /Users/huazhou/.julia/v0.6/VarianceComponentModels/src/two_variance_component.jl:893
 [11] (::VarianceComponentModels.#kw##mle_fs!)(::Array{Any,1}, ::VarianceComponentModels.#mle_fs!, ::VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}, ::VarianceComponentModels.TwoVarCompVariateRotate{Float64,Array{Float64,2},Array{Float64,2}}) at ./<missing>:0
 [12] include_string(::String, ::String) at ./loading.jl:515

From the output we can see the Fisher scoring algorithm ran into some numerical issues. Let's try the MM algorithm.


In [58]:
# reset model parameters
traitall_model = VarianceComponentModel(cg10kdata_rotated)
# estimate variance components using Fisher scoring algorithm
@time mle_mm!(traitall_model, cg10kdata_rotated; verbose=true)


     MM Algorithm
  Iter      Objective  
--------  -------------
       0  -1.311337e+05
       1  -8.002108e+04
       2  -5.806935e+04
       3  -4.926111e+04
       4  -4.611059e+04
       5  -4.511606e+04
       6  -4.482679e+04
       7  -4.474294e+04
       8  -4.471496e+04
       9  -4.470174e+04
      10  -4.469246e+04
      20  -4.462243e+04
      30  -4.456888e+04
      40  -4.452774e+04
      50  -4.449601e+04
      60  -4.447134e+04
      70  -4.445199e+04
      80  -4.443665e+04
      90  -4.442436e+04
     100  -4.441442e+04
     110  -4.440630e+04
     120  -4.439961e+04
     130  -4.439405e+04
     140  -4.438938e+04
     150  -4.438544e+04
     160  -4.438210e+04
     170  -4.437923e+04
     180  -4.437676e+04
     190  -4.437463e+04
     200  -4.437277e+04
     210  -4.437115e+04
     220  -4.436972e+04
     230  -4.436846e+04
     240  -4.436735e+04
     250  -4.436636e+04
     260  -4.436548e+04
     270  -4.436469e+04
     280  -4.436399e+04
     290  -4.436335e+04
     300  -4.436278e+04
     310  -4.436226e+04
     320  -4.436179e+04
     330  -4.436137e+04
     340  -4.436098e+04
     350  -4.436063e+04
     360  -4.436030e+04
     370  -4.436001e+04
     380  -4.435974e+04
     390  -4.435949e+04
     400  -4.435926e+04
     410  -4.435905e+04
     420  -4.435886e+04
     430  -4.435868e+04
     440  -4.435851e+04
     450  -4.435836e+04
     460  -4.435822e+04
     470  -4.435809e+04
     480  -4.435797e+04
     490  -4.435785e+04
     500  -4.435775e+04
     510  -4.435765e+04
     520  -4.435756e+04
     530  -4.435747e+04
     540  -4.435739e+04
     550  -4.435732e+04
     560  -4.435725e+04
     570  -4.435718e+04
     580  -4.435712e+04
     590  -4.435706e+04
     600  -4.435701e+04
     610  -4.435696e+04
     620  -4.435691e+04
     630  -4.435687e+04
     640  -4.435683e+04
     650  -4.435679e+04
     660  -4.435675e+04
     670  -4.435671e+04
     680  -4.435668e+04
     690  -4.435665e+04
     700  -4.435662e+04
     710  -4.435659e+04
     720  -4.435657e+04
     730  -4.435654e+04
     740  -4.435652e+04
     750  -4.435649e+04
     760  -4.435647e+04
     770  -4.435645e+04
     780  -4.435643e+04
     790  -4.435642e+04
     800  -4.435640e+04
     810  -4.435638e+04
     820  -4.435637e+04
     830  -4.435635e+04
     840  -4.435634e+04
     850  -4.435633e+04
     860  -4.435631e+04
     870  -4.435630e+04
     880  -4.435629e+04
     890  -4.435628e+04
     900  -4.435627e+04
     910  -4.435626e+04
     920  -4.435625e+04
     930  -4.435624e+04
     940  -4.435623e+04
     950  -4.435622e+04
     960  -4.435621e+04
     970  -4.435621e+04
     980  -4.435620e+04
     990  -4.435619e+04
    1000  -4.435619e+04
    1010  -4.435618e+04
    1020  -4.435617e+04
    1030  -4.435617e+04
    1040  -4.435616e+04
    1050  -4.435616e+04
    1060  -4.435615e+04
    1070  -4.435615e+04
    1080  -4.435614e+04
    1090  -4.435614e+04

  3.551301 seconds (178.42 k allocations: 70.115 MiB, 0.42% gc time)
Out[58]:
(-44356.138529861186, VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}(Array{Float64}(0,13), ([0.272384 0.190358 … -0.128222 -0.0980655; 0.190358 0.21692 … -0.0689912 -0.0444349; … ; -0.128222 -0.0689912 … 0.118227 0.0909188; -0.0980655 -0.0444349 … 0.0909188 0.107456], [0.724562 0.56992 … -0.0590518 -0.124939; 0.56992 0.782639 … 0.0238629 0.0475408; … ; -0.0590518 0.0238629 … 0.880671 0.550889; -0.124939 0.0475408 … 0.550889 0.891929]), Array{Float64}(0,0), Char[], Float64[], -Inf, Inf), ([0.0111619 0.0131088 … 0.0128956 0.0127641; 0.0131091 0.0151759 … 0.017162 0.0171466; … ; 0.0128956 0.017162 … 0.0173994 0.0182002; 0.0127643 0.0171461 … 0.0182003 0.0187848], [0.0112235 0.0133094 … 0.0130111 0.0127861; 0.01331 0.0158262 … 0.017867 0.017798; … ; 0.013011 0.0178666 … 0.0179487 0.0187579; 0.012786 0.0177975 … 0.0187578 0.0193328]), [0.000124587 7.24074e-5 … -3.35716e-7 -1.40982e-5; 7.24411e-5 0.000171849 … -2.05381e-5 -3.17975e-6; … ; -3.60221e-7 -2.05683e-5 … 0.000351859 -1.5168e-5; -1.40799e-5 -3.16738e-6 … -1.51641e-5 0.000373756], Array{Float64}(0,13), Array{Float64}(0,0))

It converges after ~1000 iterations.

Save analysis results


In [59]:
#using JLD
#@save "copd.jld"
#whos()